fit_min.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 #ifndef O2SCL_FIT_MIN_H
24 #define O2SCL_FIT_MIN_H
25 
26 /** \file fit_min.h
27  \brief File defining \ref o2scl::fit_min
28 */
29 
30 #include <o2scl/mmin.h>
31 #include <o2scl/multi_funct.h>
32 #include <o2scl/mmin_simp2.h>
33 #include <o2scl/fit_base.h>
34 #include <o2scl/fit_nonlin.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Non-linear least-squares fitting class with generic minimizer
41 
42  This minimizes a generic fitting function using any \ref
43  o2scl::mmin_base object, and then uses the GSL routines to
44  calculate the uncertainties in the parameters and the covariance
45  matrix.
46 
47  This can be useful for fitting problems which might be better
48  handled by more complex minimizers than those that are used in
49  \ref o2scl::fit_nonlin. For problems with many local minima near
50  the global minimum, using a \ref o2scl::anneal_base object with
51  this class can sometimes produce better results than \ref
52  o2scl::fit_nonlin.
53 
54  Default template arguments
55  - \c func_t - \ref gen_fit_funct<>
56  - \c vec_t - \ref boost::numeric::ublas::vector <double >
57  - \c mat_t - \ref boost::numeric::ublas::matrix <double >
58  */
59  template<class func_t=gen_fit_funct<>,
60  class vec_t=boost::numeric::ublas::vector<double>,
61  class mat_t=boost::numeric::ublas::matrix<double> > class fit_min :
62  public fit_base<func_t,vec_t,mat_t>, public fit_nonlin_b<vec_t,mat_t> {
63 
64  public:
65 
66  fit_min() {
67  min_set=false;
68  mmp=&def_mmin;
69  }
70 
71  virtual ~fit_min() {}
72 
73  /** \brief Fit the data specified in (xdat,ydat) to
74  the function \c fitfun with the parameters in \c par.
75 
76  The covariance matrix for the parameters is returned in \c
77  covar and the value of \f$ \chi^2 \f$ is returned in \c chi2.
78  */
79  virtual int fit(size_t npar, vec_t &par, mat_t &covar, double &chi2,
80  func_t &fitfun) {
81 
82  func=&fitfun;
83 
84  size_t ndata=fitfun.get_ndata();
85  fval.resize(ndata);
86 
87  // ---------------------------------------------------
88  // First minimize with the specified mmin object
89  // ---------------------------------------------------
90 
91  multi_funct11 mmf=std::bind(std::mem_fn<double(size_t, const vec_t &)>
93  this,std::placeholders::_1,
94  std::placeholders::_2);
95  //multi_funct_mfptr<fit_min>
96  //mmf(this,&fit_min::min_func);
97 
98  double dtemp;
99  mmp->mmin(npar,par,dtemp,mmf);
100 
101  // ---------------------------------------------------
102  // Now compute the Jacobian and do the QR decomposition
103  // ---------------------------------------------------
104 
105  // Allocate space
106  int signum;
107  permutation perm(npar);
108  vec_t diag, tau, norm;
109  diag.resize(npar);
110  if (ndata<npar) {
111  tau.resize(ndata);
112  } else {
113  tau.resize(npar);
114  }
115  norm.resize(npar);
116  mat_t J, r;
117  J.resize(ndata,npar);
118  r.resize(ndata,npar);
119 
120  // Evaluate function and Jacobian at optimal parameter values
121  fitfun(npar,par,ndata,fval);
122  fitfun.jac(npar,par,ndata,fval,J);
123 
124  double fnorm=o2scl_cblas::dnrm2(ndata,fval);
125 
126  // Use scaled version
127  this->compute_diag(npar,ndata,J,diag);
128 
129  double xnorm=this->scaled_enorm(diag,npar,fval);
130  double delta=this->compute_delta(diag,npar,fval);
131  matrix_copy(ndata,npar,J,r);
132 
133  o2scl_linalg::QRPT_decomp(ndata,npar,r,tau,perm,signum,norm);
134 
135  // ---------------------------------------------------
136  // Compute the covariance matrix
137  // ---------------------------------------------------
138 
139  this->covariance(ndata,npar,J,covar,norm,r,tau,perm,
140  this->tol_rel_covar);
141 
142  chi2=fnorm*fnorm;
143 
144  // Free previously allocated space
145 
146  diag.clear();
147  tau.clear();
148  norm.clear();
149  J.clear();
150  r.clear();
151  fval.clear();
152 
153  return 0;
154  }
155 
156  /** \brief Set the mmin object to use (default is
157  of type \ref o2scl::mmin_simp2)
158  */
160  mmp=&mm;
161  min_set=true;
162  return 0;
163  }
164 
165  /// The default minimizer
167 
168  /// Return string denoting type ("fit_min")
169  virtual const char *type() { return "fit_min"; }
170 
171 #ifndef DOXYGEN_INTERNAL
172 
173  protected:
174 
175  /// Storage for the function values
176  vec_t fval;
177 
178  /// Pointer to the user-specified fitting function
179  func_t *func;
180 
181  /// The minimizer
183 
184  /// True if the minimizer has been set by the user
185  bool min_set;
186 
187  /// The function to minimize, \f$ \chi^2 \f$
188  double min_func(size_t np, const vec_t &xp) {
189  double ret=0.0;
190  (*func)(np,xp,func->get_ndata(),fval);
191  for(size_t i=0;i<func->get_ndata();i++) {
192  ret+=fval[i]*fval[i];
193  }
194  return ret;
195  }
196 
197 #endif
198 
199  };
200 
201 #ifndef DOXYGEN_NO_O2NS
202 }
203 #endif
204 
205 #endif
double compute_delta(vec_t &diag2, size_t n, const vec_t &x)
Desc.
Definition: fit_nonlin.h:243
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
func_t * func
Pointer to the user-specified fitting function.
Definition: fit_min.h:179
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
virtual int mmin(size_t nvar, vec_t &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t. the array x of size nvar.
Non-linear least-squares fitting class with generic minimizer.
Definition: fit_min.h:61
A class for representing permutations.
Definition: permutation.h:70
mmin_base< multi_funct11 > * mmp
The minimizer.
Definition: fit_min.h:182
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
vec_t fval
Storage for the function values.
Definition: fit_min.h:176
Base routines for the nonlinear fitting classes.
Definition: fit_nonlin.h:67
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
Definition: cblas_base.h:156
virtual const char * type()
Return string denoting type ("fit_min")
Definition: fit_min.h:169
void matrix_copy(mat_t &src, mat2_t &dest)
Simple matrix copy.
Definition: vector.h:179
int set_mmin(mmin_base< multi_funct11 > &mm)
Set the mmin object to use (default is of type o2scl::mmin_simp2)
Definition: fit_min.h:159
bool min_set
True if the minimizer has been set by the user.
Definition: fit_min.h:185
double min_func(size_t np, const vec_t &xp)
The function to minimize, .
Definition: fit_min.h:188
virtual int fit(size_t npar, vec_t &par, 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_min.h:79
Non-linear least-squares fitting [abstract base].
Definition: fit_base.h:327
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
mmin_simp2< multi_funct11 > def_mmin
The default minimizer.
Definition: fit_min.h:166
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
double tol_rel_covar
The relative tolerance for the computation of the covariance matrix (default 0)
Definition: fit_nonlin.h:843

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