jacobian.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_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) in src/root/jacobian.h
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_funct).
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_funct
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_funct,
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_funct
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_funct,
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 
247  if (mem_size_x!=nx || mem_size_y!=ny) {
248  f.resize(ny);
249  xx.resize(nx);
250  mem_size_x=nx;
251  mem_size_y=ny;
252  }
253 
254  vector_copy(nx,x,xx);
255 
256  for (j=0;j<nx;j++) {
257 
258  // Thanks to suggestion from Conrad Curry.
259  h=epsrel*fabs(x[j]);
260  if (h<epsmin) h=epsmin;
261  if (h==0.0) h=epsrel;
262 
263  xx[j]=x[j]+h;
264  int ret=(this->func)(nx,xx,f);
265  xx[j]=x[j];
266 
267  // The function returned a non-zero value, so try a different step
268  size_t it=0;
269  while (ret!=0 && h>=epsmin && it<max_shrink_iters) {
270 
271  // First try flipping the sign
272  h=-h;
273  xx[j]=x[j]+h;
274  ret=(this->func)(nx,xx,f);
275  xx[j]=x[j];
276 
277  if (ret!=0) {
278 
279  // If that didn't work, flip to positive and try a smaller
280  // stepsize
281  h/=-shrink_fact;
282  if (h>=epsmin) {
283  xx[j]=x[j]+h;
284  ret=(this->func)(nx,xx,f);
285  xx[j]=x[j];
286  }
287 
288  }
289 
290  it++;
291  }
292 
293  if (ret!=0) {
294  O2SCL_CONV2_RET("Jacobian failed to find valid step in ",
295  "jacobian_gsl::operator().",exc_ebadfunc,
296  this->err_nonconv);
297  }
298 
299  // This is the equivalent of GSL's test of
300  // gsl_vector_isnull(&col.vector)
301 
302  bool nonzero=false;
303  for (i=0;i<ny;i++) {
304  temp=(f[i]-y[i])/h;
305  if (temp!=0.0) nonzero=true;
306  jac(i,j)=temp;
307  }
308  if (nonzero==false) {
309  O2SCL_CONV_RET((((std::string)"Row ")+o2scl::szttos(j)+
310  " of the Jacobian is zero "+
311  "in jacobian_gsl::operator().").c_str(),exc_esing,
312  this->err_nonconv);
313  }
314 
315  }
316 
317  return 0;
318  }
319 
320  };
321 
322  /** \brief A direct calculation of the jacobian using a \ref
323  deriv_base object
324 
325  By default, the stepsize, \ref deriv_gsl::h is set to \f$
326  10^{-4} \f$ in the \ref jacobian_exact constructor.
327 
328  Note that it is most often wasteful to use this Jacobian in a
329  root-finding routine and using more approximate Jacobians is
330  more efficient.
331 
332  Default template arguments
333  - \c func_t - \ref mm_funct
334  - \c vec_t - boost::numeric::ublas::vector<double>
335  - \c mat_t - boost::numeric::ublas::matrix<double>
336  */
337  template<class func_t=mm_funct,
340  public jacobian<func_t,vec_t,mat_t> {
341 
342  public:
343 
344  jacobian_exact() {
345  def_deriv.h=1.0e-4;
346  dptr=&def_deriv;
347  }
348 
349  /** \brief Parameter structure for passing information
350 
351  This class is primarily useful for specifying derivatives
352  for using the jacobian::set_deriv() function.
353 
354  \comment
355  This type needs to be publicly available so that the
356  user can properly specify a base 1-dimensional derivative
357  object.
358  \endcomment
359  */
360  typedef struct {
361  /// The number of variables
362  size_t nx;
363  /// The number of variables
364  size_t ny;
365  /// The current x value
366  size_t xj;
367  /// The current y value
368  size_t yi;
369  /// The x vector
370  vec_t *x;
371  /// The y vector
372  vec_t *y;
373  } ej_parms;
374 
375  /// The default derivative object
377 
378  /// Set the derivative object
380  dptr=&de;
381  return 0;
382  }
383 
384  /** \brief The operator()
385  */
386  virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y,
387  mat_t &jac) {
388 
389  // The function return value
390  int ret=0;
391 
392  // Temporary storage fo the derivative uncertainty
393  double err;
394 
395  ej_parms ejp;
396  ejp.nx=nx;
397  ejp.ny=ny;
398  ejp.x=&x;
399  ejp.y=&y;
400 
401  funct dfnp=std::bind(std::mem_fn<double(double,ej_parms &)>
403  this,std::placeholders::_1,std::ref(ejp));
404 
405  for (size_t j=0;j<nx;j++) {
406  ejp.xj=j;
407  for (size_t i=0;i<ny;i++) {
408  ejp.yi=i;
409  double tmp=(*ejp.x)[j];
410  int dret=dptr->deriv_err(tmp,dfnp,jac(i,j),err);
411  if (dret!=0) {
412  if (this->err_nonconv==true) {
413  O2SCL_ERR2("Derivative object tailed in jacobian_exact::",
414  "operator().",o2scl::exc_efailed);
415  }
416  ret=1;
417  }
418  (*ejp.x)[j]=tmp;
419  }
420  }
421 
422  return ret;
423  }
424 
425  /** \brief Compute the Jacobian and its uncertainty
426  from the numerical differentiation
427  */
428  virtual int jac_err(size_t nx, vec_t &x, size_t ny, vec_t &y,
429  mat_t &jac, mat_t &err) {
430 
431  // The function return value
432  int ret=0;
433 
434  ej_parms ejp;
435  ejp.nx=nx;
436  ejp.ny=ny;
437  ejp.x=&x;
438  ejp.y=&y;
439 
440  funct dfnp=std::bind(std::mem_fn<double(double,ej_parms &)>
442  this,std::placeholders::_1,std::ref(ejp));
443 
444  for (size_t j=0;j<nx;j++) {
445  ejp.xj=j;
446  for (size_t i=0;i<ny;i++) {
447  ejp.yi=i;
448  double tmp=(*ejp.x)[j];
449  int dret=dptr->deriv_err(tmp,dfnp,jac(i,j),err(i,j));
450  if (dret!=0) {
451  if (this->err_nonconv==true) {
452  O2SCL_ERR2("Derivative object tailed in jacobian_exact::",
453  "jac_err().",o2scl::exc_efailed);
454  }
455  ret=1;
456  }
457  (*ejp.x)[j]=tmp;
458  }
459  }
460 
461  return ret;
462  }
463 
464 #ifndef DOXYGEN_INTERNAL
465 
466  protected:
467 
468  /// Pointer to the derivative object
470 
471  /// Function for the derivative object
472  double dfn(double x, ej_parms &ejp) {
473  (*ejp.x)[ejp.xj]=x;
474  (this->func)(ejp.nx,*ejp.x,*ejp.y);
475  return (*ejp.y)[ejp.yi];
476  }
477 
478 #endif
479 
480  };
481 
482 #ifndef DOXYGEN_NO_O2NS
483 }
484 #endif
485 
486 #endif
boost::numeric::ublas::matrix< double >
o2scl::exc_esing
@ exc_esing
apparent singularity detected
Definition: err_hnd.h:93
o2scl::jacobian_gsl::mem_size_x
size_t mem_size_x
Size of allocated memory in x.
Definition: jacobian.h:157
o2scl::jacobian_exact
A direct calculation of the jacobian using a deriv_base object.
Definition: jacobian.h:339
O2SCL_CONV2_RET
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:303
boost::numeric::ublas::vector< double >
o2scl::jacobian_exact::ej_parms::nx
size_t nx
The number of variables.
Definition: jacobian.h:362
o2scl::jacobian::func
func_t func
A pointer to the user-specified function.
Definition: jacobian.h:93
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::jacobian_exact::ej_parms::yi
size_t yi
The current y value.
Definition: jacobian.h:368
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::jacobian_exact::ej_parms::x
vec_t * x
The x vector.
Definition: jacobian.h:370
o2scl::deriv_base
Numerical differentiation base [abstract base].
Definition: deriv.h:63
o2scl::jacobian_exact::dptr
deriv_base * dptr
Pointer to the derivative object.
Definition: jacobian.h:469
o2scl::jacobian_gsl::xx
vec_t xx
Function arguments.
Definition: jacobian.h:154
o2scl::jacobian_exact::ej_parms::xj
size_t xj
The current x value.
Definition: jacobian.h:366
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::jacobian_gsl::set_max_shrink_iters
void set_max_shrink_iters(size_t it)
Set number of times to decrease step size.
Definition: jacobian.h:234
o2scl::jacobian_exact::ej_parms::y
vec_t * y
The y vector.
Definition: jacobian.h:372
o2scl::deriv_base::deriv_err
virtual int deriv_err(fp_t x, func_t &func, fp_t &dfdx, fp_t &err)=0
Calculate the first derivative of func w.r.t. x and the uncertainty.
o2scl::jacobian_gsl::max_shrink_iters
size_t max_shrink_iters
Maximum number of times to shrink the step size.
Definition: jacobian.h:170
o2scl::jacobian_gsl::set_epsrel
void set_epsrel(double l_epsrel)
Set the relative stepsize (must be )
Definition: jacobian.h:201
o2scl::jacobian_exact::set_deriv
int set_deriv(deriv_base<> &de)
Set the derivative object.
Definition: jacobian.h:379
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::jacobian_gsl::shrink_fact
double shrink_fact
Factor to shrink stepsize by.
Definition: jacobian.h:173
o2scl::jacobian_gsl::get_epsmin
double get_epsmin()
Get the minimum stepsize (default )
Definition: jacobian.h:197
o2scl::jacobian::set_function
virtual int set_function(func_t &f)
Set the function to compute the Jacobian of.
Definition: jacobian.h:78
O2SCL_CONV_RET
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:297
o2scl::jacobian_exact::ej_parms::ny
size_t ny
The number of variables.
Definition: jacobian.h:364
o2scl::jacobian
Base for providing a numerical jacobian [abstract base].
Definition: jacobian.h:64
o2scl::jacobian_gsl
Simple automatic Jacobian.
Definition: jacobian.h:144
o2scl::jacobian_gsl::operator()
virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &jac)
The operator()
Definition: jacobian.h:241
o2scl::jacobian_gsl::epsrel
double epsrel
The relative stepsize for finite-differencing.
Definition: jacobian.h:164
o2scl::mm_funct
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef in src/base/mm_funct.h.
Definition: mm_funct.h:43
o2scl::deriv_gsl::h
fp_t h
Initial stepsize.
Definition: deriv_gsl.h:144
o2scl::deriv_gsl
Numerical differentiation (GSL)
Definition: deriv_gsl.h:124
o2scl::jacobian_exact::dfn
double dfn(double x, ej_parms &ejp)
Function for the derivative object.
Definition: jacobian.h:472
o2scl::jac_funct
std::function< int(size_t, boost::numeric::ublas::vector< double > &, size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &) > jac_funct
Jacobian function (not necessarily square) in src/root/jacobian.h.
Definition: jacobian.h:44
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::funct
std::function< double(double)> funct
One-dimensional function typedef in src/base/funct.h.
Definition: funct.h:48
o2scl::jacobian::err_nonconv
bool err_nonconv
If true, call the error handler if the routine does not converge.
Definition: jacobian.h:72
o2scl::jacobian_exact::def_deriv
deriv_gsl def_deriv
The default derivative object.
Definition: jacobian.h:376
o2scl::jacobian_exact::operator()
virtual int operator()(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &jac)
The operator()
Definition: jacobian.h:386
o2scl::jacobian_gsl::get_epsrel
double get_epsrel()
Get the relative stepsize (default )
Definition: jacobian.h:193
o2scl::jacobian::operator()
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)
o2scl::jacobian_gsl::set_epsmin
void set_epsmin(double l_epsmin)
Set the minimum stepsize (must be )
Definition: jacobian.h:212
o2scl::jacobian_gsl::f
vec_t f
Function values.
Definition: jacobian.h:151
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::jacobian_gsl::epsmin
double epsmin
The minimum stepsize.
Definition: jacobian.h:167
o2scl::jacobian_exact::ej_parms
Parameter structure for passing information.
Definition: jacobian.h:360
o2scl::jacobian_gsl::mem_size_y
size_t mem_size_y
Size of allocated memory in y.
Definition: jacobian.h:160
o2scl::jacobian_exact::jac_err
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:428
o2scl::jacobian_gsl::set_shrink_fact
void set_shrink_fact(double l_shrink_fact)
Set shrink factor for decreasing step size.
Definition: jacobian.h:223
o2scl::exc_ebadfunc
@ exc_ebadfunc
problem with user-supplied function
Definition: err_hnd.h:69

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