tforward.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
---
tforward.py (7313B)
---
1 ############################################################################
2 #
3 # This file is a part of siple.
4 #
5 # Copyright 2010 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 import numpy
15 from siple.reporting import msg
16
17 class LinearForwardProblem:
18 """
19 We are interested in solving an ill-posed problem
20
21 .. math:: T(x) = y
22
23 where :math:`T:X\\rta Y` is a linear map between Hilbert spaces. One class of methods
24 for solving such a problem involve iteratively methods that approximately minimize
25
26 .. math:: J(x) = || y - T(x) ||^2_Y
27
28 These techniques require the following maps:
29
30 The forward problem: T
31
32 the adjoint of T: T*
33
34 This class encapsulates the solution of all these maps for a single forward problem,
35 as well as defining the norms on the domain and range spaces.
36 """
37
38 def T(self,x,out=None):
39 """
40 Implements the forward linear problem. Returns T(x) in the variable :data:`out`.
41 """
42 raise NotImplementedError()
43
44 def TStar(self,y,out=None):
45 """
46 Implements the adjoint of the forward linear problem. Returns T^*(y) in the variable :data:`out`.
47 """
48 raise NotImplementedError()
49
50 def domainIP(self,u,v):
51 """
52 Implements inner product in the domain X.
53 """
54 raise NotImplementedError()
55
56 def rangeIP(self,u,v):
57 """
58 Implements inner product in the range Y.
59 """
60 raise NotImplementedError()
61
62 def testTStar(self, d, r):
63 """
64 Helper method for determining the adjoint has been coded correctly.
65
66 Returns the following inner products:
67
68 <T(d),r>_Y
69
70 <d,T*(r)>_X
71
72 If the linearization is correct and its adjoint is correct, then these should be the same number.
73 """
74 Td = self.T(d)
75 TStarR = self.TStar(r)
76
77 return (self.domainIP(d,TStarR), self.rangeIP(Td,r))
78
79
80 class NonlinearForwardProblem(LinearForwardProblem):
81 """
82 We are interested in solving an ill-posed problem
83
84 F(x) = y
85
86 where F:X->Y is a nonlinear map between Hilbert spaces. One class of methods
87 for solving such a problem involve iteratively methods that approximately minimize
88
89 J(x) = 1/2 || y - F(x) ||^2_Y
90
91 These techniques require the following maps:
92
93 the forward problem: F
94 its linearization: T
95 the adjoint of T: T*
96
97 This class encapsulates the solution of all these maps for a single forward problem,
98 as well as defining the norms on the domain and range spaces.
99 """
100
101 def F(self, x,out=None,guess=None):
102 """
103 Returns the value of the forward problem at x.
104
105 Nonlinear problems often make use of an initial guess; this can be provided in 'guess'.
106
107 Storage in 'out', if given, is used for the return value.
108 """
109 raise NotImplementedError()
110
111 def T(self,d,out=None):
112 """
113 Returns the value of the linearization, T, of F, at the point x specified previously in linearizeAt,
114 in the direction d.
115
116 Storage in 'out', if given, is used for the return value.
117 """
118 raise NotImplementedError()
119
120 def TStar(self,r,out=None):
121 """
122 Let T be the linearization of F at the point x (at the point x specified previously in linearizeAt).
123 Its adjoint is T*. This method returns the value of T* in the direction r.
124
125 Storage in 'out', if given, is used for the return value.
126 """
127 raise NotImplementedError()
128
129 def linearizeAt(self,x,guess=None):
130 """
131 Instructs the class that subsequent calls to T and TStar will be conducted for the given value of x.
132
133 Nonlinear problems often make use of an initial guess; this can be provided in 'guess'.
134 """
135 raise NotImplementedError()
136
137 def evalFandLinearize(self,x,out=None,guess=None):
138 """
139 Computes the value of F(x) and locks in a linearization. Sometimes there are efficiencies that
140 can be acheived this way.
141
142 Default implementation simply calls F, then linearizeAt.
143 """
144 Fx = self.F(x,out=out,guess=guess)
145 self.linearizeAt(x,guess=Fx)
146 return Fx
147
148 def rangeIP(self,a,b):
149 """
150 Computes the inner product of two vectors in the range.
151 """
152 raise NotImplementedError()
153
154 def domainIP(self,a,b):
155 """
156 Computes the inner product of two vectors in the domain.
157 """
158 raise NotImplementedError()
159
160 def testT(self, x, h, t=1e-4,guess=None):
161 """
162 Helper method for determining the the forward linearziation has been coded correctly.
163
164 Returns the difference quotient (F(x+th)-F(x))/t and the value T(h). If the linearization
165 has been coded correctly, these vectors should be close to each other.
166 """
167 xp = x.copy()
168 xp.axpy(t,h)
169
170 y=self.F(x,guess=guess)
171 yp = self.F(xp,guess=y)
172
173 Fp = yp.copy()
174 Fp.axpy(-1,y)
175 Fp /= t
176
177 self.linearizeAt(x,guess=y)
178 Fp2 = self.T(h)
179
180 return (Fp,Fp2)
181
182 def testTStar(self, x, d, r, guess=None):
183 """
184 Helper method for determining the adjoint has been coded correctly.
185
186 Returns the following inner products:
187
188 <T(d),r>_Y
189 <d,T*(r)>_X
190
191 If the linearization is correct and its adjoint is correct, then these should be the same number.
192 """
193 self.linearizeAt(x,guess=guess)
194
195 Td = self.T(d)
196 TStarR = self.TStar(r)
197
198 return (self.domainIP(d,TStarR), self.rangeIP(Td,r))
199
200
201
202 class ForwardProblemLineSearchAdaptor:
203 """
204 Let F:X->Y be the nonlinear map between Hilbert spaces defined in a ForwardProblem.
205
206 We are interested in approximately minimizing a functional of the form
207
208 J(x) = 0.5 * || y-F(x) ||_Y^2
209
210 Doing this frequently involves a line-search in the direction 'd', i.e. the function
211
212 Phi(t) = J(x+t*d)
213
214 is minimized over t>=0.
215
216 This class encapsulates the conversion of a ForwardProblem into the map Phi defined above.
217 """
218
219 def __init__(self,forward_problem):
220 """
221 Initialization from a ForwardProblem
222 """
223 self.forward_problem = forward_problem
224 self.z = None
225 self.r = None
226 self.Fz = None
227 self.Td = None
228
229 def eval(self,x,d,y,t):
230 """
231 Efficiently evaluate the function Phi(t) and Phi'(t) where
232
233 Phi(t) = 0.5 * || y- F(x+td) ||_Y^2
234
235 as described in the classlevel documentation.
236 """
237
238 if self.z is None:
239 self.z = x.copy()
240 self.r = y.copy()
241 # self.Fz = dolfinutils.vectorLike(y)
242 self.Td = y.vector_like()
243 else:
244 self.z.set(x)
245 self.r.set(y)
246
247 forward_problem = self.forward_problem
248
249 self.z.axpy(t,d)
250
251 try:
252 self.Fz = forward_problem.evalFandLinearize(self.z,out=self.Fz,guess=self.Fz)
253 except Exception as e:
254 msg('Exception during evaluation of linesearch at t=%g;\n%s',t,str(e))
255 return (numpy.nan,numpy.nan,None)
256 self.Td=forward_problem.T(d,out=self.Td)
257
258 self.r -= self.Fz
259 J = 0.5*forward_problem.rangeIP(self.r,self.r)
260 Jp = -self.forward_problem.rangeIP(self.Td,self.r)
261
262 # Extra return value 'None' can be used to store extra data.
263 return (J,Jp,None)