Mercurial > hg > jslab
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 } |