mmin_constr_gencan.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 /*----------------------------------------------------------------------------*
24  * Open Optimization Library - Constrained Minimization
25  *
26  * gencan/gencan.c
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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
41  *
42  * Sergio Drumond Ventura
43  * Luis Alberto D'Afonseca
44  * Ricardo Biloti
45  * since: Feb 16th, 2004
46  *
47  * $Id: gencan.c,v 1.16 2005/05/17 19:08:18 biloti Exp $
48  *----------------------------------------------------------------------------*/
49 #ifndef O2SCL_OOL_MMIN_GENCAN_H
50 #define O2SCL_OOL_MMIN_GENCAN_H
51 
52 /** \file mmin_constr_gencan.h
53  \brief File defining \ref o2scl::mmin_constr_gencan
54 */
55 
56 #include <o2scl/text_file.h>
57 #include <o2scl/multi_funct.h>
58 #include <o2scl/ool_constr_min.h>
59 #include <gsl/gsl_math.h>
60 
61 #ifndef DOXYGEN_NO_O2NS
62 namespace o2scl {
63 #endif
64 
65  /** \brief Constrained minimization by the "GENCAN" method (OOL)
66 
67  \note Not working yet
68  */
69  template<class param_t, class func_t, class dfunc_t=func_t,
70  class hfunc_t=func_t, class vec_t=boost::numeric::ublas::vector<double> >
72  public ool_constr_min<param_t,func_t,dfunc_t,hfunc_t,vec_t> {
73 
74 #ifndef DOXYGEN_INTERNAL
75 
76  protected:
77 
78  /// Desc (default 1.0)
79  double cg_src;
80 
81  /// Temporary vector
82  vec_t S;
83  /// Temporary vector
84  vec_t Y;
85  /// Temporary vector
86  vec_t D;
87  /// Temporary vector
88  vec_t cg_W;
89  /// Temporary vector
90  vec_t cg_R;
91  /// Temporary vector
92  vec_t cg_D;
93  /// Temporary vector
94  vec_t cg_Sprev;
95  /// Temporary vector
96  vec_t Xtrial;
97  /// Temporary vector
98  vec_t tnls_Xtemp;
99  /// Temporary vector
100  vec_t near_l;
101  /// Temporary vector
102  vec_t near_u;
103 
104  /// Desc
105  int *Ind;
106 
107 #ifdef NEVER_DEFINED
108 
109  /// Desc
110  int spgls() {
111 
112  gsl_vector *X = M->x;
113  gsl_vector *gradient = M->gradient;
114 
115  size_t nn = X->size;
116 
117  /* Direct access to vector data */
118  double *l = st->L->data;
119  double *u = st->U->data;
120  double *d = st->D->data;
121  double *x = X->data;
122 
123  double *xtrial = st->Xtrial->data;
124 
125  /* Internal variables */
126  size_t interp;
127  size_t imax;
128 
129  double alpha;
130  double dinfn;
131  double gtd;
132  double ftrial;
133 
134  /* Compute first trial point, spectral projected gradient direction,
135  * and directional derivative <g,d> */
136  alpha = 1;
137 
138  /* Xtrial = min{ U, max[ L, ( X-lambda G ) ] } */
139  gsl_vector_memcpy( st->Xtrial, X );
140  gsl_blas_daxpy( -(st->lambda), gradient, st->Xtrial );
141  conmin_vector_minofmax( st->n, xtrial, u, l, xtrial );
142 
143  /* D = Xtrial - X */
144  gsl_vector_memcpy( st->D, st->Xtrial );
145  gsl_vector_sub( st->D, X );
146 
147  /* Inifite norm of D and < G, D > */
148  imax = gsl_blas_idamax( st->D );
149  dinfn = fabs( gsl_vector_get( st->D, imax ) );
150  gsl_blas_ddot( gradient, st->D, &gtd );
151 
152  /* Evaluate objective function */
153  OOL_CONMIN_EVAL_F( M, st->Xtrial, ftrial );
154 
155  interp = 0;
156 
157  /* While Armijo isn't satisfied repeat */
158  while (ftrial > M->f + P->gamma*alpha*gtd ) {
159 
160  /* Test if the trial point has a function value lower than fmin */
161  if (ftrial < M->f ) {
162 
163  M->f = ftrial;
164  gsl_vector_memcpy( X, st->Xtrial );
165 
166  return OOL_UNBOUNDEDF;
167  }
168 
169  interp++;
170 
171  if (alpha < P->sigma1 ) {
172  alpha /= P->nint;
173  } else {
174  /* quadratic model */
175  double atemp = ( -gtd*alpha*alpha )/
176  ( 2.0*(ftrial-M->f-alpha*gtd) );
177 
178  if (atemp < P->sigma1 || atemp > P->sigma2*alpha ) {
179  alpha /= P->nint;
180  } else {
181  alpha = atemp;
182  }
183  }
184 
185  /* Compute new trial point
186  * Xtrial = X + alpha D */
187  gsl_vector_memcpy( st->Xtrial, X );
188  gsl_blas_daxpy( alpha, st->D, st->Xtrial );
189 
190  /* Evaluate objective function */
191  OOL_CONMIN_EVAL_F( M, st->Xtrial, ftrial );
192 
193  /* Test whether at least mininterp interpolations were made
194  * and two consecutive iterates are close enough */
195  if( interp > P->mininterp &&
196  are_close( nn, alpha, d, x, P->epsrel, P->epsabs )) {
197 
198  M->f = ftrial;
199  gsl_vector_memcpy( X, st->Xtrial );
200 
201  return OOL_FLSEARCH;
202  }
203  }
204 
205  /* Set new values of f and X */
206  M->f = ftrial;
207  gsl_vector_memcpy( X, st->Xtrial );
208 
209  return OOL_SUCCESS;
210  }
211 
212  /// Desc
213  int tnls() {
214 
215  gsl_vector *X = M->x;
216  gsl_vector *gradient = M->gradient;
217  gsl_vector *Xplus = st->Xtrial;
218 
219  /* Direct access to vector data */
220  double *x = X->data;
221  double *g = gradient->data;
222  double *d = st->D->data;
223  double *xplus = Xplus->data;
224 
225  /* Constant values */
226  const size_t nind = st->nind;
227 
228  /* Internal variables */
229  double fplus;
230  double gtd;
231  double alpha;
232  double gptd;
233 
234  /* Compute directional derivative */
235  gtd = cblas_ddot( (int)nind, g, 1, d, 1 );
236 
237  /* Compute first trial */
238  alpha = GSL_MIN( 1.0, st->tnls_amax );
239 
240  /* Xplus = X + alpha D */
241  conmin_vector_memcpy( nind, xplus, x );
242  cblas_daxpy( alpha, (int)nind,d, 1, xplus, 1 );
243 
244  /* Evaluate objective function */
245  fplus = conmin_calc_f( M, nind, st->Ind, Xplus, X );
246 
247  /* Test Armijo and beta-condition and decide for:
248  * 1 - accepting the trial point,
249  * 2 - interpolate or
250  * 3 - extrapolate. */
251  if( st->tnls_amax > 1.0 ) {
252 
253  /* X+D belongs to the interior of the feasible set (amax > 1) */
254 
255  /* Verify Armijo */
256  if( fplus <= M->f + P->gamma * alpha * gtd ) {
257 
258  /* Armijo condition holds */
259 
260  /* Evaluate the gradient of objective function */
261  conmin_calc_g( M, nind, st->Ind, Xplus, X, gradient );
262  /* Eval gptd = < g, d > */
263  gptd = cblas_ddot( (int)nind, g, 1, d, 1 );
264 
265  /* Verify directional derivative (beta condition) */
266  if ( gptd >= P->beta * gtd ) {
267 
268  /* Step = 1 was ok, finish the line search */
269 
270  M->f = fplus;
271  conmin_vector_memcpy( nind, x, xplus );
272 
273  return OOL_SUCCESS;
274  } else {
275  return tnls_extrapolation( M, st, P, alpha, fplus );
276  }
277  } else {
278  return tnls_interpolation(M, st, P, alpha, fplus, gtd);
279  }
280  } else {
281  /* x + d does not belong to the feasible set (amax <= 1) */
282  if( fplus < M->f ) {
283  return tnls_extrapolation( M, st, P, alpha, fplus );
284  } else {
285  return tnls_interpolation(M, st, P, alpha, fplus, gtd);
286  }
287  }
288  }
289 
290  /// Desc
291  int tnls_extrapolation(double alpha, double fplus) {
292 
293  gsl_vector *X = M->x;
294  gsl_vector *gradient = M->gradient;
295  gsl_vector *Xplus = st->Xtrial;
296 
297  /* Direct access to vector data */
298  double *x = X->data;
299  double *d = st->D->data;
300  double *l = st->L->data;
301  double *u = st->U->data;
302 
303  double *xplus = Xplus->data;
304  double *xtemp = st->tnls_Xtemp->data;
305 
306  /* Constant values */
307  const size_t nind = st->nind;
308 
309  /* Internal variables */
310  double atemp;
311  double ftemp;
312 
313  size_t ii, extrap;
314  short same;
315 
316  /* Iterations */
317  extrap = 0;
318  do {
319 
320  extrap++;
321 
322  /* Test if maximum number of extrapolation was exceeded */
323  if ( extrap > P->maxextrap ) {
324 
325  M->f = fplus;
326  conmin_vector_memcpy( nind, x, xplus );
327 
328  if (extrap > 0 || st->tnls_amax < 1){
329  conmin_calc_g( M, nind, st->Ind, Xplus, X, gradient );
330  }
331  return TNLS_MAXEXTRAP;
332  }
333 
334  /* Chose new step */
335  if (alpha < st->tnls_amax && st->tnls_amax < P->next*alpha ) {
336  atemp = st->tnls_amax;
337  } else {
338  atemp = P->next * alpha;
339  }
340 
341  /* Compute new trial point. Xtemp = X + atemp*D */
342  conmin_vector_memcpy( nind, xtemp, x );
343  cblas_daxpy( atemp, (int)nind, d, 1, xtemp, 1 );
344 
345  /* Project */
346  if (atemp > st->tnls_amax ) {
347  conmin_vector_maxofmin( nind, xtemp, l, xtemp, u );
348  }
349 
350  /* Test if this is not the same point as the previous one.
351  * This test is performed only when alpha >= amax. */
352  if( alpha > st->tnls_amax ) {
353 
354  same = 1;
355  for (ii = 0; ii<nind && same; ii++) {
356 
357  double aux;
358 
359  aux = P->epsrel * fabs( xplus[ii] );
360 
361  if ( fabs(xtemp[ii]-xplus[ii]) > GSL_MAX(aux,P->epsabs)) {
362  same = 0;
363  }
364  }
365 
366  if (same) {
367 
368  /* Finish the extrapolation with the current point */
369  M->f = fplus;
370 
371  conmin_vector_memcpy( nind, x, xplus );
372 
373  if (extrap > 0 || st->tnls_amax < 1){
374  conmin_calc_g( M, nind, st->Ind, Xplus, X, gradient );
375  }
376  return OOL_SUCCESS;
377  }
378  }
379 
380  ftemp = conmin_calc_f( M, nind, st->Ind, st->tnls_Xtemp, X );
381 
382  if (ftemp < fplus) {
383 
384  /* If the functional value decreases then set the current
385  * point and continue the extrapolation */
386 
387  alpha = atemp;
388  fplus = ftemp;
389  conmin_vector_memcpy( nind, xplus, xtemp );
390 
391  continue;
392 
393  } else {
394 
395  /* If the functional value does not decrease then discard the
396  * last trial and finish the extrapolation with the previous
397  * point */
398 
399  M->f = fplus;
400 
401  conmin_vector_memcpy( nind, x, xplus );
402  if (extrap > 0 || st->tnls_amax < 1) {
403  conmin_calc_g( M, nind, st->Ind, X, X, gradient );
404  }
405 
406  return OOL_SUCCESS;
407  }
408 
409  } while (1);
410 
411  /* Just to make gcc happy */
412  return OOL_SUCCESS;
413  }
414 
415  /// Desc
416  int tnls_interpolation(double alpha, double fplus, double gtd) {
417 
418  gsl_vector *X = M->x;
419  gsl_vector *gradient = M->gradient;
420  gsl_vector *Xplus = st->Xtrial;
421 
422  /* Direct access to vector data */
423  double *x = X->data;
424  double *d = st->D->data;
425  double *xplus = Xplus->data;
426 
427  /* Constant values */
428  const size_t nind = st->nind;
429 
430  /* Internal variables */
431  size_t interp;
432  double atemp;
433 
434  /* Initialization */
435  interp = 0;
436 
437  /* Iterations */
438  do {
439  interp++;
440 
441  /* Test Armijo condition */
442  if (fplus <= M->f + P->gamma * alpha * gtd ) {
443 
444  /* Finish the line search */
445  M->f = fplus;
446 
447  /* X = Xplus */
448  conmin_vector_memcpy( nind, x, xplus );
449 
450  /* Evaluate objective function gradient */
451  conmin_calc_g( M, nind, st->Ind, X, X, gradient );
452 
453  return OOL_SUCCESS;
454  }
455  /* Compute new step */
456  if (alpha < P->sigma1 ) {
457  alpha= alpha / P->nint;
458  } else {
459  /* quadratic model */
460  atemp = -gtd * alpha*alpha /
461  (2 * (fplus - M->f - alpha * gtd));
462 
463  if( atemp < P->sigma1 ||
464  atemp > P->sigma2*alpha ) {
465  alpha = alpha / P->nint;
466  } else {
467  alpha = atemp;
468  }
469  }
470 
471  /* Compute new trial point: xplus = x + alpha d */
472  conmin_vector_memcpy( nind, xplus, x );
473  cblas_daxpy( alpha, (int)nind, d, 1, xplus, 1 );
474 
475  /* Evaluate objective function */
476  fplus = conmin_calc_f( M, nind, st->Ind, Xplus, X );
477 
478  /* Test whether at least mininterp interpolations were made
479  * and the steplength is soo small */
480  if ( are_close( nind, alpha, d, x, P->epsrel, P->epsabs ) &&
481  interp > P->mininterp ){
482  return OOL_FLSEARCH;
483  }
484  }
485  while( 1 );
486 
487  /* Just to make gcc happy */
488  return OOL_SUCCESS;
489  }
490 
491  /** Truncated Newton maximum step*/
492  double tnls_maximum_step() {
493 
494  /* Direct access to vector data */
495  double *x = M->x->data;
496  double *l = st->L->data;
497  double *u = st->U->data;
498  double *d = st->D->data;
499 
500  double step = P->infabs;
501  size_t ii;
502 
503  for( ii = 0; ii < st->nind; ii++ ) {
504 
505  if( *d > 0 ) {
506  double aux = ( *u - *x ) / *d;
507  step = GSL_MIN( step, aux );
508  } else if( *d < 0 ) {
509  double aux = ( *l - *x ) / *d;
510  step = GSL_MIN( step, aux );
511  }
512 
513  x++;
514  l++;
515  u++;
516  d++;
517  }
518 
519  return step;
520  }
521 
522  /** Spectral step length */
523  void spg_steplength() {
524 
525  if (st->sty <= 0.0) {
526  st->lambda = GSL_MAX( 1.0, st->xeucn ) / sqrt( st->gpeucn2 );
527  } else {
528  double aux;
529  double ss = st->sts / st->sty;
530 
531  aux = GSL_MAX( P->lspgmi, ss );
532  st->lambda = GSL_MIN( P->lspgma, aux );
533  }
534  }
535 
536  /** Iterate */
537  int actual_iterate() {
538 
539  /* Direct access to vector data */
540  double *x = M->x->data;
541  double *l = st->L->data;
542  double *u = st->U->data;
543  /* double *d = st->D->data; */
544 
545  /* Status of internal iterations */
546  int lsflag;
547  int cgflag;
548 
549  /* Internal variables */
550  size_t ii, imax;
551 
552  /* Saving previous values */
553  gsl_vector_memcpy( st->S, M->x );
554  gsl_vector_memcpy( st->Y, M->gradient );
555 
556  /* The actual iteration */
557  if ( st->gieucn2 <= st->ometa2 * st->gpeucn2 ) {
558  /* Compute the new iterate using an SPG iteration */
559 
560  /* Perform a line search with spectral continuous
561  projected gradient */
562  lsflag = spgls( M, st, P );
563 
564  /* Compute the gradient for the new iterate */
565  OOL_CONMIN_EVAL_DF( M, M->x, M->gradient );
566 
567  } else {
568 
569  /* The new iterate will belong to the closure of the actual face */
570 
571  /* Shrink the point, its gradient and the bounds */
572  conmin_shrink( st->nind, st->Ind, M->x );
573  conmin_shrink( st->nind, st->Ind, M->gradient );
574  conmin_shrink( st->nind, st->Ind, st->L );
575  conmin_shrink( st->nind, st->Ind, st->U );
576 
577  /* Compute the descent direction solving the newtonian system */
578  cgflag = cg( M, st, P );
579 
580  /* Compute maximum step for truncated newton line search */
581  if ( cgflag == CG_BOUNDARY ) {
582  st->tnls_amax = 1.0;
583  } else {
584  st->tnls_amax = tnls_maximum_step( M, st, P );
585  }
586 
587  /* Perform the line search */
588  lsflag = tnls( M, st, P );
589 
590  /* Expand the point, its gradient and the bounds */
591  conmin_expand( st->nind, st->Ind, M->x );
592  conmin_expand( st->nind, st->Ind, M->gradient );
593  conmin_expand( st->nind, st->Ind, st->L );
594  conmin_expand( st->nind, st->Ind, st->U );
595 
596  /* If the line search in the Truncated Newton direction
597  stopped due to a very small step discard this iteration
598  and force a SPG iteration */
599  if ( lsflag == OOL_FLSEARCH ) {
600 
601  /* Perform a line search with spectral projected gradient */
602  lsflag = spgls( M, st, P );
603 
604  /* Compute the gradient for the new iterate */
605  OOL_CONMIN_EVAL_DF( M, M->x, M->gradient );
606  }
607  }
608 
609  /* Prepare for the next iteration */
610  /* Adjustment */
611  for( ii = 0; ii < st->n; ii++ ) {
612  /* In principle, the current point could be slightly changed
613  * here, requiring a new function and gradient
614  * evaluation. But, according to the algorithms authors, this
615  * is done just to account for points that are "numerically¨
616  * at faces already. Thus, no additional evaluations are
617  * performed. (May 11th, 2005).
618  */
619  if ( x[ii] <= st->near_l[ii] ) x[ii] = l[ii];
620  else if( x[ii] >= st->near_u[ii] ) x[ii] = u[ii];
621  }
622 
623  /* Compute infinite and Euclidian-norm of X */
624  imax = gsl_blas_idamax( M->x );
625  st->xsupn = fabs( gsl_vector_get( M->x, imax ) );
626  st->xeucn = gsl_blas_dnrm2 ( M->x );
627 
628  /* Until now S = X_prev, now S = X - X_prev
629  * Compute s = x_{k+1} - x_k = X - S
630  * and y = g_{k+1} - g_k = G - Y */
631  gsl_vector_sub ( st->S, M->x ); /* S = S - X */
632  gsl_vector_scale( st->S, -1.0 ); /* S = -S = X - S_prev */
633  gsl_vector_sub ( st->Y, M->gradient ); /* Y = Y - G */
634  gsl_vector_scale( st->Y, -1.0 ); /* Y = -Y = G - Y_prev */
635 
636  /* Compute sts = s dot s
637  * sty = s dot y
638  * and sinf = |s|_inf */
639  gsl_blas_ddot( st->S, st->S, &(st->sts) );
640  gsl_blas_ddot( st->S, st->Y, &(st->sty) );
641  imax = gsl_blas_idamax( st->S );
642  st->sinf = fabs( gsl_vector_get( st->S, imax ) );
643 
644  /* Compute continuous project gradient */
645  projected_gradient( st, M->x, M->gradient );
646 
647  /* Update spectral steplength */
648  spg_steplength( st, P );
649 
650  /* Update trust-region radius */
651  if ( P->trtype ) st->cg_delta = GSL_MAX( P->delmin, 10*sqrt( st->sts ) );
652  else st->cg_delta = GSL_MAX( P->delmin, 10* ( st->sinf) );
653 
654  return lsflag;
655  }
656 
657 #endif
658 
659 #endif
660 
661  public:
662 
663  mmin_constr_gencan() {
664  epsgpen=1.0e-5;
665  epsgpsn=1.0e-5;
666  fmin=-1.0e99;
667  udelta0=-1.0;
668  ucgmia=-1.0;
669  ucgmib=-1.0;
670  cg_src=1.0;
671  cg_scre=1.0;
672  cg_gpnf=1.0e-5;
673  cg_epsi=1.0e-1;
674  cg_epsf=1.0e-5;
675  cg_epsnqmp=1.0e-4;
676  cg_maxitnqmp=5;
677  nearlyq=0;
678  nint=2.0;
679  next=2.0;
680  mininterp=4;
681  maxextrap=100;
682  trtype=0;
683  eta=0.9;
684  delmin=0.1;
685  lspgmi=1.0e-10;
686  lspgma=1.0e10;
687  theta=1.0e-6;
688  gamma=1.0e-4;
689  beta=0.5;
690  sigma1=0.1;
691  sigma2=0.9;
692  epsrel=1.0e-7;
693  epsabs=1.0e-10;
694  infrel=1.0e20;
695  infabs=1.0e99;
696  }
697 
698  /// Tolerance on Euclidean norm of projected gradient (default 1.0e-5)
699  double epsgpen;
700  /// Tolerance on infinite norm of projected gradient (default 1.0e-5)
701  double epsgpsn;
702  /** \brief Minimum function value (default \f$ 10^{-99} \f$)
703 
704  If the function value is below this value, then the algorithm
705  assumes that the function is not bounded and exits.
706  */
707  double fmin;
708  /// Trust-region radius (default -1.0)
709  double udelta0;
710  /// Maximum interations per variable (default -1.0)
711  double ucgmia;
712  /// Extra maximum iterations (default -1.0)
713  double ucgmib;
714  /// Conjugate gradient condition type (default 1)
715  int cg_scre;
716  /// Projected gradient norm (default 1.0e-5)
717  double cg_gpnf;
718  /// Desc (default 1.0e-1)
719  double cg_epsi;
720  /// Desc (default 1.0e-5)
721  double cg_epsf;
722  /// Stopping fractional tolerance for conjugate gradient (default 1.0e-4)
723  double cg_epsnqmp;
724  /// Maximum iterations for conjugate gradient (default 5)
726  /// Set to 1 if the function is nearly quadratic (default 0)
727  int nearlyq;
728  /// Interpolation constant (default 2.0)
729  double nint;
730  /// Extrapolation constant (default 2.0)
731  double next;
732  /// Minimum interpolation size (default 4)
734  /// Maximum extrapolations in truncated Newton direction (default 100)
736  /// Type of trust region (default 0)
737  int trtype;
738  /// Threshold for abandoning current face (default 0.9)
739  double eta;
740  /// Minimum trust region for truncated Newton direction (default 0.1)
741  double delmin;
742  /// Minimum spectral steplength (default 1.0e-10)
743  double lspgmi;
744  /// Maximum spectral steplength (default 1.0e10)
745  double lspgma;
746  /// Constant for the angle condition (default 1.0e-6)
747  double theta;
748  /// Constant for Armijo condition (default 1.0e-4)
749  double gamma;
750  /// Constant for beta condition (default 0.5)
751  double beta;
752  /// Lower bound to the step length reduction (default 0.1)
753  double sigma1;
754  /// Upper bound to the step length reduction (default 0.9)
755  double sigma2;
756  /// Relative small number (default 1.0e-7)
757  double epsrel;
758  /// Absolute small number (default 1.0e-10)
759  double epsabs;
760  /// Relative infinite number (default 1.0e20)
761  double infrel;
762  /// Absolute infinite number (default 1.0e99)
763  double infabs;
764 
765  /// Allocate memory
766  virtual int alloc(const size_t n) {
767  if (this->dim>0) free();
768  this->ao.allocate(xx,n);
769  this->ao.allocate(d,n);
770  this->ao.allocate(s,n);
771  this->ao.allocate(y,n);
772  return ool_constr_min<param_t,func_t,dfunc_t,hfunc_t,vec_t,vec_t,
773  alloc_t>::alloc(n);
774  }
775 
776  /// Free previously allocated memory
777  virtual int free() {
778  if (this->dim>0) this->ao.free(xx);
779  return ool_constr_min<param_t,func_t,dfunc_t,hfunc_t,vec_t,vec_t,
780  alloc_t>::free();
781  }
782 
783  /// Set the function, the initial guess, and the parameters
784  virtual int set(func_t &fn, dfunc_t &dfn, hfunc_t &hfn,
785  vec_t &init, param_t &par) {
786 
787  int ret=ool_constr_min<param_t,func_t,dfunc_t,hfunc_t,vec_t,vec_t,
788  alloc_t>::set(fn,dfn,hfn,init,par);
789 
790 #ifdef NEVER_DEFINED
791  // Internal variables
792  size_t nn = M->x->size;
793 
794  // Checking dimensions
795  if( nn != st->n || nn != M->fdf->n || nn != M->con->n )
796  {
797  return OOL_EINVAL;
798  }
799 
800  // Copy boundary vectors
801  gsl_vector_memcpy( st->L, M->con->L );
802  gsl_vector_memcpy( st->U, M->con->U );
803 
804 #endif
805 
806  prepare_iteration();
807 
808  return 0;
809  }
810 
811 #ifdef NEVER_DEFINED
812 
813  int prepare_iteration {
814 
815  /* Direct access to vector data */
816  double *x = M->x->data;
817  double *l = st->L->data;
818  double *u = st->U->data;
819 
820  /* Internal variables */
821  size_t nn = M->x->size;
822  size_t ii, imax;
823 
824  /* Impose factibility */
825  conmin_vector_maxofmin( nn, x, l, u, x );
826 
827  /* Eval Euclidean and infinity norms of X */
828  st->xeucn = gsl_blas_dnrm2 ( M->x );
829  imax = gsl_blas_idamax( M->x );
830  st->xsupn = fabs( gsl_vector_get( M->x, imax ) );
831 
832  /* Evaluate objective function and its gradient */
833  OOL_CONMIN_EVAL_FDF( M, M->x, &(M->f), M->gradient );
834 
835  /* Define near_l and near_u vector */
836  for (ii=0; ii < nn; ii++){
837  st->near_l[ii] = l[ii] + GSL_MAX( P->epsrel*fabs( l[ii] ), P->epsabs );
838  st->near_u[ii] = u[ii] - GSL_MAX( P->epsrel*fabs( u[ii] ), P->epsabs );
839  }
840 
841  /* Setting constant parameters */
842  st->ometa2 = gsl_pow_2( 1.0 - P->eta );
843  st->epsgpen2 = gsl_pow_2( P->epsgpen );
844 
845  /* Compute continuous project gradient */
846  projected_gradient( st, M->x, M->gradient );
847 
848  /* Compute a linear relation between gpeucn2 and cgeps2, i.e.,
849  * scalars a and b such that
850  *
851  * a * log10(||g_P(x_0)||_2^2) + b = log10(cgeps_0^2) and
852  *
853  * a * log10(||g_P(x_f)||_2^2) + b = log10(cgeps_f^2),
854  *
855  * where cgeps_0 and cgeps_f are provided. Note that if
856  * cgeps_0 is equal to cgeps_f then cgeps will be always
857  * equal to cgeps_0 and cgeps_f.
858  *
859  * We introduce now a linear relation between gpsupn and cgeps also.
860  */
861  if (P->cg_scre == 1 ) {
862  st->acgeps = 2 *( log10( P->cg_epsf / P->cg_epsi ) /
863  log10( P->cg_gpnf * P->cg_gpnf / st->gpeucn2 ));
864 
865  st->bcgeps = 2 * log10( P->cg_epsi ) -
866  st->acgeps * log10( st->gpeucn2 );
867  } else {
868  st->acgeps = ( log10( P->cg_epsf / P->cg_epsi ) /
869  log10( P->cg_gpnf / st->gpsupn ) );
870  st->bcgeps = log10( P->cg_epsi ) - st->acgeps * log10( st->gpsupn );
871  }
872 
873  /* And it will be used for the linear relation of cgmaxit */
874  st->gpsupn0 = st->gpsupn;
875  st->gpeucn20 = st->gpeucn2;
876 
877  /* Initialize the spectral steplength */
878  if ( st->gpeucn2 != 0.0 ) {
879  st->lambda = GSL_MAX( 1.0, st->xeucn ) / sqrt( st->gpeucn2 );
880  }
881 
882  /* Initialize the trust-region radius */
883  if (P->udelta0 < 0.0 ) {
884 
885  double aux;
886  if ( P->trtype ) {
887  aux = 0.1 * GSL_MAX( 1.0, st->xeucn );
888  } else {
889  aux = 0.1 * GSL_MAX( 1.0, st->xsupn );
890  }
891 
892  st->cg_delta = GSL_MAX( P->delmin, aux );
893 
894  } else {
895  st->cg_delta = GSL_MAX( P->delmin, P->udelta0 );
896  }
897 
898  /* Export size */
899  M->size = st->gpsupn;
900 
901  return OOL_SUCCESS;
902  }
903 
904 #endif
905 
906  /// Restart the minimizer
907  virtual int restart() {
908 
909  /*
910  // Restarting dx
911  gsl_vector_set_zero( M->dx );
912 
913  return prepare_iteration( M );
914  */
915 
916  return 0;
917  }
918 
919  /// Perform an iteration
920  virtual int iterate() {
921 
922 #ifdef NEVER_DEFINED
923 
924  int status;
925 
926  status = actual_iterate( M, st, P );
927 
928  /* Export size and dx variables */
929  M->size = st->gpsupn;
930 
931  /* In the next version does dx replace st->S ? */
932  gsl_vector_memcpy( M->dx, st->S );
933 
934  return status;
935 
936 #endif
937 
938  return 0;
939  }
940 
941  /// See if we're finished
942  virtual int is_optimal() {
943 
944  //return (( st->gpeucn2 <= st->epsgpen2 ||
945  //st->gpsupn <= P->epsgpsn ||
946  //M->f <= P->fmin )? OOL_SUCCESS : OOL_CONTINUE );
947 
948  }
949 
950  /// Return string denoting type ("mmin_constr_gencan")
951  const char *type() { return "mmin_constr_gencan"; }
952 
953 #ifndef DOXYGEN_INTERNAL
954 
955  private:
956 
961 
962 #endif
963 
964  };
965 
966 #ifndef DOXYGEN_NO_O2NS
967 }
968 #endif
969 
970 #endif
971 
o2scl::mmin_constr_gencan::cg_D
vec_t cg_D
Temporary vector.
Definition: mmin_constr_gencan.h:92
o2scl::mmin_constr_gencan::set
virtual int set(func_t &fn, dfunc_t &dfn, hfunc_t &hfn, vec_t &init, param_t &par)
Set the function, the initial guess, and the parameters.
Definition: mmin_constr_gencan.h:784
o2scl::mmin_constr_gencan::gamma
double gamma
Constant for Armijo condition (default 1.0e-4)
Definition: mmin_constr_gencan.h:749
o2scl::mmin_constr_gencan
Constrained minimization by the "GENCAN" method (OOL)
Definition: mmin_constr_gencan.h:71
o2scl::mmin_constr_gencan::infabs
double infabs
Absolute infinite number (default 1.0e99)
Definition: mmin_constr_gencan.h:763
o2scl::mmin_constr_gencan::lspgma
double lspgma
Maximum spectral steplength (default 1.0e10)
Definition: mmin_constr_gencan.h:745
o2scl::mmin_constr_gencan::cg_epsi
double cg_epsi
Desc (default 1.0e-1)
Definition: mmin_constr_gencan.h:719
o2scl::mmin_constr_gencan::fmin
double fmin
Minimum function value (default )
Definition: mmin_constr_gencan.h:707
o2scl::mmin_constr_gencan::Ind
int * Ind
Desc.
Definition: mmin_constr_gencan.h:105
o2scl::mmin_constr_gencan::cg_epsf
double cg_epsf
Desc (default 1.0e-5)
Definition: mmin_constr_gencan.h:721
o2scl::mmin_constr_gencan::epsrel
double epsrel
Relative small number (default 1.0e-7)
Definition: mmin_constr_gencan.h:757
boost::numeric::ublas::vector< double >
o2scl::gradient
Class for automatically computing gradients [abstract base].
Definition: mmin.h:54
o2scl::mmin_constr_gencan::alloc
virtual int alloc(const size_t n)
Allocate memory.
Definition: mmin_constr_gencan.h:766
o2scl::mmin_constr_gencan::restart
virtual int restart()
Restart the minimizer.
Definition: mmin_constr_gencan.h:907
o2scl::mmin_constr_gencan::near_u
vec_t near_u
Temporary vector.
Definition: mmin_constr_gencan.h:102
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::mmin_constr_gencan::mininterp
int mininterp
Minimum interpolation size (default 4)
Definition: mmin_constr_gencan.h:733
o2scl::mmin_constr_gencan::near_l
vec_t near_l
Temporary vector.
Definition: mmin_constr_gencan.h:100
o2scl::mmin_constr_gencan::tnls_Xtemp
vec_t tnls_Xtemp
Temporary vector.
Definition: mmin_constr_gencan.h:98
o2scl::mmin_constr_gencan::type
const char * type()
Return string denoting type ("mmin_constr_gencan")
Definition: mmin_constr_gencan.h:951
o2scl::mmin_constr_gencan::cg_scre
int cg_scre
Conjugate gradient condition type (default 1)
Definition: mmin_constr_gencan.h:715
o2scl::mmin_constr_gencan::Xtrial
vec_t Xtrial
Temporary vector.
Definition: mmin_constr_gencan.h:96
o2scl::mmin_constr_gencan::lspgmi
double lspgmi
Minimum spectral steplength (default 1.0e-10)
Definition: mmin_constr_gencan.h:743
o2scl::mmin_constr_gencan::cg_src
double cg_src
Desc (default 1.0)
Definition: mmin_constr_gencan.h:79
o2scl::mmin_constr_gencan::cg_maxitnqmp
int cg_maxitnqmp
Maximum iterations for conjugate gradient (default 5)
Definition: mmin_constr_gencan.h:725
o2scl::mmin_constr_gencan::S
vec_t S
Temporary vector.
Definition: mmin_constr_gencan.h:82
o2scl::mmin_constr_gencan::next
double next
Extrapolation constant (default 2.0)
Definition: mmin_constr_gencan.h:731
o2scl::mmin_constr_gencan::Y
vec_t Y
Temporary vector.
Definition: mmin_constr_gencan.h:84
o2scl::mmin_constr_gencan::cg_Sprev
vec_t cg_Sprev
Temporary vector.
Definition: mmin_constr_gencan.h:94
o2scl::mmin_constr_gencan::epsgpsn
double epsgpsn
Tolerance on infinite norm of projected gradient (default 1.0e-5)
Definition: mmin_constr_gencan.h:701
o2scl::mmin_constr_gencan::cg_epsnqmp
double cg_epsnqmp
Stopping fractional tolerance for conjugate gradient (default 1.0e-4)
Definition: mmin_constr_gencan.h:723
o2scl::mmin_constr_gencan::nearlyq
int nearlyq
Set to 1 if the function is nearly quadratic (default 0)
Definition: mmin_constr_gencan.h:727
o2scl::mmin_constr_gencan::beta
double beta
Constant for beta condition (default 0.5)
Definition: mmin_constr_gencan.h:751
o2scl::mmin_constr_gencan::cg_gpnf
double cg_gpnf
Projected gradient norm (default 1.0e-5)
Definition: mmin_constr_gencan.h:717
o2scl::mmin_constr_gencan::ucgmia
double ucgmia
Maximum interations per variable (default -1.0)
Definition: mmin_constr_gencan.h:711
o2scl::mmin_constr_gencan::epsgpen
double epsgpen
Tolerance on Euclidean norm of projected gradient (default 1.0e-5)
Definition: mmin_constr_gencan.h:699
o2scl::mmin_constr_gencan::iterate
virtual int iterate()
Perform an iteration.
Definition: mmin_constr_gencan.h:920
o2scl::mmin_constr_gencan::cg_W
vec_t cg_W
Temporary vector.
Definition: mmin_constr_gencan.h:88
o2scl::mmin_constr_gencan::delmin
double delmin
Minimum trust region for truncated Newton direction (default 0.1)
Definition: mmin_constr_gencan.h:741
o2scl::mmin_constr_gencan::theta
double theta
Constant for the angle condition (default 1.0e-6)
Definition: mmin_constr_gencan.h:747
o2scl::mmin_constr_gencan::ucgmib
double ucgmib
Extra maximum iterations (default -1.0)
Definition: mmin_constr_gencan.h:713
o2scl::mmin_constr_gencan::udelta0
double udelta0
Trust-region radius (default -1.0)
Definition: mmin_constr_gencan.h:709
o2scl::mmin_constr_gencan::nint
double nint
Interpolation constant (default 2.0)
Definition: mmin_constr_gencan.h:729
o2scl::interp
Interpolation class for general vectors.
Definition: interp.h:1654
o2scl::mmin_constr_gencan::epsabs
double epsabs
Absolute small number (default 1.0e-10)
Definition: mmin_constr_gencan.h:759
o2scl::mmin_constr_gencan::free
virtual int free()
Free previously allocated memory.
Definition: mmin_constr_gencan.h:777
o2scl::mmin_constr_gencan::cg_R
vec_t cg_R
Temporary vector.
Definition: mmin_constr_gencan.h:90
o2scl::mmin_constr_gencan::is_optimal
virtual int is_optimal()
See if we're finished.
Definition: mmin_constr_gencan.h:942
o2scl::mmin_constr_gencan::maxextrap
int maxextrap
Maximum extrapolations in truncated Newton direction (default 100)
Definition: mmin_constr_gencan.h:735
o2scl::mmin_constr_gencan::infrel
double infrel
Relative infinite number (default 1.0e20)
Definition: mmin_constr_gencan.h:761
o2scl::mmin_constr_gencan::D
vec_t D
Temporary vector.
Definition: mmin_constr_gencan.h:86
o2scl::mmin_constr_gencan::sigma2
double sigma2
Upper bound to the step length reduction (default 0.9)
Definition: mmin_constr_gencan.h:755
o2scl::mmin_constr_gencan::trtype
int trtype
Type of trust region (default 0)
Definition: mmin_constr_gencan.h:737
o2scl::mmin_constr_gencan::eta
double eta
Threshold for abandoning current face (default 0.9)
Definition: mmin_constr_gencan.h:739
o2scl::mmin_constr_gencan::sigma1
double sigma1
Lower bound to the step length reduction (default 0.1)
Definition: mmin_constr_gencan.h:753

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