/[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 4100 - (hide annotations)
Tue Dec 11 06:55:20 2012 UTC (6 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 5812 byte(s)
minimizer uses now a relative change to solution to terminate iteration. Gravity still has problems
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 gross 4099 self.setTradeOffFactors()
54 gross 4074
55 gross 4099 def setTradeOffFactors(self, mu_model=1., mu_reg=1.):
56 gross 4074 """
57 gross 4099 sets the trade-off factors for the forward model and regularization
58 gross 4074 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 gross 4099 :param mu_reg: Weighting factor for the regularization (default=1.)
63     :type mu_reg: non-negative `float`
64 gross 4074 """
65     if mu_model<0:
66 gross 4099 raise ValueError("trade-off factors must be non-negative.")
67 gross 4074 self.mu_model=mu_model
68 gross 4099 self.regularization.setTradeOffFactors(mu_reg)
69    
70 gross 4074 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 caltinay 4097 returns precomputed values that are shared in the calculation of
81 gross 4074 *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 caltinay 4097 If the precomputed values are not supplied `getArguments()` is called.
94 gross 4074
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 gross 4079
103 gross 4074 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 caltinay 4097 If the precomputed values are not supplied `getArguments()` is called.
110 gross 4074
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 gross 4079 print "grad forward = ", Y
120     print "grad regularization Y = ", g_reg[0]
121     print "grad regularization X = ", g_reg[1]
122 caltinay 4097
123 gross 4074 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 caltinay 4097
131 gross 4074 :param m: level set approximation where to calculate Hessian inverse
132 caltinay 4097 :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 gross 4074 """
141 caltinay 4095 print "inverse Hessian approximation:"
142 gross 4079 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 caltinay 4097
148 gross 4074 def updateHessian(self):
149     """
150     notifies the class that the Hessian operator needs to be updated.
151     """
152 caltinay 4095 self.regularization.updateHessian()
153 gross 4100
154     def _getNorm(self, m):
155     """
156     returns the norm of ``m``
157 caltinay 4095
158 gross 4100 :param m: level set function
159     :type m: `Data`
160     :rtype: ``float``
161     """
162     return self.regularization.getNorm(m)

  ViewVC Help
Powered by ViewVC 1.1.26