fit_fix.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_FIX_H
24 #define O2SCL_FIT_FIX_H
25 
26 /** \file fit_fix.h
27  \brief File defining \ref o2scl::fit_fix_pars
28 */
29 
30 #include <boost/numeric/ublas/vector.hpp>
31 #include <boost/numeric/ublas/matrix.hpp>
32 
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 Multidimensional fitting class fixing some parameters and
41  varying others
42 
43  The number of trials used in the fit can be specified in
44  the data member of the parent class \ref fit_base::ntrial
45  associated with the fit_fix_pars object. Similarly for the
46  verbosity parameter in \ref fit_base::verbose, the absolute
47  tolerance in \ref fit_base::tol_abs, and the relative tolerance
48  in \ref fit_base::tol_abs. These values are copied to the
49  minimizer used by <tt>fit_fix_pars::mmin()</tt> during each call.
50  After the minimizer is called, the value of \ref fit_base::ntrial
51  associated with the \ref fit_fix_pars object is filled
52  with the last number of trials required for the last
53  minimization.
54  \comment
55  For some reason the reference to mmin() above doesn't work
56  in doxygen.
57  \endcomment
58 
59  Default template arguments
60  - \c func_t - \ref gen_fit_funct<>
61  - \c vec_t - \ref boost::numeric::ublas::vector <double >
62  - \c mat_t - \ref boost::numeric::ublas::matrix <double >
63  */
64  template<class bool_vec_t, class func_t=gen_fit_funct<>,
65  class vec_t=boost::numeric::ublas::vector<double>,
66  class mat_t=boost::numeric::ublas::matrix<double> >
67  class fit_fix_pars : public fit_base<func_t,vec_t,mat_t>,
68  public gen_fit_funct<vec_t,mat_t> {
69 
70  public:
71 
72  /// The generic fitter type
74  vec_t,mat_t> base_fit_t;
75 
76  /// The default fitter type
78  vec_t,mat_t> def_fit_t;
79 
80  /** \brief Specify the member function pointer
81  */
83  fitp=&def_fit;
84  expand_covar=false;
85  }
86 
87  virtual ~fit_fix_pars() {}
88 
89  /** \brief If true, expand the covariance matrix to the
90  larger space by filling with the identity matrix (default false)
91 
92  If this varable is false (the default), then the covariance
93  matrix is computed in the smaller space which enumerates only
94  the parameters which are not fixed. If this variable is true,
95  then the covariance matrix must be a full \c np by \c np
96  matrix (where \c np is the argument to \ref fit() or \ref
97  fit_fix() ) and rows and columns which correspond with
98  parameters which are fixed are replaced by elements from the
99  identity matrix.
100 
101  The optimal parameters and \f$ \chi^2 \f$ reported by the
102  fit are unchanged.
103  */
105 
106  /** \brief Fit the data specified in (xdat,ydat) to
107  the function \c fitfun with the parameters in \c par.
108 
109  The covariance matrix for the parameters is returned in \c covar
110  and the value of \f$ \chi^2 \f$ is returned in \c chi2.
111  */
112  virtual int fit(size_t np, vec_t &par, mat_t &covar, double &chi2,
113  func_t &fitfun) {
114 
115  x_par=&par;
116  funcp=&fitfun;
117  u_np=np;
118  u_np_new=np;
119  size_t nd=fitfun.get_ndata();
120 
121  u_par.resize(np);
122  J.resize(nd,np);
123 
124  fitp->verbose=this->verbose;
125  fitp->ntrial=this->ntrial;
126  fitp->tol_rel=this->tol_rel;
127  fitp->tol_abs=this->tol_abs;
128  fitp->fit(u_np_new,par,covar,chi2,*this);
129 
130  u_par.clear();
131  J.clear();
132 
133  return 0;
134  }
135 
136  /** \brief Fit function \c func while fixing some parameters as
137  specified in \c fix.
138 
139  \note When some parameters are fixed, each fixed parameter
140  leads to a row and column in the covariance matrix which
141  is unused by this function. If there are <tt>npar2<npar</tt>
142  unfixed parameters, then only the first <tt>npar2</tt> rows
143  and columns of \c covar are filled.
144  */
145  virtual int fit_fix(size_t np, vec_t &par, mat_t &covar, double &chi2,
146  func_t &fitfun, bool_vec_t &fix) {
147 
148  x_par=&par;
149  funcp=&fitfun;
150  u_np=np;
151  fix_par=&fix;
152  size_t nd=fitfun.get_ndata();
153 
154  // Count number of new parameters
155  u_np_new=0;
156  for(size_t i=0;i<np;i++) {
157  if (fix[i]==false) u_np_new++;
158  }
159  if (u_np_new==0) return 0;
160 
161  // Allocate space
162  u_par.resize(np);
163  u_par_new.resize(u_np_new);
164  J.resize(nd,np);
165 
166  // Copy to smaller vector containing only parameters to be varied
167  size_t j=0;
168  for(size_t i=0;i<np;i++) {
169  if (fix[i]==false) {
170  u_par_new[j]=par[i];
171  j++;
172  }
173  }
174 
175  fitp->verbose=this->verbose;
176  fitp->ntrial=this->ntrial;
177  fitp->tol_rel=this->tol_rel;
178  fitp->tol_abs=this->tol_abs;
179  fitp->fit(u_np_new,u_par_new,covar,chi2,*this);
180 
181  // Copy optimal parameter set back to initial guess object
182  j=0;
183  for(size_t i=0;i<np;i++) {
184  if (fix[i]==false) {
185  par[i]=u_par_new[j];
186  j++;
187  }
188  }
189 
190  // If required, expand covariance matrix
191  if (expand_covar && u_np_new<u_np) {
192  int i_new=((int)u_np_new)-1;
193  for(int i=((int)np)-1;i>=0;i--) {
194  int k_new=((int)u_np_new)-1;
195  for(int k=((int)np)-1;k>=0;k--) {
196  if (fix[i]==false && fix[k]==false) {
197  if (i_new<0 || k_new<0 ||
198  i_new>=((int)u_np_new) || k_new>=((int)u_np_new)) {
199  O2SCL_ERR2("Covariance alignment in ",
200  "fit_fix::fit_fix().",exc_esanity);
201  }
202  covar(i,k)=covar(i_new,k_new);
203  k_new--;
204  } else if (i==k) {
205  covar(i,k)=1.0;
206  } else {
207  covar(i,k)=0.0;
208  }
209  }
210  if (fix[i]==false) i_new--;
211  }
212  }
213 
214  // Free space
215  u_par_new.clear();
216  u_par.clear();
217  J.clear();
218 
219  return 0;
220  }
221 
222  /// Change the base fitter
223  int set_fit(base_fit_t &fitter) {
224  fitp=&fitter;
225  return 0;
226  }
227 
228  /// The default base fitter
229  def_fit_t def_fit;
230 
231  /// \name Reimplementation of \ref gen_fit_funct
232  //@{
233  /// The function to return the number of data points
234  virtual size_t get_ndata() {
235  return (*funcp).get_ndata();
236  }
237 
238  /** \brief The function computing deviations
239 
240  This function operates in the smaller space of size np_new.
241  */
242  virtual void operator()(size_t np_new, const vec_t &par_new,
243  size_t nd, vec_t &f) {
244 
245  // Variable 'i' runs over the user space and 'i_new' runs
246  // over the smaller space of size np_new.
247  size_t i_new=0;
248  for(size_t i=0;i<u_np;i++) {
249  if (u_np_new<u_np && (*fix_par)[i]==true) {
250  u_par[i]=(*x_par)[i];
251  } else {
252  u_par[i]=par_new[i_new];
253  i_new++;
254  }
255  }
256  if (i_new!=np_new) {
257  O2SCL_ERR("Alignment failure 1 in fit_fix::operator().",
258  exc_esanity);
259  }
260 
261  // Call the function in the larger space of size np
262  (*funcp)(u_np,u_par,nd,f);
263 
264  return;
265  }
266 
267  /** \brief The function computing the Jacobian
268  */
269  virtual void jac(size_t np_new, vec_t &par_new, size_t nd, vec_t &f,
270  mat_t &J_new) {
271 
272  size_t i_new=0;
273  for(size_t i=0;i<u_np;i++) {
274  if (u_np_new<u_np && (*fix_par)[i]==true) {
275  u_par[i]=(*x_par)[i];
276  } else {
277  u_par[i]=par_new[i_new];
278  i_new++;
279  }
280  }
281  if (i_new!=np_new) {
282  O2SCL_ERR("Alignment failure 2 in fit_fix::operator().",
283  exc_esanity);
284  }
285 
286  // Call the Jacobian in the larger space of size (nd,u_np)
287  (*funcp).jac(u_np,u_par,nd,f,J);
288 
289  // Copy the Jacobian over to the smaller space, ignoring
290  // rows corresponding to parameters which are fixed
291  for(size_t id=0;id<nd;id++) {
292  i_new=0;
293  for(size_t i=0;i<u_np;i++) {
294  if (u_np_new==u_np || (*fix_par)[i]==false) {
295  J_new(id,i_new)=J(id,i);
296  i_new++;
297  }
298  }
299  }
300 
301  return;
302  }
303  //@}
304 
305 #ifndef DOXYGEN_INTERNAL
306 
307  protected:
308 
309  /// Temporary vector to store full parameter list of size u_np
310  vec_t u_par;
311 
312  /// Vector for smaller parameter list of size u_np_new
313  vec_t u_par_new;
314 
315  /// Jacobian in the user space of size (nd,u_np)
316  mat_t J;
317 
318  /// The fitter
319  base_fit_t *fitp;
320 
321  /// The user-specified function
322  func_t *funcp;
323 
324  /// The user-specified number of variables
325  size_t u_np;
326 
327  /// The new number of variables
328  size_t u_np_new;
329 
330  /// Specify which parameters to fix (vector of size u_np)
331  bool_vec_t *fix_par;
332 
333  /// The user-specified initial vector of size u_np
334  vec_t *x_par;
335 
336  private:
337 
338  fit_fix_pars(const fit_fix_pars &);
339  fit_fix_pars& operator=(const fit_fix_pars&);
340 
341 #endif
342 
343  };
344 
345 
346 #ifndef DOXYGEN_NO_O2NS
347 }
348 #endif
349 
350 #endif
int set_fit(base_fit_t &fitter)
Change the base fitter.
Definition: fit_fix.h:223
base_fit_t * fitp
The fitter.
Definition: fit_fix.h:319
def_fit_t def_fit
The default base fitter.
Definition: fit_fix.h:229
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
mat_t J
Jacobian in the user space of size (nd,u_np)
Definition: fit_fix.h:316
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
fit_fix_pars()
Specify the member function pointer.
Definition: fit_fix.h:82
Non-linear least-squares fitting class (GSL)
Definition: fit_nonlin.h:921
virtual int fit_fix(size_t np, vec_t &par, mat_t &covar, double &chi2, func_t &fitfun, bool_vec_t &fix)
Fit function func while fixing some parameters as specified in fix.
Definition: fit_fix.h:145
vec_t u_par_new
Vector for smaller parameter list of size u_np_new.
Definition: fit_fix.h:313
vec_t * x_par
The user-specified initial vector of size u_np.
Definition: fit_fix.h:334
size_t ntrial
Maximum number of iterations (default 500)
Definition: fit_base.h:341
bool_vec_t * fix_par
Specify which parameters to fix (vector of size u_np)
Definition: fit_fix.h:331
fit_base< fit_fix_pars< bool_vec_t, func_t, vec_t, mat_t >, vec_t, mat_t > base_fit_t
The generic fitter type.
Definition: fit_fix.h:74
virtual size_t get_ndata()
The function to return the number of data points.
Definition: fit_fix.h:234
virtual int fit(size_t np, 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_fix.h:112
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
bool expand_covar
If true, expand the covariance matrix to the larger space by filling with the identity matrix (defaul...
Definition: fit_fix.h:104
func_t * funcp
The user-specified function.
Definition: fit_fix.h:322
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
double tol_rel
(default 1.0e-4)
Definition: fit_base.h:347
int verbose
An integer describing the verbosity of the output.
Definition: fit_base.h:391
Multidimensional fitting class fixing some parameters and varying others.
Definition: fit_fix.h:67
fit_nonlin< fit_fix_pars< bool_vec_t, func_t, vec_t, mat_t >, vec_t, mat_t > def_fit_t
The default fitter type.
Definition: fit_fix.h:78
vec_t u_par
Temporary vector to store full parameter list of size u_np.
Definition: fit_fix.h:310
size_t u_np_new
The new number of variables.
Definition: fit_fix.h:328
virtual void operator()(size_t np_new, const vec_t &par_new, size_t nd, vec_t &f)
The function computing deviations.
Definition: fit_fix.h:242
Non-linear least-squares fitting [abstract base].
Definition: fit_base.h:327
virtual void jac(size_t np_new, vec_t &par_new, size_t nd, vec_t &f, mat_t &J_new)
The function computing the Jacobian.
Definition: fit_fix.h:269
double tol_abs
Absolute tolerance (default 1.0e-4)
Definition: fit_base.h:344
size_t u_np
The user-specified number of variables.
Definition: fit_fix.h:325
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.

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