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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4097 - (hide annotations)
Fri Dec 7 01:18:35 2012 UTC (6 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 5636 byte(s)
Minor edits.

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

  ViewVC Help
Powered by ViewVC 1.1.26