1 |
#include <math.h>
|
2 |
#include "minimizers/OOPSEMinimizer.hpp"
|
3 |
#include "integrators/Integrator.cpp"
|
4 |
|
5 |
OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
|
6 |
MinimizerParameterSet * param) :
|
7 |
RealIntegrator(theInfo, the_ff), bShake(true), bVerbose(false) {
|
8 |
dumpOut = NULL;
|
9 |
statOut = NULL;
|
10 |
|
11 |
tStats = new Thermo(info);
|
12 |
|
13 |
|
14 |
paramSet = param;
|
15 |
|
16 |
calcDim();
|
17 |
|
18 |
curX = getCoor();
|
19 |
curG.resize(ndim);
|
20 |
|
21 |
preMove();
|
22 |
}
|
23 |
|
24 |
OOPSEMinimizer::~OOPSEMinimizer(){
|
25 |
delete tStats;
|
26 |
if(dumpOut)
|
27 |
delete dumpOut;
|
28 |
if(statOut)
|
29 |
delete statOut;
|
30 |
delete paramSet;
|
31 |
}
|
32 |
|
33 |
void OOPSEMinimizer::calcEnergyGradient(vector<double>& x, vector<double>& grad,
|
34 |
double& energy, int& status){
|
35 |
|
36 |
DirectionalAtom* dAtom;
|
37 |
int index;
|
38 |
double force[3];
|
39 |
double dAtomGrad[6];
|
40 |
int shakeStatus;
|
41 |
|
42 |
status = 1;
|
43 |
|
44 |
setCoor(x);
|
45 |
|
46 |
if (nConstrained && bShake){
|
47 |
shakeStatus = shakeR();
|
48 |
}
|
49 |
|
50 |
calcForce(1, 1);
|
51 |
|
52 |
if (nConstrained && bShake){
|
53 |
shakeStatus = shakeF();
|
54 |
}
|
55 |
|
56 |
x = getCoor();
|
57 |
|
58 |
|
59 |
index = 0;
|
60 |
|
61 |
for(int i = 0; i < integrableObjects.size(); i++){
|
62 |
|
63 |
if (integrableObjects[i]->isDirectional()) {
|
64 |
|
65 |
integrableObjects[i]->getGrad(dAtomGrad);
|
66 |
|
67 |
//gradient is equal to -f
|
68 |
grad[index++] = -dAtomGrad[0];
|
69 |
grad[index++] = -dAtomGrad[1];
|
70 |
grad[index++] = -dAtomGrad[2];
|
71 |
grad[index++] = -dAtomGrad[3];
|
72 |
grad[index++] = -dAtomGrad[4];
|
73 |
grad[index++] = -dAtomGrad[5];
|
74 |
|
75 |
}
|
76 |
else{
|
77 |
integrableObjects[i]->getFrc(force);
|
78 |
|
79 |
grad[index++] = -force[0];
|
80 |
grad[index++] = -force[1];
|
81 |
grad[index++] = -force[2];
|
82 |
|
83 |
}
|
84 |
|
85 |
}
|
86 |
|
87 |
energy = tStats->getPotential();
|
88 |
|
89 |
}
|
90 |
|
91 |
void OOPSEMinimizer::setCoor(vector<double>& x){
|
92 |
|
93 |
DirectionalAtom* dAtom;
|
94 |
int index;
|
95 |
double position[3];
|
96 |
double eulerAngle[3];
|
97 |
|
98 |
|
99 |
index = 0;
|
100 |
|
101 |
for(int i = 0; i < integrableObjects.size(); i++){
|
102 |
|
103 |
position[0] = x[index++];
|
104 |
position[1] = x[index++];
|
105 |
position[2] = x[index++];
|
106 |
|
107 |
integrableObjects[i]->setPos(position);
|
108 |
|
109 |
if (integrableObjects[i]->isDirectional()){
|
110 |
|
111 |
eulerAngle[0] = x[index++];
|
112 |
eulerAngle[1] = x[index++];
|
113 |
eulerAngle[2] = x[index++];
|
114 |
|
115 |
integrableObjects[i]->setEuler(eulerAngle[0],
|
116 |
eulerAngle[1],
|
117 |
eulerAngle[2]);
|
118 |
|
119 |
}
|
120 |
|
121 |
}
|
122 |
|
123 |
}
|
124 |
|
125 |
vector<double> OOPSEMinimizer::getCoor(){
|
126 |
|
127 |
DirectionalAtom* dAtom;
|
128 |
int index;
|
129 |
double position[3];
|
130 |
double eulerAngle[3];
|
131 |
vector<double> x;
|
132 |
|
133 |
x.resize(getDim());
|
134 |
|
135 |
index = 0;
|
136 |
|
137 |
for(int i = 0; i < integrableObjects.size(); i++){
|
138 |
position = integrableObjects[i]->getPos();
|
139 |
|
140 |
x[index++] = position[0];
|
141 |
x[index++] = position[1];
|
142 |
x[index++] = position[2];
|
143 |
|
144 |
if (integrableObjects[i]->isDirectional()){
|
145 |
|
146 |
integrableObjects[i]->getEulerAngles(eulerAngle);
|
147 |
|
148 |
x[index++] = eulerAngle[0];
|
149 |
x[index++] = eulerAngle[1];
|
150 |
x[index++] = eulerAngle[2];
|
151 |
|
152 |
}
|
153 |
|
154 |
}
|
155 |
|
156 |
return x;
|
157 |
|
158 |
}
|
159 |
|
160 |
|
161 |
int OOPSEMinimizer::shakeR(){
|
162 |
int i, j;
|
163 |
int done;
|
164 |
double posA[3], posB[3];
|
165 |
double velA[3], velB[3];
|
166 |
double pab[3];
|
167 |
double rab[3];
|
168 |
int a, b, ax, ay, az, bx, by, bz;
|
169 |
double rma, rmb;
|
170 |
double dx, dy, dz;
|
171 |
double rpab;
|
172 |
double rabsq, pabsq, rpabsq;
|
173 |
double diffsq;
|
174 |
double gab;
|
175 |
int iteration;
|
176 |
|
177 |
for (i = 0; i < nAtoms; i++){
|
178 |
moving[i] = 0;
|
179 |
moved[i] = 1;
|
180 |
}
|
181 |
|
182 |
iteration = 0;
|
183 |
done = 0;
|
184 |
while (!done && (iteration < maxIteration)){
|
185 |
done = 1;
|
186 |
for (i = 0; i < nConstrained; i++){
|
187 |
a = constrainedA[i];
|
188 |
b = constrainedB[i];
|
189 |
|
190 |
ax = (a * 3) + 0;
|
191 |
ay = (a * 3) + 1;
|
192 |
az = (a * 3) + 2;
|
193 |
|
194 |
bx = (b * 3) + 0;
|
195 |
by = (b * 3) + 1;
|
196 |
bz = (b * 3) + 2;
|
197 |
|
198 |
if (moved[a] || moved[b]){
|
199 |
posA = atoms[a]->getPos();
|
200 |
posB = atoms[b]->getPos();
|
201 |
|
202 |
for (j = 0; j < 3; j++)
|
203 |
pab[j] = posA[j] - posB[j];
|
204 |
|
205 |
//periodic boundary condition
|
206 |
|
207 |
info->wrapVector(pab);
|
208 |
|
209 |
pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
|
210 |
|
211 |
rabsq = constrainedDsqr[i];
|
212 |
diffsq = rabsq - pabsq;
|
213 |
|
214 |
// the original rattle code from alan tidesley
|
215 |
if (fabs(diffsq) > (tol * rabsq * 2)){
|
216 |
rab[0] = oldPos[ax] - oldPos[bx];
|
217 |
rab[1] = oldPos[ay] - oldPos[by];
|
218 |
rab[2] = oldPos[az] - oldPos[bz];
|
219 |
|
220 |
info->wrapVector(rab);
|
221 |
|
222 |
rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
|
223 |
|
224 |
rpabsq = rpab * rpab;
|
225 |
|
226 |
|
227 |
if (rpabsq < (rabsq * -diffsq)){
|
228 |
#ifdef IS_MPI
|
229 |
a = atoms[a]->getGlobalIndex();
|
230 |
b = atoms[b]->getGlobalIndex();
|
231 |
#endif //is_mpi
|
232 |
//cerr << "Waring: constraint failure" << endl;
|
233 |
gab = sqrt(rabsq/pabsq);
|
234 |
rab[0] = (posA[0] - posB[0])*gab;
|
235 |
rab[1]= (posA[1] - posB[1])*gab;
|
236 |
rab[2] = (posA[2] - posB[2])*gab;
|
237 |
|
238 |
info->wrapVector(rab);
|
239 |
|
240 |
rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
|
241 |
|
242 |
}
|
243 |
|
244 |
//rma = 1.0 / atoms[a]->getMass();
|
245 |
//rmb = 1.0 / atoms[b]->getMass();
|
246 |
rma = 1.0;
|
247 |
rmb =1.0;
|
248 |
|
249 |
gab = diffsq / (2.0 * (rma + rmb) * rpab);
|
250 |
|
251 |
dx = rab[0] * gab;
|
252 |
dy = rab[1] * gab;
|
253 |
dz = rab[2] * gab;
|
254 |
|
255 |
posA[0] += rma * dx;
|
256 |
posA[1] += rma * dy;
|
257 |
posA[2] += rma * dz;
|
258 |
|
259 |
atoms[a]->setPos(posA);
|
260 |
|
261 |
posB[0] -= rmb * dx;
|
262 |
posB[1] -= rmb * dy;
|
263 |
posB[2] -= rmb * dz;
|
264 |
|
265 |
atoms[b]->setPos(posB);
|
266 |
|
267 |
moving[a] = 1;
|
268 |
moving[b] = 1;
|
269 |
done = 0;
|
270 |
}
|
271 |
}
|
272 |
}
|
273 |
|
274 |
for (i = 0; i < nAtoms; i++){
|
275 |
moved[i] = moving[i];
|
276 |
moving[i] = 0;
|
277 |
}
|
278 |
|
279 |
iteration++;
|
280 |
}
|
281 |
|
282 |
if (!done){
|
283 |
cerr << "Waring: can not constraint within maxIteration" << endl;
|
284 |
return -1;
|
285 |
}
|
286 |
else
|
287 |
return 1;
|
288 |
}
|
289 |
|
290 |
|
291 |
//remove constraint force along the bond direction
|
292 |
int OOPSEMinimizer::shakeF(){
|
293 |
int i, j;
|
294 |
int done;
|
295 |
double posA[3], posB[3];
|
296 |
double frcA[3], frcB[3];
|
297 |
double rab[3], fpab[3];
|
298 |
int a, b, ax, ay, az, bx, by, bz;
|
299 |
double rma, rmb;
|
300 |
double rvab;
|
301 |
double gab;
|
302 |
double rabsq;
|
303 |
double rfab;
|
304 |
int iteration;
|
305 |
|
306 |
for (i = 0; i < nAtoms; i++){
|
307 |
moving[i] = 0;
|
308 |
moved[i] = 1;
|
309 |
}
|
310 |
|
311 |
done = 0;
|
312 |
iteration = 0;
|
313 |
while (!done && (iteration < maxIteration)){
|
314 |
done = 1;
|
315 |
|
316 |
for (i = 0; i < nConstrained; i++){
|
317 |
a = constrainedA[i];
|
318 |
b = constrainedB[i];
|
319 |
|
320 |
ax = (a * 3) + 0;
|
321 |
ay = (a * 3) + 1;
|
322 |
az = (a * 3) + 2;
|
323 |
|
324 |
bx = (b * 3) + 0;
|
325 |
by = (b * 3) + 1;
|
326 |
bz = (b * 3) + 2;
|
327 |
|
328 |
if (moved[a] || moved[b]){
|
329 |
|
330 |
posA = atoms[a]->getPos();
|
331 |
posB = atoms[b]->getPos();
|
332 |
|
333 |
for (j = 0; j < 3; j++)
|
334 |
rab[j] = posA[j] - posB[j];
|
335 |
|
336 |
info->wrapVector(rab);
|
337 |
|
338 |
atoms[a]->getFrc(frcA);
|
339 |
atoms[b]->getFrc(frcB);
|
340 |
|
341 |
//rma = 1.0 / atoms[a]->getMass();
|
342 |
//rmb = 1.0 / atoms[b]->getMass();
|
343 |
rma = 1.0;
|
344 |
rmb = 1.0;
|
345 |
|
346 |
|
347 |
fpab[0] = frcA[0] * rma - frcB[0] * rmb;
|
348 |
fpab[1] = frcA[1] * rma - frcB[1] * rmb;
|
349 |
fpab[2] = frcA[2] * rma - frcB[2] * rmb;
|
350 |
|
351 |
|
352 |
gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
|
353 |
|
354 |
if (gab < 1.0)
|
355 |
gab = 1.0;
|
356 |
|
357 |
rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
|
358 |
rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
|
359 |
|
360 |
if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001){
|
361 |
|
362 |
gab = -rfab /(rabsq*(rma + rmb));
|
363 |
|
364 |
frcA[0] = rab[0] * gab;
|
365 |
frcA[1] = rab[1] * gab;
|
366 |
frcA[2] = rab[2] * gab;
|
367 |
|
368 |
atoms[a]->addFrc(frcA);
|
369 |
|
370 |
|
371 |
frcB[0] = -rab[0] * gab;
|
372 |
frcB[1] = -rab[1] * gab;
|
373 |
frcB[2] = -rab[2] * gab;
|
374 |
|
375 |
atoms[b]->addFrc(frcB);
|
376 |
|
377 |
moving[a] = 1;
|
378 |
moving[b] = 1;
|
379 |
done = 0;
|
380 |
}
|
381 |
}
|
382 |
}
|
383 |
|
384 |
for (i = 0; i < nAtoms; i++){
|
385 |
moved[i] = moving[i];
|
386 |
moving[i] = 0;
|
387 |
}
|
388 |
|
389 |
iteration++;
|
390 |
}
|
391 |
|
392 |
if (!done){
|
393 |
cerr << "Waring: can not constraint within maxIteration" << endl;
|
394 |
return -1;
|
395 |
}
|
396 |
else
|
397 |
return 1;
|
398 |
}
|
399 |
|
400 |
|
401 |
|
402 |
//calculate the value of object function
|
403 |
void OOPSEMinimizer::calcF(){
|
404 |
calcEnergyGradient(curX, curG, curF, egEvalStatus);
|
405 |
}
|
406 |
|
407 |
void OOPSEMinimizer::calcF(vector<double>& x, double&f, int& status){
|
408 |
vector<double> tempG;
|
409 |
tempG.resize(x.size());
|
410 |
|
411 |
calcEnergyGradient(x, tempG, f, status);
|
412 |
}
|
413 |
|
414 |
//calculate the gradient
|
415 |
void OOPSEMinimizer::calcG(){
|
416 |
calcEnergyGradient(curX, curG, curF, egEvalStatus);
|
417 |
}
|
418 |
|
419 |
void OOPSEMinimizer::calcG(vector<double>& x, vector<double>& g, double& f, int& status){
|
420 |
calcEnergyGradient(x, g, f, status);
|
421 |
}
|
422 |
|
423 |
void OOPSEMinimizer::calcDim(){
|
424 |
DirectionalAtom* dAtom;
|
425 |
|
426 |
ndim = 0;
|
427 |
|
428 |
for(int i = 0; i < integrableObjects.size(); i++){
|
429 |
ndim += 3;
|
430 |
if (integrableObjects[i]->isDirectional())
|
431 |
ndim += 3;
|
432 |
}
|
433 |
}
|
434 |
|
435 |
void OOPSEMinimizer::setX(vector < double > & x){
|
436 |
|
437 |
if (x.size() != ndim && bVerbose){
|
438 |
//sprintf(painCave.errMsg,
|
439 |
// "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
|
440 |
// painCave.isFatal = 1;
|
441 |
// simError();
|
442 |
}
|
443 |
|
444 |
curX = x;
|
445 |
}
|
446 |
|
447 |
void OOPSEMinimizer::setG(vector < double > & g){
|
448 |
|
449 |
if (g.size() != ndim && bVerbose){
|
450 |
//sprintf(painCave.errMsg,
|
451 |
// "OOPSEMinimizer Error: dimesion of g and curG does not match\n");
|
452 |
// painCave.isFatal = 1;
|
453 |
//simError();
|
454 |
}
|
455 |
|
456 |
curG = g;
|
457 |
}
|
458 |
|
459 |
void OOPSEMinimizer::writeOut(vector<double>& x, double iter){
|
460 |
|
461 |
setX(x);
|
462 |
|
463 |
calcG();
|
464 |
|
465 |
dumpOut->writeDump(iter);
|
466 |
statOut->writeStat(iter);
|
467 |
}
|
468 |
|
469 |
|
470 |
void OOPSEMinimizer::printMinimizerInfo(){
|
471 |
cout << "--------------------------------------------------------------------" << endl;
|
472 |
cout << minimizerName << endl;
|
473 |
cout << "minimization parameter set" << endl;
|
474 |
cout << "function tolerance = " << paramSet->getFTol() << endl;
|
475 |
cout << "gradient tolerance = " << paramSet->getGTol() << endl;
|
476 |
cout << "step tolerance = "<< paramSet->getFTol() << endl;
|
477 |
cout << "absolute gradient tolerance = " << endl;
|
478 |
cout << "max iteration = " << paramSet->getMaxIteration() << endl;
|
479 |
cout << "max line search iteration = " << paramSet->getLineSearchMaxIteration() <<endl;
|
480 |
cout << "shake algorithm = " << bShake << endl;
|
481 |
cout << "--------------------------------------------------------------------" << endl;
|
482 |
|
483 |
}
|
484 |
|
485 |
/**
|
486 |
* In thoery, we need to find the minimum along the search direction
|
487 |
* However, function evaluation is too expensive.
|
488 |
* At the very begining of the problem, we check the search direction and make sure
|
489 |
* it is a descent direction
|
490 |
* we will compare the energy of two end points,
|
491 |
* if the right end point has lower energy, we just take it
|
492 |
*
|
493 |
*
|
494 |
*
|
495 |
*/
|
496 |
|
497 |
int OOPSEMinimizer::doLineSearch(vector<double>& direction, double stepSize){
|
498 |
vector<double> xa;
|
499 |
vector<double> xb;
|
500 |
vector<double> xc;
|
501 |
vector<double> ga;
|
502 |
vector<double> gb;
|
503 |
vector<double> gc;
|
504 |
double fa;
|
505 |
double fb;
|
506 |
double fc;
|
507 |
double a;
|
508 |
double b;
|
509 |
double c;
|
510 |
int status;
|
511 |
double initSlope;
|
512 |
double slopeA;
|
513 |
double slopeB;
|
514 |
double slopeC;
|
515 |
bool foundLower;
|
516 |
int iter;
|
517 |
int maxLSIter;
|
518 |
double mu;
|
519 |
double eta;
|
520 |
double ftol;
|
521 |
double lsTol;
|
522 |
|
523 |
xa.resize(ndim);
|
524 |
xb.resize(ndim);
|
525 |
xc.resize(ndim);
|
526 |
|
527 |
ga.resize(ndim);
|
528 |
gb.resize(ndim);
|
529 |
gc.resize(ndim);
|
530 |
|
531 |
a = 0.0;
|
532 |
fa = curF;
|
533 |
xa = curX;
|
534 |
ga = curG;
|
535 |
c = a + stepSize;
|
536 |
ftol = paramSet->getFTol();
|
537 |
lsTol = paramSet->getLineSearchTol();
|
538 |
|
539 |
//calculate the derivative at a = 0
|
540 |
slopeA = 0;
|
541 |
for (size_t i = 0; i < ndim; i++)
|
542 |
slopeA += curG[i]*direction[i];
|
543 |
|
544 |
initSlope = slopeA;
|
545 |
|
546 |
// if going uphill, use negative gradient as searching direction
|
547 |
if (slopeA > 0) {
|
548 |
|
549 |
if (bVerbose){
|
550 |
cout << "LineSearch Warning: initial searching direction is not a descent searching direction, "
|
551 |
<< " use negative gradient instead. Therefore, finding a smaller vaule of function "
|
552 |
<< " is guaranteed"
|
553 |
<< endl;
|
554 |
}
|
555 |
|
556 |
for (size_t i = 0; i < ndim; i++)
|
557 |
direction[i] = -curG[i];
|
558 |
|
559 |
for (size_t i = 0; i < ndim; i++)
|
560 |
slopeA += curG[i]*direction[i];
|
561 |
|
562 |
initSlope = slopeA;
|
563 |
}
|
564 |
|
565 |
// Take a trial step
|
566 |
for(size_t i = 0; i < ndim; i++)
|
567 |
xc[i] = curX[i] + direction[i] * c;
|
568 |
|
569 |
calcG(xc, gc, fc, status);
|
570 |
|
571 |
if (status < 0){
|
572 |
if (bVerbose)
|
573 |
cerr << "Function Evaluation Error" << endl;
|
574 |
}
|
575 |
|
576 |
//calculate the derivative at c
|
577 |
slopeC = 0;
|
578 |
for (size_t i = 0; i < ndim; i++)
|
579 |
slopeC += gc[i]*direction[i];
|
580 |
|
581 |
// found a lower point
|
582 |
if (fc < fa) {
|
583 |
curX = xc;
|
584 |
curG = gc;
|
585 |
curF = fc;
|
586 |
return LS_SUCCEED;
|
587 |
}
|
588 |
else {
|
589 |
|
590 |
if (slopeC > 0)
|
591 |
stepSize *= 0.618034;
|
592 |
}
|
593 |
|
594 |
maxLSIter = paramSet->getLineSearchMaxIteration();
|
595 |
|
596 |
iter = 0;
|
597 |
|
598 |
do {
|
599 |
// Select a new trial point.
|
600 |
// If the derivatives at points a & c have different sign we use cubic interpolate
|
601 |
//if (slopeC > 0){
|
602 |
eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
|
603 |
mu = sqrt(eta * eta - slopeA * slopeC);
|
604 |
b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));
|
605 |
|
606 |
if (b < lsTol){
|
607 |
if (bVerbose)
|
608 |
cout << "stepSize is less than line search tolerance" << endl;
|
609 |
break;
|
610 |
}
|
611 |
//}
|
612 |
|
613 |
// Take a trial step to this new point - new coords in xb
|
614 |
for(size_t i = 0; i < ndim; i++)
|
615 |
xb[i] = curX[i] + direction[i] * b;
|
616 |
|
617 |
//function evaluation
|
618 |
calcG(xb, gb, fb, status);
|
619 |
|
620 |
if (status < 0){
|
621 |
if (bVerbose)
|
622 |
cerr << "Function Evaluation Error" << endl;
|
623 |
}
|
624 |
|
625 |
//calculate the derivative at c
|
626 |
slopeB = 0;
|
627 |
for (size_t i = 0; i < ndim; i++)
|
628 |
slopeB += gb[i]*direction[i];
|
629 |
|
630 |
//Amijo Rule to stop the line search
|
631 |
if (fb <= curF + initSlope * ftol * b) {
|
632 |
curF = fb;
|
633 |
curX = xb;
|
634 |
curG = gb;
|
635 |
return LS_SUCCEED;
|
636 |
}
|
637 |
|
638 |
if (slopeB <0 && fb < fa) {
|
639 |
//replace a by b
|
640 |
fa = fb;
|
641 |
a = b;
|
642 |
slopeA = slopeB;
|
643 |
|
644 |
// swap coord a/b
|
645 |
std::swap(xa, xb);
|
646 |
std::swap(ga, gb);
|
647 |
}
|
648 |
else {
|
649 |
//replace c by b
|
650 |
fc = fb;
|
651 |
c = b;
|
652 |
slopeC = slopeB;
|
653 |
|
654 |
// swap coord b/c
|
655 |
std::swap(gb, gc);
|
656 |
std::swap(xb, xc);
|
657 |
}
|
658 |
|
659 |
|
660 |
iter++;
|
661 |
} while((fb > fa || fb > fc) && (iter < maxLSIter));
|
662 |
|
663 |
if(fb < curF || iter >= maxLSIter) {
|
664 |
//could not find a lower value, we might just go uphill.
|
665 |
return LS_ERROR;
|
666 |
}
|
667 |
|
668 |
//select the end point
|
669 |
|
670 |
if (fa <= fc) {
|
671 |
curX = xa;
|
672 |
curG = ga;
|
673 |
curF = fa;
|
674 |
}
|
675 |
else {
|
676 |
curX = xc;
|
677 |
curG = gc;
|
678 |
curF = fc;
|
679 |
}
|
680 |
|
681 |
return LS_SUCCEED;
|
682 |
|
683 |
}
|
684 |
|
685 |
void OOPSEMinimizer::minimize(){
|
686 |
|
687 |
int convgStatus;
|
688 |
int stepStatus;
|
689 |
int maxIter;
|
690 |
int writeFrq;
|
691 |
int nextWriteIter;
|
692 |
|
693 |
if (bVerbose)
|
694 |
printMinimizerInfo();
|
695 |
|
696 |
dumpOut = new DumpWriter(info);
|
697 |
statOut = new StatWriter(info);
|
698 |
|
699 |
init();
|
700 |
|
701 |
writeFrq = paramSet->getWriteFrq();
|
702 |
nextWriteIter = writeFrq;
|
703 |
|
704 |
maxIter = paramSet->getMaxIteration();
|
705 |
|
706 |
for (curIter = 1; curIter <= maxIter; curIter++){
|
707 |
|
708 |
stepStatus = step();
|
709 |
|
710 |
if (nConstrained && bShake)
|
711 |
preMove();
|
712 |
|
713 |
if (stepStatus < 0){
|
714 |
saveResult();
|
715 |
minStatus = MIN_LSERROR;
|
716 |
cerr << "OOPSEMinimizer Error: line search error, please try a small stepsize" << endl;
|
717 |
return;
|
718 |
}
|
719 |
|
720 |
if (curIter == nextWriteIter){
|
721 |
nextWriteIter += writeFrq;
|
722 |
writeOut(curX, curIter);
|
723 |
}
|
724 |
|
725 |
convgStatus = checkConvg();
|
726 |
|
727 |
if (convgStatus > 0){
|
728 |
saveResult();
|
729 |
minStatus = MIN_CONVERGE;
|
730 |
return;
|
731 |
}
|
732 |
|
733 |
prepareStep();
|
734 |
|
735 |
}
|
736 |
|
737 |
|
738 |
if (bVerbose) {
|
739 |
cout << "OOPSEMinimizer Warning: "
|
740 |
<< minimizerName << " algorithm did not converge within "
|
741 |
<< maxIter << " iteration" << endl;
|
742 |
}
|
743 |
minStatus = MIN_MAXITER;
|
744 |
saveResult();
|
745 |
|
746 |
}
|