1
package org.openscience.cdk.modeling.forcefield;
3
import javax.vecmath.GVector;
5
//import org.openscience.cdk.tools.LoggingTool;
10
*@cdk.module forcefield
13
public class LineSearchForTheWolfeConditions {
16
private GVector x = null;
17
private double linearFunctionInAlpha0;
18
private GVector dfx = null;
19
private GVector direction = null;
20
private double linearFunctionDerivativeInAlpha0;
21
private IPotentialFunction pf = null;
22
private double alphaInitialStep;
25
//line search algorithm
26
private double[] alpha = new double[3];
27
private double[] linearFunctionInAlpha = new double[3];
28
private double[] linearFunctionDerivativeInAlpha = new double[3];
30
private GVector[] dfInAlpha = new GVector[3];
31
private double[] brentStep = new double[3];
33
private final double c1 = 0.0001;
36
//private double linearFunctionGoldenAlpha;
37
private double linearFunctionAlphaInterpolation;
39
public boolean derivativeSmallEnough = true;
41
public double alphaOptimum;
42
public double linearFunctionInAlphaOptimum;
43
public GVector dfOptimum = null;
46
private double alphaj;
47
private double linearFunctionInAlphaj;
48
private double linearFunctionDerivativeInAlphaj;
49
private GVector dfInAlphaj;
50
private int functionEvaluationNumber;
52
//energy function evaluation
53
private GVector xAlpha = null;
60
private double alphaTemporal;
61
private double linearFunctionInAlphaTemporal;
62
private double linearFunctionDerivativeInAlphaTemporal;
65
private double alphaiplus1;
67
//private LoggingTool logger;
70
public LineSearchForTheWolfeConditions(IPotentialFunction pfUser, String method) {
72
if ((method == "sdm") | (method == "cgm")) {c2 = 0.07;}
76
public void initialize(GVector xUser, double fxUser, GVector dfxUser, GVector directionUser, double linearFunctionDerivativeUser, double alphaInitialStepUser) {
78
this.linearFunctionInAlpha0 = fxUser;
80
//logger.debug("derivativeSmallEnough = " + this.derivativeSmallEnough);
81
this.direction = directionUser;
82
this.linearFunctionDerivativeInAlpha0 = linearFunctionDerivativeUser;
83
//logger.debug("linearFunctionDerivativeInAlpha0 = " + linearFunctionDerivativeInAlpha0);
84
this.alphaOptimum = 0;
85
this.linearFunctionInAlphaOptimum = linearFunctionInAlpha0;
87
this.alphaInitialStep = alphaInitialStepUser;
88
this.derivativeSmallEnough = false;
89
this.xAlpha = new GVector(x.getSize());
92
/* Line Search Algorithm for the Wolfe conditions. Jorge Nocedal and Stephen J.Wright. Numerical Optimization. 1999.
93
* The algorithm has two stages. This first stage begins with a trial estimate alpha1,
94
* and keeps increasing it until it finds either an acceptable step length or an interval
95
* that brackets the desired step lengths. In the later case, the second stage is invoked
96
* by calling a function called zoom, which successively decreases the size of the interval
97
* until an acceptable step length is identified.
99
* @param alphaMax Maximum step length
101
public void lineSearchAlgorithm (double alphaMax) {
103
//logger.debug("Line search for the strong wolfe conditions");
106
linearFunctionInAlpha[0] = linearFunctionInAlpha0;
107
linearFunctionDerivativeInAlpha[0] = linearFunctionDerivativeInAlpha0; //To Analyse the possibility of eliminate linearFunctionDerivativeInAlpha[0]
108
dfInAlpha[0] = this.dfx;
110
alpha[1] = this.alphaInitialStep;
112
//logger.debug("alpha[1] = this.alphaInitialStep = " + alpha[1]);
114
brentStep[0] = alpha[0];
115
brentStep[1] = alpha[1];
119
this.functionEvaluationNumber = 0;
121
if (alpha[1] > alphaMax) {
123
//logger.debug("line search algorithm error: alphaInitialStep > alphaMax");
125
// alpha[1] = alphaMax/2;
132
System.out.println("alpha[1] == 0");
136
//logger.debug("alpha[" + i + "] = " + alpha[i]);
137
linearFunctionInAlpha[i] = evaluateEnergyFunction(alpha[i]);
138
//logger.debug("linearFunctionInAlpha[" + i + "] = " + linearFunctionInAlpha[i]);
140
if ((linearFunctionInAlpha[i] > linearFunctionInAlpha[0] + c1 * alpha[i] * linearFunctionDerivativeInAlpha[0]) |
141
((linearFunctionInAlpha[i] >= linearFunctionInAlpha[i-1]) & (i>1))) { //The interval alpha[i-1] and alpha[i] brackets the desired step lengths.
142
//logger.debug("zoom(" + alpha[i-1] + ", " + linearFunctionInAlpha[i-1] + ", " + linearFunctionDerivativeInAlpha[i-1] + ", " + dfInAlpha[i-1] + ", " + alpha[i] + ", " + linearFunctionInAlpha[i] + ")");
143
//dfInAlpha[i] = evaluateEnergyFunctionDerivative(alpha[i]);
144
//linearFunctionDerivativeInAlpha[i] = dfInAlpha[i].dot(direction);
145
//zoom(alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i-1], dfInAlpha[i-1], alpha[i], linearFunctionInAlpha[i], linearFunctionDerivativeInAlpha[i], dfInAlpha[i]);
146
zoom(alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i-1], dfInAlpha[i-1], alpha[i], linearFunctionInAlpha[i]);
150
//The first strong Wolfe condition is satisfied for alpha[i].
151
dfInAlpha[i] = evaluateEnergyFunctionDerivative(alpha[i]);
152
//logger.debug("dfOptimum = " + dfOptimum);
153
linearFunctionDerivativeInAlpha[i] = dfInAlpha[i].dot(direction);
154
//logger.debug("linearFunctionDerivativeInAlpha[" + i + "] = " + linearFunctionDerivativeInAlpha[i]);
156
if (Math.abs(linearFunctionDerivativeInAlpha[i]) <= -c2 * linearFunctionDerivativeInAlpha[0]) { //The second strong Wolfe condition is also satisfied for alpha[i]
157
//logger.debug("The second strong Wolfe condition is also satisfied for " + alpha[i]);
158
alphaOptimum = alpha[i];
159
linearFunctionInAlphaOptimum = linearFunctionInAlpha[i];
160
dfOptimum = dfInAlpha[i];
161
//logger.debug("alphaOptimun = " + alphaOptimum);
162
//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);
163
//logger.debug("dfOptimum = " + dfOptimum);
164
this.derivativeSmallEnough = true;
168
if (linearFunctionDerivativeInAlpha[i] >= 0) { //The interval alpha[i-1] and alpha[i] brackets the desired step lengths.
169
/*System.out.println("zoom(" + alpha[i-1] + ", " + linearFunctionInAlpha[i-1] + ", " + linearFunctionDerivativeInAlpha[i-1] + ", " + dfInAlpha[i-1] + ", " +
170
alpha[i] + ", " + linearFunctionInAlpha[i] + ")");*/
172
/*zoom(alpha[i], linearFunctionInAlpha[i], linearFunctionDerivativeInAlpha[i], dfInAlpha[i],
173
alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i], dfInAlpha[i]);*/
174
zoom(alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i-1], dfInAlpha[i-1],
175
alpha[i], linearFunctionInAlpha[i]);
179
if (alpha[i] == alphaMax) {
180
//logger.debug("LINE SEARCH ALGORITHM WAS TERMINATE EARLIER BECAUSE alpha[i] == alphaMax");
181
alphaOptimum = alpha[i];
182
linearFunctionInAlphaOptimum = linearFunctionInAlpha[i];
183
dfOptimum = dfInAlpha[i];
184
//logger.debug("alphaOptimun = " + alphaOptimum);
185
//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);
186
//logger.debug("dfOptimum = " + dfOptimum);
190
functionEvaluationNumber = functionEvaluationNumber + 1;
191
if (functionEvaluationNumber == 10) {
192
//logger.debug("LINE SEARCH ALGORITHM WAS TERMINATE EARLIER BECAUSE functionEvaluationNumber == 10");
193
alphaOptimum = alpha[i];
194
linearFunctionInAlphaOptimum = linearFunctionInAlpha[i];
195
dfOptimum = dfInAlpha[i];
196
//logger.debug("alphaOptimun = " + alphaOptimum);
197
//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);
198
//logger.debug("dfOptimum = " + dfOptimum);
203
brentStep[0] = brentStep[1];
204
brentStep[1] = brentStep[2];
207
linearFunctionInAlpha[1] = linearFunctionInAlpha[2];
208
linearFunctionDerivativeInAlpha[1] = linearFunctionDerivativeInAlpha[2];
209
dfInAlpha[1] = dfInAlpha[2];
212
brentStep[2] = brentStep[1] + 1.618 * (brentStep[1]-brentStep[0]);
213
//logger.debug("brentStep[2] = " + brentStep[2]);
215
if (brentStep[2] > alphaMax) {brentStep[2] = alphaMax;}
216
/*linearFunctionInBrentStep = this.evaluateEnergyFunction(brentStep[2]);
217
linearFunctionDerivativeInBrentStep = this.evaluateEnergyFunctionDerivative(brentStep[2]).dot(direction);
219
alpha[2] = brentStep[2];
220
/*alpha[2] = this.cubicInterpolation(alpha[1], linearFunctionInAlpha[1], linearFunctionDerivativeInAlpha[1],
221
brentStep[2], linearFunctionInBrentStep, linearFunctionDerivativeInBrentStep, alpha[1], brentStep[2]);
226
} while ((alpha[2] <= alphaMax) & (alpha[1] < alpha[2]) & (functionEvaluationNumber < 10));
228
} catch (Exception exception) {
229
System.out.println("Line search for the strong wolfe conditions: " + exception.getMessage());
230
System.out.println(exception);
237
/* Each iteration of zoom generates an iterate alphaj between alphaLow and alphaHigh,
238
* and then replaces one of these endpoints by alphaj in such a way that the properties
239
* (a), (b) and (c) continue to hold.
240
* (a)The interval bounded by alphaLow and alphaHigh contains step lengths that satisfy the strong Wolfe conditions.
241
* (b)alphaLow is, among all step lengths generated so far and satisfying the sufficient decrease condition,
242
* the one giving the smallest function value.
243
* (c)alphaHigh is chosen so that linearFunctionDerivativeInAlphaj * (alphaHigh-alphaLow) < 0
245
*@param alphaLow Among all step lengths generated so far and satisfying the sufficient decrease condition, the one giving the smallest function value.
246
*@param linearFunctionInAlphaLow Function value at alphaLow.
247
*@param linearFunctionDerivativeInAlphaLow Derivative value at alphaLow.
248
*@param dfInAlphaLow Gradient at alphaLow.
249
*@param alphaHigh AlphaHigh is chosen so that linearFunctionDerivativeInAlphaj * (alphaHigh-alphaLow) < 0
250
*@param linearFunctionInAlphaHigh Function value at alphaHigh.
252
private void zoom (double alphaLow, double linearFunctionInAlphaLow, double linearFunctionDerivativeInAlphaLow, GVector dfInAlphaLow,
253
double alphaHigh, double linearFunctionInAlphaHigh) {
255
//logger.debug("zoom");
257
functionEvaluationNumber = 0;
261
if (alphaLow < alphaHigh) {a = alphaLow; b = alphaHigh;}
262
else {a = alphaHigh; b = alphaLow;}
268
//alphaj = this.cubicInterpolation(alphaLow, linearFunctionInAlphaLow, linearFunctionDerivativeInAlphaLow, alphaHigh, linearFunctionInAlphaHigh, linearFunctionDerivativeInAlphaHigh, a, b);
269
/*System.out.println("interpolation(" + alphaLow + ", " + linearFunctionInAlphaLow + ", " + linearFunctionDerivativeInAlphaLow + ", "
270
+ alphaHigh + ", " + linearFunctionInAlphaHigh + ");");*/
272
alphaj = this.interpolation(alphaLow, linearFunctionInAlphaLow, linearFunctionDerivativeInAlphaLow, alphaHigh, linearFunctionInAlphaHigh);
273
//logger.debug("alphaj = " + alphaj);
274
linearFunctionInAlphaj = this.linearFunctionAlphaInterpolation;
275
//logger.debug("linearFunctionInAlphaj = " + linearFunctionInAlphaj);
277
if ((linearFunctionInAlphaj > linearFunctionInAlpha0 + c1 * alphaj * linearFunctionDerivativeInAlpha0) | //The interval 0 and alphaj brackets the desired step lengths.
278
(linearFunctionInAlphaj >= linearFunctionInAlphaLow)) {
280
//logger.debug("The minimum is between alpha1 and alphaj");
282
linearFunctionInAlphaHigh = linearFunctionInAlphaj;
283
//dfInAlphaHigh = this.evaluateEnergyFunctionDerivative(alphaHigh);
284
//linearFunctionDerivativeInAlphaHigh = dfInAlphaHigh.dot(direction);
287
dfInAlphaj = evaluateEnergyFunctionDerivative(alphaj);
288
linearFunctionDerivativeInAlphaj = dfInAlphaj.dot(direction);
289
//logger.debug("linearFunctionDerivativeInAlphaj = " + linearFunctionDerivativeInAlphaj);
290
if (Math.abs(linearFunctionDerivativeInAlphaj) <= -c2 * linearFunctionDerivativeInAlpha0) { //alphaj satisfied the second strong Wolfe condition.
291
//logger.debug("Derivative small enough : " + Math.abs(linearFunctionDerivativeInAlphaj) + " <= " + (-c2 * linearFunctionDerivativeInAlpha0));
292
this.derivativeSmallEnough = true;
293
alphaOptimum = alphaj;
294
linearFunctionInAlphaOptimum = linearFunctionInAlphaj;
295
dfOptimum = dfInAlphaj;
296
//logger.debug("alphaOptimun = " + alphaOptimum);
297
//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);
300
if (linearFunctionDerivativeInAlphaj * (alphaHigh-alphaLow) >= 0) {
301
alphaHigh = alphaLow;
302
linearFunctionInAlphaHigh = linearFunctionInAlphaLow;
303
//linearFunctionDerivativeInAlphaHigh = linearFunctionDerivativeInAlphaLow;
306
linearFunctionInAlphaLow = linearFunctionInAlphaj;
307
linearFunctionDerivativeInAlphaLow = linearFunctionDerivativeInAlphaj;
308
dfInAlphaLow = dfInAlphaj;
311
//logger.debug("AlphaLow = " + alphaLow + ", AlphaHigh = " + alphaHigh);
312
//logger.debug("linearFunctionInAlphaLow = " + linearFunctionInAlphaLow + ", linearFunctionInAlphaHigh = " + linearFunctionInAlphaHigh);
313
functionEvaluationNumber = functionEvaluationNumber + 1;
314
//logger.debug("functionEvaluationNumber = " + functionEvaluationNumber);
316
if ((functionEvaluationNumber == 10) | (Math.abs(linearFunctionInAlphaHigh - linearFunctionInAlphaLow) <= 0.000001) | (Math.abs(alphaLow - alphaHigh) <= 0.000000000001)) {
317
//logger.debug("ZOOM WAS TERMINATE EARLIER");
318
/*System.out.println("functionEvaluationNumber = " + functionEvaluationNumber +
319
", Math.abs(linearFunctionInAlphaHigh - linearFunctionInAlphaLow) = " + Math.abs(linearFunctionInAlphaHigh - linearFunctionInAlphaLow) +
320
", Math.abs(alphaLow - alphaHigh) = " + Math.abs(alphaLow - alphaHigh));*/
321
this.alphaOptimum = alphaLow;
322
this.linearFunctionInAlphaOptimum = linearFunctionInAlphaLow;
323
this.dfOptimum = dfInAlphaLow;
325
//logger.debug("(functionEvaluationNumber == 10) | (Math.abs(linearFunctionInAlphaHigh - linearFunctionInAlphaLow) <= 0.000001) | (Math.abs(alphaLow - alphaHigh) <= 0.0000001)");
326
//logger.debug("zoom end -> this.alphaOptimum = " + this.alphaOptimum);
327
//logger.debug("zoom end -> this.linearFunctionInAlphaOptimum = " + this.linearFunctionInAlphaOptimum);
332
} while ((Math.abs(linearFunctionInAlphaHigh - linearFunctionInAlphaLow) > 0.000001)
333
& (functionEvaluationNumber < 10)
334
& (Math.abs(alphaLow - alphaHigh) > 0.000000000001));
336
//logger.debug("zoom end");
342
* Cubic interpolation in the interval [a,b] known to contain desirable step length
343
* and given two previous step length estimates in this interval.
345
*@param alphai Previous step length.
346
*@param linearFunctionInAlphai Function value at the previous step length alphai.
347
*@param linearFunctionDerivativeInAlphai Derivative at the previous step length alphai.
348
*@param alphaiMinus1 Previous step length.
349
*@param linearFunctionInAlphaiMinus1 Function value at the previous step length alphaiMinus1.
350
*@param linearFunctionDerivativeInAlphaiMinus1 Derivative value at the previous step length alphaiMinus1.
351
*@param a Inferior value of the interval [a,b].
352
*@param b Superior value of the interval [a,b].
354
* @return Cubic interpolation in the interval [a,b]
356
public double cubicInterpolation(double alphai, double linearFunctionInAlphai, double linearFunctionDerivativeInAlphai,
357
double alphaiMinus1, double linearFunctionInAlphaiMinus1, double linearFunctionDerivativeInAlphaiMinus1,
358
double a, double b) {
360
//logger.debug("The interval [" + a + ", " + b + "] contains acceptable step lengths.");
362
if (alphai < alphaiMinus1) {
363
this.alphaTemporal = alphai;
364
this.linearFunctionInAlphaTemporal = linearFunctionInAlphai;
365
this.linearFunctionDerivativeInAlphaTemporal = linearFunctionDerivativeInAlphai;
366
alphai = alphaiMinus1;
367
linearFunctionInAlphai = linearFunctionInAlphaiMinus1;
368
linearFunctionDerivativeInAlphai = linearFunctionDerivativeInAlphaiMinus1;
369
alphaiMinus1 = this.alphaTemporal;
370
linearFunctionInAlphaiMinus1 = this.linearFunctionInAlphaTemporal;
371
linearFunctionDerivativeInAlphaiMinus1 = this.linearFunctionDerivativeInAlphaTemporal;
374
this.d1 = linearFunctionDerivativeInAlphaiMinus1 + linearFunctionDerivativeInAlphai - 3 * ((linearFunctionInAlphaiMinus1 - linearFunctionInAlphai)/(alphaiMinus1 - alphai));
375
//logger.debug("d1 = " + d1);
377
//logger.debug("linearFunctionDerivativeInAlphaiMinus1 = " + linearFunctionDerivativeInAlphaiMinus1);
378
//logger.debug("linearFunctionDerivativeInAlphai = " + linearFunctionDerivativeInAlphai);
380
this.d2 = Math.sqrt(Math.abs(Math.pow(d1,2) - linearFunctionDerivativeInAlphaiMinus1 * linearFunctionDerivativeInAlphai));
381
//logger.debug("d2 = " + d2);
383
this.alphaiplus1 = alphai-(alphai-alphaiMinus1) * ((linearFunctionDerivativeInAlphai + d2 - d1) / (linearFunctionDerivativeInAlphai - linearFunctionDerivativeInAlphaiMinus1 + 2 * d2));
385
//logger.debug("alphaiplus1 = " + alphaiplus1);
387
if (alphaiplus1 < a) {alphaiplus1 = a;}
388
if (alphaiplus1 > b) {alphaiplus1 = b;}
390
//logger.debug("alphaiplus1 = " + alphaiplus1);
392
if (Math.abs(alphaiplus1 - alphai) < 0.000000001) {
393
/*System.out.println("We reset alphaiplus1 = (alphaiMinus1 + alphai) / 2, because alphaiplus1 = " + alphaiplus1 + " is too close to its predecessor " +
394
"alphaiMinus1 = " + alphaiMinus1); */
395
alphaiplus1 = (alphaiMinus1 + alphai) / 2;
396
} else {if (alphaiplus1 < (alphai - 9 * (alphai-alphaiMinus1) / 10)) {
397
//logger.debug("We reset alphaiplus1 = (alphaiMinus1 + alphai) / 2, because alphaiplus1 = " + alphaiplus1 + " is too much smaller than alphai = " + alphai);
398
alphaiplus1 = (alphaiMinus1 + alphai) / 2;;
407
* The aim is to find a value of alpha that satisfies the sufficient decrease condition, without being too small.
408
* The procedures generate a value alphai such that is not too much smaller than its predecesor alphai-1.
409
* The interpolation in the first is quadratic but if the sufficient decrease condition is not satisfied
410
* then the interpolation is cubic.
412
* @param alphaLow Among all step lengths generated so far and satisfying the sufficient decrease condition, the one giving the smallest function value.
413
* @param linearFunctionInAlphaLow Energy function value at alphaLow.
414
* @param linearFunctionDerivativeInAlphaLow Derivative value at alphaLow.
415
* @param alphaHigh AlphaHigh is chosen so that linearFunctionDerivativeInAlphaj * (alphaHigh-alphaLow) < 0
416
* @param linearFunctionInAlphaHigh Energy function value at alphaHigh.
417
* @return Value of alpha that satisfies the sufficient decrease condition, without being too small.
419
private double interpolation(double alphaLow, double linearFunctionInAlphaLow, double linearFunctionDerivativeInAlphaLow,
420
double alphaHigh, double linearFunctionInAlphaHigh) {
422
double minAlpha = Math.min(alphaLow, alphaHigh);
423
double alphaDiff = Math.abs(alphaHigh - alphaLow);
424
double alphaInterpolation;
426
//logger.debug("We form a quadratic approximation to the linear function");
427
double alpha1 = -1 * ((linearFunctionDerivativeInAlphaLow * Math.pow(alphaDiff,2)) / (2 * (linearFunctionInAlphaHigh - linearFunctionInAlphaLow - linearFunctionDerivativeInAlphaLow * alphaDiff)));
429
//logger.debug("The value alpha1 = " + alpha1 + ", is the minimizer of this quadratic function");
431
if ((alpha1 > alphaDiff) | (Math.abs(alpha1 - alphaDiff) < 0.000000001)) {
432
if (alpha1 < 1E-7) {}
434
/*System.out.println("We reset alpha1 = alphaDiff / 2, because alphaInterpolation = " + alpha1 + " is too close to its predecessor " +
435
"alphaiMinus1 = " + alphaDiff); */
436
alpha1 = alphaDiff / 2;
439
if ((alpha1 < 0) & (alpha1 < (alphaDiff - 9 * alphaDiff / 10))) {
440
if (alpha1 < 1E-7) {}
442
//logger.debug("We reset alphai = alphaiMinus1 / 2, because alphaInterpolation = " + alpha1 + " is too much smaller than alphaiMinus1 = " + alphaDiff);
443
alpha1 = alphaDiff / 2;
448
//logger.debug("alpha1 = " + alpha1);
450
alphaInterpolation = minAlpha + alpha1;
451
this.linearFunctionAlphaInterpolation = this.evaluateEnergyFunction(alphaInterpolation);
452
//logger.debug("alphaInterpolation = " + alphaInterpolation);
453
//logger.debug("linearFunctionAlphaInterpolation = " + this.linearFunctionAlphaInterpolation);
454
if (this.linearFunctionAlphaInterpolation <= this.linearFunctionInAlpha0 + this.c1 * (alphaInterpolation) * this.linearFunctionDerivativeInAlpha0) {
455
//logger.debug("The sufficient decrease condition is satisfied at alpha1 and we termine the interpolation");
458
//double alphaiMinus2;
459
//double alphaiMinus1 = alphaDiff;
460
//double linearFunctionInAlphaiMinus2;
461
//double linearFunctionInAlphaiMinus1 = linearFunctionInAlphaHigh;
462
double alphai; // = alpha1;
463
//double linearFunctionInAlphai = this.linearFunctionAlphaInterpolation;
466
//alphaiMinus2 = alphaiMinus1;
467
//alphaiMinus1 = alphai;
468
//linearFunctionInAlphaiMinus2 = linearFunctionInAlphaiMinus1;
469
//linearFunctionInAlphaiMinus1 = linearFunctionInAlphai;
471
//logger.debug("We construct a cubic function that interpolates the fours pieces of information");
472
a = 1/(Math.pow(alphaDiff,2) * Math.pow(alpha1, 2) * (alpha1-alphaDiff));
474
a = a * (Math.pow(alphaDiff,2) * (this.linearFunctionAlphaInterpolation - linearFunctionInAlphaLow - linearFunctionDerivativeInAlphaLow * alpha1)
475
+ (-Math.pow(alpha1,2)) * (linearFunctionInAlphaHigh - linearFunctionInAlphaLow - linearFunctionDerivativeInAlphaLow * alphaDiff));
476
b = b * (- Math.pow(alphaDiff,3) * (this.linearFunctionAlphaInterpolation - linearFunctionInAlphaLow - linearFunctionDerivativeInAlphaLow * alpha1)
477
+ Math.pow(alpha1,3) * (linearFunctionInAlphaHigh - linearFunctionInAlphaLow - linearFunctionDerivativeInAlphaLow * alphaDiff));
479
//logger.debug("a = " + a);
480
//logger.debug("b = " + b);
482
alphai = (-b + Math.sqrt(Math.pow(b,2) - 3 * a * linearFunctionDerivativeInAlphaLow)) / (3 * a);
483
//logger.debug("alphai = " + alphai);
485
if (Math.abs(alphai - alpha1) < 0.000000001) {
486
/*System.out.println("We reset alphai = alpha1 / 2, because alphaInterpolation = " + alphai + " is too close to its predecessor " +
487
"alpha1 = " + alpha1); */
489
} else {if (alphai < (alpha1 - 9 * alpha1 / 10)) {
490
//logger.debug("We reset alphai = alpha1 / 2, because alphaInterpolation = " + alphai + " is too much smaller than alpha1 = " + alpha1);
495
alphaInterpolation = minAlpha + alphai;
496
this.linearFunctionAlphaInterpolation = this.evaluateEnergyFunction(alphaInterpolation);
497
//logger.debug("alphaInterpolation = " + alphaInterpolation);
498
//logger.debug("linearFunctionAlphaInterpolation = " + this.linearFunctionAlphaInterpolation);
499
//functionEvaluationNumber = functionEvaluationNumber + 1;
501
/*} while (((linearFunctionInAlphai > linearFunctionInAlphaLow + this.c1 * (alphaLow + alphai) * linearFunctionDerivativeInAlphaLow) & (functionEvaluationNumber < 5))
502
| ((linearFunctionInAlphai - this.linearFunctionAlphaInterpolation) < 0.00000001) | ((alphai - alpha1) < 0.00000001));*/
507
return alphaInterpolation;
511
/**Evaluate the energy function from an alpha value, using the current coordinates and the current direction.
514
* @return Energy function value.
516
private double evaluateEnergyFunction(double alpha) {
517
//logger.debug("alpha= " + alpha);
518
this.xAlpha.set(this.x);
519
//logger.debug("xAlpha = " + xAlpha);
520
GVector directionStep = direction;
521
//logger.debug("directionStep = " + directionStep);
522
xAlpha.scaleAdd(alpha, directionStep, xAlpha);
523
//logger.debug("xAlpha = " + xAlpha);
524
double fxAlpha = pf.energyFunction(xAlpha);
525
//logger.debug("fxAlpha = " + fxAlpha);
530
/**Evaluate the gradient of the energy function from an alpha value,
531
* using the current coordinates and the current direction.
533
* @param alpha Alpha value for the one-dimensional problem generate from the current coordinates and the current direction.
534
* @return Gradient of the energy function at alpha.
536
private GVector evaluateEnergyFunctionDerivative(double alpha) {
537
//logger.debug("alpha= " + alpha);
538
this.xAlpha.set(this.x);
539
//logger.debug("xAlpha = " + xAlpha);
540
GVector directionStep = direction;
541
//logger.debug("directionStep = " + directionStep);
542
xAlpha.scaleAdd(alpha, directionStep, xAlpha);
543
//logger.debug("xAlpha = " + xAlpha);
544
pf.setEnergyGradient(xAlpha);
545
GVector dfxAlpha = pf.getEnergyGradient();
546
//logger.debug("dfxAlpha = " + dfxAlpha);
552
* From the interval [a, b] that bracket the minimum, evaluates the energy function at an intermediate point x
553
* and obtain a new, smaller bracketing interval, either (a,x) or (x,b).
556
* @param flambdaMin Energy function at a.
558
* @param flambdaMax Energy function at b.
559
* @return An intermediate point x
561
/*private double goldenSectionMethod(double lambdaMin, double flambdaMin, double lambdaMax, double flambdaMax) {
563
//logger.debug("Golden Section Search");
566
double lambda1 = lambdaMin;
567
double lambda4 = lambdaMax;
568
double lambda2 = lambda1 + 0.3819660112 * (lambda4 - lambda1);
569
double lambda3 = lambda1 + 0.6180339887498948482 * (lambda4 - lambda1);
571
//double flambda1 = flambdaMin;
572
double flambda2 = evaluateEnergyFunction(lambda2);
573
double flambda3 = evaluateEnergyFunction(lambda3);
574
//double flambda4 = flambdaMax;
576
//logger.debug("lambda1 = " + lambda1 + ", flambda1 = " + flambda1);
577
//logger.debug("lambda2 = " + lambda2 + ", flambda2 = " + flambda2);
578
//logger.debug("lambda3 = " + lambda3 + ", flambda3 = " + flambda3);
579
//logger.debug("lambda4 = " + lambda4 + ", flambda4 = " + flambda4);
581
if (flambda2 < flambda3) { //we can bracket the minimum by the interval [lambda1, lambda3]
582
goldenLambda = lambda2;
583
linearFunctionGoldenAlpha = flambda2;
585
else { //we can bracket the minimum by the interval [lambda2, lambda4]
586
goldenLambda = lambda3;
587
linearFunctionGoldenAlpha = flambda3;
added
removed
src/org/openscience/cdk/modeling/forcefield/LineSearchForTheWolfeConditions.java
