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

Revision 4074 - (hide annotations)
Thu Nov 15 03:30:59 2012 UTC (6 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 5753 byte(s)
```this will break the downundermodule: major revision of the classes.
```
 1 gross 4074 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2012 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 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 23 __url__= 24 25 __all__ = ['SimpleInversionCostFunction'] 26 27 from costfunctions import MeteredCostFunction 28 from esys.escript.pdetools import ArithmeticTuple 29 from esys.escript import Data 30 31 32 class SimpleInversionCostFunction(MeteredCostFunction): 33 """ 34 This is a simple cost function with a single continuous (mapped) variable. 35 It is the sum of two weighted terms, a single forward model and a single 36 regularization term. This cost function is used in the gravity inversion. 37 """ 38 def __init__(self, regularization, mapping, forwardmodel): 39 """ 40 constructor stores the supplied object references and sets default 41 weights. 42 43 :param regularization: The regularization part of the cost function 44 :param mapping: Parametrization object 45 :param forwardmodel: The forward model part of the cost function 46 """ 47 super(SimpleInversionCostFunction, self).__init__() 48 self.forwardmodel=forwardmodel 49 self.regularization=regularization 50 self.mapping=mapping 51 self.setWeights() 52 53 def setWeights(self, mu_model=1., mu_reg_0=None,mu_reg_1=None): 54 """ 55 sets the weighting factors for the forward model and regularization 56 terms. 57 58 :param mu_model: Weighting factor for the forward model (default=1.) 59 :type mu_model: non-negative `float` 60 :param mu_reg_0: Weighting factor for the regularization (default=1.) 61 :type mu_reg_0: non-negative `float` 62 """ 63 if mu_model<0: 64 raise ValueError("weighting factors must be non-negative.") 65 self.mu_model=mu_model 66 self.regularization.setWeightsForS0(mu_reg_0) 67 self.regularization.setWeightsForS1(mu_reg_1) 68 69 def _getDualProduct(self, x, r): 70 """ 71 returns ``regularization.getDualProduct(x, r)`` 72 73 :rtype: `float` 74 """ 75 return self.regularization.getDualProduct(x, r) 76 77 def _getArguments(self, m): 78 """ 79 returns precalculated values that are shared in the calculation of 80 *f(x)* and *grad f(x)*. In this implementation returns a tuple with the 81 mapped value of ``m``, the arguments from the forward model and the 82 arguments from the regularization. 83 84 :rtype: `tuple` 85 """ 86 p=self.mapping(m) 87 return p, self.forwardmodel.getArguments(p), self.regularization.getArguments(m) 88 89 def _getValue(self, m, *args): 90 """ 91 returns the function value at m. 92 If the precalculated values are not supplied `getArguments()` is called. 93 94 :rtype: `float` 95 """ 96 # if there is more than one forward_model and/or regularization their 97 # contributions need to be added up. But this implementation allows 98 # only one of each... 99 if len(args)==0: 100 args=self.getArguments(m) 101 return self.mu_model * self.forwardmodel.getValue(args[0],*args[1]) \ 102 + self.regularization.getValue(m, *args[2]) 103 104 def _getGradient(self, m, *args): 105 """ 106 returns the gradient of *f* at *m*. 107 If the precalculated values are not supplied `getArguments()` is called. 108 109 :rtype: `esys.escript.Data` 110 """ 111 dpdm = self.mapping.getDerivative(m) 112 if len(args)==0: 113 args = self.getArguments(m) 114 115 Y = self.forwardmodel.getGradient(args[0],*args[1]) * dpdm 116 g_reg = self.regularization.getGradient(m, *args[2]) 117 118 return self.mu_model * ArithmeticTuple(Y, Data()) + g_reg 119 120 121 def _getDirectionalDerivative(self, m, d, *args): 122 """ 123 returns the directional derivative at *m* in direction *d*. 124 125 :rtype: `float` 126 """ 127 dpdm = self.mapping.getDerivative(m) 128 return self.mu_model* self.forwardmodel.getDirectionalDerivative(args[0], dpdm*d,*args[1]) + \ 129 getDirectionalDerivative(self, m, d, *args[2]) 130 131 def _getInverseHessianApproximation(self, m, r, *args): 132 """ 133 returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function 134 for a given gradient type *r* at a given location *m*: *H(m) p = r* 135 136 :param m: level set approximation where to calculate Hessian inverse 137 :type m: ``Data`` 138 :param r: a given gradient 139 :type r: ``ArithmeticTuple`` 140 :param args: pre-calculated values for ``m`` from ``getArguments()`` 141 :rtype: ``Data`` 142 :note: in the current implementation only the regularization term is considered in the 143 inverse Hessian approximation. 144 """ 145 return self.regularization.getInverseHessianApproximation(m, r, *args[2]) 146 147 def updateHessian(self): 148 """ 149 notifies the class that the Hessian operator needs to be updated. 150 """ 151 self.regularization.updateHessian()