/[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 4122 by gross, Thu Dec 20 05:42:35 2012 UTC revision 4386 by jfenwick, Tue Apr 30 02:58:02 2013 UTC
# Line 1  Line 1 
1    
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # http://www.uq.edu.au  # http://www.uq.edu.au
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
# Line 13  Line 13 
13  #  #
14  ##############################################################################  ##############################################################################
15    
16  """Collection of cost functions for inversion"""  """Cost functions for inversions with one or more forward models"""
17    
18  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19  http://www.uq.edu.au  http://www.uq.edu.au
20  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
21  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
# Line 24  __url__="https://launchpad.net/escript-f Line 24  __url__="https://launchpad.net/escript-f
24    
25  __all__ = [ 'InversionCostFunction']  __all__ = [ 'InversionCostFunction']
26    
27  from costfunctions import MeteredCostFunction  from .costfunctions import MeteredCostFunction
28  from mappings import Mapping  from .mappings import Mapping
29  from forwardmodels import ForwardModel  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  import numpy as np
33    
34    
35  class InversionCostFunction(MeteredCostFunction):  class InversionCostFunction(MeteredCostFunction):
36      """      """
37      Class to define cost function *J(m)* for inversion with one or more forward model      Class to define cost function *J(m)* for inversion with one or more
38      based on a multi-valued level set function *m*:      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)*      *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 cost function applied      where *J_reg(m)* is the regularization and cross gradient component of the
43      to a level set function *m*, *J_f(p)* are the data defect cost function involving a      cost function applied to a level set function *m*, *J_f(p)* are the data
44      physical forward model using the physical parameter(s) *p* and mu_f is trade-off factors for model f.      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 constructed from  
48      components of the level set *m* function via mappings.      A forward model depends on a set of physical parameters *p* which are
49        constructed from components of the level set function *m* via mappings.
50    
51      Example 1 (single forward model):      Example 1 (single forward model):
52           m=Mapping()           m=Mapping()
53           f=ForwardModel()           f=ForwardModel()
54           J=InversionCostFunction(Regularization(), m, f)           J=InversionCostFunction(Regularization(), m, f)
55        
56      Example 2 (two forward models on a single valued level set)      Example 2 (two forward models on a single valued level set)
57           m0=Mapping()           m0=Mapping()
58           m1=Mapping()           m1=Mapping()
# Line 58  class InversionCostFunction(MeteredCostF Line 60  class InversionCostFunction(MeteredCostF
60           f1=ForwardModel()           f1=ForwardModel()
61    
62           J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])           J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])
63            
64       Example 2 (two forward models on 2-valued level set)       Example 2 (two forward models on 2-valued level set)
65           m0=Mapping()           m0=Mapping()
66           m1=Mapping()           m1=Mapping()
67           f0=ForwardModel()           f0=ForwardModel()
68           f1=ForwardModel()           f1=ForwardModel()
69    
70           J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])             J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
71        
72      :var provides_inverse_Hessian_approximation: if true the class provides an approximative inverse of the      :cvar provides_inverse_Hessian_approximation: if true the class provides an
73                                                   Hessian operator.            approximative inverse of the Hessian operator.
74      """      """
75      provides_inverse_Hessian_approximation=True      provides_inverse_Hessian_approximation=True
76    
77      def __init__(self, regularization, mappings, forward_models):      def __init__(self, regularization, mappings, forward_models):
78          """          """
79          constructor for the cost function.          constructor for the cost function.
80          stores the supplied object references and sets default          Stores the supplied object references and sets default weights.
         weights.  
