/[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 4122 - (hide annotations)
Thu Dec 20 05:42:35 2012 UTC (6 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 15664 byte(s)
More work towards joint inversion. There is now one class for inversion cost function which can handle all relevant cases 
for a single model inversion, cross gradient case and functional dependence of physical parameters.


1 gross 4074
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 gross 4122 __all__ = [ 'InversionCostFunction']
26 gross 4074
27     from costfunctions import MeteredCostFunction
28 gross 4122 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 gross 4122 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 gross 4074
49 gross 4122 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 gross 4074 """
73 caltinay 4095 provides_inverse_Hessian_approximation=True
74    
75 gross 4122 def __init__(self, regularization, mappings, forward_models):
76 gross 4074 """
77 gross 4122 constructor for the cost function.
78     stores the supplied object references and sets default
79 gross 4074 weights.
80    
81 gross 4122 :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 gross 4074 """
97 gross 4122 super(InversionCostFunction, self).__init__()
98 gross 4074 self.regularization=regularization
99 gross 4122
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 gross 4074
128 gross 4122 :rtype: ``int``
129 gross 4074 """
130 gross 4122 return self.__num_tradeoff_factors
131     def getForwardModels(self):
132     """
133     returns the forward models as a list
134     """
135     return self.forward_models
136    
137     def getRegularization(self):
138     """
139     returns the regularization
140     """
141     return self.regularization
142 caltinay 4097
143 gross 4122
144     def setTradeOffFactorsModels(self,mu=None):
145 gross 4074 """
146 gross 4122 sets the trade-off factors for the forward model components.
147 gross 4099
148 gross 4122 :param mu: list of the trade-off factors. If not present ones are used.
149     :type mu: `float` in case of a single model or a `list` of `float` with the length of the number of models.
150     """
151     if mu==None:
152     self.mu_model=np.ones((self.numModels, ))
153     else:
154     if self.numModels > 1:
155     mu=np.asarray(mu)
156     if min(mu) > 0:
157     self.mu_model= mu
158     else:
159     raise ValueError("All value for trade-off factor mu must be positive.")
160     else:
161     mu=float(mu)
162     if mu > 0:
163     self.mu_model= [mu, ]
164     else:
165     raise ValueError("Trade-off factor must be positive.")
166    
167     def setTradeOffFactorsRegularization(self,mu=None, mu_c=None):
168     """
169     sets the trade of factors for the regularization component of the cost function, see
170     `Regularization` for details.
171    
172     :param m: trade-off factors for the level-set variation part
173     :param mu_c: trade-off factors for the cross gradient variation part
174     """
175     self.regularization.setTradeOffFactorsForVariation(mu)
176     self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
177    
178     def setTradeOffFactors(self, mu=None):
179     """
180     sets the trade-off factors for the forward model and regularization
181     terms.
182    
183     :param mu_model: Weighting factor for the forward model (default=1.)
184     :type mu_model: non-negative `float`
185     :param mu: list of trade-off factors.
186     :type mu: `list` of `float`
187     """
188     if mu==None:
189     mu=mp.ones((self.__num_tradeoff_factors,))
190     self.setTradeOffFactorsModels(mu[:self.numModels])
191     self.regularization.setTradeOffFactors(mu[self.numModels:])
192    
193     def createLevelSetFunction(self, *props):
194     """
195     return an instance of an object used to represent a level set function initialed
196     with zeros. Components can be overwritten by physical properties 'props'. If present
197     entries must correspond to the a `mappings` arguments in the constructor. Use `None`
198     for properties for which no value is given.
199     """
200     m=self.regularization.getPDE().createSolution()
201     if len(props) > 0:
202     for i in xrange(self.numMappings):
203     if props[i]:
204     mm=self.mappings[i]
205     if isinstance(mm, Mapping):
206     m=mm.getInverse(props[i])
207     elif len(mm) == 1:
208     m=mm[0].getInverse(props[i])
209     else:
210     m[mm[1]]=mm[0].getInverse(props[i])
211     return m
212    
213     def getProperties(self, m, return_list=False):
214     """
215     returns a list of the physical properties from a given level set function *m* using the
216     mappings of the cost function.
217    
218     :param m: level set function
219     :type m: `Data`
220     :param return_list: if True a list is returned.
221     :type return_list: `bool`
222     :rtype m: `list` of `Data`
223     """
224     props=[]
225     for i in xrange(self.numMappings):
226     mm=self.mappings[i]
227     if isinstance(mm, Mapping):
228     p=mm.getValue(m)
229     elif len(mm) == 1:
230     p=mm[0].getValue(m)
231     else:
232     p=mm[0].getValue(m[mm[1]])
233     props.append(p)
234     if self.numMappings > 1 or return_list:
235     return props
236     else:
237     return props[0]
238    
239 gross 4074 def _getDualProduct(self, x, r):
240     """
241 gross 4122 Returns the dual product, see `Regularization.getDualProduct`
242 gross 4074
243 gross 4122 :type x: `Data`
244     :type r: `ArithmeticTuple`
245 gross 4074 :rtype: `float`
246     """
247     return self.regularization.getDualProduct(x, r)
248    
249     def _getArguments(self, m):
250     """
251 gross 4122 returns pre-computed values that are shared in the calculation of
252     *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
253 gross 4074 mapped value of ``m``, the arguments from the forward model and the
254     arguments from the regularization.
255 gross 4122
256     :param m: current approximation of the level set function
257     :type m: `Data`
258     :return: tuple of of values of the parameters, pre-computed values for the forward model and
259     pre-computed values for the regularization
260 gross 4074 :rtype: `tuple`
261     """
262 gross 4122 args_reg=self.regularization.getArguments(m)
263     # cache for physical parameters:
264     props=self.getProperties(m, return_list=True)
265     args_f=[]
266     for i in xrange(self.numModels):
267     f=self.forward_models[i]
268     if isinstance(f, ForwardModel):
269     aa=f.getArguments(props[0])
270     elif len(f) == 1:
271     aa=f[0].getArguments(props[0])
272     else:
273     idx = f[1]
274     f=f[0]
275     if isinstance(idx, int):
276     aa=f.getArguments(props[idx])
277     else:
278     pp=tuple( [ props[i] for i in idx] )
279     aa=f.getArguments(*pp)
280     args_f.append(aa)
281    
282     return props, args_f, args_reg
283 gross 4074
284     def _getValue(self, m, *args):
285     """
286 gross 4122 Returns the value *J(m)* of the cost function at *m*.
287     If the pre-computed values are not supplied `getArguments()` is called.
288 gross 4074
289 gross 4122 :param m: current approximation of the level set function
290     :type m: `Data`
291     :param args: tuple of of values of the parameters, pre-computed values for the forward model and
292     pre-computed values for the regularization
293 gross 4074 :rtype: `float`
294     """
295     # if there is more than one forward_model and/or regularization their
296     # contributions need to be added up. But this implementation allows
297     # only one of each...
298     if len(args)==0:
299     args=self.getArguments(m)
300 gross 4102
301 gross 4122 props=args[0]
302     args_f=args[1]
303     args_reg=args[2]
304    
305     J = self.regularization.getValue(m, *args_reg)
306     print "J_reg = %e"%J
307    
308     for i in xrange(self.numModels):
309    
310     f=self.forward_models[i]
311     if isinstance(f, ForwardModel):
312     J_f = f.getValue(props[0],*args_f[i])
313     elif len(f) == 1:
314     J_f=f[0].getValue(props[0],*args_f[i])
315     else:
316     idx = f[1]
317     f=f[0]
318     if isinstance(idx, int):
319     J_f = f.getValue(props[idx],*args_f[i])
320     else:
321     args=tuple( [ props[j] for j in idx] + args_f[i])
322     J_f = f.getValue(*args)
323     print "J_f[%d] = %e"%(i, J_f)
324     print "mu_model[%d] = %e"%(i, self.mu_model[i])
325     J += self.mu_model[i] * J_f
326    
327     return J
328 gross 4079
329 gross 4074 def _getGradient(self, m, *args):
330     """
331 gross 4122 returns the gradient of the cost function at *m*.
332     If the pre-computed values are not supplied `getArguments()` is called.
333 gross 4074
334 gross 4122 :param m: current approximation of the level set function
335     :type m: `Data`
336     :param args: tuple of of values of the parameters, pre-computed values for the forward model and
337     pre-computed values for the regularization
338    
339     :rtype: `ArithmeticTuple`
340 gross 4074 """
341     if len(args)==0:
342     args = self.getArguments(m)
343 gross 4122
344     props=args[0]
345     args_f=args[1]
346     args_reg=args[2]
347    
348     g_J = self.regularization.getGradient(m, *args_reg)
349     p_diffs=[]
350     for i in xrange(self.numMappings):
351     mm=self.mappings[i]
352     if isinstance(mm, Mapping):
353     dpdm = mm.getDerivative(m)
354     elif len(mm) == 1:
355     dpdm = mm[0].getDerivative(m)
356     else:
357     dpdm = mm[0].getDerivative(m[mm[1]])
358     p_diffs.append(dpdm)
359    
360     Y=g_J[0]
361     for i in xrange(self.numModels):
362     mu=self.mu_model[i]
363     f=self.forward_models[i]
364     if isinstance(f, ForwardModel):
365     Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
366     if self.numLevelSets == 1 :
367     Y +=Ys
368     else:
369     Y[0] +=Ys
370     elif len(f) == 1:
371     Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
372     if self.numLevelSets == 1 :
373     Y +=Ys
374     else:
375     Y[0] +=Ys
376     else:
377     idx = f[1]
378     f=f[0]
379     if isinstance(idx, int):
380     Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
381     if self.numLevelSets == 1 :
382     if idx == 0:
383     Y+=Ys
384     else:
385     raise IndexError("Illegal mapping index.")
386     else:
387     Y[idx] += Ys
388     else:
389     args=tuple( [ props[j] for j in idx] + args_f[i])
390     Ys = f.getGradient(*args)
391     for ii in xrange(len(idx)):
392     Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
393 gross 4074
394 gross 4122 return g_J
395 gross 4074
396    
397     def _getInverseHessianApproximation(self, m, r, *args):
398     """
399     returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
400     for a given gradient type *r* at a given location *m*: *H(m) p = r*
401 caltinay 4097
402 gross 4074 :param m: level set approximation where to calculate Hessian inverse
403 caltinay 4097 :type m: `Data`
404     :param r: a given gradient
405     :type r: `ArithmeticTuple`
406 gross 4122 :param args: tuple of of values of the parameters, pre-computed values for the forward model and
407     pre-computed values for the regularization
408 caltinay 4097 :rtype: `Data`
409     :note: in the current implementation only the regularization term is
410     considered in the inverse Hessian approximation.
411    
412 gross 4074 """
413 gross 4079 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
414     return m
415 caltinay 4097
416 gross 4074 def updateHessian(self):
417     """
418     notifies the class that the Hessian operator needs to be updated.
419     """
420 caltinay 4095 self.regularization.updateHessian()
421 gross 4100
422     def _getNorm(self, m):
423     """
424     returns the norm of ``m``
425 caltinay 4095
426 gross 4100 :param m: level set function
427     :type m: `Data`
428     :rtype: ``float``
429     """
430 gross 4122 return self.regularization.getNorm(m)

  ViewVC Help
Powered by ViewVC 1.1.26