|
samer@0
|
1 /*
|
|
samer@0
|
2 * Copyright (c) 2000, Samer Abdallah, King's College London.
|
|
samer@0
|
3 * All rights reserved.
|
|
samer@0
|
4 *
|
|
samer@0
|
5 * This software is provided AS iS and WITHOUT ANY WARRANTY;
|
|
samer@0
|
6 * without even the implied warranty of MERCHANTABILITY or
|
|
samer@0
|
7 * FITNESS FOR A PARTICULAR PURPOSE.
|
|
samer@0
|
8 */
|
|
samer@0
|
9
|
|
samer@0
|
10 package samer.maths.opt;
|
|
samer@0
|
11 import samer.maths.*;
|
|
samer@0
|
12 import samer.core.*;
|
|
samer@0
|
13
|
|
samer@0
|
14 /**
|
|
samer@0
|
15 Line search from <b>Numerical Recipes in C</b>, pretty
|
|
samer@0
|
16 much as supplied.
|
|
samer@0
|
17
|
|
samer@0
|
18 I removed the max step and slope checking and made it the
|
|
samer@0
|
19 responsibility of the calling code to check these.
|
|
samer@0
|
20 Also, this works with ONE-DIMENSIONAL functions: the stuff
|
|
samer@0
|
21 about vectors can be handled by a <code>LineFunction</code> object.
|
|
samer@0
|
22 */
|
|
samer@0
|
23
|
|
samer@0
|
24 public class PolynomialLineSearch
|
|
samer@0
|
25 {
|
|
samer@0
|
26 State S;
|
|
samer@0
|
27
|
|
samer@0
|
28 public PolynomialLineSearch(State s) { S=s; }
|
|
samer@0
|
29
|
|
samer@0
|
30 /**
|
|
samer@0
|
31 This assumes that the initial step has already been taken
|
|
samer@0
|
32 @param alamin minimum step size
|
|
samer@0
|
33 */
|
|
samer@0
|
34
|
|
samer@0
|
35 public void run(Condition termination)
|
|
samer@0
|
36 {
|
|
samer@0
|
37 double alam, alam2, f3;
|
|
samer@0
|
38 double tmplam;
|
|
samer@0
|
39
|
|
samer@0
|
40 // if initial step is ok...
|
|
samer@0
|
41 if (termination.test()) return;
|
|
samer@0
|
42
|
|
samer@0
|
43 // backtrack using quadratic interpolation
|
|
samer@0
|
44 tmplam = -S.P1.s/(2*(S.P1.f-S.P2.f-S.P1.s));
|
|
samer@0
|
45 alam = S.alpha;
|
|
samer@0
|
46
|
|
samer@0
|
47 for (;;) {
|
|
samer@0
|
48
|
|
samer@0
|
49 // remember failed parameters
|
|
samer@0
|
50 alam2=alam; f3=S.P2.f;
|
|
samer@0
|
51 alam=Math.max(tmplam,0.1*alam);
|
|
samer@0
|
52
|
|
samer@0
|
53 S.step(alam);
|
|
samer@0
|
54
|
|
samer@0
|
55 // if step was ok...
|
|
samer@0
|
56 if (termination.test()) return;
|
|
samer@0
|
57
|
|
samer@0
|
58 // backtrack using cubic interpolation
|
|
samer@0
|
59 double rhs1 = S.P2.f - S.P1.f - alam*S.P1.s;
|
|
samer@0
|
60 double rhs2 = f3 - S.P1.f - alam2*S.P1.s;
|
|
samer@0
|
61 double a = (rhs1/(alam*alam) - rhs2/(alam2*alam2)) / (alam-alam2);
|
|
samer@0
|
62 double b = (-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
|
|
samer@0
|
63
|
|
samer@0
|
64 if (a==0.0) tmplam = -S.P1.s/(2*b);
|
|
samer@0
|
65 else {
|
|
samer@0
|
66 double disc=b*b-2*a*S.P1.s;
|
|
samer@0
|
67 if (disc<0) tmplam=alam/2;
|
|
samer@0
|
68 else if (b<=0) tmplam=(-b+Math.sqrt(disc))/(3*a);
|
|
samer@0
|
69 else tmplam=-S.P1.s/(b+Math.sqrt(disc));
|
|
samer@0
|
70 }
|
|
samer@0
|
71 if (tmplam>alam/2) tmplam=alam/2;
|
|
samer@0
|
72 }
|
|
samer@0
|
73 }
|
|
samer@0
|
74 }
|
