fit_nonlin.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /*
24  * Based on the multifit routines in GSL
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Brian Gough
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 #ifndef O2SCL_FIT_NONLIN_H
43 #define O2SCL_FIT_NONLIN_H
44 
45 /** \file fit_nonlin.h
46  \brief File defining \ref o2scl::fit_nonlin
47 */
48 
49 #include <boost/numeric/ublas/vector.hpp>
50 #include <boost/numeric/ublas/matrix.hpp>
51 
52 #include <o2scl/vector.h>
53 #include <o2scl/fit_base.h>
54 #include <o2scl/permutation.h>
55 #include <o2scl/cblas.h>
56 #include <o2scl/qr.h>
57 #include <o2scl/qrpt.h>
58 #include <o2scl/columnify.h>
59 
60 #ifndef DOXYGEN_NO_O2NS
61 namespace o2scl {
62 #endif
63 
64  /** \brief Base routines for the nonlinear fitting classes
65  */
66  template<class vec_t=boost::numeric::ublas::vector<double>,
67  class mat_t=boost::numeric::ublas::matrix<double> > class fit_nonlin_b {
68 
69  protected:
70 
71  /** \brief Desc
72  */
73  double compute_actual_reduction(double fnorm0, double fnorm1) {
74  double actred;
75 
76  if (0.1*fnorm1<fnorm0) {
77  double u=fnorm1/fnorm0;
78  actred=1-u*u;
79  } else {
80  actred=-1;
81  }
82 
83  return actred;
84  }
85 
86  /** \brief Desc
87  */
88  size_t count_nsing(const size_t ncols, const mat_t &r2) {
89 
90  size_t i;
91  for (i=0;i<ncols;i++) {
92  double rii=r2(i,i);
93  if (rii == 0) {
94  break;
95  }
96  }
97 
98  return i;
99  }
100 
101  /** \brief Desc
102  */
104  (size_t n, const mat_t &r2, const permutation &perm2,
105  const vec_t &qtf2, vec_t &x) {
106 
107  size_t i, j, nsing;
108 
109  for (i=0;i<n;i++) {
110  double qtfi=qtf2[i];
111  x[i]=qtfi;
112  }
113 
114  nsing=count_nsing(n,r2);
115 
116  for (i=nsing;i<n;i++) {
117  x[i]=0.0;
118  }
119 
120  if (nsing > 0) {
121  for (j=nsing;j>0 && j--;) {
122  double rjj=r2(j,j);
123  double temp=x[j]/rjj;
124 
125  x[j]=temp;
126 
127  for (i=0;i<j;i++) {
128  double rij=r2(i,j);
129  double xi=x[i];
130  x[i]=xi-rij*temp;
131  }
132  }
133  }
134 
135  perm2.apply_inverse(x);
136  return;
137  }
138 
139  /** \brief Desc
140  */
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) {
144 
145  size_t i, j;
146 
147  size_t nsing=count_nsing(np,r2);
148 
149  if (nsing<np) {
150  for(i=0;i<nd;i++) {
151  w[i]=0.0;
152  }
153  return;
154  }
155 
156  for (i=0;i<np;i++) {
157  size_t pi=perm[i];
158  double dpi=diag[pi];
159  double xpi=x[pi];
160  w[i]=dpi*(dpi*xpi/dxnorm);
161  }
162 
163  for (j=0;j<np;j++) {
164  double sum=0.0;
165 
166  for (i=0;i<j;i++) {
167  sum +=r2(i,j)*w[i];
168  }
169 
170  {
171  double rjj=r2(j,j);
172  double wj=w[j];
173  w[j]=(wj-sum)/rjj;
174  }
175  }
176  return;
177  }
178 
179  /** \brief Desc
180  */
182  (size_t n, const mat_t &r, const permutation &p, const vec_t &qtf2,
183  const vec_t &diag, vec_t &g) {
184 
185  size_t i, j;
186 
187  for (j=0;j<n;j++) {
188  double sum=0;
189 
190  for (i=0;i <= j;i++) {
191  sum +=r(i,j)*qtf2[i];
192  }
193 
194  {
195  size_t pj=p[j];
196  double dpj=diag[pj];
197 
198  g[j]=sum/dpj;
199  }
200  }
201  return;
202  }
203 
204  /** \brief Desc
205  */
206  void update_diag(size_t n, const mat_t &J, vec_t &diag2) {
207 
208  for (size_t j=0;j<n;j++) {
209  double cnorm,diagj,sum=0;
210  for (size_t i=0;i<n;i++) {
211  double Jij=J(i,j);
212  sum+=Jij*Jij;
213  }
214  if (sum == 0) {
215  sum=1.0;
216  }
217 
218  cnorm=sqrt(sum);
219  diagj=diag2[j];
220 
221  if (cnorm>diagj) {
222  diag2[j]=cnorm;
223  }
224  }
225  return;
226  }
227 
228  /** \brief Euclidean norm of vector \c f of length \c n, scaled by
229  vector \c d
230  */
231  double scaled_enorm(const vec_t &d, size_t n, const vec_t &f) {
232  double e2=0;
233  for (size_t i=0;i<n;i++) {
234  double di=d[i];
235  double u=di*f[i];
236  e2+=u*u;
237  }
238  return sqrt(e2);
239  }
240 
241  /** \brief Desc
242  */
243  double compute_delta(vec_t &diag2, size_t n, const vec_t &x) {
244  double Dx=scaled_enorm(diag2,n,x);
245  // [GSL] The generally recommended value from MINPACK is 100
246  double factor=100;
247  return (Dx>0) ? factor*Dx : factor;
248  }
249 
250  /** \brief Desc
251  */
252  void compute_rptdx(const mat_t &r2, const permutation &p,
253  size_t N, vec_t &dx, vec_t &rptdx2) {
254 
255  for (size_t i=0;i<N;i++) {
256  double sum=0.0;
257  for (size_t j=i;j<N;j++) {
258  sum+=r2(i,j)*dx[p[j]];
259  }
260  rptdx2[i]=sum;
261  }
262 
263  return;
264  }
265 
266  /** \brief Compute the solution to a least squares system
267 
268  \verbatim
269  This function computes the solution to the least squares system
270 
271  phi=[ A x=b ,lambda D x=0 ]^2
272 
273  where A is an M by N matrix,D is an N by N diagonal matrix,lambda
274  is a scalar parameter and b is a vector of length M.
275 
276  The function requires the factorization of A into A=Q R P^T,
277  where Q is an orthogonal matrix,R is an upper triangular matrix
278  with diagonal elements of non-increasing magnitude and P is a
279  permuation matrix. The system above is then equivalent to
280 
281  [ R z=Q^T b,P^T (lambda D) P z=0 ]
282 
283  where x=P z. If this system does not have full rank then a least
284  squares solution is obtained. On output the function also provides
285  an upper triangular matrix S such that
286 
287  P^T (A^T A+lambda^2 D^T D) P=S^T S
288 
289  Parameters,
290 
291  r: On input,contains the full upper triangle of R. On output the
292  strict lower triangle contains the transpose of the strict upper
293  triangle of S,and the diagonal of S is stored in sdiag. The full
294  upper triangle of R is not modified.
295 
296  p: the encoded form of the permutation matrix P. column j of P is
297  column p[j] of the identity matrix.
298 
299  lambda,diag: contains the scalar lambda and the diagonal elements
300  of the matrix D
301 
302  qtb: contains the product Q^T b
303 
304  x: on output contains the least squares solution of the system
305 
306  wa: is a workspace of length N
307  \endverbatim
308 
309  */
310  int qrsolv(size_t n, mat_t &r2, const permutation &p,
311  const double lambda, const vec_t &diag2,
312  const vec_t &qtb, vec_t &x, vec_t &sdiag2, vec_t &wa) {
313 
314  size_t i, j, k, nsing;
315 
316  // [GSL] Copy r and qtb to preserve input and initialise s. In
317  // particular, save the diagonal elements of r in x
318 
319  for (j=0;j<n;j++) {
320  double rjj=r2(j,j);
321  double qtbj=qtb[j];
322 
323  for (i=j+1;i<n;i++) {
324  r2(i,j)=r2(j,i);
325  }
326 
327  x[j]=rjj;
328  wa[j]=qtbj;
329  }
330 
331  // [GSL] Eliminate the diagonal matrix d using a Givens rotation
332 
333  for (j=0;j<n;j++) {
334 
335  double qtbpj;
336 
337  size_t pj=p[j];
338  double diagpj=lambda*diag2[pj];
339  if (diagpj == 0) {
340  continue;
341  }
342 
343  sdiag2[j]=diagpj;
344 
345  for (k=j+1;k<n;k++) {
346  sdiag2[k]=0.0;
347  }
348 
349  // [GSL] The transformations to eliminate the row of d modify only a
350  // single element of qtb beyond the first n, which is initially
351  // zero
352 
353  qtbpj=0;
354 
355  for (k=j;k<n;k++) {
356 
357  // [GSL] Determine a Givens rotation which eliminates the
358  // appropriate element in the current row of d
359 
360  double sine, cosine;
361 
362  double wak=wa[k];
363  double rkk=r2(k,k);
364  double sdiagk=sdiag2[k];
365 
366  if (sdiagk == 0) {
367  continue;
368  }
369 
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;
374  } else {
375  double tangent=sdiagk/rkk;
376  cosine=0.5/sqrt(0.25+0.25*tangent*tangent);
377  sine=cosine*tangent;
378  }
379 
380  // [GSL] Compute the modified diagonal element of r and the
381  // modified element of [qtb,0]
382 
383  {
384  double new_rkk=cosine*rkk+sine*sdiagk;
385  double new_wak=cosine*wak+sine*qtbpj;
386 
387  qtbpj=-sine*wak+cosine*qtbpj;
388 
389  r2(k,k)=new_rkk;
390  wa[k]=new_wak;
391  }
392 
393  // [GSL] Accumulate the transformation in the row of s
394 
395  for (i=k+1;i<n;i++) {
396  double rik=r2(i,k);
397  double sdiagi=sdiag2[i];
398 
399  double new_rik=cosine*rik+sine*sdiagi;
400  double new_sdiagi=-sine*rik+cosine*sdiagi;
401 
402  r2(i,k)=new_rik;
403  sdiag2[i]=new_sdiagi;
404  }
405  }
406 
407  // [GSL] Store the corresponding diagonal element of s and
408  // restore the corresponding diagonal element of r
409  {
410  double rjj=r2(j,j);
411  double xj=x[j];
412 
413  sdiag2[j]=rjj;
414  r2(j,j)=xj;
415  }
416 
417  }
418 
419  // [GSL] Solve the triangular system for z. If the system is singular
420  // then obtain a least squares solution
421 
422  nsing=n;
423 
424  for (j=0;j<n;j++) {
425  double sdiagj=sdiag2[j];
426 
427  if (sdiagj == 0) {
428  nsing=j;
429  break;
430  }
431  }
432 
433  for (j=nsing;j<n;j++) {
434  wa[j]=0.0;
435  }
436 
437  for (k=0;k<nsing;k++) {
438  double sum=0;
439 
440  j=(nsing-1)-k;
441 
442  for (i=j+1;i<nsing;i++) {
443  sum+=r2(i,j)*wa[i];
444  }
445 
446  {
447  double waj=wa[j];
448  double sdiagj=sdiag2[j];
449 
450  wa[j]=(waj-sum)/sdiagj;
451  }
452  }
453 
454  // [GSL] Permute the components of z back to the components of x
455 
456  for (j=0;j<n;j++) {
457  size_t pj=p[j];
458  double waj=wa[j];
459 
460  x[pj]=waj;
461  }
462 
463  return success;
464  }
465 
466  /** \brief Desc
467  */
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) {
471 
472  for (size_t i=0;i<n;i++) {
473  size_t pi=p[i];
474  double dpi=diag2[pi];
475  double xpi=x[pi];
476 
477  w2[i]=dpi*(dpi*xpi)/dxnorm;
478  }
479 
480  for (size_t j=0;j<n;j++) {
481  double sj=sdiag2[j];
482  double wj=w2[j];
483  double tj=wj/sj;
484 
485  w2[j]=tj;
486 
487  for (size_t i=j+1;i<n;i++) {
488  double rij=r2(i,j);
489  double wi=w2[i];
490 
491  w2[i]=wi-rij*tj;
492  }
493  }
494  return;
495  }
496 
497  /** \brief Determine Levenburg-Marquardt parameter
498  */
499  void lmpar(mat_t &r2, const permutation &perm2, const vec_t &qtf2,
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) {
503 
504  double dxnorm, gnorm, fp, fp_old, par_lower, par_upper, par_c;
505 
506  double par2=*par_inout;
507 
508  size_t iter2=0;
509 
510  this->compute_newton_direction(nparm,r2,perm2,qtf2,newton2);
511 
512  // [GSL] Evaluate the function at the origin and test for acceptance of
513  // the Gauss-Newton direction.
514 
515  dxnorm=this->scaled_enorm(diag2,nparm,newton2);
516 
517  fp=dxnorm-delta2;
518 
519  if (fp<=0.1*delta2) {
520 
521  for(size_t i=0;i<nparm;i++) {
522  x[i]=newton2[i];
523  }
524 
525  *par_inout=0;
526 
527  return;
528  }
529 
530  this->compute_newton_bound(ndata,nparm,r2,newton2,dxnorm,perm2,diag2,w2);
531 
532  {
533  double wnorm=o2scl_cblas::dnrm2(ndata,w2);
534  double phider=wnorm*wnorm;
535 
536  // [GSL] w == zero if r rank-deficient,
537  // then set lower bound to zero form MINPACK
538  // Hans E. Plesser 2002-02-25 (hans.plesser@itf.nlh.no)
539  if (wnorm>0) {
540  par_lower=fp/(delta2*phider);
541  } else {
542  par_lower=0.0;
543  }
544  }
545 
546  this->compute_gradient_direction(nparm,r2,perm2,qtf2,diag2,gradient2);
547 
548  gnorm=o2scl_cblas::dnrm2(nparm,gradient2);
549 
550  par_upper=gnorm/delta2;
551 
552  if (par_upper == 0) {
553  double mint;
554  if (delta2<0.1) mint=delta2;
555  else mint=0.1;
556  double dbl_min=std::numeric_limits<double>::min();
557  par_upper=dbl_min/mint;
558  }
559 
560  if (par2>par_upper) {
561  par2=par_upper;
562  } else if (par2<par_lower) {
563  par2=par_lower;
564  }
565 
566  if (par2 == 0) {
567  par2=gnorm/dxnorm;
568  }
569 
570  // [GSL] Beginning of iteration
571 
572  while (true) {
573 
574  iter2++;
575 
576 #ifdef O2SCL_NEVER_DEFINED
577 #ifdef BRIANSFIX
578  // [GSL] Seems like this is described in the paper but
579  // not in the MINPACK code
580 
581  if (par2<par_lower || par2>par_upper) {
582  par2=GSL_MAX_DBL(0.001*par_upper,sqrt(par_lower*par_upper));
583  }
584 #endif
585 #endif
586 
587  // [GSL] Evaluate the function at the current value of par
588 
589  if (par2 == 0) {
590  par2=GSL_MAX_DBL(0.001*par_upper,GSL_DBL_MIN);
591  }
592 
593  // [GSL] Compute the least squares solution of
594  // [ R P x-Q^T f, sqrt(par) D x] for A=Q R P^T
595 
596  {
597  double sqrt_par=sqrt(par2);
598 
599  qrsolv(nparm,r2,perm2,sqrt_par,diag2,qtf2,x,sdiag2,w2);
600  }
601 
602  dxnorm=scaled_enorm(diag2,nparm,x);
603 
604  fp_old=fp;
605 
606  fp=dxnorm-delta2;
607 
608  // [GSL] If the function is small enough, accept the current
609  // value of par
610 
611  if (fabs(fp)<=0.1*delta2) {
612  *par_inout=par2;
613  return;
614  }
615  if (par_lower == 0 && fp<=fp_old && fp_old<0) {
616  *par_inout=par2;
617  return;
618  }
619 
620  // [GSL] Check for maximum number of iterations
621 
622  if (iter2 == 10) {
623  *par_inout=par2;
624  return;
625  }
626 
627  // [GSL] Compute the Newton correction
628 
629  this->compute_newton_correction(nparm,r2,sdiag2,perm2,x,
630  dxnorm,diag2,w2);
631  {
632  double wnorm=o2scl_cblas::dnrm2(ndata,w2);
633  par_c=fp/(delta2*wnorm*wnorm);
634  }
635 
636  // [GSL] Depending on the sign of the function,
637  // update par_lower or par_upper
638 
639  if (fp>0) {
640  if (par2>par_lower) {
641  par_lower=par2;
642  }
643  } else if (fp<0) {
644  if (par2<par_upper) {
645  par_upper=par2;
646  }
647  }
648 
649  // [GSL] Compute an improved estimate for par
650 
651  par2=GSL_MAX_DBL(par_lower,par2+par_c);
652 
653  }
654 
655  O2SCL_ERR("Sanity check failed in fit_nonlin_b::lmpar().",
656  exc_esanity);
657  return;
658  }
659 
660  /// Compute trial step, \f$ \mathrm{trial}=\mathrm{x}+\mathrm{dx} \f$
661  void compute_trial_step(size_t N, vec_t &x, vec_t &dx,
662  vec_t &trial) {
663  for (size_t i=0;i<N;i++) {
664  trial[i]=x[i]+dx[i];
665  }
666  return;
667  }
668 
669  /** \brief Compute the root of the sum of the squares of
670  the columns of \c J
671 
672  This computes
673  \f[
674  \mathrm{diag\_vec}_j = \sqrt{\sum_{i=0}^{\mathrm{ndata}-1} J_{ij}}
675  \f]
676  for \f$ 0\leq j \leq \mathrm{nparm}-1 \f$ . If any of the
677  columns of \c J is all zero, then the corresponding entry in
678  \c diag_vec is set to one instead.
679  */
680  int compute_diag(size_t nparm, size_t ndata,
681  const mat_t &J, vec_t &diag_vec) {
682 
683  for (size_t j=0;j<nparm;j++) {
684  double sum=0.0;
685 
686  for (size_t i=0;i<ndata;i++) {
687  sum+=J(i,j)*J(i,j);
688  }
689  if (sum == 0) {
690  sum=1.0;
691  }
692 
693  diag_vec[j]=sqrt(sum);
694  }
695  return 0;
696  }
697 
698  /** \brief Compute the covarance matrix \c covar given
699  the Jacobian \c J
700 
701  Given a \c m by \c n Jacobian matrix \c J (where \c m must not
702  be less than \c n), and a relative tolerance \c epsrel, this
703  function computes the entries of the \c n by \c n covariance
704  matrix \c covar. The allocation for \c covar must be performed
705  beforehand.
706 
707  This function is basically the equivalent of the function
708  <tt>gsl_multifit_covar()</tt>, but rewritten for generic
709  vector and matrix types.
710 
711  The workspace \c work1 is used here.
712 
713  \comment
714  Note that in the GSL example for non-linear fitting, the value
715  of \c epsrel is set to zero, so this class sets \ref
716  tol_rel_covar to zero by default. We could remove the \c epsrel
717  parameter in this function, but it may be one day useful to
718  separate this function so we leave \c epsrel as a parameter.
719  \endcomment
720  */
721  int covariance(size_t m, size_t n, const mat_t &J,
722  mat_t &covar, vec_t &norm, mat_t &r,
723  vec_t &tau, permutation &perm, double epsrel) {
724 
725  if (m<n) {
726  O2SCL_ERR2("Jacobian must have m>=n in ",
727  "fit_nonlin_b::covariance().",exc_efailed);
728  }
729 
730  double tolr;
731 
732  size_t kmax=0;
733 
734  int signum=0;
735 
736  for(size_t i=0;i<m;i++) {
737  for(size_t j=0;j<n;j++) {
738  r(i,j)=J(i,j);
739  }
740  }
741  o2scl_linalg::QRPT_decomp(m,n,r,tau,perm,signum,norm);
742 
743  // [GSL] Form the inverse of R in the full upper triangle of R
744 
745  tolr=epsrel*fabs(r(0,0));
746 
747  for (size_t k=0;k<n;k++) {
748  double rkk=r(k,k);
749 
750  if (fabs(rkk)<=tolr) {
751  break;
752  }
753 
754  r(k,k)=1.0/rkk;
755 
756  for (size_t j=0;j<k;j++) {
757  double t=r(j,k)/rkk;
758  r(j,k)=0.0;
759 
760  for (size_t i=0;i<=j;i++) {
761  double rik=r(i,k);
762  double rij=r(i,j);
763  r(i,k)=rik-t*rij;
764  }
765  }
766  kmax=k;
767  }
768 
769  // [GSL] Form the full upper triangle of the inverse of R^T R in
770  // the full upper triangle of R
771 
772  for (size_t k=0;k<=kmax;k++) {
773  for (size_t j=0;j<k;j++) {
774  double rjk=r(j,k);
775 
776  for (size_t i=0;i<=j;i++) {
777  double rij=r(i,j);
778  double rik=r(i,k);
779  r(i,j)=rij+rjk*rik;
780  }
781  }
782 
783  double t=r(k,k);
784 
785  for (size_t i=0;i<=k;i++) {
786  r(i,k)*=t;
787  }
788  }
789 
790  // [GSL] Form the full lower triangle of the covariance matrix in
791  // the strict lower triangle of R and in w
792 
793  for (size_t j=0;j<n;j++) {
794  size_t pj=perm[j];
795 
796  for (size_t i=0;i<=j;i++) {
797  size_t pi=perm[i];
798 
799  double rij;
800 
801  if (j>kmax) {
802  r(i,j)=0.0;
803  rij=0.0;
804  } else {
805  rij=r(i,j);
806  }
807 
808  if (pi>pj) {
809  r(pi,pj)=rij;
810  } else if (pi<pj) {
811  r(pj,pi)=rij;
812  }
813 
814  }
815 
816  double rjj=r(j,j);
817  covar(pj,pj)=rjj;
818  }
819 
820 
821  // [GSL] symmetrize the covariance matrix
822 
823  for (size_t j=0;j<n;j++) {
824  for (size_t i=0;i<j;i++) {
825  double rji=r(j,i);
826  covar(j,i)=rji;
827  covar(i,j)=rji;
828  }
829  }
830 
831  return success;
832  }
833 
834  public:
835 
836  fit_nonlin_b() {
837  tol_rel_covar=0.0;
838  }
839 
840  /** \brief The relative tolerance for the computation of
841  the covariance matrix (default 0)
842  */
844 
845  /// Test if converged
846  int test_delta_f(size_t nparm, vec_t &dx, vec_t &x,
847  double l_epsabs, double l_epsrel) {
848 
849  if (l_epsrel<0.0) {
850  O2SCL_ERR2("Relative tolerance less than zero ",
851  "in fit_nonlin_b::test_delta_f().",exc_einval);
852  }
853 
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) {
857  return gsl_continue;
858  }
859  }
860  return success;
861  }
862 
863  /// Test if converged
864  int test_gradient_f(size_t nparm, vec_t &g, double l_epsabs) {
865 
866  double residual=0.0;
867 
868  if (l_epsabs<0.0) {
869  O2SCL_ERR2("Absolute tolerance less than zero ",
870  "in fit_nonlin_b::test_gradient_f().",exc_einval);
871  }
872 
873  for(size_t i=0;i<nparm;i++) {
874  residual+=fabs(g[i]);
875  }
876  if (residual<l_epsabs) {
877  return success;
878  }
879  return gsl_continue;
880  }
881 
882  };
883 
884  /** \brief Non-linear least-squares fitting class (GSL)
885 
886  The GSL-based fitting class using a Levenberg-Marquardt type
887  algorithm. The algorithm stops when
888  \f[
889  |dx_i| < \mathrm{tol\_abs}+\mathrm{tol\_rel}\times|x_i|
890  \f]
891  where \f$dx\f$ is the last step and \f$x\f$ is the current
892  position. If test_gradient is true, then additionally fit()
893  requires that
894  \f[
895  \sum_i |g_i| < \mathrm{tol\_abs}
896  \f]
897  where \f$g_i\f$ is the \f$i\f$-th component of the gradient of
898  the function \f$\Phi(x)\f$ where
899  \f[
900  \Phi(x) = || F(x) ||^2
901  \f]
902 
903  Default template arguments
904  - \c func_t - \ref gen_fit_funct<>
905  - \c vec_t - \ref boost::numeric::ublas::vector <double >
906  - \c mat_t - \ref boost::numeric::ublas::matrix <double >
907 
908  \todo Allow the user to specify the derivatives
909  \todo Fix so that the user can specify automatic
910  scaling of the fitting parameters, where the initial
911  guess are used for scaling so that the fitting parameters
912  are near unity.
913 
914  \future Some of these member functions (like
915  update_diag()) don't depend on member data and could be
916  possibly be moved to a parent class?
917  */
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> >
921  class fit_nonlin : public fit_nonlin_b<vec_t,mat_t>,
922  public fit_base<func_t,vec_t,mat_t> {
923 
924  protected:
925 
926  /// Function to fit
927  func_t *cff;
928 
929  /// Trial step
930  vec_t x_trial;
931 
932  /// Trial function value
933  vec_t f_trial;
934 
935  /// Number of iterations
936  size_t iter;
937 
938  /** \brief Desc
939  */
940  double xnorm;
941 
942  /** \brief Desc */
943  double fnorm;
944 
945  /** \brief Desc */
946  double delta;
947 
948  /** \brief Desc */
949  double par;
950 
951  /** \brief Desc */
952  mat_t r;
953 
954  /** \brief Desc */
955  vec_t tau;
956 
957  /** \brief Desc */
958  vec_t diag;
959 
960  /** \brief Desc */
961  vec_t qtf;
962 
963  /** \brief Desc */
964  vec_t df;
965 
966  /** \brief Desc */
967  vec_t rptdx;
968 
969  /** \brief Desc */
970  vec_t newton;
971 
972  /** \brief Desc */
973  vec_t gradient;
974 
975  /** \brief Desc */
976  vec_t sdiag;
977 
978  /** \brief Desc */
979  vec_t w;
980 
981  /** \brief Desc */
982  vec_t work1;
983 
984  /** \brief Desc */
986 
987  /// Number of data points
988  size_t ndata;
989 
990  /// Number of parameters
991  size_t nparm;
992 
993  /// Desc
994  vec_t g_;
995 
996  /// Desc
997  mat_t J_;
998 
999  /// Desc
1000  vec_t *x_;
1001 
1002  /// Free allocated memory
1003  void free() {
1004  if (ndata>0) {
1005  perm.free();
1006  f_trial.clear();
1007  x_trial.clear();
1008  f_.clear();
1009  dx_.clear();
1010  g_.clear();
1011  J_.clear();
1012  r.clear();
1013  rptdx.clear();
1014  w.clear();
1015  work1.clear();
1016  newton.clear();
1017  gradient.clear();
1018  sdiag.clear();
1019  qtf.clear();
1020  df.clear();
1021  diag.clear();
1022  tau.clear();
1023  ndata=0;
1024  nparm=0;
1025  }
1026  }
1027 
1028  public:
1029 
1030  fit_nonlin() {
1031  use_scaled=true;
1032  test_gradient=false;
1033  ndata=0;
1034  nparm=0;
1035  err_nonconv=true;
1036  }
1037 
1038  virtual ~fit_nonlin() {
1039  free();
1040  }
1041 
1042  /** \brief If true, call the error handler if fit()
1043  does not converge (default true)
1044  */
1046 
1047  /// Print the progress in the current iteration
1048  virtual int print_iter(size_t nv, vec_t &x, vec_t &dx, int iter2,
1049  double l_epsabs, double l_epsrel) {
1050 
1051  if (this->verbose<=0) return 0;
1052 
1053  // Variables 'val' and 'lim' set to zero to avoid
1054  // uninitialized variable errors. In reality, they're
1055  // always set in the loop below.
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) {
1060  val=fabs(dx[i]);
1061  lim=tolerance;
1062  i=nv;
1063  } else {
1064  double diff=tolerance-fabs(dx[i]);
1065  if (diff>diff_max) {
1066  diff_max=diff;
1067  val=fabs(dx[i]);
1068  lim=tolerance;
1069  }
1070  }
1071  }
1072 
1073  std::cout << "Iteration: " << iter2 << std::endl;
1074  std::cout << "x: ";
1075  for(size_t i=0;i<nv;i++) {
1076  std::cout << x[i] << " ";
1077  }
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. ";
1082  char ch;
1083  std::cin >> ch;
1084  }
1085 
1086  return 0;
1087  }
1088 
1089  /// Allocate memory with \c n data points and \c p parameters
1090  void resize(size_t n, size_t p) {
1091 
1092  free();
1093 
1094  ndata=n;
1095  nparm=p;
1096 
1097  perm.allocate(nparm);
1098  r.resize(ndata,nparm);
1099  if (ndata<nparm) {
1100  tau.resize(ndata);
1101  for(size_t i=0;i<ndata;i++) tau[i]=0.0;
1102  } else {
1103  tau.resize(nparm);
1104  for(size_t i=0;i<nparm;i++) tau[i]=0.0;
1105  }
1106  df.resize(ndata);
1107  rptdx.resize(ndata);
1108  qtf.resize(ndata);
1109  w.resize(ndata);
1110  diag.resize(nparm);
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++) {
1117  r(i,j)=0.0;
1118  }
1119  rptdx[i]=0.0;
1120  w[i]=0.0;
1121  qtf[i]=0.0;
1122  df[i]=0.0;
1123  }
1124  for(size_t i=0;i<nparm;i++) {
1125  work1[i]=0.0;
1126  newton[i]=0.0;
1127  gradient[i]=0.0;
1128  sdiag[i]=0.0;
1129  diag[i]=0.0;
1130  }
1131 
1132  x_trial.resize(nparm);
1133  f_trial.resize(ndata);
1134  f_.resize(ndata);
1135  dx_.resize(nparm);
1136  g_.resize(ndata);
1137  J_.resize(ndata,nparm);
1138 
1139  return;
1140  }
1141 
1142  /// The last step taken in parameter space
1143  vec_t dx_;
1144 
1145  /// Desc
1146  vec_t f_;
1147 
1148  /** \brief Set the initial values of the parameters and
1149  the fitting function to use for the next call to
1150  iterate()
1151  */
1152  int set(size_t npar, vec_t &parms, func_t &fitfun) {
1153 
1154  cff=&fitfun;
1155 
1156  if (fitfun.get_ndata()==0 || npar==0) {
1157  O2SCL_ERR2("Either zero data or zero parameters in ",
1158  "fit_nonlin::fit().",exc_einval);
1159  }
1160 
1161  if (ndata!=fitfun.get_ndata() || nparm!=npar) {
1162  resize(fitfun.get_ndata(),npar);
1163  }
1164 
1165  x_=&parms;
1166 
1167  int signum;
1168 
1169  // Evaluate function and Jacobian at x
1170  (*cff)(nparm,*x_,ndata,f_);
1171  cff->jac(nparm,*x_,ndata,f_,J_);
1172 
1173  par=0;
1174  iter=1;
1175  fnorm=o2scl_cblas::dnrm2(ndata,f_);
1176 
1177  for(size_t i=0;i<nparm;i++) {
1178  dx_[i]=0.0;
1179  }
1180 
1181  // [GSL] store column norms in diag
1182 
1183  if (use_scaled) {
1184  this->compute_diag(nparm,ndata,J_,diag);
1185  } else {
1186  for(size_t i=0;i<nparm;i++) {
1187  diag[i]=1.0;
1188  }
1189  }
1190 
1191  // [GSL] set delta to 100 |D x| or to 100 if |D x| is zero
1192 
1193  xnorm=this->scaled_enorm(diag,nparm,*x_);
1194  delta=this->compute_delta(diag,nparm,*x_);
1195 
1196  // [GSL] Factorize J into QR decomposition
1197 
1198  for(size_t i=0;i<ndata;i++) {
1199  for(size_t j=0;j<nparm;j++) {
1200  r(i,j)=J_(i,j);
1201  }
1202  }
1203 
1204  o2scl_linalg::QRPT_decomp(ndata,nparm,r,tau,perm,signum,work1);
1205 
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;
1208 
1209  // [GSL] Zero the trial vector, as in the alloc function
1210 
1211  for(size_t i=0;i<ndata;i++) f_trial[i]=0.0;
1212 
1213  return success;
1214  }
1215 
1216  /** \brief Perform an iteration
1217  */
1218  int iterate() {
1219 
1220  vec_t &f=f_;
1221  mat_t &J=J_;
1222  vec_t &dx=dx_;
1223  vec_t &x=*x_;
1224 
1225  double prered, actred;
1226  double pnorm, fnorm1, fnorm1p, gnorm;
1227  double ratio;
1228  double dirder;
1229 
1230  int niter=0;
1231 
1232  double p1=0.1, p25=0.25, p5=0.5, p75=0.75, p0001=0.0001;
1233 
1234  if (fnorm == 0.0) {
1235  return success;
1236  }
1237 
1238  // [GSL] Compute qtf=Q^T f
1239 
1240  for(size_t i=0;i<ndata;i++) {
1241  qtf[i]=f[i];
1242  }
1243  o2scl_linalg::QR_QTvec(ndata,nparm,r,tau,qtf);
1244 
1245  // [GSL] Compute norm of scaled gradient
1246 
1247  this->compute_gradient_direction(nparm,r,perm,qtf,diag,gradient);
1248  size_t iamax=vector_max_index<vec_t,double>(nparm,gradient);
1249  gnorm=fabs(gradient[iamax]/fnorm);
1250 
1251  // [GSL] Determine the Levenberg-Marquardt parameter
1252 
1253  bool lm_iteration=true;
1254  while (lm_iteration) {
1255 
1256  niter++;
1257 
1258  this->lmpar(r,perm,qtf,diag,delta,&par,newton,
1259  gradient,sdiag,dx,w,nparm,ndata);
1260 
1261  // [GSL] Take a trial step
1262 
1263  // [GSL] reverse the step to go downhill
1264  for(size_t i=0;i<nparm;i++) dx[i]*=-1.0;
1265 
1266  this->compute_trial_step(nparm,x,dx,x_trial);
1267 
1268  pnorm=this->scaled_enorm(diag,nparm,dx);
1269 
1270  if (iter == 1) {
1271  if (pnorm<delta) {
1272  delta=pnorm;
1273  }
1274  }
1275 
1276  // [GSL] Evaluate function at x+p
1277  (*cff)(nparm,x_trial,ndata,f_trial);
1278 
1279  fnorm1=o2scl_cblas::dnrm2(ndata,f_trial);
1280 
1281  // [GSL] Compute the scaled actual reduction
1282 
1283  actred=this->compute_actual_reduction(fnorm,fnorm1);
1284 
1285  // [GSL] Compute rptdx=R P^T dx,noting that |J dx|=|R P^T dx|
1286 
1287  this->compute_rptdx(r,perm,nparm,dx,rptdx);
1288 
1289  fnorm1p=o2scl_cblas::dnrm2(ndata,rptdx);
1290 
1291  // [GSL] Compute the scaled predicted reduction,
1292  // |J dx|^2+2 par |D dx|^2
1293 
1294  double t1=fnorm1p/fnorm;
1295  double t2=(sqrt(par)*pnorm)/fnorm;
1296 
1297  prered=t1*t1+t2*t2/p5;
1298  dirder=-(t1*t1+t2*t2);
1299 
1300  // [GSL] compute the ratio of the actual to predicted reduction
1301 
1302  if (prered>0) {
1303  ratio=actred/prered;
1304  } else {
1305  ratio=0;
1306  }
1307 
1308  // [GSL] update the step bound
1309 
1310  if (ratio>p25) {
1311  if (par == 0 || ratio >= p75) {
1312  delta=pnorm/p5;
1313  par *= p5;
1314  }
1315  } else {
1316  double temp=(actred >= 0) ? p5 : p5*dirder/
1317  (dirder+p5*actred);
1318 
1319  if (p1*fnorm1 >= fnorm || temp<p1 ) {
1320  temp=p1;
1321  }
1322 
1323  if (delta<pnorm/p1) delta=temp*delta;
1324  else delta=temp*pnorm/p1;
1325 
1326  par /= temp;
1327  }
1328 
1329 
1330  // Test for successful iteration, termination and
1331  // stringent tolerances.
1332 
1333  lm_iteration=false;
1334 
1335  double dbl_eps=std::numeric_limits<double>::epsilon();
1336 
1337  if (ratio >= p0001) {
1338 
1339  vector_copy(nparm,x_trial,x);
1340  vector_copy(ndata,f_trial,f);
1341 
1342  cff->jac(nparm,x_trial,ndata,f_trial,J_);
1343 
1344  // [GSL] wa2_j =diag_j*x_j
1345  xnorm=this->scaled_enorm(diag,nparm,x);
1346  fnorm=fnorm1;
1347  iter++;
1348 
1349  // [GSL] Rescale if necessary
1350 
1351  if (use_scaled) {
1352  this->update_diag(nparm,J,diag);
1353  }
1354 
1355  {
1356  int signum;
1357  for(size_t i=0;i<ndata;i++) {
1358  for(size_t j=0;j<nparm;j++) {
1359  r(i,j)=J(i,j);
1360  }
1361  }
1362  o2scl_linalg::QRPT_decomp(ndata,nparm,r,tau,perm,signum,work1);
1363  }
1364 
1365  return success;
1366  } else if (fabs(actred)<=dbl_eps &&
1367  prered<=dbl_eps && p5*ratio<=1.0) {
1368  return exc_etolf;
1369  } else if (delta<=dbl_eps*xnorm) {
1370  return exc_etolx;
1371  } else if (gnorm<=dbl_eps) {
1372  return exc_etolg;
1373  } else if (niter<10) {
1374  lm_iteration=true;
1375  }
1376 
1377  }
1378 
1379  return gsl_continue;
1380  }
1381 
1382  /** \brief Fit the data specified in (xdat,ydat) to
1383  the function \c fitfun with the parameters in \c par.
1384 
1385  The covariance matrix for the parameters is returned in \c covar
1386  and the value of \f$ \chi^2 \f$ is returned in \c chi2.
1387  */
1388  virtual int fit(size_t npar, vec_t &parms, mat_t &covar,
1389  double &chi2, func_t &fitfun) {
1390 
1391  set(npar,parms,fitfun);
1392 
1393  int status;
1394  size_t niter=0;
1395  do {
1396 
1397  niter++;
1398  status=iterate();
1399 
1400  if (status) {
1401  break;
1402  }
1403 
1404  status=this->test_delta_f(nparm,dx_,parms,this->tol_abs,
1405  this->tol_rel);
1406  if (test_gradient && status==gsl_continue) {
1407  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
1408  o2scl_cblas::o2cblas_NoTrans,
1409  fitfun.get_ndata(),npar,1.0,J_,f_,0.0,g_);
1410  status=this->test_gradient_f(nparm,g_,this->tol_abs);
1411  }
1412 
1413  this->print_iter(npar,parms,dx_,niter,this->tol_abs,this->tol_rel);
1414 
1415  } while (status == gsl_continue && niter<this->ntrial);
1416 
1417  if (niter>=this->ntrial) {
1418  O2SCL_CONV2_RET("Function fit_nonlin::fit() ",
1419  "exceeded max. number of iterations.",
1420  exc_emaxiter,err_nonconv);
1421  }
1422 
1423  // Compute covariance matrix
1424  this->covariance(fitfun.get_ndata(),npar,J_,covar,work1,
1425  r,tau,perm,this->tol_rel_covar);
1426 
1427  // Compute chi squared
1428  chi2=o2scl_cblas::dnrm2(fitfun.get_ndata(),f_);
1429  chi2*=chi2;
1430 
1431  return 0;
1432  }
1433 
1434  /// If true, test the gradient also (default false)
1436 
1437  /** \brief Use the scaled routine if true (default true)
1438  */
1440 
1441  /// Return string denoting type ("fit_nonlin")
1442  virtual const char *type() { return "fit_nonlin"; }
1443 
1444  };
1445 
1446 #ifndef DOXYGEN_NO_O2NS
1447 }
1448 #endif
1449 
1450 #endif
bool test_gradient
If true, test the gradient also (default false)
Definition: fit_nonlin.h:1435
double compute_delta(vec_t &diag2, size_t n, const vec_t &x)
Desc.
Definition: fit_nonlin.h:243
int test_delta_f(size_t nparm, vec_t &dx, vec_t &x, double l_epsabs, double l_epsrel)
Test if converged.
Definition: fit_nonlin.h:846
vec_t newton
Desc.
Definition: fit_nonlin.h:970
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
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.
Definition: fit_nonlin.h:499
vec_t diag
Desc.
Definition: fit_nonlin.h:958
bool err_nonconv
If true, call the error handler if fit() does not converge (default true)
Definition: fit_nonlin.h:1045
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.
Definition: fit_nonlin.h:182
const double pi
Definition: constants.h:62
mat_t r
Desc.
Definition: fit_nonlin.h:952
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.
Definition: fit_nonlin.h:680
vec_t df
Desc.
Definition: fit_nonlin.h:964
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
Non-linear least-squares fitting class (GSL)
Definition: fit_nonlin.h:921
invalid argument supplied by user
Definition: err_hnd.h:59
A class for representing permutations.
Definition: permutation.h:70
exceeded max number of iterations
Definition: err_hnd.h:73
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.
Definition: fit_nonlin.h:1048
vec_t dx_
The last step taken in parameter space.
Definition: fit_nonlin.h:1143
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
size_t nparm
Number of parameters.
Definition: fit_nonlin.h:991
void compute_trial_step(size_t N, vec_t &x, vec_t &dx, vec_t &trial)
Compute trial step, .
Definition: fit_nonlin.h:661
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.
Definition: qrpt_base.h:54
iteration has not converged
Definition: err_hnd.h:51
int free()
Free the memory.
Definition: permutation.h:233
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.
Definition: fit_nonlin.h:469
int apply_inverse(vec_t &v) const
Apply the inverse permutation to a vector.
Definition: permutation.h:315
double par
Desc.
Definition: fit_nonlin.h:949
size_t ndata
Number of data points.
Definition: fit_nonlin.h:988
vec_t w
Desc.
Definition: fit_nonlin.h:979
generic failure
Definition: err_hnd.h:61
vec_t * x_
Desc.
Definition: fit_nonlin.h:1000
vec_t work1
Desc.
Definition: fit_nonlin.h:982
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.
Definition: fit_nonlin.h:310
Base routines for the nonlinear fitting classes.
Definition: fit_nonlin.h:67
vec_t tau
Desc.
Definition: fit_nonlin.h:955
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
Definition: cblas_base.h:156
permutation perm
Desc.
Definition: fit_nonlin.h:985
vec_t g_
Desc.
Definition: fit_nonlin.h:994
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 .
Definition: cblas_base.h:218
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
vec_t f_trial
Trial function value.
Definition: fit_nonlin.h:933
size_t iter
Number of iterations.
Definition: fit_nonlin.h:936
double xnorm
Desc.
Definition: fit_nonlin.h:940
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
vec_t rptdx
Desc.
Definition: fit_nonlin.h:967
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...
Definition: fit_nonlin.h:1388
vec_t gradient
Desc.
Definition: fit_nonlin.h:973
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
double delta
Desc.
Definition: fit_nonlin.h:946
void compute_rptdx(const mat_t &r2, const permutation &p, size_t N, vec_t &dx, vec_t &rptdx2)
Desc.
Definition: fit_nonlin.h:252
mat_t J_
Desc.
Definition: fit_nonlin.h:997
cannot reach the specified tolerance in x
Definition: err_hnd.h:111
cannot reach the specified tolerance in gradient gradient
Definition: err_hnd.h:113
size_t count_nsing(const size_t ncols, const mat_t &r2)
Desc.
Definition: fit_nonlin.h:88
void compute_newton_direction(size_t n, const mat_t &r2, const permutation &perm2, const vec_t &qtf2, vec_t &x)
Desc.
Definition: fit_nonlin.h:104
void resize(size_t n, size_t p)
Allocate memory with n data points and p parameters.
Definition: fit_nonlin.h:1090
void free()
Free allocated memory.
Definition: fit_nonlin.h:1003
double compute_actual_reduction(double fnorm0, double fnorm1)
Desc.
Definition: fit_nonlin.h:73
double fnorm
Desc.
Definition: fit_nonlin.h:943
int iterate()
Perform an iteration.
Definition: fit_nonlin.h:1218
vec_t sdiag
Desc.
Definition: fit_nonlin.h:976
Non-linear least-squares fitting [abstract base].
Definition: fit_base.h:327
bool use_scaled
Use the scaled routine if true (default true)
Definition: fit_nonlin.h:1439
void update_diag(size_t n, const mat_t &J, vec_t &diag2)
Desc.
Definition: fit_nonlin.h:206
vec_t f_
Desc.
Definition: fit_nonlin.h:1146
vec_t x_trial
Trial step.
Definition: fit_nonlin.h:930
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.
Definition: fit_nonlin.h:142
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.
Definition: qr_base.h:71
int allocate(size_t dim)
Allocate memory for a permutation of size dim.
Definition: permutation.h:219
vec_t qtf
Desc.
Definition: fit_nonlin.h:961
Success.
Definition: err_hnd.h:47
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.
Definition: fit_nonlin.h:721
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.
Definition: fit_nonlin.h:231
virtual const char * type()
Return string denoting type ("fit_nonlin")
Definition: fit_nonlin.h:1442
int test_gradient_f(size_t nparm, vec_t &g, double l_epsabs)
Test if converged.
Definition: fit_nonlin.h:864
func_t * cff
Function to fit.
Definition: fit_nonlin.h:927
double tol_rel_covar
The relative tolerance for the computation of the covariance matrix (default 0)
Definition: fit_nonlin.h:843
cannot reach the specified tolerance in f
Definition: err_hnd.h:109

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).