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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4121 by gross, Wed Dec 12 06:17:03 2012 UTC revision 4122 by gross, Thu Dec 20 05:42:35 2012 UTC
# Line 22  __license__="""Licensed under the Open S Line 22  __license__="""Licensed under the Open S
22  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
23  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
24    
25  __all__ = ['SimpleInversionCostFunction']  __all__ = [ 'InversionCostFunction']
26    
27  from costfunctions import MeteredCostFunction  from costfunctions import MeteredCostFunction
28    from mappings import Mapping
29    from forwardmodels import ForwardModel
30  from esys.escript.pdetools import ArithmeticTuple  from esys.escript.pdetools import ArithmeticTuple
31  from esys.escript import Data  from esys.escript import Data
32    import numpy as np
33    
34    class InversionCostFunction(MeteredCostFunction):
 class SimpleInversionCostFunction(MeteredCostFunction):  
35      """      """
36      This is a simple cost function with a single continuous (mapped) variable.      Class to define cost function *J(m)* for inversion with one or more forward model
37      It is the sum of two weighted terms, a single forward model and a single      based on a multi-valued level set function *m*:
38      regularization term. This cost function is used in the gravity inversion.      
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    
49        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      """      """
73      provides_inverse_Hessian_approximation=True      provides_inverse_Hessian_approximation=True
74    
75      def __init__(self, regularization, mapping, forwardmodel):      def __init__(self, regularization, mappings, forward_models):
76          """          """
77          constructor stores the supplied object references and sets default          constructor for the cost function.
78            stores the supplied object references and sets default
79          weights.          weights.
80    
81          :param regularization: The regularization part of the cost function          :param regularization: the regularization part of the cost function
82          :param mapping: Parametrization object          :type regularization: `Regularization`
83          :param forwardmodel: The forward model part of the cost function          :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          """          """
97          super(SimpleInversionCostFunction, self).__init__()          super(InversionCostFunction, self).__init__()
         self.forwardmodel=forwardmodel  
98          self.regularization=regularization          self.regularization=regularization
99          self.mapping=mapping          
100          self.setTradeOffFactors()          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    
128      def setTradeOffFactors(self, mu_model=1., mu_reg=1.):          :rtype: ``int``
129            """
130            return self.__num_tradeoff_factors
131        def getForwardModels(self):
132          """          """
133          sets the trade-off factors for the forward model and regularization          returns the forward models as a list
134          terms.          """
135            return self.forward_models
136            
137        def getRegularization(self):
138            """
139            returns the regularization
140            """
141            return self.regularization
142    
         :param mu_model: Weighting factor for the forward model (default=1.)  
         :type mu_model: non-negative `float`  
         :param mu_reg: Weighting factor for the regularization (default=1.)  
         :type mu_reg: non-negative `float`  
         """  
         if mu_model<0:  
             raise ValueError("trade-off factors must be non-negative.")  
         self.mu_model=mu_model  
         self.regularization.setTradeOffFactors(mu_reg)  
143                    
144        def setTradeOffFactorsModels(self,mu=None):
145            """
146            sets the trade-off factors for the forward model components.
147            
148            :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      def _getDualProduct(self, x, r):      def _getDualProduct(self, x, r):
240          """          """
241          returns ``regularization.getDualProduct(x, r)``          Returns the dual product, see `Regularization.getDualProduct`
242    
243            :type x: `Data`
244            :type r: `ArithmeticTuple`            
245          :rtype: `float`          :rtype: `float`
246          """          """
247          return self.regularization.getDualProduct(x, r)          return self.regularization.getDualProduct(x, r)
248    
249      def _getArguments(self, m):      def _getArguments(self, m):
250          """          """
251          returns precomputed values that are shared in the calculation of          returns pre-computed values that are shared in the calculation of
252          *f(x)* and *grad f(x)*. In this implementation returns a tuple with the          *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
253          mapped value of ``m``, the arguments from the forward model and the          mapped value of ``m``, the arguments from the forward model and the
254          arguments from the regularization.          arguments from the regularization.
255            
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          :rtype: `tuple`          :rtype: `tuple`
261          """          """
262          p=self.mapping(m)          args_reg=self.regularization.getArguments(m)
263          return p, self.forwardmodel.getArguments(p), self.regularization.getArguments(m)          # 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    
284      def _getValue(self, m, *args):      def _getValue(self, m, *args):
285          """          """
286          returns the function value at m.          Returns the value *J(m)* of the cost function at *m*.
287          If the precomputed values are not supplied `getArguments()` is called.          If the pre-computed values are not supplied `getArguments()` is called.
288    
289            :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          :rtype: `float`          :rtype: `float`
294          """          """
295          # if there is more than one forward_model and/or regularization their          # if there is more than one forward_model and/or regularization their
# Line 100  class SimpleInversionCostFunction(Metere Line 298  class SimpleInversionCostFunction(Metere
298          if len(args)==0:          if len(args)==0:
299              args=self.getArguments(m)              args=self.getArguments(m)
300                    
301          A = self.forwardmodel.getValue(args[0],*args[1])          props=args[0]
302          B = self.regularization.getValue(m, *args[2])          args_f=args[1]
303          print "J_f=",A          args_reg=args[2]
304          print "J_reg = ",B          
305          print "mu_model =",self.mu_model          J = self.regularization.getValue(m, *args_reg)
306          print "J = ", A * self.mu_model + B          print "J_reg = %e"%J
307          return   A * self.mu_model + B                  
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    
329      def _getGradient(self, m, *args):      def _getGradient(self, m, *args):
330          """          """
331          returns the gradient of *f* at *m*.          returns the gradient of the cost function  at *m*.
332          If the precomputed values are not supplied `getArguments()` is called.          If the pre-computed values are not supplied `getArguments()` is called.
333    
334          :rtype: `esys.escript.Data`          :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          """          """
         dpdm = self.mapping.getDerivative(m)  
341          if len(args)==0:          if len(args)==0:
342              args = self.getArguments(m)              args = self.getArguments(m)
343            
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    
394          Y = self.forwardmodel.getGradient(args[0],*args[1]) * dpdm          return g_J
         g_reg = self.regularization.getGradient(m, *args[2])  
         return self.mu_model * ArithmeticTuple(Y, Data()) + g_reg  
395    
396    
397      def _getInverseHessianApproximation(self, m, r, *args):      def _getInverseHessianApproximation(self, m, r, *args):
# Line 133  class SimpleInversionCostFunction(Metere Line 403  class SimpleInversionCostFunction(Metere
403          :type m: `Data`          :type m: `Data`
404          :param r: a given gradient          :param r: a given gradient
405          :type r: `ArithmeticTuple`          :type r: `ArithmeticTuple`
406          :param args: pre-calculated values for ``m`` from `getArguments()`          :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          :rtype: `Data`          :rtype: `Data`
409          :note: in the current implementation only the regularization term is          :note: in the current implementation only the regularization term is
410                 considered in the inverse Hessian approximation.                 considered in the inverse Hessian approximation.
# Line 156  class SimpleInversionCostFunction(Metere Line 427  class SimpleInversionCostFunction(Metere
427          :type m: `Data`          :type m: `Data`
428          :rtype: ``float``          :rtype: ``float``
429          """          """
430          return self.regularization.getNorm(m)          return self.regularization.getNorm(m)

Legend:
Removed from v.4121  
changed lines
  Added in v.4122

  ViewVC Help
Powered by ViewVC 1.1.26