fit_min.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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_funct 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
o2scl::mmin_simp2< multi_funct >
o2scl::fit_min::func
func_t * func
Pointer to the user-specified fitting function.
Definition: fit_min.h:179
o2scl::permutation
A class for representing permutations.
Definition: permutation.h:70
o2scl::matrix_copy
void matrix_copy(mat_t &src, mat2_t &dest)
Simple matrix copy.
Definition: vector.h:176
o2scl::fit_min::fval
vec_t fval
Storage for the function values.
Definition: fit_min.h:176
o2scl::fit_min::min_set
bool min_set
True if the minimizer has been set by the user.
Definition: fit_min.h:185
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl_linalg::QRPT_decomp
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
o2scl::fit_min::set_mmin
int set_mmin(mmin_base< multi_funct > &mm)
Set the mmin object to use (default is of type o2scl::mmin_simp2)
Definition: fit_min.h:159
o2scl::multi_funct
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef in src/base/multi_funct.h.
Definition: multi_funct.h:45
o2scl::fit_nonlin_b
Base routines for the nonlinear fitting classes.
Definition: fit_nonlin.h:67
o2scl::fit_min::fit
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
o2scl::fit_nonlin_b< boost::numeric::ublas::vector< double >, boost::numeric::ublas::matrix< double > >::tol_rel_covar
double tol_rel_covar
The relative tolerance for the computation of the covariance matrix (default 0)
Definition: fit_nonlin.h:843
o2scl::fit_base
Non-linear least-squares fitting [abstract base].
Definition: fit_base.h:329
o2scl::fit_min
Non-linear least-squares fitting class with generic minimizer.
Definition: fit_min.h:61
o2scl::mmin_base< multi_funct >
o2scl::fit_min::type
virtual const char * type()
Return string denoting type ("fit_min")
Definition: fit_min.h:169
o2scl::fit_min::def_mmin
mmin_simp2< multi_funct > def_mmin
The default minimizer.
Definition: fit_min.h:166
o2scl::fit_min::mmp
mmin_base< multi_funct > * mmp
The minimizer.
Definition: fit_min.h:182
o2scl::fit_nonlin_b< boost::numeric::ublas::vector< double >, boost::numeric::ublas::matrix< double > >::compute_delta
double compute_delta(boost::numeric::ublas::vector< double > &diag2, size_t n, const boost::numeric::ublas::vector< double > &x)
Desc.
Definition: fit_nonlin.h:243
o2scl::fit_nonlin_b< boost::numeric::ublas::vector< double >, boost::numeric::ublas::matrix< double > >::compute_diag
int compute_diag(size_t nparm, size_t ndata, const boost::numeric::ublas::matrix< double > &J, boost::numeric::ublas::vector< double > &diag_vec)
Compute the root of the sum of the squares of the columns of J.
Definition: fit_nonlin.h:680
o2scl::mmin_base::mmin
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.
o2scl_cblas::dnrm2
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
Definition: cblas_base.h:183
o2scl::fit_nonlin_b< boost::numeric::ublas::vector< double >, boost::numeric::ublas::matrix< double > >::covariance
int covariance(size_t m, size_t n, const boost::numeric::ublas::matrix< double > &J, boost::numeric::ublas::matrix< double > &covar, boost::numeric::ublas::vector< double > &norm, boost::numeric::ublas::matrix< double > &r, boost::numeric::ublas::vector< double > &tau, permutation &perm, double epsrel)
Compute the covarance matrix covar given the Jacobian J.
Definition: fit_nonlin.h:721
o2scl::fit_min::min_func
double min_func(size_t np, const vec_t &xp)
The function to minimize, .
Definition: fit_min.h:188
o2scl::fit_nonlin_b< boost::numeric::ublas::vector< double >, boost::numeric::ublas::matrix< double > >::scaled_enorm
double scaled_enorm(const boost::numeric::ublas::vector< double > &d, size_t n, const boost::numeric::ublas::vector< double > &f)
Euclidean norm of vector f of length n, scaled by vector d.
Definition: fit_nonlin.h:231

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