jacobian.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_JACOBIAN_H
24 #define O2SCL_JACOBIAN_H
25 
26 /** \file jacobian.h
27  \brief File for Jacobian evaluation and function classes
28 */
29 
30 #include <string>
31 #include <o2scl/mm_funct.h>
32 #include <o2scl/deriv_gsl.h>
33 #include <o2scl/columnify.h>
34 #include <o2scl/vector.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /// Jacobian function (not necessarily square)
41  typedef std::function<
45 
46  /** \brief Base for providing a numerical jacobian [abstract base]
47 
48  This is provides a Jacobian which is numerically determined
49  by differentiating a user-specified function (typically
50  of the form of \ref mm_funct11).
51 
52  By convention, the Jacobian is stored in the order
53  <tt>J[i][j]</tt> (or <tt>J(i,j)</tt>) where the rows have index
54  \c i which runs from 0 to <tt>ny-1</tt> and the columns have
55  index \c j with runs from 0 to <tt>nx-1</tt>.
56 
57  Default template arguments
58  - \c func_t - \ref mm_funct11
59  - \c vec_t - boost::numeric::ublas::vector<double>
60  - \c mat_t - boost::numeric::ublas::matrix<double>
61  */
62  template<class func_t=mm_funct11,
65 
66  public:
67 
68  jacobian() {
69  err_nonconv=true;
70  };
71 
72  virtual ~jacobian() {};
73 
74  /// If true, call the error handler if the routine does not converge
75  bool err_nonconv;
76 
77  /// Set the function to compute the Jacobian of
78  virtual int set_function(func_t &f) {
79  func=f;
80  return 0;
81  }
82 
83  /** \brief Evaluate the Jacobian \c j at point \c y(x)
84  */
85  virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y,
86  mat_t &j)=0;
87 
88 #ifndef DOXYGEN_INTERNAL
89 
90  protected:
91 
92  /// A pointer to the user-specified function
93  func_t func;
94 
95  private:
96 
97  jacobian(const jacobian &);
98  jacobian& operator=(const jacobian&);
99 
100 #endif
101 
102  };
103 
104  /** \brief Simple automatic Jacobian
105 
106  This class computes a numerical Jacobian by finite differencing.
107  The stepsize is initially chosen to be
108  \f$ h_j = \mathrm{max}(\mathrm{epsrel}~|x_j|,\mathrm{epsmin}) \f$.
109  Then if \f$ h_j = 0 \f$, the value of \f$ h_j \f$ is set to
110  \f$ \mathrm{epsrel}) \f$ .
111 
112  Values of \c epsmin which are non-zero are useful, for example,
113  in \ref mroot_hybrids when one of the variables is either
114  very small or zero, so that the step size doesn't become too
115  small.
116 
117  If the function evaluation leads to a non-zero return value,
118  then the step size is alternately flipped in sign or decreased
119  by a fixed factor (default \f$ 10^2 \f$, set in \ref
120  set_shrink_fact() ) in order to obtain a valid result. This
121  process is repeated a fixed number of times (default 10, set in
122  \ref set_max_shrink_iters() ).
123 
124  This is equivalent to the GSL method for computing Jacobians as
125  in \c multiroots/fdjac.c if one calls \ref
126  set_max_shrink_iters() with a parameter value of zero.
127 
128  If one row of the Jacobian is all zero, or if there was no
129  step-size found which would give a zero return value from
130  the user-specified function, then the error handler is called
131  depending on the value of \ref err_nonconv.
132 
133  This class does not separately check the vector and matrix sizes
134  to ensure they are commensurate.
135 
136  Default template arguments
137  - \c func_t - \ref mm_funct11
138  - \c vec_t - boost::numeric::ublas::vector<double>
139  - \c mat_t - boost::numeric::ublas::matrix<double>
140  */
141  template<class func_t=mm_funct11,
144  class jacobian_gsl : public jacobian<func_t,vec_t,mat_t> {
145 
146 #ifndef DOXYGEN_INTERNAL
147 
148  protected:
149 
150  /// Function values
151  vec_t f;
152 
153  /// Function arguments
154  vec_t xx;
155 
156  /// Size of allocated memory in x
157  size_t mem_size_x;
158 
159  /// Size of allocated memory in y
160  size_t mem_size_y;
161 
162  /** \brief The relative stepsize for finite-differencing
163  */
164  double epsrel;
165 
166  /// The minimum stepsize
167  double epsmin;
168 
169  /// Maximum number of times to shrink the step size
171 
172  /// Factor to shrink stepsize by
173  double shrink_fact;
174 
175 #endif
176 
177  public:
178 
179  jacobian_gsl() {
180  epsrel=sqrt(std::numeric_limits<double>::epsilon());
181  epsmin=0.0;
182  mem_size_x=0;
183  mem_size_y=0;
184  max_shrink_iters=10;
185  shrink_fact=1.0e2;
186  }
187 
188  virtual ~jacobian_gsl() {
189  }
190 
191  /** \brief Get the relative stepsize (default \f$ 10^{-4} \f$ )
192  */
193  double get_epsrel() { return epsrel; }
194 
195  /** \brief Get the minimum stepsize (default \f$ 10^{-15} \f$)
196  */
197  double get_epsmin() { return epsmin; }
198 
199  /** \brief Set the relative stepsize (must be \f$ > 0 \f$)
200  */
201  void set_epsrel(double l_epsrel) {
202  if (l_epsrel<=0.0) {
203  O2SCL_ERR2("Negative or zero value specified in ",
204  "jacobian_gsl::set_epsrel().",exc_einval);
205  }
206  epsrel=l_epsrel;
207  return;
208  }
209 
210  /** \brief Set the minimum stepsize (must be \f$ \geq 0 \f$)
211  */
212  void set_epsmin(double l_epsmin) {
213  if (l_epsmin<0.0) {
214  O2SCL_ERR("Negative value specified in jacobian_gsl::set_epsmin().",
215  exc_einval);
216  }
217  epsmin=l_epsmin;
218  return;
219  }
220 
221  /** \brief Set shrink factor for decreasing step size
222  */
223  void set_shrink_fact(double l_shrink_fact) {
224  if (l_shrink_fact<0.0) {
225  O2SCL_ERR("Negative value specified in jacobian_gsl::set_shrink_fact().",
226  exc_einval);
227  }
228  shrink_fact=l_shrink_fact;
229  return;
230  }
231 
232  /** \brief Set number of times to decrease step size
233  */
234  void set_max_shrink_iters(size_t it) {
235  max_shrink_iters=it;
236  return;
237  }
238 
239  /** \brief The operator()
240  */
241  virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y,
242  mat_t &jac) {
243 
244  size_t i,j;
245  double h,temp;
246  bool success=true;
247 
248  if (mem_size_x!=nx || mem_size_y!=ny) {
249  f.resize(ny);
250  xx.resize(nx);
251  mem_size_x=nx;
252  mem_size_y=ny;
253  }
254 
255  vector_copy(nx,x,xx);
256 
257  for (j=0;j<nx;j++) {
258 
259  // Thanks to suggestion from Conrad Curry.
260  h=epsrel*fabs(x[j]);
261  if (h<epsmin) h=epsmin;
262  if (h==0.0) h=epsrel;
263 
264  xx[j]=x[j]+h;
265  int ret=(this->func)(nx,xx,f);
266  xx[j]=x[j];
267 
268  // The function returned a non-zero value, so try a different step
269  size_t it=0;
270  while (ret!=0 && h>=epsmin && it<max_shrink_iters) {
271 
272  // First try flipping the sign
273  h=-h;
274  xx[j]=x[j]+h;
275  ret=(this->func)(nx,xx,f);
276  xx[j]=x[j];
277 
278  if (ret!=0) {
279 
280  // If that didn't work, flip to positive and try a smaller
281  // stepsize
282  h/=-shrink_fact;
283  if (h>=epsmin) {
284  xx[j]=x[j]+h;
285  ret=(this->func)(nx,xx,f);
286  xx[j]=x[j];
287  }
288 
289  }
290 
291  it++;
292  }
293 
294  if (ret!=0) {
295  O2SCL_CONV2_RET("Jacobian failed to find valid step in ",
296  "jacobian_gsl::operator().",exc_ebadfunc,
297  this->err_nonconv);
298  }
299 
300  // This is the equivalent of GSL's test of
301  // gsl_vector_isnull(&col.vector)
302 
303  bool nonzero=false;
304  for (i=0;i<ny;i++) {
305  temp=(f[i]-y[i])/h;
306  if (temp!=0.0) nonzero=true;
307  jac(i,j)=temp;
308  }
309  if (nonzero==false) success=false;
310 
311  }
312 
313  if (success==false) {
314  O2SCL_CONV2_RET("At least one row of the Jacobian is zero ",
315  "in jacobian_gsl::operator().",exc_esing,
316  this->err_nonconv);
317  }
318  return 0;
319  }
320 
321  };
322 
323  /** \brief A direct calculation of the jacobian using a \ref
324  deriv_base object
325 
326  By default, the stepsize, \ref deriv_gsl::h is set to \f$
327  10^{-4} \f$ in the \ref jacobian_exact constructor.
328 
329  Note that it is most often wasteful to use this Jacobian in a
330  root-finding routine and using more approximate Jacobians is
331  more efficient.
332 
333  Default template arguments
334  - \c func_t - \ref mm_funct11
335  - \c vec_t - boost::numeric::ublas::vector<double>
336  - \c mat_t - boost::numeric::ublas::matrix<double>
337  */
338  template<class func_t=mm_funct11,
341  public jacobian<func_t,vec_t,mat_t> {
342 
343  public:
344 
345  jacobian_exact() {
346  def_deriv.h=1.0e-4;
347  dptr=&def_deriv;
348  }
349 
350  /** \brief Parameter structure for passing information
351 
352  This class is primarily useful for specifying derivatives
353  for using the jacobian::set_deriv() function.
354 
355  \comment
356  This type needs to be publicly available so that the
357  user can properly specify a base 1-dimensional derivative
358  object.
359  \endcomment
360  */
361  typedef struct {
362  /// The number of variables
363  size_t nx;
364  /// The number of variables
365  size_t ny;
366  /// The current x value
367  size_t xj;
368  /// The current y value
369  size_t yi;
370  /// The x vector
371  vec_t *x;
372  /// The y vector
373  vec_t *y;
374  } ej_parms;
375 
376  /// The default derivative object
378 
379  /// Set the derivative object
381  dptr=&de;
382  return 0;
383  }
384 
385  /** \brief The operator()
386  */
387  virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y,
388  mat_t &jac) {
389 
390  // The function return value
391  int ret=0;
392 
393  // Temporary storage fo the derivative uncertainty
394  double err;
395 
396  ej_parms ejp;
397  ejp.nx=nx;
398  ejp.ny=ny;
399  ejp.x=&x;
400  ejp.y=&y;
401 
402  funct11 dfnp=std::bind(std::mem_fn<double(double,ej_parms &)>
404  this,std::placeholders::_1,std::ref(ejp));
405 
406  for (size_t j=0;j<nx;j++) {
407  ejp.xj=j;
408  for (size_t i=0;i<ny;i++) {
409  ejp.yi=i;
410  double tmp=(*ejp.x)[j];
411  int dret=dptr->deriv_err(tmp,dfnp,jac(i,j),err);
412  if (dret!=0) {
413  if (this->err_nonconv==true) {
414  O2SCL_ERR2("Derivative object tailed in jacobian_exact::",
415  "operator().",o2scl::exc_efailed);
416  }
417  ret=1;
418  }
419  (*ejp.x)[j]=tmp;
420  }
421  }
422 
423  return ret;
424  }
425 
426  /** \brief Compute the Jacobian and its uncertainty
427  from the numerical differentiation
428  */
429  virtual int jac_err(size_t nx, vec_t &x, size_t ny, vec_t &y,
430  mat_t &jac, mat_t &err) {
431 
432  // The function return value
433  int ret=0;
434 
435  ej_parms ejp;
436  ejp.nx=nx;
437  ejp.ny=ny;
438  ejp.x=&x;
439  ejp.y=&y;
440 
441  funct11 dfnp=std::bind(std::mem_fn<double(double,ej_parms &)>
443  this,std::placeholders::_1,std::ref(ejp));
444 
445  for (size_t j=0;j<nx;j++) {
446  ejp.xj=j;
447  for (size_t i=0;i<ny;i++) {
448  ejp.yi=i;
449  double tmp=(*ejp.x)[j];
450  int dret=dptr->deriv_err(tmp,dfnp,jac(i,j),err(i,j));
451  if (dret!=0) {
452  if (this->err_nonconv==true) {
453  O2SCL_ERR2("Derivative object tailed in jacobian_exact::",
454  "jac_err().",o2scl::exc_efailed);
455  }
456  ret=1;
457  }
458  (*ejp.x)[j]=tmp;
459  }
460  }
461 
462  return ret;
463  }
464 
465 #ifndef DOXYGEN_INTERNAL
466 
467  protected:
468 
469  /// Pointer to the derivative object
471 
472  /// Function for the derivative object
473  double dfn(double x, ej_parms &ejp) {
474  (*ejp.x)[ejp.xj]=x;
475  (this->func)(ejp.nx,*ejp.x,*ejp.y);
476  return (*ejp.y)[ejp.yi];
477  }
478 
479 #endif
480 
481  };
482 
483 #ifndef DOXYGEN_NO_O2NS
484 }
485 #endif
486 
487 #endif
vec_t f
Function values.
Definition: jacobian.h:151
virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &j)=0
Evaluate the Jacobian j at point y(x)
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
virtual int set_function(func_t &f)
Set the function to compute the Jacobian of.
Definition: jacobian.h:78
size_t xj
The current x value.
Definition: jacobian.h:367
vec_t * y
The y vector.
Definition: jacobian.h:373
std::function< double(double)> funct11
One-dimensional function typedef.
Definition: funct.h:44
virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &jac)
The operator()
Definition: jacobian.h:387
invalid argument supplied by user
Definition: err_hnd.h:59
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct11
Array of multi-dimensional functions typedef.
Definition: mm_funct.h:43
apparent singularity detected
Definition: err_hnd.h:93
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
size_t mem_size_x
Size of allocated memory in x.
Definition: jacobian.h:157
double get_epsrel()
Get the relative stepsize (default )
Definition: jacobian.h:193
size_t yi
The current y value.
Definition: jacobian.h:369
generic failure
Definition: err_hnd.h:61
double epsrel
The relative stepsize for finite-differencing.
Definition: jacobian.h:164
double dfn(double x, ej_parms &ejp)
Function for the derivative object.
Definition: jacobian.h:473
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &jac)
The operator()
Definition: jacobian.h:241
bool err_nonconv
If true, call the error handler if the routine does not converge.
Definition: jacobian.h:72
void set_epsmin(double l_epsmin)
Set the minimum stepsize (must be )
Definition: jacobian.h:212
Numerical differentiation base [abstract base].
Definition: deriv.h:63
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
size_t ny
The number of variables.
Definition: jacobian.h:365
double get_epsmin()
Get the minimum stepsize (default )
Definition: jacobian.h:197
std::function< int(size_t, boost::numeric::ublas::vector< double > &, size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &) > jac_funct11
Jacobian function (not necessarily square)
Definition: jacobian.h:44
vec_t xx
Function arguments.
Definition: jacobian.h:154
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void set_epsrel(double l_epsrel)
Set the relative stepsize (must be )
Definition: jacobian.h:201
deriv_base * dptr
Pointer to the derivative object.
Definition: jacobian.h:470
size_t max_shrink_iters
Maximum number of times to shrink the step size.
Definition: jacobian.h:170
Numerical differentiation (GSL)
Definition: deriv_gsl.h:123
int set_deriv(deriv_base<> &de)
Set the derivative object.
Definition: jacobian.h:380
void set_shrink_fact(double l_shrink_fact)
Set shrink factor for decreasing step size.
Definition: jacobian.h:223
void set_max_shrink_iters(size_t it)
Set number of times to decrease step size.
Definition: jacobian.h:234
func_t func
A pointer to the user-specified function.
Definition: jacobian.h:93
problem with user-supplied function
Definition: err_hnd.h:69
size_t nx
The number of variables.
Definition: jacobian.h:363
double epsmin
The minimum stepsize.
Definition: jacobian.h:167
deriv_gsl def_deriv
The default derivative object.
Definition: jacobian.h:377
double shrink_fact
Factor to shrink stepsize by.
Definition: jacobian.h:173
virtual int jac_err(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &jac, mat_t &err)
Compute the Jacobian and its uncertainty from the numerical differentiation.
Definition: jacobian.h:429
size_t mem_size_y
Size of allocated memory in y.
Definition: jacobian.h:160
vec_t * x
The x vector.
Definition: jacobian.h:371
Parameter structure for passing information.
Definition: jacobian.h:361
Success.
Definition: err_hnd.h:47
Simple automatic Jacobian.
Definition: jacobian.h:144
Base for providing a numerical jacobian [abstract base].
Definition: jacobian.h:64
A direct calculation of the jacobian using a deriv_base object.
Definition: jacobian.h:340

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