42 #ifndef O2SCL_FIT_NONLIN_H 43 #define O2SCL_FIT_NONLIN_H 49 #include <boost/numeric/ublas/vector.hpp> 50 #include <boost/numeric/ublas/matrix.hpp> 52 #include <o2scl/vector.h> 53 #include <o2scl/fit_base.h> 54 #include <o2scl/permutation.h> 55 #include <o2scl/cblas.h> 57 #include <o2scl/qrpt.h> 58 #include <o2scl/columnify.h> 60 #ifndef DOXYGEN_NO_O2NS 66 template<
class vec_t=boost::numeric::ublas::vector<
double>,
67 class mat_t=boost::numeric::ublas::matrix<
double> >
class fit_nonlin_b {
76 if (0.1*fnorm1<fnorm0) {
77 double u=fnorm1/fnorm0;
91 for (i=0;i<ncols;i++) {
105 const vec_t &qtf2, vec_t &x) {
116 for (i=nsing;i<n;i++) {
121 for (j=nsing;j>0 && j--;) {
123 double temp=x[j]/rjj;
142 (
size_t nd,
size_t np,
const mat_t &r2,
const vec_t &x,
double dxnorm,
143 const permutation &perm,
const vec_t &diag, vec_t &w) {
160 w[i]=dpi*(dpi*xpi/dxnorm);
182 (
size_t n,
const mat_t &r,
const permutation &p,
const vec_t &qtf2,
183 const vec_t &diag, vec_t &g) {
190 for (i=0;i <= j;i++) {
191 sum +=r(i,j)*qtf2[i];
208 for (
size_t j=0;j<n;j++) {
209 double cnorm,diagj,sum=0;
210 for (
size_t i=0;i<n;i++) {
233 for (
size_t i=0;i<n;i++) {
247 return (Dx>0) ? factor*Dx : factor;
253 size_t N, vec_t &dx, vec_t &rptdx2) {
255 for (
size_t i=0;i<N;i++) {
257 for (
size_t j=i;j<N;j++) {
258 sum+=r2(i,j)*dx[p[j]];
311 const double lambda,
const vec_t &diag2,
312 const vec_t &qtb, vec_t &x, vec_t &sdiag2, vec_t &wa) {
314 size_t i, j, k, nsing;
323 for (i=j+1;i<n;i++) {
338 double diagpj=lambda*diag2[pj];
345 for (k=j+1;k<n;k++) {
364 double sdiagk=sdiag2[k];
370 if (fabs(rkk)<fabs(sdiagk)) {
371 double cotangent=rkk/sdiagk;
372 sine=0.5/sqrt(0.25+0.25*cotangent*cotangent);
373 cosine=sine*cotangent;
375 double tangent=sdiagk/rkk;
376 cosine=0.5/sqrt(0.25+0.25*tangent*tangent);
384 double new_rkk=cosine*rkk+sine*sdiagk;
385 double new_wak=cosine*wak+sine*qtbpj;
387 qtbpj=-sine*wak+cosine*qtbpj;
395 for (i=k+1;i<n;i++) {
397 double sdiagi=sdiag2[i];
399 double new_rik=cosine*rik+sine*sdiagi;
400 double new_sdiagi=-sine*rik+cosine*sdiagi;
403 sdiag2[i]=new_sdiagi;
425 double sdiagj=sdiag2[j];
433 for (j=nsing;j<n;j++) {
437 for (k=0;k<nsing;k++) {
442 for (i=j+1;i<nsing;i++) {
448 double sdiagj=sdiag2[j];
450 wa[j]=(waj-sum)/sdiagj;
469 (
size_t n,
const mat_t &r2,
const vec_t &sdiag2,
const permutation &p,
470 vec_t &x,
double dxnorm,
const vec_t &diag2, vec_t &w2) {
472 for (
size_t i=0;i<n;i++) {
474 double dpi=diag2[
pi];
477 w2[i]=dpi*(dpi*xpi)/dxnorm;
480 for (
size_t j=0;j<n;j++) {
487 for (
size_t i=j+1;i<n;i++) {
500 const vec_t &diag2,
double delta2,
double *par_inout,
501 vec_t &newton2, vec_t &gradient2, vec_t &sdiag2,
502 vec_t &x, vec_t &w2,
size_t nparm,
size_t ndata) {
504 double dxnorm, gnorm, fp, fp_old, par_lower, par_upper, par_c;
506 double par2=*par_inout;
519 if (fp<=0.1*delta2) {
521 for(
size_t i=0;i<nparm;i++) {
534 double phider=wnorm*wnorm;
540 par_lower=fp/(delta2*phider);
550 par_upper=gnorm/delta2;
552 if (par_upper == 0) {
554 if (delta2<0.1) mint=delta2;
556 double dbl_min=std::numeric_limits<double>::min();
557 par_upper=dbl_min/mint;
560 if (par2>par_upper) {
562 }
else if (par2<par_lower) {
576 #ifdef O2SCL_NEVER_DEFINED 581 if (par2<par_lower || par2>par_upper) {
582 par2=GSL_MAX_DBL(0.001*par_upper,sqrt(par_lower*par_upper));
590 par2=GSL_MAX_DBL(0.001*par_upper,GSL_DBL_MIN);
597 double sqrt_par=sqrt(par2);
599 qrsolv(nparm,r2,perm2,sqrt_par,diag2,qtf2,x,sdiag2,w2);
611 if (fabs(fp)<=0.1*delta2) {
615 if (par_lower == 0 && fp<=fp_old && fp_old<0) {
633 par_c=fp/(delta2*wnorm*wnorm);
640 if (par2>par_lower) {
644 if (par2<par_upper) {
651 par2=GSL_MAX_DBL(par_lower,par2+par_c);
655 O2SCL_ERR(
"Sanity check failed in fit_nonlin_b::lmpar().",
663 for (
size_t i=0;i<N;i++) {
681 const mat_t &J, vec_t &diag_vec) {
683 for (
size_t j=0;j<nparm;j++) {
686 for (
size_t i=0;i<ndata;i++) {
693 diag_vec[j]=sqrt(sum);
722 mat_t &covar, vec_t &norm, mat_t &r,
736 for(
size_t i=0;i<m;i++) {
737 for(
size_t j=0;j<n;j++) {
745 tolr=epsrel*fabs(r(0,0));
747 for (
size_t k=0;k<n;k++) {
750 if (fabs(rkk)<=tolr) {
756 for (
size_t j=0;j<k;j++) {
760 for (
size_t i=0;i<=j;i++) {
772 for (
size_t k=0;k<=kmax;k++) {
773 for (
size_t j=0;j<k;j++) {
776 for (
size_t i=0;i<=j;i++) {
785 for (
size_t i=0;i<=k;i++) {
793 for (
size_t j=0;j<n;j++) {
796 for (
size_t i=0;i<=j;i++) {
823 for (
size_t j=0;j<n;j++) {
824 for (
size_t i=0;i<j;i++) {
847 double l_epsabs,
double l_epsrel) {
850 O2SCL_ERR2(
"Relative tolerance less than zero ",
851 "in fit_nonlin_b::test_delta_f().",
exc_einval);
854 for(
size_t i=0;i<nparm;i++) {
855 double tolerance=l_epsabs+l_epsrel*fabs(x[i]);
856 if (fabs(dx[i])>=tolerance) {
869 O2SCL_ERR2(
"Absolute tolerance less than zero ",
870 "in fit_nonlin_b::test_gradient_f().",
exc_einval);
873 for(
size_t i=0;i<nparm;i++) {
874 residual+=fabs(g[i]);
876 if (residual<l_epsabs) {
918 template<
class func_t=gen_fit_funct<>,
919 class vec_t=boost::numeric::ublas::vector<
double>,
920 class mat_t=boost::numeric::ublas::matrix<
double> >
922 public fit_base<func_t,vec_t,mat_t> {
1032 test_gradient=
false;
1048 virtual int print_iter(
size_t nv, vec_t &x, vec_t &dx,
int iter2,
1049 double l_epsabs,
double l_epsrel) {
1051 if (this->verbose<=0)
return 0;
1056 double val=0.0, lim=0.0, diff_max=0.0;
1057 for(
size_t i=0;i<nv;i++) {
1058 double tolerance=l_epsabs+l_epsrel*fabs(x[i]);
1059 if (fabs(dx[i])>=tolerance) {
1064 double diff=tolerance-fabs(dx[i]);
1065 if (diff>diff_max) {
1073 std::cout <<
"Iteration: " << iter2 << std::endl;
1075 for(
size_t i=0;i<nv;i++) {
1076 std::cout << x[i] <<
" ";
1078 std::cout << std::endl;
1079 std::cout <<
" Val: " << val <<
" Lim: " << lim << std::endl;
1080 if (this->verbose>1) {
1081 std::cout <<
"Press a key and type enter to continue. ";
1098 r.resize(ndata,nparm);
1101 for(
size_t i=0;i<ndata;i++) tau[i]=0.0;
1104 for(
size_t i=0;i<nparm;i++) tau[i]=0.0;
1107 rptdx.resize(ndata);
1111 work1.resize(nparm);
1112 newton.resize(nparm);
1113 gradient.resize(nparm);
1114 sdiag.resize(nparm);
1115 for(
size_t i=0;i<ndata;i++) {
1116 for(
size_t j=0;j<nparm;j++) {
1124 for(
size_t i=0;i<nparm;i++) {
1132 x_trial.resize(nparm);
1133 f_trial.resize(ndata);
1137 J_.resize(ndata,nparm);
1152 int set(
size_t npar, vec_t &parms, func_t &fitfun) {
1156 if (fitfun.get_ndata()==0 || npar==0) {
1157 O2SCL_ERR2(
"Either zero data or zero parameters in ",
1161 if (ndata!=fitfun.get_ndata() || nparm!=npar) {
1162 resize(fitfun.get_ndata(),npar);
1170 (*cff)(nparm,*x_,ndata,f_);
1171 cff->jac(nparm,*x_,ndata,f_,J_);
1177 for(
size_t i=0;i<nparm;i++) {
1186 for(
size_t i=0;i<nparm;i++) {
1198 for(
size_t i=0;i<ndata;i++) {
1199 for(
size_t j=0;j<nparm;j++) {
1206 for(
size_t ii=0;ii<rptdx.size();ii++) rptdx[ii]=0.0;
1207 for(
size_t ii=0;ii<w.size();ii++) w[ii]=0.0;
1211 for(
size_t i=0;i<ndata;i++) f_trial[i]=0.0;
1225 double prered, actred;
1226 double pnorm, fnorm1, fnorm1p, gnorm;
1232 double p1=0.1, p25=0.25, p5=0.5, p75=0.75, p0001=0.0001;
1240 for(
size_t i=0;i<ndata;i++) {
1248 size_t iamax=vector_max_index<vec_t,double>(nparm,gradient);
1249 gnorm=fabs(gradient[iamax]/fnorm);
1253 bool lm_iteration=
true;
1254 while (lm_iteration) {
1258 this->
lmpar(r,perm,qtf,diag,delta,&par,newton,
1259 gradient,sdiag,dx,w,nparm,ndata);
1264 for(
size_t i=0;i<nparm;i++) dx[i]*=-1.0;
1277 (*cff)(nparm,x_trial,ndata,f_trial);
1294 double t1=fnorm1p/fnorm;
1295 double t2=(sqrt(par)*pnorm)/fnorm;
1297 prered=t1*t1+t2*t2/p5;
1298 dirder=-(t1*t1+t2*t2);
1303 ratio=actred/prered;
1311 if (par == 0 || ratio >= p75) {
1316 double temp=(actred >= 0) ? p5 : p5*dirder/
1319 if (p1*fnorm1 >= fnorm || temp<p1 ) {
1323 if (delta<pnorm/p1) delta=temp*delta;
1324 else delta=temp*pnorm/p1;
1335 double dbl_eps=std::numeric_limits<double>::epsilon();
1337 if (ratio >= p0001) {
1342 cff->jac(nparm,x_trial,ndata,f_trial,J_);
1357 for(
size_t i=0;i<ndata;i++) {
1358 for(
size_t j=0;j<nparm;j++) {
1366 }
else if (fabs(actred)<=dbl_eps &&
1367 prered<=dbl_eps && p5*ratio<=1.0) {
1369 }
else if (delta<=dbl_eps*xnorm) {
1371 }
else if (gnorm<=dbl_eps) {
1373 }
else if (niter<10) {
1388 virtual int fit(
size_t npar, vec_t &parms, mat_t &covar,
1389 double &chi2, func_t &fitfun) {
1391 set(npar,parms,fitfun);
1404 status=this->
test_delta_f(nparm,dx_,parms,this->tol_abs,
1408 o2scl_cblas::o2cblas_NoTrans,
1409 fitfun.get_ndata(),npar,1.0,J_,f_,0.0,g_);
1413 this->print_iter(npar,parms,dx_,niter,this->tol_abs,this->tol_rel);
1415 }
while (status ==
gsl_continue && niter<this->ntrial);
1417 if (niter>=this->ntrial) {
1419 "exceeded max. number of iterations.",
1424 this->
covariance(fitfun.get_ndata(),npar,J_,covar,work1,
1442 virtual const char *
type() {
return "fit_nonlin"; }
1446 #ifndef DOXYGEN_NO_O2NS bool test_gradient
If true, test the gradient also (default false)
double compute_delta(vec_t &diag2, size_t n, const vec_t &x)
Desc.
int test_delta_f(size_t nparm, vec_t &dx, vec_t &x, double l_epsabs, double l_epsrel)
Test if converged.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
void lmpar(mat_t &r2, const permutation &perm2, const vec_t &qtf2, const vec_t &diag2, double delta2, double *par_inout, vec_t &newton2, vec_t &gradient2, vec_t &sdiag2, vec_t &x, vec_t &w2, size_t nparm, size_t ndata)
Determine Levenburg-Marquardt parameter.
bool err_nonconv
If true, call the error handler if fit() does not converge (default true)
void compute_gradient_direction(size_t n, const mat_t &r, const permutation &p, const vec_t &qtf2, const vec_t &diag, vec_t &g)
Desc.
int compute_diag(size_t nparm, size_t ndata, const mat_t &J, vec_t &diag_vec)
Compute the root of the sum of the squares of the columns of J.
sanity check failed - shouldn't happen
Non-linear least-squares fitting class (GSL)
invalid argument supplied by user
A class for representing permutations.
exceeded max number of iterations
virtual int print_iter(size_t nv, vec_t &x, vec_t &dx, int iter2, double l_epsabs, double l_epsrel)
Print the progress in the current iteration.
vec_t dx_
The last step taken in parameter space.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
size_t nparm
Number of parameters.
void compute_trial_step(size_t N, vec_t &x, vec_t &dx, vec_t &trial)
Compute trial step, .
void QRPT_decomp(size_t M, size_t N, mat_t &A, vec_t &tau, o2scl::permutation &p, int &signum, vec2_t &norm)
Compute the QR decomposition of matrix A.
iteration has not converged
int free()
Free the memory.
void compute_newton_correction(size_t n, const mat_t &r2, const vec_t &sdiag2, const permutation &p, vec_t &x, double dxnorm, const vec_t &diag2, vec_t &w2)
Desc.
int apply_inverse(vec_t &v) const
Apply the inverse permutation to a vector.
size_t ndata
Number of data points.
int qrsolv(size_t n, mat_t &r2, const permutation &p, const double lambda, const vec_t &diag2, const vec_t &qtb, vec_t &x, vec_t &sdiag2, vec_t &wa)
Compute the solution to a least squares system.
Base routines for the nonlinear fitting classes.
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
void dgemv(const enum o2cblas_order order, const enum o2cblas_transpose TransA, const size_t M, const size_t N, const double alpha, const mat_t &A, const vec_t &X, const double beta, vec2_t &Y)
Compute .
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
vec_t f_trial
Trial function value.
size_t iter
Number of iterations.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
virtual int fit(size_t npar, vec_t &parms, mat_t &covar, double &chi2, func_t &fitfun)
Fit the data specified in (xdat,ydat) to the function fitfun with the parameters in par...
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void compute_rptdx(const mat_t &r2, const permutation &p, size_t N, vec_t &dx, vec_t &rptdx2)
Desc.
cannot reach the specified tolerance in x
cannot reach the specified tolerance in gradient gradient
size_t count_nsing(const size_t ncols, const mat_t &r2)
Desc.
void compute_newton_direction(size_t n, const mat_t &r2, const permutation &perm2, const vec_t &qtf2, vec_t &x)
Desc.
void resize(size_t n, size_t p)
Allocate memory with n data points and p parameters.
void free()
Free allocated memory.
double compute_actual_reduction(double fnorm0, double fnorm1)
Desc.
int iterate()
Perform an iteration.
Non-linear least-squares fitting [abstract base].
bool use_scaled
Use the scaled routine if true (default true)
void update_diag(size_t n, const mat_t &J, vec_t &diag2)
Desc.
void compute_newton_bound(size_t nd, size_t np, const mat_t &r2, const vec_t &x, double dxnorm, const permutation &perm, const vec_t &diag, vec_t &w)
Desc.
void QR_QTvec(const size_t M, const size_t N, const mat_t &QR, const vec_t &tau, vec2_t &v)
Form the product Q^T v from a QR factorized matrix.
int allocate(size_t dim)
Allocate memory for a permutation of size dim.
int covariance(size_t m, size_t n, const mat_t &J, mat_t &covar, vec_t &norm, mat_t &r, vec_t &tau, permutation &perm, double epsrel)
Compute the covarance matrix covar given the Jacobian J.
double scaled_enorm(const vec_t &d, size_t n, const vec_t &f)
Euclidean norm of vector f of length n, scaled by vector d.
virtual const char * type()
Return string denoting type ("fit_nonlin")
int test_gradient_f(size_t nparm, vec_t &g, double l_epsabs)
Test if converged.
func_t * cff
Function to fit.
double tol_rel_covar
The relative tolerance for the computation of the covariance matrix (default 0)
cannot reach the specified tolerance in f