comparison src/samer/maths/opt/PolynomialLineSearch.java @ 0:bf79fb79ee13

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