mmin_constr_gencan.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 /*----------------------------------------------------------------------------*
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 
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 
960  (const mmin_constr_gencan<func_t,dfunc_t,hfunc_t,vec_t>&);
961 
962 #endif
963 
964  };
965 
966 #ifndef DOXYGEN_NO_O2NS
967 }
968 #endif
969 
970 #endif
971 
vec_t cg_W
Temporary vector.
Class for automatically computing gradients [abstract base].
Definition: mmin.h:54
double epsrel
Relative small number (default 1.0e-7)
double infabs
Absolute infinite number (default 1.0e99)
virtual int is_optimal()
See if we&#39;re finished.
double cg_epsf
Desc (default 1.0e-5)
int maxextrap
Maximum extrapolations in truncated Newton direction (default 100)
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
double sigma1
Lower bound to the step length reduction (default 0.1)
vec_t tnls_Xtemp
Temporary vector.
vec_t near_u
Temporary vector.
double cg_epsi
Desc (default 1.0e-1)
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.
virtual int alloc(const size_t n)
Allocate memory.
const char * type()
Return string denoting type ("mmin_constr_gencan")
virtual int restart()
Restart the minimizer.
double ucgmia
Maximum interations per variable (default -1.0)
double cg_epsnqmp
Stopping fractional tolerance for conjugate gradient (default 1.0e-4)
double gamma
Constant for Armijo condition (default 1.0e-4)
double eta
Threshold for abandoning current face (default 0.9)
vec_t cg_D
Temporary vector.
vec_t Xtrial
Temporary vector.
vec_t Y
Temporary vector.
double cg_gpnf
Projected gradient norm (default 1.0e-5)
double sigma2
Upper bound to the step length reduction (default 0.9)
int cg_maxitnqmp
Maximum iterations for conjugate gradient (default 5)
double theta
Constant for the angle condition (default 1.0e-6)
virtual int iterate()
Perform an iteration.
double fmin
Minimum function value (default )
int mininterp
Minimum interpolation size (default 4)
double delmin
Minimum trust region for truncated Newton direction (default 0.1)
double epsgpsn
Tolerance on infinite norm of projected gradient (default 1.0e-5)
double cg_src
Desc (default 1.0)
double nint
Interpolation constant (default 2.0)
double next
Extrapolation constant (default 2.0)
int trtype
Type of trust region (default 0)
virtual int free()
Free previously allocated memory.
vec_t cg_Sprev
Temporary vector.
int nearlyq
Set to 1 if the function is nearly quadratic (default 0)
double epsabs
Absolute small number (default 1.0e-10)
vec_t D
Temporary vector.
double ucgmib
Extra maximum iterations (default -1.0)
Constrained minimization by the "GENCAN" method (OOL)
double udelta0
Trust-region radius (default -1.0)
vec_t S
Temporary vector.
vec_t cg_R
Temporary vector.
double lspgma
Maximum spectral steplength (default 1.0e10)
double epsgpen
Tolerance on Euclidean norm of projected gradient (default 1.0e-5)
double infrel
Relative infinite number (default 1.0e20)
vec_t near_l
Temporary vector.
double beta
Constant for beta condition (default 0.5)
int cg_scre
Conjugate gradient condition type (default 1)
Interpolation class for general vectors.
Definition: interp.h:1558
double lspgmi
Minimum spectral steplength (default 1.0e-10)

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