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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4099 - (show annotations)
Tue Dec 11 04:04:47 2012 UTC (6 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 5585 byte(s)
Clarification of terminology and scaling.
Joint inversion is not working yet.

1
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 __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 provides_inverse_Hessian_approximation=True
39
40 def __init__(self, regularization, mapping, forwardmodel):
41 """
42 constructor stores the supplied object references and sets default
43 weights.
44
45 :param regularization: The regularization part of the cost function
46 :param mapping: Parametrization object
47 :param forwardmodel: The forward model part of the cost function
48 """
49 super(SimpleInversionCostFunction, self).__init__()
50 self.forwardmodel=forwardmodel
51 self.regularization=regularization
52 self.mapping=mapping
53 self.setTradeOffFactors()
54
55 def setTradeOffFactors(self, mu_model=1., mu_reg=1.):
56 """
57 sets the trade-off factors for the forward model and regularization
58 terms.
59
60 :param mu_model: Weighting factor for the forward model (default=1.)
61 :type mu_model: non-negative `float`
62 :param mu_reg: Weighting factor for the regularization (default=1.)
63 :type mu_reg: non-negative `float`
64 """
65 if mu_model<0:
66 raise ValueError("trade-off factors must be non-negative.")
67 self.mu_model=mu_model
68 self.regularization.setTradeOffFactors(mu_reg)
69
70 def _getDualProduct(self, x, r):
71 """
72 returns ``regularization.getDualProduct(x, r)``
73
74 :rtype: `float`
75 """
76 return self.regularization.getDualProduct(x, r)
77
78 def _getArguments(self, m):
79 """
80 returns precomputed values that are shared in the calculation of
81 *f(x)* and *grad f(x)*. In this implementation returns a tuple with the
82 mapped value of ``m``, the arguments from the forward model and the
83 arguments from the regularization.
84
85 :rtype: `tuple`
86 """
87 p=self.mapping(m)
88 return p, self.forwardmodel.getArguments(p), self.regularization.getArguments(m)
89
90 def _getValue(self, m, *args):
91 """
92 returns the function value at m.
93 If the precomputed values are not supplied `getArguments()` is called.
94
95 :rtype: `float`
96 """
97 # if there is more than one forward_model and/or regularization their
98 # contributions need to be added up. But this implementation allows
99 # only one of each...
100 if len(args)==0:
101 args=self.getArguments(m)
102
103 return self.mu_model * self.forwardmodel.getValue(args[0],*args[1]) \
104 + self.regularization.getValue(m, *args[2])
105
106 def _getGradient(self, m, *args):
107 """
108 returns the gradient of *f* at *m*.
109 If the precomputed values are not supplied `getArguments()` is called.
110
111 :rtype: `esys.escript.Data`
112 """
113 dpdm = self.mapping.getDerivative(m)
114 if len(args)==0:
115 args = self.getArguments(m)
116
117 Y = self.forwardmodel.getGradient(args[0],*args[1]) * dpdm
118 g_reg = self.regularization.getGradient(m, *args[2])
119 print "grad forward = ", Y
120 print "grad regularization Y = ", g_reg[0]
121 print "grad regularization X = ", g_reg[1]
122
123 return self.mu_model * ArithmeticTuple(Y, Data()) + g_reg
124
125
126 def _getInverseHessianApproximation(self, m, r, *args):
127 """
128 returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
129 for a given gradient type *r* at a given location *m*: *H(m) p = r*
130
131 :param m: level set approximation where to calculate Hessian inverse
132 :type m: `Data`
133 :param r: a given gradient
134 :type r: `ArithmeticTuple`
135 :param args: pre-calculated values for ``m`` from `getArguments()`
136 :rtype: `Data`
137 :note: in the current implementation only the regularization term is
138 considered in the inverse Hessian approximation.
139
140 """
141 print "inverse Hessian approximation:"
142 print "Y = ",r[0]
143 print "X = ",r[1]
144 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
145 print "m = ",m
146 return m
147
148 def updateHessian(self):
149 """
150 notifies the class that the Hessian operator needs to be updated.
151 """
152 self.regularization.updateHessian()
153

  ViewVC Help
Powered by ViewVC 1.1.26