mmin_simp2.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 /* multimin/simplex2.c
24  *
25  * Copyright (C) 2007 Brian Gough
26  * Copyright (C) 2002 Tuomo Keskitalo, Ivo Alxneit
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 #ifndef O2SCL_MMIN_SIMP2_H
44 #define O2SCL_MMIN_SIMP2_H
45 
46 /** \file mmin_simp2.h
47  \brief File defining \ref o2scl::mmin_simp2
48 */
49 #include <string>
50 
51 #include <boost/numeric/ublas/vector.hpp>
52 
53 #include <gsl/gsl_multimin.h>
54 
55 #include <o2scl/mmin.h>
56 #include <o2scl/cblas.h>
57 #include <o2scl/vec_stats.h>
58 
59 #ifndef DOXYGEN_NO_O2NS
60 namespace o2scl {
61 #endif
62 
63  /** \brief Multidimensional minimization by the Simplex method (v2) (GSL)
64 
65  This class mins a function using Nelder and Mead's Simplex
66  algorithm. A simplex in a N-dimensional space is defined as a
67  set of N+1 points which describe an N-dimensional volume
68  surrounding the minimum. The algorithm proceeds by shifting the
69  simplex points until the simplex is sufficiently small and thus
70  the minimum is known with sufficient accuracy.
71 
72  This class has a high-level interface using mmin(),
73  mmin_twovec() or mmin_simplex() which automatically performs the
74  memory allocation and minimization, or a GSL-like interface
75  using allocate(), free(), interate() and set() or set_simplex().
76 
77  The simplex can be completely specified by the user (see
78  mmin_simplex() and set_simplex()). Alternatively, the simplex is
79  automatically specified given initial guess \f$ x_j \f$ and a
80  step size vector \f$ s_k \f$ for \f$ 0\leq k<n_s \f$. The
81  simplex \f$ p_{ij} \f$ with \f$ 0\leq i\leq n \f$ and \f$ 0\leq
82  j<n \f$ is chosen with \f$ p_{0j} = x_j \f$ and
83  \f{eqnarray*}
84  p_{i+1,j} &=& x_j\quad\mathrm{for}\quad i\neq j \\
85  p_{i+1,j} &=& x_j+s_{j~\mathrm{mod}~n_s}\quad\mathrm{for}\quad i=j
86  \f}
87  for \f$ 0<i<n \f$. The step size vector \f$ s \f$ is set by the
88  set_step() member function. The prsence of \f$ \mathrm{mod} \f$
89  in the recipe above just indicates that elements of the step
90  size vector are automatically re-used if there are less step
91  sizes than dimensions in the minimization problem.
92 
93  \note It is important that the initial simplex contains
94  sufficient variation in every direction of the parameter space
95  over which one is minimizing. For example, if all three points
96  in a simplex for minimizing over a two-dimensional space contain
97  nearly the same value for the second parameter, then the
98  minimizer may only minimize the function with respect to the
99  first parameter.
100 
101  \note The algorithm used to estimate the simplex size does not
102  work well any of the parameters in the minimization problem has
103  a scale which is not close to 1.
104 
105  See an example for the usage of this class
106  in \ref ex_mmin_sect .
107 
108  Default template arguments
109  - \c param_t - no default
110  - \c func_t - \ref multi_funct11
111  - \c vec_t - \ref boost::numeric::ublas::vector < double >
112 
113  Based on \ref Nelder65 .
114 
115  A variable <tt>count</tt> originally defined in the GSL simplex
116  state is not present here, because it was unused.
117 
118  \future Double check that the updates in gsl-1.13 are included
119  here, and also add support for the nmsimplex2rand algorithm
120  in GSL.
121  */
122  template<class func_t=multi_funct11,
124  public mmin_base<func_t,func_t,vec_t> {
125 
126  public:
127 
129 
130 #ifndef DOXYGEN_INTERNAL
131 
132  protected:
133 
134  /// An array of n+1 vectors containing the simplex
135  vec_t *x1;
136  /** \brief The n+1 function values at the simplex points
137 
138  \comment
139  This can't be the template type because it has to be
140  a different size (n+1) rather than n.
141  \endcomment
142  */
143  ubvector y1;
144  /// Workspace vector 1
145  vec_t ws1;
146  /// Workspace vector 2
147  vec_t ws2;
148  /// Workspace vector 3
149  vec_t ws3;
150  /// Center of simplex
151  vec_t center;
152  /// Desc
153  vec_t delta;
154  /// Distance of vector from center
155  vec_t xmc;
156  /// Squared simplex size
157  double S2;
158 
159  /// Compute the center of the simplex
161  size_t i, j;
162  double val;
163 
164  for (j = 0; j < dim; j++) {
165  val = 0.0;
166  for (i = 0; i < dim+1; i++) {
167  val += x1[i][j];
168  }
169  val /= ((double)dim)+1;
170  center[j]=val;
171  }
172 
173  return 0;
174  }
175 
176  /** \brief Compute the size of the simplex
177 
178  Calculates simplex size as average sum of length of vectors
179  from simplex center to corner points:
180 
181  \f$ (1/n) \sum || y - y_{\mathrm{center}} || \f$
182  */
183  double compute_size() {
184 
185  size_t P=dim+1;
186  double ss=0.0;
187 
188  for(size_t i=0;i<P;i++) {
189  for(size_t j=0;j<dim;j++) ws1[j]=x1[i][j];
190  o2scl_cblas::daxpy(-1.0,dim,center,ws1);
191  double t=o2scl_cblas::dnrm2(dim,ws1);
192  ss += t*t;
193  }
194 
195  /* Store squared size in the state */
196  S2=ss/((double)P);
197  return sqrt(ss/(double)(P));
198  }
199 
200  /** \brief Move a corner of a simplex
201 
202  Moves a simplex corner scaled by coeff (negative value
203  represents mirroring by the middle point of the "other" corner
204  points) and gives new corner in xc and function value at xc
205  in \c newval.
206  */
207  virtual int try_corner_move(const double coeff, size_t corner,
208  vec_t &xc, func_t &f, size_t nvar,
209  double &newval) {
210 
211  size_t P=dim+1;
212 
213  /*
214  GSL: xc = (1-coeff)*(P/(P-1)) * center(all) +
215  ((P*coeff-1)/(P-1))*x_corner
216  */
217  double alpha=(1.0-coeff)*((double)P)/((double)dim);
218  double beta=(((double)P)*coeff-1.0)/((double)dim);
219 
220  for(size_t j=0;j<dim;j++) {
221  xc[j]=center[j]*alpha;
222  xc[j]+=x1[corner][j]*beta;
223  }
224 
225  newval=f(nvar,xc);
226 
227  return 0;
228  }
229 
230  /// Update point \c i in the simplex with values \c xx
231  virtual int update_point(size_t i, vec_t &xx, double val) {
232 
233  const size_t P=dim+1;
234 
235  /* GSL: Compute delta = x - x_orig */
236  for(size_t j=0;j<dim;j++) {
237  delta[j]=xx[j];
238  delta[j]-=x1[i][j];
239  }
240 
241  /* GSL: Compute xmc = x_orig - c */
242  for(size_t j=0;j<dim;j++) {
243  xmc[j]=x1[i][j];
244  xmc[j]-=center[j];
245  }
246 
247  /* GSL: Update size: S2' = S2 + (2/P) *
248  (x_orig - c).delta + (P-1)*(delta/P)^2
249  */
250  double d=o2scl_cblas::dnrm2(dim,delta);
251  double xmcd=o2scl_cblas::ddot(dim,xmc,delta);
252  S2 += (2.0 / ((double)P)) * xmcd +
253  ((((double)P) - 1.0) / ((double)P)) * (d * d / ((double)P));
254 
255  /* GSL: Update center: c' = c + (x - x_orig) / P */
256  double alpha=1.0/((double)P);
257  for(size_t j=0;j<dim;j++) {
258  center[j]-=alpha*x1[i][j];
259  center[j]+=alpha*xx[j];
260  }
261 
262  for(size_t j=0;j<dim;j++) {
263  x1[i][j]=xx[j];
264  }
265  y1[i]=val;
266 
267  return 0;
268  }
269 
270  /** \brief Contract the simplex towards the best point
271 
272  Function contracts the simplex in respect to best valued
273  corner. All corners besides the best corner are moved.
274 
275  \comment
276  The GSL version requires a vector for workspace
277  named 'xc' which is not required here.
278  \endcomment
279  */
280  virtual int contract_by_best(size_t best, func_t &f,
281  size_t nvar) {
282 
283  /* GSL: Function contracts the simplex in respect to best valued
284  corner. That is, all corners besides the best corner are
285  moved. (This function is rarely called in practice, since
286  it is the last choice, hence not optimised - BJG)
287  */
288  size_t i,j, it;
289  double newval;
290  bool failed;
291 
292  int status=success;
293 
294  for (i=0;i<dim+1;i++) {
295 
296  if (i!=best) {
297 
298  for (j=0;j<dim;j++) {
299  x1[i][j]=0.5*(x1[i][j]+x1[best][j]);
300  }
301  y1[i]=f(nvar,x1[i]);
302  if (!std::isfinite(y1[i])) {
303  std::string err=((std::string)"Function not finite (returned ")+
304  dtos(y1[i])+" in mmin_simp2::contract_by_best().";
305  O2SCL_ERR(err.c_str(),exc_ebadfunc);
306  }
307 
308  /*
309  if (avoid_nonzero==true && ret!=0) {
310  std::cout << "Found problem in contract." << std::endl;
311 
312  // copy the old point to ws3
313  for (j=0;j<dim;j++) {
314  ws3[j]=2.0*x1[i][j]-x1[best][j];
315  }
316 
317  // Try (21*best+20*xold)/41, (22*best+19*xold)/41, ...,
318  // (40*best+xold)/41 to see if we can find a contraction
319  // that works
320  for(it=0;ret!=0 && it<20;it++) {
321  for (j=0;j<dim;j++) {
322  x1[i][j]=((20-it)*ws3[j]+(it+21)*x1[best][j])/41.0;
323  }
324  y1[i]=f(nvar,x1[i]);
325  std::cout << "it, ret, x: " << it << " " << ret << " "
326  << x << std::endl;
327  }
328  if (ret!=0) {
329  O2SCL_CONV2_RET("Failed to find suitable contraction ",
330  "in mmin_simp2::contract_by_best().",
331  exc_efailed,this->err_nonconv);
332  }
333  }
334  */
335  }
336  }
337 
338  compute_center();
339  size=compute_size();
340 
341  return success;
342  }
343 
344  /// Number of variables to be mind over
345  size_t dim;
346 
347  /// Function
348  func_t *func;
349 
350  /// True if set() has been called
352 
353  /// Vector of step sizes
354  ubvector step_vec;
355 
356  /** \brief If true, try to automatically avoid regions where the
357  function returns a non-zero value (default false)
358 
359  \note This option doesn't work yet, so I've made the variable
360  protected to prevent the casual user from changing it.
361  */
363 
364 #endif
365 
366  public:
367 
368  mmin_simp2() {
369  dim=0;
370  print_simplex=0;
371  step_vec.resize(1);
372  step_vec[0]=1.0;
373  avoid_nonzero=false;
374  }
375 
376  virtual ~mmin_simp2() {
377  }
378 
379  /// Set the step sizes for each independent variable
380  template<class vec2_t> int set_step(size_t nv, vec2_t &step) {
381  if (nv>0) {
382  step_vec.resize(nv);
383  for(size_t i=0;i<nv;i++) step_vec[i]=step[i];
384  }
385  return 0;
386  }
387 
388  /// Size of current simplex computed by iterate()
389  double size;
390 
391  /// Present minimum vector computed by iterate()
392  vec_t x;
393 
394  /// Function value at minimum computed by iterate()
395  double fval;
396 
397  /** \brief Print simplex information in print_iter() (default 0)
398 
399  If this is 1 and \ref verbose is greater than 0,
400  then print_iter() will print the function values at all
401  the simplex points. If this is 2 and \ref verbose is
402  greater than 0, then print_iter() will print the
403  simplex coordinates in addition to the function values.
404  */
406 
407  /** \brief Calculate the minimum \c min of \c func w.r.t the
408  array \c x of size \c nvar.
409  */
410  virtual int mmin(size_t nn, vec_t &xx, double &fmin,
411  func_t &ufunc) {
412 
413  if (nn==0) {
414  O2SCL_ERR2("Tried to min over zero variables ",
415  " in mmin_simp2::mmin().",exc_einval);
416  }
417 
418  int ret=0,status,iter=0;
419 
420  allocate(nn);
421 
422  vec_t ss(nn);
423  for (size_t is=0;is<nn;is++) ss[is]=step_vec[is % step_vec.size()];
424  ret=set(ufunc,nn,xx,ss);
425 
426  if(ret!=0) {
427  return ret;
428  }
429 
430  do {
431  iter++;
432 
433  status=iterate();
434  if(status) break;
435 
436  if(this->verbose>0) {
437  print_iter(nn,x,x1,fval,iter,size,this->tol_abs,
438  "mmin_simp2");
439  }
440 
441  status=gsl_multimin_test_size(size,this->tol_abs);
442 
443  } while(status == GSL_CONTINUE && iter<this->ntrial);
444 
445  for (size_t i=0;i<nn;i++) xx[i]=x[i];
446  fmin=fval;
447 
448  this->last_ntrial=iter;
449 
450  if(iter>=this->ntrial) {
451  std::string str="Exceeded maximum number of iterations ("+
452  itos(this->ntrial)+") in mmin_simp2::mmin().";
453  O2SCL_CONV_RET(str.c_str(),exc_emaxiter,this->err_nonconv);
454  }
455 
456  return status;
457  }
458 
459  /** \brief Calculate the minimum \c min of \c func w.r.t the
460  array \c x of size \c nvar, using \c xx and \c xx2 to specify
461  the simplex
462  */
463  virtual int mmin_twovec(size_t nn, vec_t &xx, vec_t &xx2, double &fmin,
464  func_t &ufunc) {
465 
466  int ret=0,i,status,iter=0;
467 
468  allocate(nn);
469 
470  vec_t ss(nn);
471  for (size_t is=0;is<nn;is++) ss[is]=xx2[is]-xx[is];
472  ret=set(ufunc,nn,xx,ss);
473 
474  if(ret!=0) {
475  return ret;
476  }
477 
478  do {
479  iter++;
480 
481  status=iterate();
482  if(status) break;
483 
484  if(this->verbose>0) {
485  print_iter(nn,x,x1,fval,iter,size,this->tol_abs,
486  "mmin_simp2");
487  }
488 
489  status=gsl_multimin_test_size(size,this->tol_abs);
490 
491  } while(status == GSL_CONTINUE && iter<this->ntrial);
492 
493  for (i=0;i<((int)nn);i++) xx[i]=x[i];
494  fmin=fval;
495 
496  this->last_ntrial=iter;
497 
498  if(iter>=this->ntrial) {
499  std::string str="Exceeded maximum number of iterations ("+
500  itos(this->ntrial)+") in mmin_simp2::mmin_twovec().";
501  O2SCL_CONV_RET(str.c_str(),exc_emaxiter,this->err_nonconv);
502  }
503 
504  return status;
505  }
506 
507  /** \brief Calculate the minimum \c min of \c func w.r.t the
508  array \c x of size \c nvar, given an initial simplex
509  */
510  template<class mat_t>
511  int mmin_simplex(size_t nn, mat_t &sx, double &fmin,
512  func_t &ufunc) {
513 
514  int ret=0,i,status,iter=0;
515 
516  allocate(nn);
517 
518  ret=set_simplex(ufunc,sx);
519  if(ret!=0) {
520  return ret;
521  }
522 
523  do {
524  iter++;
525 
526  status=iterate();
527  if(status) break;
528 
529  if(this->verbose>0) {
530  print_iter(nn,x,x1,fval,iter,size,this->tol_abs,
531  "mmin_simp2");
532  }
533 
534  status=gsl_multimin_test_size(size,this->tol_abs);
535 
536  } while(status == GSL_CONTINUE && iter<this->ntrial);
537 
538  for (i=0;i<((int)nn);i++) sx(0,i)=x[i];
539  fmin=fval;
540 
541  this->last_ntrial=iter;
542 
543  if (iter>=this->ntrial) {
544  std::string str="Exceeded maximum number of iterations ("+
545  itos(this->ntrial)+") in mmin_simp2::mmin_simplex().";
546  O2SCL_CONV_RET(str.c_str(),exc_emaxiter,this->err_nonconv);
547  }
548 
549  return status;
550  }
551 
552  /// Allocate the memory
553  virtual int allocate(size_t n) {
554 
555  set_called=false;
556  dim=n;
557 
558  x.resize(n);
559  x1=new vec_t[n+1];
560  for(size_t i=0;i<n+1;i++) x1[i].resize(n);
561  y1.resize(n+1);
562  ws1.resize(n);
563  ws2.resize(n);
564  ws3.resize(n);
565  center.resize(n);
566  delta.resize(n);
567  xmc.resize(n);
568 
569  return success;
570  }
571 
572  /// Set the function and initial guess
573  virtual int set(func_t &ufunc, size_t n, vec_t &ax,
574  vec_t &step_size) {
575  size_t i;
576 
577  if(dim!=n) allocate(n);
578  func=&ufunc;
579 
580  // Copy initial guess to x
581  for (i=0;i<dim;i++) x[i]=ax[i];
582 
583  // first point is the original x0
584 
585  y1[0]=ufunc(dim,ax);
586  if (!std::isfinite(y1[0])) {
587  std::string err=((std::string)"Function not finite (returned ")+
588  dtos(y1[0])+" in mmin_simp2::set().";
589  O2SCL_ERR(err.c_str(),exc_ebadfunc);
590  }
591  for(i=0;i<dim;i++) x1[0][i]=ax[i];
592 
593  /* following points are initialized to x0+step_size */
594 
595  for (i=1;i<dim+1;i++) {
596  for(size_t j=0;j<dim;j++) x1[i][j]=x[j];
597  x1[i][i-1]=x1[i][i-1]+step_size[i-1];
598  y1[i]=ufunc(dim,x1[i]);
599  }
600 
601  /* Initialize simplex size */
602 
603  compute_center();
604  size=compute_size();
605 
606  set_called=true;
607 
608  return success;
609  }
610 
611  /// Set the function and initial simplex
612  template<class mat_t>
613  int set_simplex(func_t &ufunc, mat_t &sx) {
614 
615  if(dim==0) {
616  O2SCL_ERR2("Memory not allocated in ",
617  "mmin_simp2::set_simplex().",exc_ebadlen);
618  }
619 
620  func=&ufunc;
621 
622  for(size_t i=0;i<dim+1;i++) {
623  for(size_t j=0;j<dim;j++) {
624  x1[i][j]=sx(i,j);
625  }
626  y1[i]=ufunc(dim,x1[i]);
627  if (!std::isfinite(y1[i])) {
628  std::string err=((std::string)"Function not finite (returned ")+
629  dtos(y1[i])+" in mmin_simp2::set_simplex().";
630  O2SCL_ERR(err.c_str(),exc_ebadfunc);
631  }
632  }
633 
634  /* Initialize simplex size */
635 
636  compute_center();
637  size=compute_size();
638 
639  set_called=true;
640 
641  return success;
642  }
643 
644  /// Perform an iteration
645  virtual int iterate() {
646 
647  size_t n=dim+1;
648  size_t i;
649  size_t lo, hi, s_hi;
650  double dhi,ds_hi,dlo;
651  int status;
652  double val,val2;
653 
654  /* get index of highest, second highest and lowest point */
655 
656  lo=hi=0;
657  dlo=dhi=y1[0];
658  s_hi=1;
659  ds_hi=y1[1];
660 
661  for (i=1;i<n;i++) {
662  val=y1[i];
663  if(val<dlo) {
664  dlo=val;
665  lo=i;
666  } else if (val > dhi) {
667  ds_hi=dhi;
668  s_hi=hi;
669  dhi=val;
670  hi=i;
671  } else if(val > ds_hi) {
672  ds_hi=val;
673  s_hi=i;
674  }
675  }
676 
677  /* reflect the highest value */
678 
679  int ret1=try_corner_move(-1.0,hi,ws1,*func,dim,val);
680 
681  if (avoid_nonzero && ret1!=0) {
682  std::cout << "Found problem move1: " << std::endl;
683  for (size_t it=0;it<20 && ret1!=0;it++) {
684  ret1=try_corner_move(-1.0+((double)it)/10.0,hi,ws1,
685  *func,dim,val);
686  std::cout << "it,ret: " << it << " " << ret1 << std::endl;
687  }
688  if (ret1!=0) {
689  O2SCL_ERR2("Failed to move corner (1) in ",
690  "mmin_simp2::iterate().",exc_efailed);
691  }
692  }
693 
694  if (std::isfinite(val) && val<y1[lo]) {
695 
696  /* reflected point becomes lowest point, try expansion */
697 
698  int ret2=try_corner_move(-2.0,hi,ws2,*func,dim,val2);
699 
700  if (avoid_nonzero && ret2!=0) {
701  std::cout << "Found problem move2: " << std::endl;
702  for (size_t it=0;it<20 && ret2!=0;it++) {
703  ret2=try_corner_move(-2.0+((double)it)/6.7,hi,ws2,
704  *func,dim,val2);
705  std::cout << "it,ret: " << it << " " << ret2 << std::endl;
706  }
707  if (ret2!=0) {
708  O2SCL_ERR2("Failed to move corner (2) in ",
709  "mmin_simp2::iterate().",exc_efailed);
710  }
711  }
712 
713  if (std::isfinite(val2) && val2<y1[lo]) {
714  update_point(hi,ws2,val2);
715  } else {
716  update_point(hi,ws1,val);
717  }
718 
719  } else if (!std::isfinite(val) || val > y1[s_hi]) {
720 
721  /* reflection does not improve things enough */
722 
723  if (std::isfinite(val) && val <= y1[hi]) {
724 
725  /* if trial point is better than highest point, replace
726  highest point */
727 
728  update_point(hi,ws1,val);
729  }
730 
731  /* try one-dimensional contraction */
732 
733  int ret3=try_corner_move(0.5,hi,ws2,*func,dim,val2);
734 
735  if (avoid_nonzero && ret3!=0) {
736  std::cout << "Found problem move3: " << std::endl;
737  for (size_t it=0;it<20 && ret3!=0;it++) {
738  ret3=try_corner_move(0.025*((double)(19-it)),hi,ws2,
739  *func,dim,val2);
740  std::cout << "it,ret: " << it << " " << ret3 << std::endl;
741  }
742  if (ret3!=0) {
743  O2SCL_ERR2("Failed to move corner (2) in ",
744  "mmin_simp2::iterate().",exc_efailed);
745  }
746  }
747 
748  if (std::isfinite(val2) && val2 <= y1[hi]) {
749 
750  update_point(hi,ws2,val2);
751 
752  } else {
753 
754  /* contract the whole simplex in respect to the best point */
755  status=contract_by_best(lo,*func,dim);
756  if(status != 0) {
757  O2SCL_ERR("Function contract_by_best() failed in iterate().",
758  exc_efailed);
759  }
760 
761  }
762 
763  } else {
764 
765  /* trial point is better than second highest point.
766  Replace highest point by it */
767 
768  update_point(hi,ws1,val);
769  }
770 
771  /* return lowest point of simplex as x */
772 
773  vector_min(dim+1,y1,lo,val);
774  for(i=0;i<dim;i++) x[i]=x1[lo][i];
775  fval=y1[lo];
776 
777  /* Update simplex size */
778 
779  if (S2 > 0) {
780  size=sqrt(S2);
781  } else {
782  /* recompute if accumulated error has made size invalid */
783  size=compute_size();
784  }
785 
786  return success;
787  }
788 
789  /** \brief Print out iteration information.
790 
791  Depending on the value of the variable verbose, this prints out
792  the iteration information. If verbose=0, then no information is
793  printed, while if verbose>1, then after each iteration, the
794  present values of x and y are output to std::cout along with the
795  iteration number. If verbose>=2 then each iteration waits for a
796  character.
797  */
798  virtual int print_iter(size_t nv, vec_t &xx, vec_t *simp,
799  double y, int iter,
800  double value, double limit,
801  std::string comment) {
802 
803  if (this->verbose<=0) return 0;
804 
805  size_t i;
806  char ch;
807 
808  (*this->outs) << comment << " Iteration: " << iter << std::endl;
809  (*this->outs) << "x: ";
810  for(i=0;i<nv;i++) (*this->outs) << xx[i] << " ";
811  (*this->outs) << std::endl;
812  if (print_simplex>0) {
813  if (print_simplex==2) {
814  (*this->outs) << "Simplex Values:" << std::endl;
815  for(i=0;i<nv+1;i++) {
816  (*this->outs) << i << ": ";
817  for(size_t j=0;j<nv;j++) {
818  (*this->outs) << simp[i][j] << " ";
819  }
820  (*this->outs) << ": " << y1[i] << std::endl;
821  }
822  } else {
823  (*this->outs) << "Simplex Values:" << std::endl;
824  for(i=0;i<nv+1;i++) {
825  (*this->outs) << y1[i] << " ";
826  }
827  (*this->outs) << std::endl;
828  }
829  }
830  (*this->outs) << "y: " << y << " Val: " << value << " Lim: "
831  << limit << std::endl;
832  if (this->verbose>1) {
833  (*this->outs) << "Press a key and type enter to continue. ";
834  (*this->ins) >> ch;
835  }
836 
837  return 0;
838  }
839 
840  /// Return string denoting type("mmin_simp2")
841  virtual const char *type() { return "mmin_simp2";}
842 
843 #ifndef DOXYGEN_INTERNAL
844 
845  private:
846 
848  (const mmin_simp2<func_t,vec_t> &);
849  mmin_simp2<func_t,vec_t>& operator=
850  (const mmin_simp2<func_t,vec_t>&);
851 
852 #endif
853 
854  };
855 
856 #ifndef DOXYGEN_NO_O2NS
857 }
858 #endif
859 
860 #endif
vec_t ws1
Workspace vector 1.
Definition: mmin_simp2.h:145
double S2
Squared simplex size.
Definition: mmin_simp2.h:157
virtual int update_point(size_t i, vec_t &xx, double val)
Update point i in the simplex with values xx.
Definition: mmin_simp2.h:231
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
std::istream * ins
Stream for verbose input.
Definition: mmin.h:174
ubvector step_vec
Vector of step sizes.
Definition: mmin_simp2.h:354
func_t * func
Function.
Definition: mmin_simp2.h:348
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
std::ostream * outs
Stream for verbose output.
Definition: mmin.h:171
virtual int allocate(size_t n)
Allocate the memory.
Definition: mmin_simp2.h:553
vec_t delta
Desc.
Definition: mmin_simp2.h:153
invalid argument supplied by user
Definition: err_hnd.h:59
virtual int mmin_twovec(size_t nn, vec_t &xx, vec_t &xx2, double &fmin, func_t &ufunc)
Calculate the minimum min of func w.r.t the array x of size nvar, using xx and xx2 to specify the sim...
Definition: mmin_simp2.h:463
exceeded max number of iterations
Definition: err_hnd.h:73
Multidimensional minimization by the Simplex method (v2) (GSL)
Definition: mmin_simp2.h:123
int mmin_simplex(size_t nn, mat_t &sx, double &fmin, func_t &ufunc)
Calculate the minimum min of func w.r.t the array x of size nvar, given an initial simplex...
Definition: mmin_simp2.h:511
void vector_min(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the minimum of the first n elements of a vector.
Definition: vector.h:1208
generic failure
Definition: err_hnd.h:61
double size
Size of current simplex computed by iterate()
Definition: mmin_simp2.h:389
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
Definition: cblas_base.h:156
virtual int iterate()
Perform an iteration.
Definition: mmin_simp2.h:645
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: mmin.h:209
matrix, vector lengths are not conformant
Definition: err_hnd.h:89
Multidimensional minimization [abstract base].
Definition: mmin.h:164
vec_t xmc
Distance of vector from center.
Definition: mmin_simp2.h:155
bool avoid_nonzero
If true, try to automatically avoid regions where the function returns a non-zero value (default fals...
Definition: mmin_simp2.h:362
vec_t x
Present minimum vector computed by iterate()
Definition: mmin_simp2.h:392
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
size_t dim
Number of variables to be mind over.
Definition: mmin_simp2.h:345
virtual int try_corner_move(const double coeff, size_t corner, vec_t &xc, func_t &f, size_t nvar, double &newval)
Move a corner of a simplex.
Definition: mmin_simp2.h:207
double ddot(const size_t N, const vec_t &X, const vec2_t &Y)
Compute .
Definition: cblas_base.h:123
vec_t * x1
An array of n+1 vectors containing the simplex.
Definition: mmin_simp2.h:135
double compute_size()
Compute the size of the simplex.
Definition: mmin_simp2.h:183
virtual int mmin(size_t nn, vec_t &xx, double &fmin, func_t &ufunc)
Calculate the minimum min of func w.r.t the array x of size nvar.
Definition: mmin_simp2.h:410
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int last_ntrial
The number of iterations for in the most recent minimization.
Definition: mmin.h:206
double tol_abs
The independent variable tolerance.
Definition: mmin.h:203
int set_simplex(func_t &ufunc, mat_t &sx)
Set the function and initial simplex.
Definition: mmin_simp2.h:613
void daxpy(const double alpha, const size_t N, const vec_t &X, vec2_t &Y)
Compute .
Definition: cblas_base.h:98
int set_step(size_t nv, vec2_t &step)
Set the step sizes for each independent variable.
Definition: mmin_simp2.h:380
vec_t center
Center of simplex.
Definition: mmin_simp2.h:151
vec_t ws3
Workspace vector 3.
Definition: mmin_simp2.h:149
problem with user-supplied function
Definition: err_hnd.h:69
double fval
Function value at minimum computed by iterate()
Definition: mmin_simp2.h:395
ubvector y1
The n+1 function values at the simplex points.
Definition: mmin_simp2.h:143
virtual int contract_by_best(size_t best, func_t &f, size_t nvar)
Contract the simplex towards the best point.
Definition: mmin_simp2.h:280
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
vec_t ws2
Workspace vector 2.
Definition: mmin_simp2.h:147
std::string itos(int x)
Convert an integer to a string.
int print_simplex
Print simplex information in print_iter() (default 0)
Definition: mmin_simp2.h:405
int compute_center()
Compute the center of the simplex.
Definition: mmin_simp2.h:160
bool set_called
True if set() has been called.
Definition: mmin_simp2.h:351
Success.
Definition: err_hnd.h:47
virtual int print_iter(size_t nv, vec_t &xx, vec_t *simp, double y, int iter, double value, double limit, std::string comment)
Print out iteration information.
Definition: mmin_simp2.h:798
virtual const char * type()
Return string denoting type("mmin_simp2")
Definition: mmin_simp2.h:841
int ntrial
Maximum number of iterations.
Definition: mmin.h:197

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