fit_base.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_BASE_H
24 #define O2SCL_FIT_BASE_H
25 
26 /** \file fit_base.h
27  \brief File defining \ref o2scl::fit_base and fitting functions
28 */
29 
30 #include <string>
31 
32 #include <o2scl/jacobian.h>
33 #include <o2scl/mm_funct.h>
34 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39  /** \brief Array of multi-dimensional functions typedef (C++11
40  version) in src/fit/fit_base.h
41  */
42  typedef std::function<
43  double(size_t,const boost::numeric::ublas::vector<double> &,
44  double)> fit_funct;
45 
46  /** \brief String fitting function
47 
48  Default template arguments
49  - \c vec_t - \ref boost::numeric::ublas::vector < double >
50  */
52 
53  public:
54 
55  /** \brief Specify a fitting function through a string
56  */
57  template<class vec_string_t=std::vector<std::string> >
58  fit_funct_strings(std::string expr, vec_string_t &parms,
59  std::string var) {
60  calc.compile(expr.c_str(),&vars);
61  st_form=expr;
62  int np=parms.size();
63  st_parms.resize(np);
64  for (int i=0;i<np;i++) {
65  st_parms[i]=parms[i];
66  }
67  st_var=var;
68  }
69 
70  /** \brief Set the values of the auxilliary parameters that were
71  specified in \c auxp in the constructor
72  */
73  int set_aux_parm(std::string name, double val) {
74  vars[name]=val;
75  return 0;
76  }
77 
78  /** \brief Using parameters in \c p, predict \c y given \c x
79  */
80  template<class vec_t=boost::numeric::ublas::vector<double> >
81  double operator()(size_t np, const vec_t &p, double x) {
82 
83  for(size_t i=0;i<np;i++) {
84  vars[st_parms[i]]=p[i];
85  }
86  vars[st_var]=x;
87  double y=calc.eval(&vars);
88  //std::cout << "Debug: " << calc.RPN_to_string() << std::endl;
89  //std::cout << "Here: " << y << " " << p[0]*exp(x)+p[1]*sqrt(x)
90  //<< std::endl;
91  //exit(-1);
92  return y;
93  }
94 
95 #ifndef DOXYGEN_INTERNAL
96 
97  protected:
98 
99  /// The function parser
101 
102  /// Desc
103  std::map<std::string,double> vars;
104 
105  /// The expression
106  std::string st_form;
107 
108  /// The parameters
109  std::vector<std::string> st_parms;
110 
111  /// The variable
112  std::string st_var;
113 
114  fit_funct_strings() {};
115 
116  /// Specify the strings which define the fitting function
117  int set_function(std::string expr, std::string parms,
118  std::string var, int nauxp=0, std::string auxp="");
119 
120  private:
121 
123  fit_funct_strings& operator=(const fit_funct_strings&);
124 
125 #endif
126 
127  };
128 
129  /** \brief Generalized fitting function [abstract base]
130 
131  Default template arguments
132  - \c vec_t - \ref boost::numeric::ublas::vector < double >
133  - \c mat_t - \ref boost::numeric::ublas::matrix < double >
134  */
135  template<class vec_t=boost::numeric::ublas::vector<double>,
136  class mat_t=boost::numeric::ublas::matrix<double> >
138 
139  public:
140 
141  gen_fit_funct() {}
142 
143  virtual ~gen_fit_funct() {}
144 
145  /** \brief Using parameters in \c p, compute the
146  relative deviations in \c f
147  */
148  virtual void operator()(size_t np, const vec_t &p, size_t nd,
149  vec_t &f)=0;
150 
151  /** \brief Using parameters in \c p, compute the Jacobian
152  in \c J
153  */
154  virtual void jac(size_t np, vec_t &p, size_t nd, vec_t &f,
155  mat_t &J)=0;
156 
157  /// Return the number of data points
158  virtual size_t get_ndata()=0;
159 
160 #ifndef DOXYGEN_INTERNAL
161 
162  private:
163 
164  gen_fit_funct(const gen_fit_funct &);
165  gen_fit_funct& operator=(const gen_fit_funct&);
166 
167 #endif
168 
169  };
170 
171  /** \brief Standard fitting function based on one-dimensional
172  data with a numerical Jacobian
173 
174  This class specifies the deviations (in <tt>operator()</tt> ) and
175  Jacobian (in \ref jac()) for a fitting class like \ref fit_nonlin.
176  It assumes a one-dimensional data set with no uncertainty in the
177  abcissae and a fitting function specified in a form similar to
178  \ref fit_funct.
179  \comment
180  For some reason the reference to operator() above doesn't work
181  in doxygen.
182  \endcomment
183 
184  The default method for numerically computing the Jacobian is
185  from \ref jacobian_gsl. This default is identical to the GSL
186  approach, except that the default value of \ref
187  jacobian_gsl::epsmin is non-zero. See \ref jacobian_gsl for more
188  details.
189 
190  Default template arguments
191  - \c vec_t - \ref boost::numeric::ublas::vector < double >
192  - \c mat_t - \ref boost::numeric::ublas::matrix < double >
193  - \c func_t - \ref fit_funct
194 
195  \future Allow a user-specified Jacobian or make that into
196  a separate class?
197  \future Default constructor?
198  */
199  template<class vec_t=boost::numeric::ublas::vector<double>,
200  class mat_t=boost::numeric::ublas::matrix<double>,
201  class fit_func_t=fit_funct> class chi_fit_funct :
202  public gen_fit_funct<vec_t,mat_t> {
203 
204  public:
205 
206  /** \brief Create an object with specified data and specified
207  fitting function
208  */
209  chi_fit_funct(size_t ndat, const vec_t &xdat, const vec_t &ydat,
210  const vec_t &yerr, fit_func_t &fun) {
211  ndat_=ndat;
212  xdat_=&xdat;
213  ydat_=&ydat;
214  yerr_=&yerr;
215  fun_=&fun;
216 
217  mfm=std::bind(std::mem_fn<int(size_t,const vec_t &,vec_t &)>
219  std::placeholders::_1,std::placeholders::_2,
220  std::placeholders::_3);
221 
222  auto_jac.set_function(mfm);
223  //double sqrt_dbl_eps=sqrt(std::numeric_limits<double>::epsilon());
224  //auto_jac.set_epsrel(sqrt_dbl_eps);
225  }
226 
227  /** \brief Set the data to be fit
228  */
229  void set_data(size_t ndat, const vec_t &xdat, const vec_t &ydat,
230  const vec_t &yerr) {
231  ndat_=ndat;
232  xdat_=&xdat;
233  ydat_=&ydat;
234  yerr_=&yerr;
235  return;
236  }
237 
238  /** \brief Set the fitting function
239  */
240  void set_func(fit_func_t &fun) {
241  fun_=&fun;
242  return;
243  }
244 
245  /** \brief Return \f$ \chi^2 \f$
246  */
247  virtual double chi2(size_t np, const vec_t &p) {
248  double ret=0.0;
249  for(size_t i=0;i<ndat_;i++) {
250  double yi=((*fun_)(np,p,(*xdat_)[i])-(*ydat_)[i])/((*yerr_)[i]);
251  ret+=yi*yi;
252  }
253  return ret;
254  }
255 
256  /** \brief Using parameters in \c p, compute the
257  relative deviations in \c f
258  */
259  virtual void operator()(size_t np, const vec_t &p, size_t nd, vec_t &f) {
260 
261  for(size_t i=0;i<nd;i++) {
262  double yi=(*fun_)(np,p,(*xdat_)[i]);
263  f[i]=(yi-(*ydat_)[i])/((*yerr_)[i]);
264  }
265  return;
266  }
267 
268  /** \brief Using parameters in \c p, compute the Jacobian
269  in \c J
270  */
271  virtual void jac(size_t np, vec_t &p, size_t nd, vec_t &f,
272  mat_t &J) {
273 
274  auto_jac(np,p,nd,f,J);
275 
276  return;
277  }
278 
279  /// Return the number of data points
280  virtual size_t get_ndata() {
281  return ndat_;
282  }
283 
284  /// Automatic Jacobian object
285  jacobian_gsl<std::function<int(size_t,const vec_t &,vec_t &)>,
286  vec_t,mat_t> auto_jac;
287 
288 #ifndef DOXYGEN_INTERNAL
289 
290  protected:
291 
292  /// Reformulate <tt>operator()</tt> into a \ref mm_funct object
293  int jac_mm_funct(size_t np, const vec_t &p, vec_t &f) {
294  operator()(np,p,ndat_,f);
295  return 0;
296  }
297 
298  /// Function object for Jacobian object
299  std::function<int(size_t,const vec_t &,vec_t &)> mfm;
300 
301  /// \name Data and uncertainties
302  //@{
303  size_t ndat_;
304  const vec_t *xdat_;
305  const vec_t *ydat_;
306  const vec_t *yerr_;
307  //@}
308 
309  /// Fitting function
310  fit_func_t *fun_;
311 
312  private:
313 
314  chi_fit_funct(const chi_fit_funct &);
315  chi_fit_funct& operator=(const chi_fit_funct&);
316 
317 #endif
318 
319  };
320 
321  // ----------------------------------------------------------------
322  // Fitter base
323  // ----------------------------------------------------------------
324 
325  /** \brief Non-linear least-squares fitting [abstract base]
326  */
327  template<class func_t=fit_funct,
330 
331  public:
332 
333  fit_base() {
334  ntrial=500;
335  tol_abs=1.0e-4;
336  tol_rel=1.0e-4;
337  verbose=0;
338  }
339 
340  virtual ~fit_base() {}
341 
342  /// Maximum number of iterations (default 500)
343  size_t ntrial;
344 
345  /// Absolute tolerance (default 1.0e-4)
346  double tol_abs;
347 
348  /// (default 1.0e-4)
349  double tol_rel;
350 
351  /** \brief Print out iteration information.
352 
353  Depending on the value of the variable verbose, this prints out
354  the iteration information. If verbose=0, then no information is
355  printed, while if verbose>1, then after each iteration, the
356  present values of x and y are output to std::cout along with the
357  iteration number. If verbose>=2 then each iteration waits for a
358  character.
359  */
360  virtual int print_iter(size_t nv, vec_t &x, double y, int iter,
361  double value=0.0, double limit=0.0) {
362  if (verbose<=0) return 0;
363 
364  size_t i;
365  char ch;
366 
367  std::cout << "Iteration: " << iter << std::endl;
368  std::cout << "x: ";
369  for(i=0;i<nv;i++) std::cout << x[i] << " ";
370  std::cout << std::endl;
371  std::cout << "y: " << y << " Val: " << value << " Lim: " << limit
372  << std::endl;
373  if (verbose>1) {
374  std::cout << "Press a key and type enter to continue. ";
375  std::cin >> ch;
376  }
377 
378  return 0;
379  }
380 
381  /** \brief Fit function \c fitfun using parameters in \c parms as
382  initial guesses
383 
384  The covariance matrix for the parameters is returned in \c covar
385  and the value of \f$ \chi^2 \f$ is returned in \c chi2.
386 
387  */
388  virtual int fit(size_t npar, vec_t &parms, mat_t &covar,
389  double &chi2, func_t &fitfun)=0;
390 
391  /** \brief An integer describing the verbosity of the output
392  */
393  int verbose;
394 
395  /// Return string denoting type ("fit_base")
396  virtual const char *type() { return "fit_base"; }
397 
398  /// The number of data points
399  size_t n_dat;
400 
401  /// The number of parameters
402  size_t n_par;
403 
404 
405  };
406 
407 #ifndef DOXYGEN_NO_O2NS
408 }
409 #endif
410 
411 #endif
boost::numeric::ublas::matrix< double >
o2scl::chi_fit_funct
Standard fitting function based on one-dimensional data with a numerical Jacobian.
Definition: fit_base.h:201
o2scl::chi_fit_funct::auto_jac
jacobian_gsl< std::function< int(size_t, const vec_t &, vec_t &)>, vec_t, mat_t > auto_jac
Automatic Jacobian object.
Definition: fit_base.h:286
o2scl::fit_base::tol_rel
double tol_rel
(default 1.0e-4)
Definition: fit_base.h:349
boost::numeric::ublas::vector< double >
o2scl::fit_funct_strings::calc
calculator calc
The function parser.
Definition: fit_base.h:100
o2scl::chi_fit_funct::jac
virtual void jac(size_t np, vec_t &p, size_t nd, vec_t &f, mat_t &J)
Using parameters in p, compute the Jacobian in J.
Definition: fit_base.h:271
o2scl::gen_fit_funct::operator()
virtual void operator()(size_t np, const vec_t &p, size_t nd, vec_t &f)=0
Using parameters in p, compute the relative deviations in f.
o2scl::fit_funct_strings::set_function
int set_function(std::string expr, std::string parms, std::string var, int nauxp=0, std::string auxp="")
Specify the strings which define the fitting function.
o2scl::fit_base::type
virtual const char * type()
Return string denoting type ("fit_base")
Definition: fit_base.h:396
o2scl::fit_funct_strings::st_var
std::string st_var
The variable.
Definition: fit_base.h:112
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::fit_funct_strings::st_form
std::string st_form
The expression.
Definition: fit_base.h:106
o2scl::fit_base::n_dat
size_t n_dat
The number of data points.
Definition: fit_base.h:399
o2scl::fit_base::tol_abs
double tol_abs
Absolute tolerance (default 1.0e-4)
Definition: fit_base.h:346
o2scl::chi_fit_funct::chi2
virtual double chi2(size_t np, const vec_t &p)
Return .
Definition: fit_base.h:247
o2scl::calculator::compile
void compile(const char *expr, std::map< std::string, double > *vars=0, bool debug=false, std::map< std::string, int > opPrec=opPrecedence)
Compile expression expr using variables specified in vars and return an integer to indicate success o...
o2scl::fit_base::fit
virtual int fit(size_t npar, vec_t &parms, mat_t &covar, double &chi2, func_t &fitfun)=0
Fit function fitfun using parameters in parms as initial guesses.
o2scl::fit_base::print_iter
virtual int print_iter(size_t nv, vec_t &x, double y, int iter, double value=0.0, double limit=0.0)
Print out iteration information.
Definition: fit_base.h:360
o2scl::fit_funct_strings::st_parms
std::vector< std::string > st_parms
The parameters.
Definition: fit_base.h:109
o2scl::fit_funct_strings::set_aux_parm
int set_aux_parm(std::string name, double val)
Set the values of the auxilliary parameters that were specified in auxp in the constructor.
Definition: fit_base.h:73
o2scl::gen_fit_funct
Generalized fitting function [abstract base].
Definition: fit_base.h:137
o2scl::fit_funct_strings::operator()
double operator()(size_t np, const vec_t &p, double x)
Using parameters in p, predict y given x.
Definition: fit_base.h:81
o2scl::calculator
Evaluate a mathematical expression in a string.
Definition: shunting_yard.h:119
o2scl::chi_fit_funct::get_ndata
virtual size_t get_ndata()
Return the number of data points.
Definition: fit_base.h:280
o2scl::fit_funct_strings
String fitting function.
Definition: fit_base.h:51
o2scl::chi_fit_funct::chi_fit_funct
chi_fit_funct(size_t ndat, const vec_t &xdat, const vec_t &ydat, const vec_t &yerr, fit_func_t &fun)
Create an object with specified data and specified fitting function.
Definition: fit_base.h:209
o2scl::fit_funct
std::function< double(size_t, const boost::numeric::ublas::vector< double > &, double)> fit_funct
Array of multi-dimensional functions typedef (C++11 version) in src/fit/fit_base.h.
Definition: fit_base.h:44
o2scl::fit_funct_strings::vars
std::map< std::string, double > vars
Desc.
Definition: fit_base.h:103
o2scl::fit_base::n_par
size_t n_par
The number of parameters.
Definition: fit_base.h:402
o2scl::jacobian_gsl
Simple automatic Jacobian.
Definition: jacobian.h:144
o2scl::fit_base
Non-linear least-squares fitting [abstract base].
Definition: fit_base.h:329
o2scl::fit_base::ntrial
size_t ntrial
Maximum number of iterations (default 500)
Definition: fit_base.h:343
o2scl::chi_fit_funct::fun_
fit_func_t * fun_
Fitting function.
Definition: fit_base.h:310
o2scl::chi_fit_funct::jac_mm_funct
int jac_mm_funct(size_t np, const vec_t &p, vec_t &f)
Reformulate operator() into a mm_funct object.
Definition: fit_base.h:293
o2scl::fit_funct_strings::fit_funct_strings
fit_funct_strings(std::string expr, vec_string_t &parms, std::string var)
Specify a fitting function through a string.
Definition: fit_base.h:58
o2scl::chi_fit_funct::operator()
virtual void operator()(size_t np, const vec_t &p, size_t nd, vec_t &f)
Using parameters in p, compute the relative deviations in f.
Definition: fit_base.h:259
o2scl::chi_fit_funct::set_func
void set_func(fit_func_t &fun)
Set the fitting function.
Definition: fit_base.h:240
o2scl::calculator::eval
double eval(std::map< std::string, double > *vars=0)
Evalate the previously compiled expression using variables specified in vars.
o2scl::fit_base::verbose
int verbose
An integer describing the verbosity of the output.
Definition: fit_base.h:393
o2scl::gen_fit_funct::get_ndata
virtual size_t get_ndata()=0
Return the number of data points.
o2scl::chi_fit_funct::set_data
void set_data(size_t ndat, const vec_t &xdat, const vec_t &ydat, const vec_t &yerr)
Set the data to be fit.
Definition: fit_base.h:229
o2scl::gen_fit_funct::jac
virtual void jac(size_t np, vec_t &p, size_t nd, vec_t &f, mat_t &J)=0
Using parameters in p, compute the Jacobian in J.
o2scl::chi_fit_funct::mfm
std::function< int(size_t, const vec_t &, vec_t &)> mfm
Function object for Jacobian object.
Definition: fit_base.h:299

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