tlinesearchCR.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tlinesearchCR.py (7403B)
---
1 ############################################################################
2 #
3 # This file is a part of siple.
4 #
5 # Copyright 2010, 2014 David Maxwell
6 #
7 # siple is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 2 of the License, or
10 # (at your option) any later version.
11 #
12 ############################################################################
13
14 from numpy import isinf, isnan, sqrt, any, isreal, real, nan, inf
15 from siple.reporting import msg
16 from siple.params import Bunch, Parameters
17
18 #This program is distributed WITHOUT ANY WARRANTY; without even the implied
19 #warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
20 #
21 # This is a modified version of:
22 #
23 #
24 #This file contains a Python version of Carl Rasmussen's Matlab-function
25 #minimize.m
26 #
27 #minimize.m is copyright (C) 1999 - 2006, Carl Edward Rasmussen.
28 #Python adaptation by Roland Memisevic 2008.
29 #
30 #
31 #The following is the original copyright notice that comes with the
32 #function minimize.m
33 #(from http://www.kyb.tuebingen.mpg.de/bs/people/carl/code/minimize/Copyright):
34 #
35 #
36 #"(C) Copyright 1999 - 2006, Carl Edward Rasmussen
37 #
38 #Permission is granted for anyone to copy, use, or modify these
39 #programs and accompanying documents for purposes of research or
40 #education, provided this copyright notice is retained, and note is
41 #made of any changes that have been made.
42 #
43 #These programs and documents are distributed without any warranty,
44 #express or implied. As the programs were written for research
45 #purposes only, they have not been tested to the degree that would be
46 #advisable in any important application. All use of these programs is
47 #entirely at the user's own risk."
48
49 class LinesearchCR:
50
51 @staticmethod
52 def defaultParameters():
53 return Parameters('linesearchCR',
54 INT= 0.1, # don't reevaluate within 0.1 of the limit of the current bracket
55 EXT = 3.0, # extrapolate maximum 3 times the current step-size
56 MAX = 20, # max 20 function evaluations per line search
57 RATIO = 10, # maximum allowed slope ratio
58 SIG = 0.1,
59 verbose=False) # RHO = SIG/2) # SIG and RHO are the constants controlling the Wolfe-
60
61 def __init__(self,params=None):
62 self.params = self.defaultParameters()
63 if not (params is None): self.params.update(params)
64
65 def error(self):
66 return self.code > 0
67
68 def search(self,f,f0,df0,t0=None):
69 INT = self.params.INT; # don't reevaluate within 0.1 of the limit of the current bracket
70 EXT = self.params.EXT; # extrapolate maximum 3 times the current step-size
71 MAX = self.params.MAX; # max 20 function evaluations per line search
72 RATIO = self.params.RATIO; # maximum allowed slope ratio
73 SIG = self.params.SIG;
74 RHO = SIG/2;
75 SMALL = 10.**-16 #minimize.m uses matlab's realmin
76
77 # SIG and RHO are the constants controlling the Wolfe-
78 #Powell conditions. SIG is the maximum allowed absolute ratio between
79 #previous and new slopes (derivatives in the search direction), thus setting
80 #SIG to low (positive) values forces higher precision in the line-searches.
81 #RHO is the minimum allowed fraction of the expected (from the slope at the
82 #initial point in the linesearch). Constants must satisfy 0 < RHO < SIG < 1.
83 #Tuning of SIG (depending on the nature of the function to be optimized) may
84 #speed up the minimization; it is probably not worth playing much with RHO.
85
86 d0 = df0;
87 fdata0 = None
88 if t0 is None:
89 t0 = 1/(1.-d0)
90 x3 = t0 # initial step is red/(|s|+1)
91 X0 = 0; F0 = f0; dF0 = df0; Fdata0=fdata0 # make a copy of current values
92 M = MAX
93
94 x2 = 0; f2 = f0; d2 = d0;
95 while True: # keep extrapolating as long as necessary
96 # x2 = 0; f2 = f0; d2 = d0;
97 f3 = f0; df3 = df0; fdata3 = fdata0;
98 success = 0
99 while (not success) and (M > 0):
100 try:
101 M = M - 1
102 (f3, df3, fdata3) = f(x3)
103 if isnan(f3) or isinf(f3) or any(isnan(df3)+isinf(df3)):
104 raise Exception('nan')
105 success = 1
106 except: # catch any error which occured in f
107 if self.params.verbose: msg('error on extrapolate. shrinking %g to %g', x3, (x2+x3)/2)
108 x3 = (x2+x3)/2 # bisect and try again
109
110 if f3 < F0:
111 X0 = x3; F0 = f3; dF0 = df3; Fdata0=fdata3 # keep best values
112 d3 = df3 # new slope
113 if d3 > SIG*d0 or f3 > f0+x3*RHO*d0 or M == 0: break # are we done extrapolati
114
115 x1 = x2; f1 = f2; d1 = d2 # move point 2 to point 1
116 x2 = x3; f2 = f3; d2 = d3 # move point 3 to point 2
117 A = 6*(f1-f2)+3*(d2+d1)*(x2-x1) # make cubic extrapolation
118 B = 3*(f2-f1)-(2*d1+d2)*(x2-x1)
119 Z = B+sqrt(complex(B*B-A*d1*(x2-x1)))
120 if Z != 0.0:
121 x3 = x1-d1*(x2-x1)**2/Z # num. error possible, ok!
122 else:
123 x3 = inf
124 if (not isreal(x3)) or isnan(x3) or isinf(x3) or (x3 < 0):
125 # num prob | wrong sign?
126 x3 = x2*EXT # extrapolate maximum amount
127 elif x3 > x2*EXT: # new point beyond extrapolation limit?
128 x3 = x2*EXT # extrapolate maximum amount
129 elif x3 < x2+INT*(x2-x1): # new point too close to previous point?
130 x3 = x2+INT*(x2-x1)
131 x3 = real(x3)
132 msg('extrapolating: x1 %g d1 %g x2 %g d2 %g x3 %g',x1,d1,x2,d3,x3)
133
134
135 while (abs(d3) > -SIG*d0 or f3 > f0+x3*RHO*d0) and M > 0:
136 if (d3 > 0) or (f3 > f0+x3*RHO*d0): # choose subinterval
137 x4 = x3; f4 = f3; d4 = d3 # move point 3 to point 4
138 else:
139 x2 = x3; f2 = f3; d2 = d3 # move point 3 to point 2
140 if self.params.verbose: msg('interpolating x2 %g x4 %g f2 %g f4 %g wolfef %g d2 %g d4 %g wolfed %g ',x2,x4,f2,f4,f0+x3*RHO*d0,d2,d4,-SIG*d0)
141
142 if f4 > f0:
143 x3 = x2-(0.5*d2*(x4-x2)**2)/(f4-f2-d2*(x4-x2)) # quadratic interpolation
144 else:
145 A = 6*(f2-f4)/(x4-x2)+3*(d4+d2) # cubic interpolation
146 B = 3*(f4-f2)-(2*d2+d4)*(x4-x2)
147 if A != 0:
148 x3=x2+(sqrt(B*B-A*d2*(x4-x2)**2)-B)/A # num. error possible, ok!
149 else:
150 x3 = inf
151 if isnan(x3) or isinf(x3):
152 x3 = (x2+x4)/2 # if we had a numerical problem then bisect
153 x3 = max(min(x3, x4-INT*(x4-x2)),x2+INT*(x4-x2)) # don't accept too close
154 (f3, df3, fdata3) = f(x3);
155 d3 =df3;
156 M = M - 1; # count epochs
157 if f3 < F0:
158 X0 = x3; F0 = f3; dF0 = df3; Fdata0 = fdata3 # keep best values
159
160 if (abs(d3) < -SIG*d0) and (f3 < f0+x3*RHO*d0): # if line search succeeded
161 self.code = 0
162 self.value = Bunch(F=f3,Fp=d3,t=x3,data=fdata3)
163 self.errMsg = ""
164 else:
165 self.code = 1
166 if M == 0:
167 self.errMsg = 'Too many function evaluations (>%d)' % MAX
168 else:
169 self.errMsg = 'unknown error';
170 self.value = Bunch(F=f0,Fp=dF0,t=X0,data=Fdata0)