78 rhoStart=
mmax(rhoStart,rhoEnd);
79 int dim=of->
dim(), info,
k, t, nerror;
80 double rho=rhoStart,
delta=rhoStart, rhoNew,
81 lambda1, normD=rhoEnd+1.0, modelStep, reduction, r, valueOF, valueFk, bound, noise;
83 bool improvement, forceTRStep=
true, evalNeeded;
96 if (poly.NewtonBasis==NULL)
98 printf(
"cannot construct lagrange basis.\n");
106 valueFk=poly.valuesF[
k];
170 poly.vBase+=poly.NewtonPoints[
k];
173 of->
xBest=poly.vBase;
178 poly.gradientHessian(poly.NewtonPoints[k],vG,mH);
194 d+=poly.NewtonPoints[
k];
198 reduction=-poly.shiftedEval(d,valueFk);
201 if ((normD<0.5*rho)&&(!forceTRStep)) { evalNeeded=
true;
break; }
207 if ((reduction<noise)&&(!forceTRStep)) { evalNeeded=
true;
break; }
208 forceTRStep=
false; evalNeeded=
false;
211 tmp=poly.vBase+
d; nerror=0; valueOF=of->
eval(tmp, &nerror);
217 if (normD>=2*rho)
continue;
222 printf(
"violation: %e\n",r);
226 r=(valueFk-valueOF)/reduction;
227 if (r<=0.1) delta=0.5*normD;
228 else if (r<0.7) delta=
mmax(0.5*delta, normD);
229 else delta=
mmax(rho+ normD,
mmax(1.25*normD, delta));
231 if (delta<1.5*rho) delta=
rho;
235 t=poly.findAGoodPointToReplace(-1, rho, d,&modelStep);
236 k=t; valueFk=valueOF;
241 t=poly.findAGoodPointToReplace(k, rho, d,&modelStep);
246 if (t<0) { poly.updateM(d, valueOF);
break; }
253 if ((!of->
isConstrained)||(improvement)||(reduction>0.0)||(normD<rho)) poly.replace(t, d, valueOF);
255 if (improvement)
continue;
257 if (modelStep>2*rho)
continue;
258 if (normD>=2*rho)
continue;
268 bound=0.5*
sqr(rho)*lambda1;
269 if (poly.nUpdateOfM<10) bound=0.0;
275 t=poly.checkIfValidityIsInBound(d, k, bound, rho );
279 tmp=poly.vBase+
d; nerror=0; valueOF=of->
eval(tmp, &nerror);
284 poly.NewtonBasis[t].gradientHessian(poly.NewtonPoints[k],GXk,H);
285 double rhot=
rho,vmax;
291 d+=poly.NewtonPoints[
k];
292 tmp=poly.vBase+
d; nerror=0; valueOF=of->
eval(tmp, &nerror);
296 poly.replace(t, d, valueOF);
297 if ((valueOF<valueFk)&&
298 (of->
isFeasible(tmp))) { k=t; valueFk=valueOF; };
305 if ((normD<=rho)||(reduction<0.0))
break;
311 if (rho<=rhoEnd)
break;
315 if (rho<16*rhoEnd) rhoNew=rhoEnd;
316 else if (rho<250*rhoEnd) rhoNew=
sqrt(rho*rhoEnd);
318 delta=
mmax(0.5*rho,rhoNew);
324 poly.translate(poly.NewtonPoints[k]);
335 tmp=poly.vBase+
d; nerror=0; valueOF=of->
eval(tmp,&nerror);
337 if ((nerror)||(valueOF<valueFk))
339 poly.vBase+=poly.NewtonPoints[
k];
340 poly.gradientHessian(poly.NewtonPoints[k],vG,mH);
344 valueFk=valueOF; poly.vBase=tmp;
345 poly.gradientHessian(d,vG,mH);
349 poly.vBase+=poly.NewtonPoints[
k];
350 poly.gradientHessian(poly.NewtonPoints[k],vG,mH);
358 of->
xBest=poly.vBase;
virtual void saveValue(Vector tmp, double valueOF, int nerror)
virtual void finalize(Vector vG, Matrix mH, Vector vLambda)
MultIndCache cacheMultInd
void fullUpdateOfM(double rho, Vector Base, Matrix data, InterPolynomial poly)
void sqrt(Image< double > &op)
void parallelImprove(InterPolynomial *p, int *_k, double _rho, double *_valueFk, Vector _Base)
MultInd * get(unsigned _dim, unsigned _deg)
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
virtual double eval(Vector v, int *nerror)=0
void initConstrainedStep(ObjectiveFunction *of)
char isFeasible(Vector vx, double *d=NULL)
Vector LAGMAXModified(Vector G, Matrix H, double rho, double &VMAX)
char checkForTermination(Vector d, Vector Base, double rhoEnd)
void(* quickHack)(InterPolynomial, int)
void startParallelThread()
void parallelInit(int _nnode, int _dim, ObjectiveFunction *_of)
static Vector emptyVector
Vector ConstrainedL2NormMinimizer(InterPolynomial poly, int k, double delta, int *info, int iterMax, double *lambda1, Vector vOBase, ObjectiveFunction *of)