/[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 4263 - (show annotations)
Thu Feb 28 05:59:04 2013 UTC (6 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 17435 byte(s)
No more xrange in downunder.

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

  ViewVC Help
Powered by ViewVC 1.1.26