/[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 4079 - (show annotations)
Fri Nov 16 07:59:01 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 5641 byte(s)
some modifications to scaling in downunder. still not perfect.
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 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
102 return self.mu_model * self.forwardmodel.getValue(args[0],*args[1]) \
103 + self.regularization.getValue(m, *args[2])
104
105 def _getGradient(self, m, *args):
106 """
107 returns the gradient of *f* at *m*.
108 If the precalculated values are not supplied `getArguments()` is called.
109
110 :rtype: `esys.escript.Data`
111 """
112 dpdm = self.mapping.getDerivative(m)
113 if len(args)==0:
114 args = self.getArguments(m)
115
116 Y = self.forwardmodel.getGradient(args[0],*args[1]) * dpdm
117 g_reg = self.regularization.getGradient(m, *args[2])
118 print "grad forward = ", Y
119 print "grad regularization Y = ", g_reg[0]
120 print "grad regularization X = ", g_reg[1]
121
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 considered in the
138 inverse Hessian approximation.
139 """
140 print "nverseHessianApproximation:"
141 print "Y = ",r[0]
142 print "X = ",r[1]
143 m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
144 print "m = ",m
145 return m
146
147 def updateHessian(self):
148 """
149 notifies the class that the Hessian operator needs to be updated.
150 """
151 self.regularization.updateHessian()

  ViewVC Help
Powered by ViewVC 1.1.26