Go to the documentation of this file.
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> {
1038 virtual ~fit_nonlin() {
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;
1081 std::cout <<
"Press a key and type enter to continue. ";
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;
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++) {
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);
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;
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);
1253 bool lm_iteration=
true;
1254 while (lm_iteration) {
1264 for(
size_t i=0;i<
nparm;i++) dx[i]*=-1.0;
1294 double t1=fnorm1p/
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 ) {
1324 else delta=temp*pnorm/p1;
1335 double dbl_eps=std::numeric_limits<double>::epsilon();
1337 if (ratio >= p0001) {
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) {
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);
1408 o2scl_cblas::o2cblas_NoTrans,
1409 fitfun.get_ndata(),npar,1.0,
J_,
f_,0.0,
g_);
1417 if (niter>=this->ntrial) {
1419 "exceeded max. number of iterations.",
1442 virtual const char *
type() {
return "fit_nonlin"; }
1446 #ifndef DOXYGEN_NO_O2NS
size_t iter
Number of iterations.
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 .
size_t nparm
Number of parameters.
A class for representing permutations.
void compute_newton_direction(size_t n, const mat_t &r2, const permutation &perm2, const vec_t &qtf2, vec_t &x)
Desc.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
double tol_rel
(default 1.0e-4)
Class for automatically computing gradients [abstract base].
int test_gradient_f(size_t nparm, vec_t &g, double l_epsabs)
Test if converged.
bool test_gradient
If true, test the gradient also (default false)
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.
@ exc_efailed
generic failure
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_ERR2(d, d2, n)
Set an error, two-string version.
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.
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.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
@ exc_etolf
cannot reach the specified tolerance in f
vec_t f_trial
Trial function value.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
virtual const char * type()
Return string denoting type ("fit_nonlin")
@ exc_emaxiter
exceeded max number of iterations
int allocate(size_t dim)
Allocate memory for a permutation of size dim.
int test_delta_f(size_t nparm, vec_t &dx, vec_t &x, double l_epsabs, double l_epsrel)
Test if converged.
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.
double tol_abs
Absolute tolerance (default 1.0e-4)
size_t ndata
Number of data points.
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.
void free()
Free allocated 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.
func_t * cff
Function to fit.
size_t count_nsing(const size_t ncols, const mat_t &r2)
Desc.
@ exc_einval
invalid argument supplied by user
int set(size_t npar, vec_t &parms, func_t &fitfun)
Set the initial values of the parameters and the fitting function to use for the next call to iterate...
@ gsl_continue
iteration has not converged
@ exc_etolg
cannot reach the specified tolerance in gradient gradient
Base routines for the nonlinear fitting classes.
int apply_inverse(vec_t &v) const
Apply the inverse permutation to a vector.
double compute_actual_reduction(double fnorm0, double fnorm1)
Desc.
int iterate()
Perform an iteration.
int free()
Free the memory.
double tol_rel_covar
The relative tolerance for the computation of the covariance matrix (default 0)
Non-linear least-squares fitting [abstract base].
@ exc_etolx
cannot reach the specified tolerance in x
vec_t dx_
The last step taken in parameter space.
size_t ntrial
Maximum number of iterations (default 500)
@ exc_esanity
sanity check failed - shouldn't happen
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
double compute_delta(vec_t &diag2, size_t n, const vec_t &x)
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.
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.
int verbose
An integer describing the verbosity of the output.
Non-linear least-squares fitting class (GSL)
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.
void resize(size_t n, size_t p)
Allocate memory with n data points and p parameters.
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
bool err_nonconv
If true, call the error handler if fit() does not converge (default true)
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.
void compute_trial_step(size_t N, vec_t &x, vec_t &dx, vec_t &trial)
Compute trial step, .
void compute_rptdx(const mat_t &r2, const permutation &p, size_t N, vec_t &dx, vec_t &rptdx2)
Desc.
bool use_scaled
Use the scaled routine if true (default true)
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.
void update_diag(size_t n, const mat_t &J, vec_t &diag2)
Desc.
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).