downunder/py_src/inversioncostfunctions.py


##############################################################################
#
# Copyright (c) 2003-2012 by University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development since 2012 by School of Earth Sciences
#
##############################################################################

"""Collection of cost functions for inversion"""

__copyright__="""Copyright (c) 2003-2012 by University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

__all__ = ['SimpleInversionCostFunction']

from costfunctions import MeteredCostFunction
from esys.escript.pdetools import ArithmeticTuple
from esys.escript import Data


class SimpleInversionCostFunction(MeteredCostFunction):
    """
    This is a simple cost function with a single continuous (mapped) variable.
    It is the sum of two weighted terms, a single forward model and a single
    regularization term. This cost function is used in the gravity inversion.
    """
    def __init__(self, regularization, mapping, forwardmodel):
        """
        constructor stores the supplied object references and sets default
        weights.

        :param regularization: The regularization part of the cost function
        :param mapping: Parametrization object
        :param forwardmodel: The forward model part of the cost function
        """
        super(SimpleInversionCostFunction, self).__init__()
        self.forwardmodel=forwardmodel
        self.regularization=regularization
        self.mapping=mapping
        self.setWeights()

    def setWeights(self, mu_model=1., mu_reg_0=None,mu_reg_1=None):
        """
        sets the weighting factors for the forward model and regularization
        terms.
        
        :param mu_model: Weighting factor for the forward model (default=1.)
        :type mu_model: non-negative `float`
        :param mu_reg_0: Weighting factor for the regularization (default=1.)
        :type mu_reg_0: non-negative `float`
        """
        if mu_model<0:
            raise ValueError("weighting factors must be non-negative.")
        self.mu_model=mu_model
        self.regularization.setWeightsForS0(mu_reg_0)
        self.regularization.setWeightsForS1(mu_reg_1)
        
    def _getDualProduct(self, x, r):
        """
        returns ``regularization.getDualProduct(x, r)``

        :rtype: `float`
        """
        return self.regularization.getDualProduct(x, r)

    def _getArguments(self, m):
        """
        returns precalculated values that are shared in the calculation of
        *f(x)* and *grad f(x)*. In this implementation returns a tuple with the
        mapped value of ``m``, the arguments from the forward model and the
        arguments from the regularization.

        :rtype: `tuple`
        """
        p=self.mapping(m)
        return p, self.forwardmodel.getArguments(p), self.regularization.getArguments(m)

    def _getValue(self, m, *args):
        """
        returns the function value at m.
        If the precalculated values are not supplied `getArguments()` is called.

        :rtype: `float`
        """
        # if there is more than one forward_model and/or regularization their
        # contributions need to be added up. But this implementation allows
        # only one of each...
        if len(args)==0:
            args=self.getArguments(m)

        return  self.mu_model * self.forwardmodel.getValue(args[0],*args[1]) \
               +  self.regularization.getValue(m, *args[2])

    def _getGradient(self, m, *args):
        """
        returns the gradient of *f* at *m*.
        If the precalculated values are not supplied `getArguments()` is called.

        :rtype: `esys.escript.Data`
        """
        dpdm = self.mapping.getDerivative(m)
        if len(args)==0:
            args = self.getArguments(m)

        Y = self.forwardmodel.getGradient(args[0],*args[1]) * dpdm
        g_reg = self.regularization.getGradient(m, *args[2])
        print "grad forward = ", Y
        print "grad regularization Y  = ", g_reg[0]
        print "grad regularization X = ", g_reg[1]
        
         
        return self.mu_model * ArithmeticTuple(Y, Data()) + g_reg


    def _getInverseHessianApproximation(self, m, r, *args):
        """
        returns an approximative evaluation *p* of the inverse of the Hessian operator of the cost function
        for a given gradient type *r* at a given location *m*: *H(m) p = r*
        
        :param m: level set approximation where to calculate Hessian inverse
        :type m: ``Data``
        :param r: a given gradient 
        :type r: ``ArithmeticTuple``
        :param args: pre-calculated values for ``m`` from ``getArguments()``
        :rtype: ``Data``
        :note: in the current implementation only the regularization term is considered in the 
        inverse Hessian approximation. 
        """
        print "nverseHessianApproximation:"
        print "Y  = ",r[0]
        print "X  = ",r[1]
        m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
        print "m  = ",m
        return m
        
    def updateHessian(self):
        """
        notifies the class that the Hessian operator needs to be updated.
        """
        self.regularization.updateHessian()
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()