81    
82          :param regularization: the regularization part of the cost function          :param regularization: the regularization part of the cost function
83          :type regularization: `Regularization`          :type regularization: `Regularization`
84          :param mappings: the mappings to calculate physical parameters from the regularization.          :param mappings: the mappings to calculate physical parameters from the
85                           This is a list of 2-tuples *(map, i)* where the first component map defines a `Mapping`                           regularization. This is a list of 2-tuples *(map, i)*
86                           and the second component *i* defines the index of the component of level set function                           where the first component map defines a `Mapping` and
87                           to be used to calculate the mapping. An item in the list can be just a `Mapping` object                           the second component *i* defines the index of the
88                           then the entire level set function function is fed into the `Mapping` (typically used for                           component of level set function to be used to
89                             calculate the mapping. An item in the list can be just
90                             a `Mapping` object then the entire level set function
91                             function is fed into the `Mapping` (typically used for
92                           a single-component level set function.                           a single-component level set function.
93          :type mappings: `list` where each item is a `tuple` of `Mapping` and `int` or just a `Mapping`.          :type mappings: ``list`` where each item is a ``tuple`` of `Mapping`
94          :param forward_models: the forward models involved in the calculation of the cost function.                          and ``int`` or just a `Mapping`
95                                This is a list of 2-tuples *(f, ii)* where the first component map defines a `ForwardModel`          :param forward_models: the forward models involved in the calculation
96                                and the second component *ii* a list of indexes referring to the physical parameters                                 of the cost function. This is a list of 2-tuples
97                                in the `mappings` list. The  2-tuple can be replaced by a `ForwardModel` if                                 *(f, ii)* where the first component f defines a
98                                a `mappings` list as a single entry.                                 `ForwardModel` and the second component *ii* a
99          :param forward_models: `list` where each item is 'tuple` of `ForwardModel` and `list` of `int' or is  `ForwardModel`.                                 list of indexes referring to the physical
100                                   parameters in the `mappings` list. The 2-tuple
101                                   can be replaced by a `ForwardModel` if a
102                                   `mappings` list as a single entry.
103            :param forward_models: ``list`` where each item is ``tuple`` of
104                                   `ForwardModel` and ``list`` of ``int`` or is
105                                   `ForwardModel`.
106          """          """
107          super(InversionCostFunction, self).__init__()          super(InversionCostFunction, self).__init__()
108          self.regularization=regularization          self.regularization=regularization
109            
110          if isinstance(mappings, Mapping):          if isinstance(mappings, Mapping):
111           self.mappings = [mappings ]               self.mappings = [mappings ]
112      else:          else:
113           self.mappings = mappings               self.mappings = mappings
114        
115      if  isinstance(forward_models, ForwardModel):          if  isinstance(forward_models, ForwardModel):
116          self.forward_models = [ forward_models ]              self.forward_models = [ forward_models ]
117      else:              else:
118              self.forward_models=forward_models              self.forward_models=forward_models
119            
120            trafo =  regularization.getCoordinateTransformation()
121            for m in self.forward_models :
122             if not m[0].getCoordinateTransformation() == trafo:
123               raise ValueError("Coordinate transformation for regularization and model %m don't match.")
124    
125          self.numMappings=len(self.mappings)          self.numMappings=len(self.mappings)
126          self.numModels=len(self.forward_models)          self.numModels=len(self.forward_models)
127          self.numLevelSets = self.regularization.getNumLevelSets()          self.numLevelSets = self.regularization.getNumLevelSets()
128          self.__num_tradeoff_factors = self.regularization.getNumTradeOffFactors() + self.numModels          self.__num_tradeoff_factors = self.regularization.getNumTradeOffFactors() + self.numModels
129          self.setTradeOffFactorsModels()          self.setTradeOffFactorsModels()
130            
131      def getDomain(self):      def getDomain(self):
132          """          """
133          returns the domain of the cost function          returns the domain of the cost function
134          :rtype: 'Domain`          :rtype: 'Domain`
135          """          """
136          self.regularization.getDomain()          self.regularization.getDomain()
137            
138      def getNumTradeOffFactors(self):      def getNumTradeOffFactors(self):
139          """          """
140          returns the number of trade-off factors being used including the trade-off factors used in          returns the number of trade-off factors being used including the
141          the regularization component.          trade-off factors used in the regularization component.
142    
143          :rtype: ``int``          :rtype: ``int``
144          """          """
145          return self.__num_tradeoff_factors          return self.__num_tradeoff_factors
146      def getForwardModels(self):  
147          """      def getForwardModel(self, idx=None):
         returns the forward models as a list  
148          """          """
149          return self.forward_models          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):      def getRegularization(self):
164          """          """
165          returns the regularization          returns the regularization
166    
167            :rtype: `Regularization`
168          """          """
169          return self.regularization          return self.regularization
170    
171                def setTradeOffFactorsModels(self, mu=None):
     def setTradeOffFactorsModels(self,mu=None):  
172          """          """
173          sets the trade-off factors for the forward model components.          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.          :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 `float` with the length of the number of models.          :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:          if mu==None:
180              self.mu_model=np.ones((self.numModels, ))              self.mu_model=np.ones((self.numModels, ))
181          else:          else:
182          if self.numModels > 1:              if self.numModels > 1:
183                 mu=np.asarray(mu)                  mu=np.asarray(mu)
184                 if min(mu) > 0:                  if min(mu) > 0:
185            self.mu_model= mu                      self.mu_model= mu
186                 else:                  else:
187                raise ValueError("All value for trade-off factor mu must be positive.")                      raise ValueError("All values for trade-off factor mu must be positive.")
188          else:              else:
189            mu=float(mu)                  mu=float(mu)
190            if mu > 0:                  if mu > 0:
191            self.mu_model= [mu, ]                      self.mu_model= [mu, ]
192                else:                  else:
193                raise ValueError("Trade-off factor must be positive.")                      raise ValueError("Trade-off factor must be positive.")
194          
195      def setTradeOffFactorsRegularization(self,mu=None, mu_c=None):      def getTradeOffFactorsModels(self):
196          """          """
197          sets the trade of factors for the regularization component of the cost function, see          returns the trade-off factors for the forward models
198          `Regularization` for details.  
199                    :rtype: ``float`` or ``list`` of ``float``
200          :param m:  trade-off factors for the level-set variation part          """
201          :param  mu_c:  trade-off factors for the cross gradient variation part          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)          self.regularization.setTradeOffFactorsForVariation(mu)
215          self.regularization.setTradeOffFactorsForCrossGradient(mu_c)          self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
216            
217      def setTradeOffFactors(self, mu=None):      def setTradeOffFactors(self, mu=None):
218      """          """
219      sets the trade-off factors for the forward model and regularization          sets the trade-off factors for the forward model and regularization
220      terms.          terms.
221    
222      :param mu_model: Weighting factor for the forward model (default=1.)          :param mu: list of trade-off factors.
223      :type mu_model: non-negative `float`          :type mu: ``list`` of ``float``
224      :param mu: list of trade-off factors.          """
225      :type mu: `list` of `float`          if mu is None:
226      """              mu=np.ones((self.__num_tradeoff_factors,))
227      if mu==None:          self.setTradeOffFactorsModels(mu[:self.numModels])
228          mu=mp.ones((self.__num_tradeoff_factors,))          self.regularization.setTradeOffFactors(mu[self.numModels:])
229      self.setTradeOffFactorsModels(mu[:self.numModels])  
230      self.regularization.setTradeOffFactors(mu[self.numModels:])      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):      def createLevelSetFunction(self, *props):
241          """          """
242          return an instance of an object used to represent a level set function initialed          returns an instance of an object used to represent a level set function
243          with zeros. Components can be overwritten by physical properties 'props'. If present          initialized with zeros. Components can be overwritten by physical
244          entries must correspond to the a `mappings` arguments in the constructor. Use `None`          properties `props`. If present entries must correspond to the
245          for properties for which no value is given.          `mappings` arguments in the constructor. Use ``None`` for properties
246            for which no value is given.
247          """          """
248          m=self.regularization.getPDE().createSolution()          m=self.regularization.getPDE().createSolution()
249          if len(props) > 0:          if len(props) > 0:
250             for i in xrange(self.numMappings):              for i in range(self.numMappings):
251                if props[i]:                  if props[i]:
252            mm=self.mappings[i]                      mm=self.mappings[i]
253            if isinstance(mm, Mapping):                      if isinstance(mm, Mapping):
254                m=mm.getInverse(props[i])                          m=mm.getInverse(props[i])
255            elif len(mm) == 1:                      elif len(mm) == 1:
256                m=mm[0].getInverse(props[i])                          m=mm[0].getInverse(props[i])
257            else:                      else:
258                m[mm[1]]=mm[0].getInverse(props[i])                          m[mm[1]]=mm[0].getInverse(props[i])
259          return m          return m
260        
261      def getProperties(self, m, return_list=False):      def getProperties(self, m, return_list=False):
262          """          """
263          returns a list of the physical properties from a given level set function *m* using the          returns a list of the physical properties from a given level set
264          mappings of the cost function.          function *m* using the mappings of the cost function.
265            
266          :param m: level set function          :param m: level set function
267          :type m: `Data`          :type m: `Data`
268          :param return_list: if True a list is returned.          :param return_list: if ``True`` a list is returned.
269          :type return_list: `bool`          :type return_list: ``bool``
270          :rtype m: `list` of `Data`          :rtype: ``list`` of `Data`
271          """          """
272          props=[]          props=[]
273          for i in xrange(self.numMappings):          for i in range(self.numMappings):
274             mm=self.mappings[i]              mm=self.mappings[i]
275             if isinstance(mm, Mapping):              if isinstance(mm, Mapping):
276             p=mm.getValue(m)                  p=mm.getValue(m)
277         elif len(mm) == 1:              elif len(mm) == 1:
278             p=mm[0].getValue(m)                  p=mm[0].getValue(m)
279         else:              else:
280             p=mm[0].getValue(m[mm[1]])                  p=mm[0].getValue(m[mm[1]])
281             props.append(p)              props.append(p)
282          if self.numMappings > 1 or return_list:          if self.numMappings > 1 or return_list:
283         return props              return props
284      else:          else:
285         return props[0]              return props[0]
286              
287      def _getDualProduct(self, x, r):      def _getDualProduct(self, x, r):
288          """          """
289          Returns the dual product, see `Regularization.getDualProduct`          Returns the dual product, see `Regularization.getDualProduct`
290    
291          :type x: `Data`          :type x: `Data`
292          :type r: `ArithmeticTuple`                      :type r: `ArithmeticTuple`
293          :rtype: `float`          :rtype: ``float``
294          """          """
295          return self.regularization.getDualProduct(x, r)          return self.regularization.getDualProduct(x, r)
296    
# Line 252  class InversionCostFunction(MeteredCostF Line 300  class InversionCostFunction(MeteredCostF
300          *J(m)* and *grad J(m)*. In this implementation returns a tuple with the          *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          mapped value of ``m``, the arguments from the forward model and the
302          arguments from the regularization.          arguments from the regularization.
303            
304          :param m: current approximation of the level set function          :param m: current approximation of the level set function
305          :type m: `Data`          :type m: `Data`
306          :return: tuple of of values of the parameters, pre-computed values for the forward model and          :return: tuple of of values of the parameters, pre-computed values
307                   pre-computed values for the regularization                   for the forward model and pre-computed values for the
308          :rtype: `tuple`                   regularization
309            :rtype: ``tuple``
310          """          """
311          args_reg=self.regularization.getArguments(m)          args_reg=self.regularization.getArguments(m)
312          # cache for physical parameters:          # cache for physical parameters:
313          props=self.getProperties(m, return_list=True)          props=self.getProperties(m, return_list=True)
314          args_f=[]          args_f=[]
315          for i in xrange(self.numModels):          for i in range(self.numModels):
316         f=self.forward_models[i]              f=self.forward_models[i]
317         if isinstance(f, ForwardModel):              if isinstance(f, ForwardModel):
318            aa=f.getArguments(props[0])                  aa=f.getArguments(props[0])
319         elif len(f) == 1:              elif len(f) == 1:
320            aa=f[0].getArguments(props[0])                  aa=f[0].getArguments(props[0])
321         else:              else:
322            idx = f[1]                  idx = f[1]
323            f=f[0]                  f=f[0]
324            if isinstance(idx, int):                  if isinstance(idx, int):
325           aa=f.getArguments(props[idx])                      aa=f.getArguments(props[idx])
326            else:                  else:
327           pp=tuple( [ props[i] for i in idx] )                      pp=tuple( [ props[i] for i in idx] )
328           aa=f.getArguments(*pp)                      aa=f.getArguments(*pp)
329         args_f.append(aa)              args_f.append(aa)
330          
331          return props, args_f, args_reg          return props, args_f, args_reg
332    
333      def _getValue(self, m, *args):      def _getValue(self, m, *args):
# Line 288  class InversionCostFunction(MeteredCostF Line 337  class InversionCostFunction(MeteredCostF
337    
338          :param m: current approximation of the level set function          :param m: current approximation of the level set function
339          :type m: `Data`          :type m: `Data`
340          :param args: tuple of of values of the parameters, pre-computed values for the forward model and          :param args: tuple of of values of the parameters, pre-computed values
341                   pre-computed values for the regularization                       for the forward model and pre-computed values for the
342          :rtype: `float`                       regularization
343          """          :rtype: ``float``
344          # if there is more than one forward_model and/or regularization their          """
         # contributions need to be added up. But this implementation allows  
         # only one of each...  
345          if len(args)==0:          if len(args)==0:
346              args=self.getArguments(m)              args=self.getArguments(m)
347            
348          props=args[0]          props=args[0]
349          args_f=args[1]          args_f=args[1]
350          args_reg=args[2]          args_reg=args[2]
351            
352          J = self.regularization.getValue(m, *args_reg)          J = self.regularization.getValue(m, *args_reg)
353          print "J_reg = %e"%J  
354                            for i in range(self.numModels):
355          for i in xrange(self.numModels):              f=self.forward_models[i]
356                        if isinstance(f, ForwardModel):
357         f=self.forward_models[i]                  J_f = f.getValue(props[0],*args_f[i])
358         if isinstance(f, ForwardModel):              elif len(f) == 1:
359                J_f = f.getValue(props[0],*args_f[i])                  J_f=f[0].getValue(props[0],*args_f[i])
360             elif len(f) == 1:              else:
361            J_f=f[0].getValue(props[0],*args_f[i])                  idx = f[1]
362         else:                  f=f[0]
363            idx = f[1]                  if isinstance(idx, int):
364            f=f[0]                      J_f = f.getValue(props[idx],*args_f[i])
365            if isinstance(idx, int):                  else:
366           J_f = f.getValue(props[idx],*args_f[i])                      args=tuple( [ props[j] for j in idx] + args_f[i])
367            else:                      J_f = f.getValue(*args)
368           args=tuple( [ props[j] for j in idx] + args_f[i])              self.logger.debug("J_f[%d] = %e"%(i, J_f))
369           J_f = f.getValue(*args)              self.logger.debug("mu_model[%d] = %e"%(i, self.mu_model[i]))
370             print "J_f[%d] = %e"%(i, J_f)              J += self.mu_model[i] * J_f
371             print "mu_model[%d] = %e"%(i, self.mu_model[i])  
372             J += self.mu_model[i] * J_f          return J
373              
374          return   J      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):      def _getGradient(self, m, *args):
418          """          """
419          returns the gradient of the cost function  at *m*.          returns the gradient of the cost function at *m*.
420          If the pre-computed values are not supplied `getArguments()` is called.          If the pre-computed values are not supplied `getArguments()` is called.
421    
422          :param m: current approximation of the level set function          :param m: current approximation of the level set function
423          :type m: `Data`          :type m: `Data`
424          :param args: tuple of of values of the parameters, pre-computed values for the forward model and          :param args: tuple of values of the parameters, pre-computed values
425                   pre-computed values for the regularization                       for the forward model and pre-computed values for the
426                                         regularization
427    
428          :rtype: `ArithmeticTuple`          :rtype: `ArithmeticTuple`
429          """          """
430          if len(args)==0:          if len(args)==0:
431              args = self.getArguments(m)              args = self.getArguments(m)
432            
433          props=args[0]          props=args[0]
434          args_f=args[1]          args_f=args[1]
435          args_reg=args[2]          args_reg=args[2]
436            
437          g_J = self.regularization.getGradient(m, *args_reg)          g_J = self.regularization.getGradient(m, *args_reg)
438          p_diffs=[]          p_diffs=[]
439          for i in xrange(self.numMappings):          for i in range(self.numMappings):
440             mm=self.mappings[i]              mm=self.mappings[i]
441             if isinstance(mm, Mapping):              if isinstance(mm, Mapping):
442             dpdm = mm.getDerivative(m)                  dpdm = mm.getDerivative(m)
443         elif len(mm) == 1:              elif len(mm) == 1:
444             dpdm = mm[0].getDerivative(m)                  dpdm = mm[0].getDerivative(m)
445         else:              else:
446             dpdm = mm[0].getDerivative(m[mm[1]])                  dpdm = mm[0].getDerivative(m[mm[1]])
447             p_diffs.append(dpdm)              p_diffs.append(dpdm)
448              
449          Y=g_J[0]            Y=g_J[0]
450          for i in xrange(self.numModels):          for i in range(self.numModels):
451         mu=self.mu_model[i]              mu=self.mu_model[i]
452         f=self.forward_models[i]              f=self.forward_models[i]
453         if isinstance(f, ForwardModel):              if isinstance(f, ForwardModel):
454            Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu                  Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
455            if self.numLevelSets == 1 :                  if self.numLevelSets == 1 :
456               Y +=Ys                      Y +=Ys
457            else:                  else:
458                    Y[0] +=Ys                      Y[0] +=Ys
459             elif len(f) == 1:              elif len(f) == 1:
460            Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0]  * mu                  Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0]  * mu
461                if self.numLevelSets == 1 :                  if self.numLevelSets == 1 :
462               Y +=Ys                      Y +=Ys
463            else:                  else:
464                    Y[0] +=Ys                      Y[0] +=Ys
465         else:              else:
466            idx = f[1]                  idx = f[1]
467            f=f[0]                  f = f[0]
468            if isinstance(idx, int):                  if isinstance(idx, int):
469           Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu                      Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
470           if self.numLevelSets == 1 :                      if self.numLevelSets == 1 :
471               if idx == 0:                          if idx == 0:
472                   Y+=Ys                              Y+=Ys
473               else:                          else:
474                   raise IndexError("Illegal mapping index.")                              raise IndexError("Illegal mapping index.")
475           else:                      else:
476               Y[idx] += Ys                          Y[idx] += Ys
477            else:                  else:
478           args=tuple( [ props[j] for j in idx] + args_f[i])                      args = tuple( [ props[j] for j in idx] + args_f[i])
479           Ys = f.getGradient(*args)                      Ys = f.getGradient(*args)
480           for ii in xrange(len(idx)):                      for ii in range(len(idx)):
481               Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]]  * mu                          Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
482    
483          return g_J          return g_J
484    
   
485      def _getInverseHessianApproximation(self, m, r, *args):      def _getInverseHessianApproximation(self, m, r, *args):
486          """          """
487          returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function          returns an approximative evaluation *p* of the inverse of the Hessian
488          for a given gradient type *r* at a given location *m*: *H(m) p = r*          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          :param m: level set approximation where to calculate Hessian inverse
492          :type m: `Data`          :type m: `Data`
493          :param r: a given gradient          :param r: a given gradient
494          :type r: `ArithmeticTuple`          :type r: `ArithmeticTuple`
495          :param args: tuple of of values of the parameters, pre-computed values for the forward model and          :param args: tuple of values of the parameters, pre-computed values
496                   pre-computed values for the regularization                       for the forward model and pre-computed values for the
497                         regularization
498          :rtype: `Data`          :rtype: `Data`
499          :note: in the current implementation only the regularization term is          :note: in the current implementation only the regularization term is
500                 considered in the inverse Hessian approximation.                 considered in the inverse Hessian approximation.
   
501          """          """
502          m=self.regularization.getInverseHessianApproximation(m, r, *args[2])          m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
503          return m          return m
# Line 418  class InversionCostFunction(MeteredCostF Line 507  class InversionCostFunction(MeteredCostF
507          notifies the class that the Hessian operator needs to be updated.          notifies the class that the Hessian operator needs to be updated.
508          """          """
509          self.regularization.updateHessian()          self.regularization.updateHessian()
510            
511      def _getNorm(self, m):      def _getNorm(self, m):
512          """          """
513          returns the norm of ``m``          returns the norm of `m`
514    
515          :param m: level set function          :param m: level set function
516          :type m: `Data`          :type m: `Data`
517          :rtype: ``float``          :rtype: ``float``
518          """          """
         return self.regularization.getNorm(m)  
519            return self.regularization.getNorm(m)
520    

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

  ViewVC Help
Powered by ViewVC 1.1.26