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

  ViewVC Help
Powered by ViewVC 1.1.26