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

revision 4079 by gross, Fri Nov 16 07:59:01 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  #  #
# 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
19  http://www.uq.edu.au  http://www.uq.edu.au
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
35  class SimpleInversionCostFunction(MeteredCostFunction):  class InversionCostFunction(MeteredCostFunction):
36      """      """
37      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
38      It is the sum of two weighted terms, a single forward model and a single      forward models based on a multi-valued level set function *m*:
39      regularization term. This cost function is used in the gravity inversion.
40        *J(m) = J_reg(m) + sum_f mu_f * J_f(p)*
41
42        where *J_reg(m)* is the regularization and cross gradient component of the
43        cost function applied to a level set function *m*, *J_f(p)* are the data
44        defect cost functions involving a physical forward model using the
45        physical parameter(s) *p* and *mu_f* is the trade-off factor for model f.
46
47
48        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):
52             m=Mapping()
53             f=ForwardModel()
54             J=InversionCostFunction(Regularization(), m, f)
55
56        Example 2 (two forward models on a single valued level set)
57             m0=Mapping()
58             m1=Mapping()
59             f0=ForwardModel()
60             f1=ForwardModel()
61
62             J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])
63
64         Example 2 (two forward models on 2-valued level set)
65             m0=Mapping()
66             m1=Mapping()
67             f0=ForwardModel()
68             f1=ForwardModel()
69
70             J=InversionCostFunction(Regularization(numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
71
72        :cvar provides_inverse_Hessian_approximation: if true the class provides an
73              approximative inverse of the Hessian operator.
74      """      """
75      def __init__(self, regularization, mapping, forwardmodel):      provides_inverse_Hessian_approximation=True
76
77        def __init__(self, regularization, mappings, forward_models):
78          """          """
79          constructor stores the supplied object references and sets default          constructor for the cost function.
80          weights.          Stores the supplied object references and sets default weights.
81
82          :param regularization: The regularization part of the cost function          :param regularization: the regularization part of the cost function
83          :param mapping: Parametrization object          :type regularization: `Regularization`
84          :param forwardmodel: The forward model part of the cost function          :param mappings: the mappings to calculate physical parameters from the
85                             regularization. This is a list of 2-tuples *(map, i)*
86                             where the first component map defines a `Mapping` and
87                             the second component *i* defines the index of the
88                             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.
93            :type mappings: ``list`` where each item is a ``tuple`` of `Mapping`
94                            and ``int`` or just a `Mapping`
95            :param forward_models: the forward models involved in the calculation
96                                   of the cost function. This is a list of 2-tuples
97                                   *(f, ii)* where the first component f defines a
98                                   `ForwardModel` and the second component *ii* a
99                                   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(SimpleInversionCostFunction, self).__init__()          super(InversionCostFunction, self).__init__()
self.forwardmodel=forwardmodel
108          self.regularization=regularization          self.regularization=regularization
self.mapping=mapping
self.setWeights()
109
110      def setWeights(self, mu_model=1., mu_reg_0=None,mu_reg_1=None):          if isinstance(mappings, Mapping):
111                 self.mappings = [mappings ]
112            else:
113                 self.mappings = mappings
114
115            if  isinstance(forward_models, ForwardModel):
116                self.forward_models = [ forward_models ]
117            else:
118                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)
126            self.numModels=len(self.forward_models)
127            self.numLevelSets = self.regularization.getNumLevelSets()
130
131        def getDomain(self):
132            """
133            returns the domain of the cost function
134            :rtype: 'Domain`
135            """
136            self.regularization.getDomain()
137
139            """
140            returns the number of trade-off factors being used including the
141            trade-off factors used in the regularization component.
142
143            :rtype: ``int``
144            """
146
147        def getForwardModel(self, idx=None):
148            """
149            returns the *idx*-th forward model.
150
151            :param idx: model index. If cost function contains one model only `idx`
152                        can be omitted.
153            :type idx: ``int``
154            """
155            if idx==None: idx=0
156            f=self.forward_models[idx]
157            if isinstance(f, ForwardModel):
158                F=f
159            else:
160                F=f[0]
161            return F
162
163        def getRegularization(self):
164            """
165            returns the regularization
166
167            :rtype: `Regularization`
168            """
169            return self.regularization
170
172            """
173            sets the trade-off factors for the forward model components.
174
175            :param mu: list of the trade-off factors. If not present ones are used.
176            :type mu: ``float`` in case of a single model or a ``list`` of
177                      ``float`` with the length of the number of models.
178            """
179            if mu==None:
180                self.mu_model=np.ones((self.numModels, ))
181            else:
182                if self.numModels > 1:
183                    mu=np.asarray(mu)
184                    if min(mu) > 0:
185                        self.mu_model= mu
186                    else:
187                        raise ValueError("All values for trade-off factor mu must be positive.")
188                else:
189                    mu=float(mu)
190                    if mu > 0:
191                        self.mu_model= [mu, ]
192                    else:
193                        raise ValueError("Trade-off factor must be positive.")
194
196          """          """
197          sets the weighting factors for the forward model and regularization          returns the trade-off factors for the forward models
198
199            :rtype: ``float`` or ``list`` of ``float``
200            """
201            if self.numModels>1:
202                return self.mu_model
203            else:
204                return self.mu_model[0]
205
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
213            """
216
218            """
219            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_reg_0: Weighting factor for the regularization (default=1.)          """
225          :type mu_reg_0: non-negative `float`          if mu is None:
228              raise ValueError("weighting factors must be non-negative.")          self.regularization.setTradeOffFactors(mu[self.numModels:])
229          self.mu_model=mu_model
231          self.regularization.setWeightsForS1(mu_reg_1)          """
232                    returns a list of the trade-off factors.
233
234            :rtype: ``list`` of ``float``
235            """
238            return [ m for m in mu1] + [ m for m in mu2]
239
240        def createLevelSetFunction(self, *props):
241            """
242            returns an instance of an object used to represent a level set function
243            initialized with zeros. Components can be overwritten by physical
244            properties `props`. If present entries must correspond to the
245            `mappings` arguments in the constructor. Use ``None`` for properties
246            for which no value is given.
247            """
248            m=self.regularization.getPDE().createSolution()
249            if len(props) > 0:
250                for i in range(self.numMappings):
251                    if props[i]:
252                        mm=self.mappings[i]
253                        if isinstance(mm, Mapping):
254                            m=mm.getInverse(props[i])
255                        elif len(mm) == 1:
256                            m=mm[0].getInverse(props[i])
257                        else:
258                            m[mm[1]]=mm[0].getInverse(props[i])
259            return m
260
261        def getProperties(self, m, return_list=False):
262            """
263            returns a list of the physical properties from a given level set
264            function *m* using the mappings of the cost function.
265
266            :param m: level set function
267            :type m: `Data`
268            :param return_list: if ``True`` a list is returned.
269            :type return_list: ``bool``
270            :rtype: ``list`` of `Data`
271            """
272            props=[]
273            for i in range(self.numMappings):
274                mm=self.mappings[i]
275                if isinstance(mm, Mapping):
276                    p=mm.getValue(m)
277                elif len(mm) == 1:
278                    p=mm[0].getValue(m)
279                else:
280                    p=mm[0].getValue(m[mm[1]])
281                props.append(p)
282            if self.numMappings > 1 or return_list:
283                return props
284            else:
285                return props[0]
286
287      def _getDualProduct(self, x, r):      def _getDualProduct(self, x, r):
288          """          """
289          returns ``regularization.getDualProduct(x, r)``          Returns the dual product, see `Regularization.getDualProduct`
290
291          :rtype: `float`          :type x: `Data`
292            :type r: `ArithmeticTuple`
293            :rtype: ``float``
294          """          """
295          return self.regularization.getDualProduct(x, r)          return self.regularization.getDualProduct(x, r)
296
297      def _getArguments(self, m):      def _getArguments(self, m):
298          """          """
299          returns precalculated values that are shared in the calculation of          returns pre-computed values that are shared in the calculation of
300          *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
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          :rtype: `tuple`          :param m: current approximation of the level set function
305          """          :type m: `Data`
306          p=self.mapping(m)          :return: tuple of of values of the parameters, pre-computed values
307          return p, self.forwardmodel.getArguments(p), self.regularization.getArguments(m)                   for the forward model and pre-computed values for the
308                     regularization
309            :rtype: ``tuple``
310            """
311            args_reg=self.regularization.getArguments(m)
312            # cache for physical parameters:
313            props=self.getProperties(m, return_list=True)
314            args_f=[]
315            for i in range(self.numModels):
316                f=self.forward_models[i]
317                if isinstance(f, ForwardModel):
318                    aa=f.getArguments(props[0])
319                elif len(f) == 1:
320                    aa=f[0].getArguments(props[0])
321                else:
322                    idx = f[1]
323                    f=f[0]
324                    if isinstance(idx, int):
325                        aa=f.getArguments(props[idx])
326                    else:
327                        pp=tuple( [ props[i] for i in idx] )
328                        aa=f.getArguments(*pp)
329                args_f.append(aa)
330
331            return props, args_f, args_reg
332
333      def _getValue(self, m, *args):      def _getValue(self, m, *args):
334          """          """
335          returns the function value at m.          Returns the value *J(m)* of the cost function at *m*.
336          If the precalculated values are not supplied `getArguments()` is called.          If the pre-computed values are not supplied `getArguments()` is called.
337
338            :param m: current approximation of the level set function
339            :type m: `Data`
340            :param args: tuple of of values of the parameters, pre-computed values
341                         for the forward model and pre-computed values for the
342                         regularization
343            :rtype: ``float``
344            """
345            if len(args)==0:
346                args=self.getArguments(m)
347
348          :rtype: `float`          props=args[0]
349            args_f=args[1]
350            args_reg=args[2]
351
352            J = self.regularization.getValue(m, *args_reg)
353
354            for i in range(self.numModels):
355                f=self.forward_models[i]
356                if isinstance(f, ForwardModel):
357                    J_f = f.getValue(props[0],*args_f[i])
358                elif len(f) == 1:
359                    J_f=f[0].getValue(props[0],*args_f[i])
360                else:
361                    idx = f[1]
362                    f=f[0]
363                    if isinstance(idx, int):
364                        J_f = f.getValue(props[idx],*args_f[i])
365                    else:
366                        args=tuple( [ props[j] for j in idx] + args_f[i])
367                        J_f = f.getValue(*args)
368                self.logger.debug("J_f[%d] = %e"%(i, J_f))
369                self.logger.debug("mu_model[%d] = %e"%(i, self.mu_model[i]))
370                J += self.mu_model[i] * J_f
371
372            return J
373
374        def getComponentValues(self, m, *args):
375            return self._getComponentValues(m, *args)
376
377        def _getComponentValues(self, m, *args):
378            """
379            returns the values of the individual cost functions that make up *f(x)*
380            using the precalculated values for *x*.
381
382            :param x: a solution approximation
383            :type x: x-type
384            :rtype: ``list<<float>>``
385          """          """
# 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...
386          if len(args)==0:          if len(args)==0:
387              args=self.getArguments(m)              args=self.getArguments(m)
388
389          return  self.mu_model * self.forwardmodel.getValue(args[0],*args[1]) \          props=args[0]
390                 +  self.regularization.getValue(m, *args[2])          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
418          """          """
419          returns the gradient of *f* at *m*.          returns the gradient of the cost function at *m*.
420          If the precalculated 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
423            :type m: `Data`
424            :param args: tuple of values of the parameters, pre-computed values
425                         for the forward model and pre-computed values for the
426                         regularization
427
428          :rtype: `esys.escript.Data`          :rtype: `ArithmeticTuple`
429          """          """
dpdm = self.mapping.getDerivative(m)
430          if len(args)==0:          if len(args)==0:
431              args = self.getArguments(m)              args = self.getArguments(m)
432
433          Y = self.forwardmodel.getGradient(args[0],*args[1]) * dpdm          props=args[0]
434          g_reg = self.regularization.getGradient(m, *args[2])          args_f=args[1]
435          print "grad forward = ", Y          args_reg=args[2]
436          print "grad regularization Y  = ", g_reg[0]
437          print "grad regularization X = ", g_reg[1]          g_J = self.regularization.getGradient(m, *args_reg)
438                    p_diffs=[]
439                    for i in range(self.numMappings):
440          return self.mu_model * ArithmeticTuple(Y, Data()) + g_reg              mm=self.mappings[i]
441                if isinstance(mm, Mapping):
442                    dpdm = mm.getDerivative(m)
443                elif len(mm) == 1:
444                    dpdm = mm[0].getDerivative(m)
445                else:
446                    dpdm = mm[0].getDerivative(m[mm[1]])
447                p_diffs.append(dpdm)
448
449            Y=g_J[0]
450            for i in range(self.numModels):
451                mu=self.mu_model[i]
452                f=self.forward_models[i]
453                if isinstance(f, ForwardModel):
454                    Ys= f.getGradient(props[0],*args_f[i]) * p_diffs[0] * mu
455                    if self.numLevelSets == 1 :
456                        Y +=Ys
457                    else:
458                        Y[0] +=Ys
459                elif len(f) == 1:
460                    Ys=f[0].getGradient(props[0],*args_f[i]) * p_diffs[0]  * mu
461                    if self.numLevelSets == 1 :
462                        Y +=Ys
463                    else:
464                        Y[0] +=Ys
465                else:
466                    idx = f[1]
467                    f = f[0]
468                    if isinstance(idx, int):
469                        Ys = f.getGradient(props[idx],*args_f[i]) * p_diffs[idx] * mu
470                        if self.numLevelSets == 1 :
471                            if idx == 0:
472                                Y+=Ys
473                            else:
474                                raise IndexError("Illegal mapping index.")
475                        else:
476                            Y[idx] += Ys
477                    else:
478                        args = tuple( [ props[j] for j in idx] + args_f[i])
480                        for ii in range(len(idx)):
481                            Y[idx[ii]]+=Ys[ii]* p_diffs[idx[ii]] * mu
482
483            return g_J
484
485      def _getInverseHessianApproximation(self, m, r, *args):      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`
494          :type r: ``ArithmeticTuple``          :type r: `ArithmeticTuple`
495          :param args: pre-calculated values for ``m`` from ``getArguments()``          :param args: tuple of values of the parameters, pre-computed values
496          :rtype: ``Data``                       for the forward model and pre-computed values for the
497          :note: in the current implementation only the regularization term is considered in the                       regularization
498          inverse Hessian approximation.          :rtype: `Data`
499          """          :note: in the current implementation only the regularization term is
500          print "nverseHessianApproximation:"                 considered in the inverse Hessian approximation.
501          print "Y  = ",r[0]          """
print "X  = ",r[1]
502          m=self.regularization.getInverseHessianApproximation(m, r, *args[2])          m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
print "m  = ",m
503          return m          return m
504
505      def updateHessian(self):      def updateHessian(self):
506          """          """
507          notifies the class that the Hessian operator needs to be updated.          notifies the class that the Hessian operator needs to be updated.
508          """          """
self.regularization.updateHessian()
509            self.regularization.updateHessian()
510
511        def _getNorm(self, m):
512            """
513            returns the norm of `m`
514
515            :param m: level set function
516            :type m: `Data`
517            :rtype: ``float``
518            """
519            return self.regularization.getNorm(m)
520

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