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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4393 - (hide 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 gross 4074
2     ##############################################################################
3     #
4 jfenwick 4154 # Copyright (c) 2003-2013 by University of Queensland
5 gross 4074 # 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 caltinay 4178 """Cost functions for inversions with one or more forward models"""
17 gross 4074
18 jfenwick 4154 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 gross 4074 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 gross 4122 __all__ = [ 'InversionCostFunction']
26 gross 4074
27 jfenwick 4238 from .costfunctions import MeteredCostFunction
28     from .mappings import Mapping
29     from .forwardmodels import ForwardModel
30 gross 4074 from esys.escript.pdetools import ArithmeticTuple
31     from esys.escript import Data
32 gross 4122 import numpy as np
33 gross 4074
34 jfenwick 4238
35 gross 4122 class InversionCostFunction(MeteredCostFunction):
36     """
37 caltinay 4178 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 gross 4122 *J(m) = J_reg(m) + sum_f mu_f * J_f(p)*
41 gross 4074
42 caltinay 4178 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 gross 4122 Example 1 (single forward model):
51     m=Mapping()
52     f=ForwardModel()
53     J=InversionCostFunction(Regularization(), m, f)
54 caltinay 4178
55 gross 4122 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 caltinay 4178
63 caltinay 4393 Example 3 (two forward models on 2-valued level set)
64 gross 4122 m0=Mapping()
65     m1=Mapping()
66     f0=ForwardModel()
67     f1=ForwardModel()
68    
69 caltinay 4178 J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
70    
71 caltinay 4283 :cvar provides_inverse_Hessian_approximation: if true the class provides an
72     approximative inverse of the Hessian operator.
73 gross 4074 """
74 caltinay 4095 provides_inverse_Hessian_approximation=True
75    
76 gross 4122 def __init__(self, regularization, mappings, forward_models):
77 gross 4074 """
78 caltinay 4178 constructor for the cost function.
79 caltinay 4132 Stores the supplied object references and sets default weights.
80 gross 4074
81 gross 4122 :param regularization: the regularization part of the cost function
82     :type regularization: `Regularization`
83 caltinay 4132 :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 caltinay 4393 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 caltinay 4132 :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 caltinay 4393 can be replaced by a `ForwardModel` if the
100     `mappings` list has a single entry.
101     :param forward_models: `ForwardModel` or ``list``
102 gross 4074 """
103 gross 4122 super(InversionCostFunction, self).__init__()
104 gross 4074 self.regularization=regularization
105 caltinay 4178
106 gross 4122 if isinstance(mappings, Mapping):
107 caltinay 4132 self.mappings = [mappings ]
108     else:
109     self.mappings = mappings
110 caltinay 4178
111 caltinay 4393 if isinstance(forward_models, ForwardModel):
112 caltinay 4132 self.forward_models = [ forward_models ]
113 caltinay 4178 else:
114 gross 4122 self.forward_models=forward_models
115 gross 4373
116 caltinay 4393 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 caltinay 4178
125 gross 4122 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 caltinay 4178
131 gross 4122 def getDomain(self):
132     """
133     returns the domain of the cost function
134     :rtype: 'Domain`
135     """
136     self.regularization.getDomain()
137 caltinay 4178
138 gross 4122 def getNumTradeOffFactors(self):
139     """
140 caltinay 4132 returns the number of trade-off factors being used including the
141     trade-off factors used in the regularization component.
142 gross 4074
143 gross 4122 :rtype: ``int``
144 gross 4074 """
145 gross 4122 return self.__num_tradeoff_factors
146 caltinay 4132
147 gross 4123 def getForwardModel(self, idx=None):
148 gross 4122 """
149 gross 4123 returns the *idx*-th forward model.
150 caltinay 4132
151     :param idx: model index. If cost function contains one model only `idx`
152     can be omitted.
153     :type idx: ``int``
154 gross 4122 """
155 gross 4123 if idx==None: idx=0
156     f=self.forward_models[idx]
157 caltinay 4178 if isinstance(f, ForwardModel):
158     F=f
159 gross 4123 else:
160 caltinay 4178 F=f[0]
161 gross 4123 return F
162 caltinay 4178
163 gross 4122 def getRegularization(self):
164     """
165     returns the regularization
166 caltinay 4178
167 gross 4142 :rtype: `Regularization`
168 gross 4122 """
169     return self.regularization
170 caltinay 4097
171 caltinay 4178 def setTradeOffFactorsModels(self, mu=None):
172 gross 4074 """
173 gross 4122 sets the trade-off factors for the forward model components.
174 caltinay 4178
175 gross 4122 :param mu: list of the trade-off factors. If not present ones are used.
176 caltinay 4178 :type mu: ``float`` in case of a single model or a ``list`` of
177     ``float`` with the length of the number of models.
178 gross 4122 """
179     if mu==None:
180 caltinay 4178 self.mu_model=np.ones((self.numModels, ))
181 gross 4122 else:
182 caltinay 4132 if self.numModels > 1:
183 caltinay 4178 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 caltinay 4132 else:
189 caltinay 4178 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 gross 4142 def getTradeOffFactorsModels(self):
196     """
197     returns the trade-off factors for the forward models
198 caltinay 4178
199 gross 4142 :rtype: ``float`` or ``list`` of ``float``
200     """
201     if self.numModels>1:
202 caltinay 4178 return self.mu_model
203     else:
204     return self.mu_model[0]
205    
206     def setTradeOffFactorsRegularization(self, mu=None, mu_c=None):
207 gross 4122 """
208 caltinay 4178 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 gross 4122 """
214     self.regularization.setTradeOffFactorsForVariation(mu)
215     self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
216 caltinay 4178
217 gross 4122 def setTradeOffFactors(self, mu=None):
218 caltinay 4132 """
219     sets the trade-off factors for the forward model and regularization
220 caltinay 4178 terms.
221 gross 4122
222 caltinay 4178 :param mu: list of trade-off factors.
223 caltinay 4132 :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 caltinay 4178
230 gross 4142 def getTradeOffFactors(self, mu=None):
231     """
232 caltinay 4178 returns a list of the trade-off factors.
233    
234 gross 4142 :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 gross 4122
240     def createLevelSetFunction(self, *props):
241     """
242 caltinay 4132 returns an instance of an object used to represent a level set function
243     initialized with zeros. Components can be overwritten by physical
244 caltinay 4178 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 gross 4122 """
248     m=self.regularization.getPDE().createSolution()
249     if len(props) > 0:
250 caltinay 4263 for i in range(self.numMappings):
251 caltinay 4178 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 gross 4122 return m
260 caltinay 4178
261 gross 4122 def getProperties(self, m, return_list=False):
262     """
263 caltinay 4132 returns a list of the physical properties from a given level set
264     function *m* using the mappings of the cost function.
265 caltinay 4178
266 gross 4122 :param m: level set function
267     :type m: `Data`
268 caltinay 4178 :param return_list: if ``True`` a list is returned.
269     :type return_list: ``bool``
270     :rtype: ``list`` of `Data`
271 gross 4122 """
272     props=[]
273 caltinay 4263 for i in range(self.numMappings):
274 caltinay 4178 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 gross 4122 if self.numMappings > 1 or return_list:
283 caltinay 4178 return props
284 caltinay 4132 else:
285 caltinay 4178 return props[0]
286    
287 gross 4074 def _getDualProduct(self, x, r):
288     """
289 gross 4122 Returns the dual product, see `Regularization.getDualProduct`
290 gross 4074
291 gross 4122 :type x: `Data`
292 caltinay 4178 :type r: `ArithmeticTuple`
293     :rtype: ``float``
294 gross 4074 """
295     return self.regularization.getDualProduct(x, r)
296    
297     def _getArguments(self, m):
298     """
299 gross 4122 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 gross 4074 mapped value of ``m``, the arguments from the forward model and the
302     arguments from the regularization.
303 caltinay 4178
304 gross 4122 :param m: current approximation of the level set function
305     :type m: `Data`
306 caltinay 4178 :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 gross 4074 """
311 gross 4122 args_reg=self.regularization.getArguments(m)
312     # cache for physical parameters:
313     props=self.getProperties(m, return_list=True)
314     args_f=[]
315 caltinay 4263 for i in range(self.numModels):
316 caltinay 4178 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 gross 4122 return props, args_f, args_reg
332 gross 4074
333     def _getValue(self, m, *args):
334     """
335 gross 4122 Returns the value *J(m)* of the cost function at *m*.
336     If the pre-computed values are not supplied `getArguments()` is called.
337 gross 4074
338 gross 4122 :param m: current approximation of the level set function
339     :type m: `Data`
340 caltinay 4178 :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 gross 4074 """
345     if len(args)==0:
346     args=self.getArguments(m)
347 caltinay 4178
348 gross 4122 props=args[0]
349     args_f=args[1]
350     args_reg=args[2]
351 caltinay 4178
352 gross 4122 J = self.regularization.getValue(m, *args_reg)
353 caltinay 4178
354 caltinay 4263 for i in range(self.numModels):
355 caltinay 4178 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 gross 4079
372 caltinay 4178 return J
373    
374 plaub 4367 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 gross 4074 def _getGradient(self, m, *args):
418     """
419 caltinay 4178 returns the gradient of the cost function at *m*.
420 gross 4122 If the pre-computed values are not supplied `getArguments()` is called.
421 gross 4074
422 gross 4122 :param m: current approximation of the level set function
423     :type m: `Data`
424 caltinay 4178 :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 gross 4122 :rtype: `ArithmeticTuple`
429 gross 4074 """
430     if len(args)==0:
431     args = self.getArguments(m)
432 caltinay 4178
433 gross 4122 props=args[0]
434     args_f=args[1]
435     args_reg=args[2]
436 caltinay 4178
437     g_J = self.regularization.getGradient(m, *args_reg)
438 gross 4122 p_diffs=[]
439 caltinay 4263 for i in range(self.numMappings):
440 caltinay 4178 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 caltinay 4263 for i in range(self.numModels):
451 caltinay 4178 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 caltinay 4263 for ii in range(len(idx)):
481 caltinay 4178 Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
482 gross 4074
483 gross 4122 return g_J
484 gross 4074
485     def _getInverseHessianApproximation(self, m, r, *args):
486     """
487 caltinay 4178 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 caltinay 4097
491 gross 4074 :param m: level set approximation where to calculate Hessian inverse
492 caltinay 4097 :type m: `Data`
493     :param r: a given gradient
494     :type r: `ArithmeticTuple`
495 caltinay 4178 :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 caltinay 4097 :rtype: `Data`
499     :note: in the current implementation only the regularization term is
500     considered in the inverse Hessian approximation.
501 gross 4074 """
502 gross 4079 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
503     return m
504 caltinay 4097
505 gross 4074 def updateHessian(self):
506     """
507     notifies the class that the Hessian operator needs to be updated.
508     """
509 caltinay 4095 self.regularization.updateHessian()
510 caltinay 4178
511 gross 4100 def _getNorm(self, m):
512     """
513 caltinay 4178 returns the norm of `m`
514 caltinay 4095
515 gross 4100 :param m: level set function
516     :type m: `Data`
517     :rtype: ``float``
518     """
519 caltinay 4132 return self.regularization.getNorm(m)
520 caltinay 4178

  ViewVC Help
Powered by ViewVC 1.1.26