/[escript]/trunk/downunder/py_src/inversioncostfunctions.py
ViewVC logotype

Contents of /trunk/downunder/py_src/inversioncostfunctions.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4123 - (show annotations)
Thu Dec 20 23:47:55 2012 UTC (6 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 15918 byte(s)
small fix
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 """Collection of cost functions for inversion"""
17
18 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 __all__ = [ 'InversionCostFunction']
26
27 from costfunctions import MeteredCostFunction
28 from mappings import Mapping
29 from forwardmodels import ForwardModel
30 from esys.escript.pdetools import ArithmeticTuple
31 from esys.escript import Data
32 import numpy as np
33
34 class InversionCostFunction(MeteredCostFunction):
35 """
36 Class to define cost function *J(m)* for inversion with one or more forward model
37 based on a multi-valued level set function *m*:
38
39 *J(m) = J_reg(m) + sum_f mu_f * J_f(p)*
40
41 where *J_reg(m)* is the regularization and cross gradient component of the cost function applied
42 to a level set function *m*, *J_f(p)* are the data defect cost function involving a
43 physical forward model using the physical parameter(s) *p* and mu_f is trade-off factors for model f.
44
45
46 A forward model depends on a set of physical parameters *p* which are constructed from
47 components of the level set *m* function via mappings.
48
49 Example 1 (single forward model):
50 m=Mapping()
51 f=ForwardModel()
52 J=InversionCostFunction(Regularization(), m, f)
53
54 Example 2 (two forward models on a single valued level set)
55 m0=Mapping()
56 m1=Mapping()
57 f0=ForwardModel()
58 f1=ForwardModel()
59
60 J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])
61
62 Example 2 (two forward models on 2-valued level set)
63 m0=Mapping()
64 m1=Mapping()
65 f0=ForwardModel()
66 f1=ForwardModel()
67
68 J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
69
70 :var provides_inverse_Hessian_approximation: if true the class provides an approximative inverse of the
71 Hessian operator.
72 """
73 provides_inverse_Hessian_approximation=True
74
75 def __init__(self, regularization, mappings, forward_models):
76 """
77 constructor for the cost function.
78 stores the supplied object references and sets default
79 weights.
80
81 :param regularization: the regularization part of the cost function
82 :type regularization: `Regularization`
83 :param mappings: the mappings to calculate physical parameters from the regularization.
84 This is a list of 2-tuples *(map, i)* where the first component map defines a `Mapping`
85 and the second component *i* defines the index of the component of level set function
86 to be used to calculate the mapping. An item in the list can be just a `Mapping` object
87 then the entire level set function function is fed into the `Mapping` (typically used for
88 a single-component level set function.
89 :type mappings: `list` where each item is a `tuple` of `Mapping` and `int` or just a `Mapping`.
90 :param forward_models: the forward models involved in the calculation of the cost function.
91 This is a list of 2-tuples *(f, ii)* where the first component map defines a `ForwardModel`
92 and the second component *ii* a list of indexes referring to the physical parameters
93 in the `mappings` list. The 2-tuple can be replaced by a `ForwardModel` if
94 a `mappings` list as a single entry.
95 :param forward_models: `list` where each item is 'tuple` of `ForwardModel` and `list` of `int' or is `ForwardModel`.
96 """
97 super(InversionCostFunction, self).__init__()
98 self.regularization=regularization
99
100 if isinstance(mappings, Mapping):
101 self.mappings = [mappings ]
102 else:
103 self.mappings = mappings
104
105 if isinstance(forward_models, ForwardModel):
106 self.forward_models = [ forward_models ]
107 else:
108 self.forward_models=forward_models
109
110 self.numMappings=len(self.mappings)
111 self.numModels=len(self.forward_models)
112 self.numLevelSets = self.regularization.getNumLevelSets()
113 self.__num_tradeoff_factors = self.regularization.getNumTradeOffFactors() + self.numModels
114 self.setTradeOffFactorsModels()
115
116 def getDomain(self):
117 """
118 returns the domain of the cost function
119 :rtype: 'Domain`
120 """
121 self.regularization.getDomain()
122
123 def getNumTradeOffFactors(self):
124 """
125 returns the number of trade-off factors being used including the trade-off factors used in
126 the regularization component.
127
128 :rtype: ``int``
129 """
130 return self.__num_tradeoff_factors
131 def getForwardModel(self, idx=None):
132 """
133 returns the *idx*-th forward model.
134 :param idx: model index. If cost function contains one model only `idx` can be omitted.
135 :type idx: `int`
136 """
137 if idx==None: idx=0
138 f=self.forward_models[idx]
139 if isinstance(f, ForwardModel):
140 F=f
141 else:
142 F=f[0]
143 return F
144
145 def getRegularization(self):
146 """
147 returns the regularization
148 """
149 return self.regularization
150
151
152 def setTradeOffFactorsModels(self,mu=None):
153 """
154 sets the trade-off factors for the forward model components.
155
156 :param mu: list of the trade-off factors. If not present ones are used.
157 :type mu: `float` in case of a single model or a `list` of `float` with the length of the number of models.
158 """
159 if mu==None:
160 self.mu_model=np.ones((self.numModels, ))
161 else:
162 if self.numModels > 1:
163 mu=np.asarray(mu)
164 if min(mu) > 0:
165 self.mu_model= mu
166 else:
167 raise ValueError("All value for trade-off factor mu must be positive.")
168 else:
169 mu=float(mu)
170 if mu > 0:
171 self.mu_model= [mu, ]
172 else:
173 raise ValueError("Trade-off factor must be positive.")
174
175 def setTradeOffFactorsRegularization(self,mu=None, mu_c=None):
176 """
177 sets the trade of factors for the regularization component of the cost function, see
178 `Regularization` for details.
179
180 :param m: trade-off factors for the level-set variation part
181 :param mu_c: trade-off factors for the cross gradient variation part
182 """
183 self.regularization.setTradeOffFactorsForVariation(mu)
184 self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
185
186 def setTradeOffFactors(self, mu=None):
187 """
188 sets the trade-off factors for the forward model and regularization
189 terms.
190
191 :param mu_model: Weighting factor for the forward model (default=1.)
192 :type mu_model: non-negative `float`
193 :param mu: list of trade-off factors.
194 :type mu: `list` of `float`
195 """
196 if mu==None:
197 mu=mp.ones((self.__num_tradeoff_factors,))
198 self.setTradeOffFactorsModels(mu[:self.numModels])
199 self.regularization.setTradeOffFactors(mu[self.numModels:])
200
201 def createLevelSetFunction(self, *props):
202 """
203 return an instance of an object used to represent a level set function initialed
204 with zeros. Components can be overwritten by physical properties 'props'. If present
205 entries must correspond to the a `mappings` arguments in the constructor. Use `None`
206 for properties for which no value is given.
207 """
208 m=self.regularization.getPDE().createSolution()
209 if len(props) > 0:
210 for i in xrange(self.numMappings):
211 if props[i]:
212 mm=self.mappings[i]
213 if isinstance(mm, Mapping):
214 m=mm.getInverse(props[i])
215 elif len(mm) == 1:
216 m=mm[0].getInverse(props[i])
217 else:
218 m[mm[1]]=mm[0].getInverse(props[i])
219 return m
220
221 def getProperties(self, m, return_list=False):
222 """
223 returns a list of the physical properties from a given level set function *m* using the
224 mappings of the cost function.
225
226 :param m: level set function
227 :type m: `Data`
228 :param return_list: if True a list is returned.
229 :type return_list: `bool`
230 :rtype m: `list` of `Data`
231 """
232 props=[]
233 for i in xrange(self.numMappings):
234 mm=self.mappings[i]
235 if isinstance(mm, Mapping):
236 p=mm.getValue(m)
237 elif len(mm) == 1:
238 p=mm[0].getValue(m)
239 else:
240 p=mm[0].getValue(m[mm[1]])
241 props.append(p)
242 if self.numMappings > 1 or return_list:
243 return props
244 else:
245 return props[0]
246
247 def _getDualProduct(self, x, r):
248 """
249 Returns the dual product, see `Regularization.getDualProduct`
250
251 :type x: `Data`
252 :type r: `ArithmeticTuple`
253 :rtype: `float`
254 """
255 return self.regularization.getDualProduct(x, r)
256
257 def _getArguments(self, m):
258 """
259 returns pre-computed values that are shared in the calculation of
260 *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
261 mapped value of ``m``, the arguments from the forward model and the
262 arguments from the regularization.
263
264 :param m: current approximation of the level set function
265 :type m: `Data`
266 :return: tuple of of values of the parameters, pre-computed values for the forward model and
267 pre-computed values for the regularization
268 :rtype: `tuple`
269 """
270 args_reg=self.regularization.getArguments(m)
271 # cache for physical parameters:
272 props=self.getProperties(m, return_list=True)
273 args_f=[]
274 for i in xrange(self.numModels):
275 f=self.forward_models[i]
276 if isinstance(f, ForwardModel):
277 aa=f.getArguments(props[0])
278 elif len(f) == 1:
279 aa=f[0].getArguments(props[0])
280 else:
281 idx = f[1]
282 f=f[0]
283 if isinstance(idx, int):
284 aa=f.getArguments(props[idx])
285 else:
286 pp=tuple( [ props[i] for i in idx] )
287 aa=f.getArguments(*pp)
288 args_f.append(aa)
289
290 return props, args_f, args_reg
291
292 def _getValue(self, m, *args):
293 """
294 Returns the value *J(m)* of the cost function at *m*.
295 If the pre-computed values are not supplied `getArguments()` is called.
296
297 :param m: current approximation of the level set function
298 :type m: `Data`
299 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
300 pre-computed values for the regularization
301 :rtype: `float`
302 """
303 # if there is more than one forward_model and/or regularization their
304 # contributions need to be added up. But this implementation allows
305 # only one of each...
306 if len(args)==0:
307 args=self.getArguments(m)
308
309 props=args[0]
310 args_f=args[1]
311 args_reg=args[2]
312
313 J = self.regularization.getValue(m, *args_reg)
314 print "J_reg = %e"%J
315
316 for i in xrange(self.numModels):
317
318 f=self.forward_models[i]
319 if isinstance(f, ForwardModel):
320 J_f = f.getValue(props[0],*args_f[i])
321 elif len(f) == 1:
322 J_f=f[0].getValue(props[0],*args_f[i])
323 else:
324 idx = f[1]
325 f=f[0]
326 if isinstance(idx, int):
327 J_f = f.getValue(props[idx],*args_f[i])
328 else:
329 args=tuple( [ props[j] for j in idx] + args_f[i])
330 J_f = f.getValue(*args)
331 print "J_f[%d] = %e"%(i, J_f)
332 print "mu_model[%d] = %e"%(i, self.mu_model[i])
333 J += self.mu_model[i] * J_f
334
335 return J
336
337 def _getGradient(self, m, *args):
338 """
339 returns the gradient of the cost function at *m*.
340 If the pre-computed values are not supplied `getArguments()` is called.
341
342 :param m: current approximation of the level set function
343 :type m: `Data`
344 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
345 pre-computed values for the regularization
346
347 :rtype: `ArithmeticTuple`
348 """
349 if len(args)==0:
350 args = self.getArguments(m)
351
352 props=args[0]
353 args_f=args[1]
354 args_reg=args[2]
355
356 g_J = self.regularization.getGradient(m, *args_reg)
357 p_diffs=[]
358 for i in xrange(self.numMappings):
359 mm=self.mappings[i]
360 if isinstance(mm, Mapping):
361 dpdm = mm.getDerivative(m)
362 elif len(mm) == 1:
363 dpdm = mm[0].getDerivative(m)
364 else:
365 dpdm = mm[0].getDerivative(m[mm[1]])
366 p_diffs.append(dpdm)
367
368 Y=g_J[0]
369 for i in xrange(self.numModels):
370 mu=self.mu_model[i]
371 f=self.forward_models[i]
372 if isinstance(f, ForwardModel):
373 Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
374 if self.numLevelSets == 1 :
375 Y +=Ys
376 else:
377 Y[0] +=Ys
378 elif len(f) == 1:
379 Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
380 if self.numLevelSets == 1 :
381 Y +=Ys
382 else:
383 Y[0] +=Ys
384 else:
385 idx = f[1]
386 f=f[0]
387 if isinstance(idx, int):
388 Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
389 if self.numLevelSets == 1 :
390 if idx == 0:
391 Y+=Ys
392 else:
393 raise IndexError("Illegal mapping index.")
394 else:
395 Y[idx] += Ys
396 else:
397 args=tuple( [ props[j] for j in idx] + args_f[i])
398 Ys = f.getGradient(*args)
399 for ii in xrange(len(idx)):
400 Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
401
402 return g_J
403
404
405 def _getInverseHessianApproximation(self, m, r, *args):
406 """
407 returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
408 for a given gradient type *r* at a given location *m*: *H(m) p = r*
409
410 :param m: level set approximation where to calculate Hessian inverse
411 :type m: `Data`
412 :param r: a given gradient
413 :type r: `ArithmeticTuple`
414 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
415 pre-computed values for the regularization
416 :rtype: `Data`
417 :note: in the current implementation only the regularization term is
418 considered in the inverse Hessian approximation.
419
420 """
421 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
422 return m
423
424 def updateHessian(self):
425 """
426 notifies the class that the Hessian operator needs to be updated.
427 """
428 self.regularization.updateHessian()
429
430 def _getNorm(self, m):
431 """
432 returns the norm of ``m``
433
434 :param m: level set function
435 :type m: `Data`
436 :rtype: ``float``
437 """
438 return self.regularization.getNorm(m)

  ViewVC Help
Powered by ViewVC 1.1.26