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

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