/[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 4393 - (show annotations)
Mon May 6 03:35:48 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 18976 byte(s)
Correct indexing of forward models.

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

  ViewVC Help
Powered by ViewVC 1.1.26