/[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 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     # 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     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()

  ViewVC Help
Powered by ViewVC 1.1.26