1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.numerics.optimization;
39
40 import ffx.numerics.OptimizationInterface;
41
42 import static ffx.numerics.optimization.LBFGS.DEFAULT_ANGLEMAX;
43 import static ffx.numerics.optimization.LBFGS.DEFAULT_CAPPA;
44 import static ffx.numerics.optimization.LBFGS.DEFAULT_INTMAX;
45 import static ffx.numerics.optimization.LBFGS.DEFAULT_SLOPEMAX;
46 import static ffx.numerics.optimization.LBFGS.DEFAULT_STEPMAX;
47 import static ffx.numerics.optimization.LBFGS.DEFAULT_STEPMIN;
48 import static ffx.numerics.optimization.LBFGS.aV1PlusV2;
49 import static ffx.numerics.optimization.LBFGS.v1DotV2;
50 import static java.lang.System.arraycopy;
51 import static org.apache.commons.math3.util.FastMath.abs;
52 import static org.apache.commons.math3.util.FastMath.acos;
53 import static org.apache.commons.math3.util.FastMath.max;
54 import static org.apache.commons.math3.util.FastMath.min;
55 import static org.apache.commons.math3.util.FastMath.sqrt;
56 import static org.apache.commons.math3.util.FastMath.toDegrees;
57
58
59
60
61
62
63
64
65 public class LineSearch {
66
67
68
69
70 private final int n;
71
72
73
74 private final double[] s;
75
76
77
78 private final double[] x0;
79
80
81
82 private OptimizationInterface optimizationSystem;
83
84
85
86 private int[] functionEvaluations;
87
88
89
90 private LineSearchResult[] info;
91
92
93
94 private double[] x;
95
96
97
98 private double[] g;
99
100
101
102 private double step;
103
104
105
106 private int interpolation;
107
108
109
110 private double f0, fA, fB, fC;
111
112
113
114 private double sg0, sgA, sgB, sgC;
115
116
117
118 private boolean restart;
119
120
121
122
123
124
125
126 LineSearch(int n) {
127 s = new double[n];
128 x0 = new double[n];
129 this.n = n;
130 }
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152 public double search(int n, double[] x, double f, double[] g, double[] p, double[] angle,
153 double fMove, LineSearchResult[] info, int[] functionEvaluations,
154 OptimizationInterface optimizationSystem) {
155
156 assert (n > 0);
157
158
159 this.x = x;
160 this.g = g;
161 this.optimizationSystem = optimizationSystem;
162 this.functionEvaluations = functionEvaluations;
163 this.info = info;
164 fA = 0.0;
165 fB = 0.0;
166 fC = 0.0;
167 sgA = 0.0;
168 sgB = 0.0;
169 sgC = 0.0;
170
171
172 info[0] = null;
173
174
175 arraycopy(p, 0, s, 0, n);
176
177
178 double gNorm = sqrt(v1DotV2(n, g, 0, 1, g, 0, 1));
179 double sNorm = sqrt(v1DotV2(n, s, 0, 1, s, 0, 1));
180
181
182
183
184
185 f0 = f;
186 arraycopy(x, 0, x0, 0, n);
187 for (int i = 0; i < n; i++) {
188 s[i] /= sNorm;
189 }
190 sg0 = v1DotV2(n, s, 0, 1, g, 0, 1);
191
192
193
194
195
196 double cosang = -sg0 / gNorm;
197 cosang = min(1.0, max(-1.0, cosang));
198 angle[0] = toDegrees(acos(cosang));
199 if (angle[0] > DEFAULT_ANGLEMAX) {
200 info[0] = LineSearchResult.WideAngle;
201 return f;
202 }
203
204
205
206
207
208 step = 2.0 * abs(fMove / sg0);
209 step = min(step, sNorm);
210 if (step > DEFAULT_STEPMAX) {
211 step = DEFAULT_STEPMAX;
212 }
213 if (step < DEFAULT_STEPMIN) {
214 step = DEFAULT_STEPMIN;
215 }
216
217 return begin();
218 }
219
220
221
222
223 private double begin() {
224 restart = true;
225 interpolation = 0;
226 fB = f0;
227 sgB = sg0;
228 return step();
229 }
230
231
232
233
234 private double step() {
235 fA = fB;
236 sgA = sgB;
237 aV1PlusV2(n, step, s, 0, 1, x, 0, 1);
238
239
240 functionEvaluations[0]++;
241 fB = optimizationSystem.energyAndGradient(x, g);
242 sgB = v1DotV2(n, s, 0, 1, g, 0, 1);
243
244
245 if (abs(sgB / sgA) >= DEFAULT_SLOPEMAX && restart) {
246 arraycopy(x0, 0, x, 0, n);
247 step /= 10.0;
248 info[0] = LineSearchResult.ScaleStep;
249 begin();
250 }
251 restart = false;
252
253
254
255
256
257 if (abs(sgB / sg0) <= DEFAULT_CAPPA && fB < fA) {
258 if (info[0] == null) {
259 info[0] = LineSearchResult.Success;
260 }
261 f0 = fB;
262 sg0 = sgB;
263 return f0;
264 }
265
266
267 if (sgB * sgA < 0.0 || fB > fA) {
268 return cubic();
269 }
270
271
272
273
274
275
276 step = 2.0 * step;
277 if (sgB > sgA) {
278 double parab = (fA - fB) / (sgB - sgA);
279 if (parab > 2.0 * step) {
280 parab = 2.0 * step;
281 }
282 if (parab < 0.5 * step) {
283 parab = 0.5 * step;
284 }
285 step = parab;
286 }
287 if (step > DEFAULT_STEPMAX) {
288 step = DEFAULT_STEPMAX;
289 }
290 return step();
291 }
292
293
294
295
296 private double cubic() {
297 interpolation++;
298 double sss = 3.0 * (fB - fA) / step - sgA - sgB;
299 double ttt = sss * sss - sgA * sgB;
300 if (ttt < 0.0) {
301 info[0] = LineSearchResult.IntplnErr;
302 f0 = fB;
303 sg0 = sgB;
304 return f0;
305 }
306 ttt = sqrt(ttt);
307 double cube = step * (sgB + ttt + sss) / (sgB - sgA + 2.0 * ttt);
308 if (cube < 0 || cube > step) {
309 info[0] = LineSearchResult.IntplnErr;
310 f0 = fB;
311 sg0 = sgB;
312 return f0;
313 }
314 aV1PlusV2(n, -cube, s, 0, 1, x, 0, 1);
315
316
317 functionEvaluations[0]++;
318 fC = optimizationSystem.energyAndGradient(x, g);
319 sgC = v1DotV2(n, s, 0, 1, g, 0, 1);
320 if (abs(sgC / sg0) <= DEFAULT_CAPPA) {
321 if (info[0] == null) {
322 info[0] = LineSearchResult.Success;
323 }
324 f0 = fC;
325 sg0 = sgC;
326 return f0;
327 }
328
329
330
331
332 if (fC <= fA || fC <= fB) {
333 double cubstp = min(abs(cube), abs(step - cube));
334 if (cubstp >= DEFAULT_STEPMIN && interpolation < DEFAULT_INTMAX) {
335 if (sgA * sgB < 0.0) {
336
337
338
339
340
341 if (sgA * sgC < 0.0) {
342 fB = fC;
343 sgB = sgC;
344 step = step - cube;
345 } else {
346 fA = fC;
347 sgA = sgC;
348 step = cube;
349 aV1PlusV2(n, cube, s, 0, 1, x, 0, 1);
350 }
351 } else {
352
353
354
355
356
357
358 if (sgA * sgC < 0.0 || fA <= fC) {
359 fB = fC;
360 sgB = sgC;
361 step -= cube;
362 } else {
363 fA = fC;
364 sgA = sgC;
365 step = cube;
366 aV1PlusV2(n, cube, s, 0, 1, x, 0, 1);
367 }
368 }
369 return cubic();
370 }
371 }
372
373
374 double f1 = min(fA, min(fB, fC));
375 double sg1;
376 if (f1 == fA) {
377 sg1 = sgA;
378 aV1PlusV2(n, cube - step, s, 0, 1, x, 0, 1);
379 } else if (f1 == fB) {
380 sg1 = sgB;
381 aV1PlusV2(n, cube, s, 0, 1, x, 0, 1);
382 } else {
383 sg1 = sgC;
384 }
385
386
387 if (f1 > f0) {
388 functionEvaluations[0]++;
389 f0 = optimizationSystem.energyAndGradient(x, g);
390 sg0 = v1DotV2(n, s, 0, 1, g, 0, 1);
391 info[0] = LineSearchResult.IntplnErr;
392 return f0;
393 }
394 f0 = f1;
395 sg0 = sg1;
396 if (sg1 > 0.0) {
397 for (int i = 0; i < n; i++) {
398 s[i] *= -1.0;
399 }
400 sg0 = -sg1;
401 }
402 step = max(cube, step - cube) / 10.0;
403 if (step < DEFAULT_STEPMIN) {
404 step = DEFAULT_STEPMIN;
405 }
406
407
408 if (info[0] == LineSearchResult.ReSearch) {
409 functionEvaluations[0]++;
410 f0 = optimizationSystem.energyAndGradient(x, g);
411 sg0 = v1DotV2(n, s, 0, 1, g, 0, 1);
412 info[0] = LineSearchResult.BadIntpln;
413 return f0;
414 } else {
415
416 info[0] = LineSearchResult.ReSearch;
417 return begin();
418 }
419 }
420
421
422
423
424
425
426 public enum LineSearchResult {
427
428
429
430 Success,
431
432
433
434 WideAngle,
435
436
437
438 ScaleStep,
439
440
441
442 IntplnErr,
443
444
445
446 ReSearch,
447
448
449
450 BadIntpln
451 }
452 }