49 #ifndef O2SCL_OOL_MMIN_GENCAN_H 50 #define O2SCL_OOL_MMIN_GENCAN_H 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> 61 #ifndef DOXYGEN_NO_O2NS 69 template<
class param_t,
class func_t,
class dfunc_t=func_t,
72 public ool_constr_min<param_t,func_t,dfunc_t,hfunc_t,vec_t> {
74 #ifndef DOXYGEN_INTERNAL 112 gsl_vector *X = M->x;
118 double *l = st->L->data;
119 double *u = st->U->data;
120 double *d = st->D->data;
123 double *xtrial = st->Xtrial->data;
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 );
144 gsl_vector_memcpy( st->D, st->Xtrial );
145 gsl_vector_sub( st->D, X );
148 imax = gsl_blas_idamax( st->D );
149 dinfn = fabs( gsl_vector_get( st->D, imax ) );
150 gsl_blas_ddot( gradient, st->D, >d );
153 OOL_CONMIN_EVAL_F( M, st->Xtrial, ftrial );
158 while (ftrial > M->f + P->gamma*alpha*gtd ) {
161 if (ftrial < M->f ) {
164 gsl_vector_memcpy( X, st->Xtrial );
166 return OOL_UNBOUNDEDF;
175 double atemp = ( -gtd*alpha*alpha )/
176 ( 2.0*(ftrial-M->f-alpha*gtd) );
178 if (atemp < P->sigma1 || atemp > P->sigma2*alpha ) {
187 gsl_vector_memcpy( st->Xtrial, X );
188 gsl_blas_daxpy( alpha, st->D, st->Xtrial );
191 OOL_CONMIN_EVAL_F( M, st->Xtrial, ftrial );
195 if( interp > P->mininterp &&
196 are_close( nn, alpha, d, x, P->epsrel, P->epsabs )) {
199 gsl_vector_memcpy( X, st->Xtrial );
207 gsl_vector_memcpy( X, st->Xtrial );
215 gsl_vector *X = M->x;
217 gsl_vector *Xplus = st->Xtrial;
221 double *g = gradient->data;
222 double *d = st->D->data;
223 double *xplus = Xplus->data;
226 const size_t nind = st->nind;
235 gtd = cblas_ddot( (
int)nind, g, 1, d, 1 );
238 alpha = GSL_MIN( 1.0, st->tnls_amax );
241 conmin_vector_memcpy( nind, xplus, x );
242 cblas_daxpy( alpha, (
int)nind,d, 1, xplus, 1 );
245 fplus = conmin_calc_f( M, nind, st->Ind, Xplus, X );
251 if( st->tnls_amax > 1.0 ) {
256 if( fplus <= M->f + P->gamma * alpha * gtd ) {
261 conmin_calc_g( M, nind, st->Ind, Xplus, X, gradient );
263 gptd = cblas_ddot( (
int)nind, g, 1, d, 1 );
266 if ( gptd >= P->beta * gtd ) {
271 conmin_vector_memcpy( nind, x, xplus );
275 return tnls_extrapolation( M, st, P, alpha, fplus );
278 return tnls_interpolation(M, st, P, alpha, fplus, gtd);
283 return tnls_extrapolation( M, st, P, alpha, fplus );
285 return tnls_interpolation(M, st, P, alpha, fplus, gtd);
291 int tnls_extrapolation(
double alpha,
double fplus) {
293 gsl_vector *X = M->x;
295 gsl_vector *Xplus = st->Xtrial;
299 double *d = st->D->data;
300 double *l = st->L->data;
301 double *u = st->U->data;
303 double *xplus = Xplus->data;
304 double *xtemp = st->tnls_Xtemp->data;
307 const size_t nind = st->nind;
323 if ( extrap > P->maxextrap ) {
326 conmin_vector_memcpy( nind, x, xplus );
328 if (extrap > 0 || st->tnls_amax < 1){
329 conmin_calc_g( M, nind, st->Ind, Xplus, X, gradient );
331 return TNLS_MAXEXTRAP;
335 if (alpha < st->tnls_amax && st->tnls_amax < P->next*alpha ) {
336 atemp = st->tnls_amax;
338 atemp = P->next * alpha;
342 conmin_vector_memcpy( nind, xtemp, x );
343 cblas_daxpy( atemp, (
int)nind, d, 1, xtemp, 1 );
346 if (atemp > st->tnls_amax ) {
347 conmin_vector_maxofmin( nind, xtemp, l, xtemp, u );
352 if( alpha > st->tnls_amax ) {
355 for (ii = 0; ii<nind && same; ii++) {
359 aux = P->epsrel * fabs( xplus[ii] );
361 if ( fabs(xtemp[ii]-xplus[ii]) > GSL_MAX(aux,P->epsabs)) {
371 conmin_vector_memcpy( nind, x, xplus );
373 if (extrap > 0 || st->tnls_amax < 1){
374 conmin_calc_g( M, nind, st->Ind, Xplus, X, gradient );
380 ftemp = conmin_calc_f( M, nind, st->Ind, st->tnls_Xtemp, X );
389 conmin_vector_memcpy( nind, xplus, xtemp );
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 );
416 int tnls_interpolation(
double alpha,
double fplus,
double gtd) {
418 gsl_vector *X = M->x;
420 gsl_vector *Xplus = st->Xtrial;
424 double *d = st->D->data;
425 double *xplus = Xplus->data;
428 const size_t nind = st->nind;
442 if (fplus <= M->f + P->gamma * alpha * gtd ) {
448 conmin_vector_memcpy( nind, x, xplus );
451 conmin_calc_g( M, nind, st->Ind, X, X, gradient );
457 alpha= alpha / P->nint;
460 atemp = -gtd * alpha*alpha /
461 (2 * (fplus - M->f - alpha * gtd));
463 if( atemp < P->sigma1 ||
464 atemp > P->sigma2*alpha ) {
465 alpha = alpha / P->nint;
472 conmin_vector_memcpy( nind, xplus, x );
473 cblas_daxpy( alpha, (
int)nind, d, 1, xplus, 1 );
476 fplus = conmin_calc_f( M, nind, st->Ind, Xplus, X );
480 if ( are_close( nind, alpha, d, x, P->epsrel, P->epsabs ) &&
481 interp > P->mininterp ){
492 double tnls_maximum_step() {
495 double *x = M->x->data;
496 double *l = st->L->data;
497 double *u = st->U->data;
498 double *d = st->D->data;
500 double step = P->infabs;
503 for( ii = 0; ii < st->nind; ii++ ) {
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 );
523 void spg_steplength() {
525 if (st->sty <= 0.0) {
526 st->lambda = GSL_MAX( 1.0, st->xeucn ) / sqrt( st->gpeucn2 );
529 double ss = st->sts / st->sty;
531 aux = GSL_MAX( P->lspgmi, ss );
532 st->lambda = GSL_MIN( P->lspgma, aux );
537 int actual_iterate() {
540 double *x = M->x->data;
541 double *l = st->L->data;
542 double *u = st->U->data;
553 gsl_vector_memcpy( st->S, M->x );
554 gsl_vector_memcpy( st->Y, M->gradient );
557 if ( st->gieucn2 <= st->ometa2 * st->gpeucn2 ) {
562 lsflag = spgls( M, st, P );
565 OOL_CONMIN_EVAL_DF( M, M->x, M->gradient );
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 );
578 cgflag = cg( M, st, P );
581 if ( cgflag == CG_BOUNDARY ) {
584 st->tnls_amax = tnls_maximum_step( M, st, P );
588 lsflag = tnls( M, st, P );
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 );
599 if ( lsflag == OOL_FLSEARCH ) {
602 lsflag = spgls( M, st, P );
605 OOL_CONMIN_EVAL_DF( M, M->x, M->gradient );
611 for( ii = 0; ii < st->n; ii++ ) {
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];
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 );
631 gsl_vector_sub ( st->S, M->x );
632 gsl_vector_scale( st->S, -1.0 );
633 gsl_vector_sub ( st->Y, M->gradient );
634 gsl_vector_scale( st->Y, -1.0 );
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 ) );
645 projected_gradient( st, M->x, M->gradient );
648 spg_steplength( st, P );
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) );
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,
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,
784 virtual int set(func_t &fn, dfunc_t &dfn, hfunc_t &hfn,
785 vec_t &init, param_t &par) {
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);
792 size_t nn = M->x->size;
795 if( nn != st->n || nn != M->fdf->n || nn != M->con->n )
801 gsl_vector_memcpy( st->L, M->con->L );
802 gsl_vector_memcpy( st->U, M->con->U );
813 int prepare_iteration {
816 double *x = M->x->data;
817 double *l = st->L->data;
818 double *u = st->U->data;
821 size_t nn = M->x->size;
825 conmin_vector_maxofmin( nn, x, l, u, 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 ) );
833 OOL_CONMIN_EVAL_FDF( M, M->x, &(M->f), M->gradient );
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 );
842 st->ometa2 = gsl_pow_2( 1.0 - P->eta );
843 st->epsgpen2 = gsl_pow_2( P->epsgpen );
846 projected_gradient( st, M->x, M->gradient );
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 ));
865 st->bcgeps = 2 * log10( P->cg_epsi ) -
866 st->acgeps * log10( st->gpeucn2 );
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 );
874 st->gpsupn0 = st->gpsupn;
875 st->gpeucn20 = st->gpeucn2;
878 if ( st->gpeucn2 != 0.0 ) {
879 st->lambda = GSL_MAX( 1.0, st->xeucn ) / sqrt( st->gpeucn2 );
883 if (P->udelta0 < 0.0 ) {
887 aux = 0.1 * GSL_MAX( 1.0, st->xeucn );
889 aux = 0.1 * GSL_MAX( 1.0, st->xsupn );
892 st->cg_delta = GSL_MAX( P->delmin, aux );
895 st->cg_delta = GSL_MAX( P->delmin, P->udelta0 );
899 M->size = st->gpsupn;
926 status = actual_iterate( M, st, P );
929 M->size = st->gpsupn;
932 gsl_vector_memcpy( M->dx, st->S );
951 const char *
type() {
return "mmin_constr_gencan"; }
953 #ifndef DOXYGEN_INTERNAL 960 (
const mmin_constr_gencan<func_t,dfunc_t,hfunc_t,vec_t>&);
966 #ifndef DOXYGEN_NO_O2NS vec_t cg_W
Temporary vector.
Class for automatically computing gradients [abstract base].
double epsrel
Relative small number (default 1.0e-7)
double infabs
Absolute infinite number (default 1.0e99)
virtual int is_optimal()
See if we'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='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
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.
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)
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 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.
double lspgmi
Minimum spectral steplength (default 1.0e-10)