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;
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;
294 gsl_vector *gradient = M->gradient;
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;
419 gsl_vector *gradient = M->gradient;
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));
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) );
663 mmin_constr_gencan() {
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
966 #ifndef DOXYGEN_NO_O2NS