downunder/py_src/forwardmodels.py


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

__copyright__="""Copyright (c) 2003-2012 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
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__ = ['ForwardModel','GravityModel']

from esys.escript import unitsSI as U
from esys.escript import *
from esys.escript.linearPDEs import LinearSinglePDE

PI = 3.14159265358979323846
G = 6.6742e-11*U.m**3/(U.kg*U.sec**2)

class ForwardModel(object):
    """
    An abstract forward model that can be plugged into a cost function.
    Subclasses should implement getValue() and getGradient().
    """
    def __init__(self):
        pass

    def getArguments(self, x):
        return None

    def getValue(self, x, *args):
        pass

    def getGradient(self, x, *args):
        pass

class GravityModel(ForwardModel):
    """
    Forward Model for gravity inversion as described in the inversion cookbook.
    """
    def __init__(self, domain, chi, g, constrain_top=True, gravity_constant=G, tol=1e-8):
        """
        Creates a new gravity model on the given domain with one or more
        surveys (chi, g).
        """
        self.__domain = domain
        try:
            n=len(chi)
            m=len(g)
            if m != n:
                raise ValueError("Length of chi and g must be the same.")
            self.__chi = chi
            self.__g = g
        except TypeError:
            self.__chi = [chi]
            self.__g = [g]

        A=0
        for s in xrange(len(self.__chi)):
            A2 = integrate(inner(self.__chi[s], self.__g[s]**2))
            if A2 < 0:
                raise ValueError("Negative weighting factor for survey %s"%s)
            A=max(A2, A)
        if not A > 0:
            raise ValueError("No reference data set.")
        self.__G = gravity_constant
        BX = boundingBox(domain)
        DIM = domain.getDim()
        x = domain.getX()
        self.__pde=LinearSinglePDE(domain)
        self.__pde.getSolverOptions().setTolerance(tol)
        if constrain_top:
            constraint=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0])
        else:
            constraint=whereZero(x[DIM-1]-BX[DIM-1][0])
        for i in xrange(DIM-1):
            constraint=constraint+whereZero(x[i]-BX[i][1])+whereZero(x[i]-BX[i][0])
        self.__pde.setValue(A=kronecker(DIM), q=constraint)

    def getArguments(self, rho):
        """
        Returns precomputed values shared by getValue() and getGradient().
        """
        phi = self.getPotential(rho)
        return phi, -grad(phi)

    def getPotential(self, rho):
        self.__pde.setValue(Y=(-4*PI*self.__G) * rho, X=Data())
        phi=self.__pde.getSolution()
        return phi

    def getValue(self, rho, phi, gravity_force):
        A=0.
        for s in xrange(len(self.__chi)):
            A = inner(self.__chi[s], (gravity_force-self.__g[s])**2) + A
        return 0.5*integrate(A)

    def getGradient(self, rho, phi, gravity_force):
        Z=0.
        for s in xrange(len(self.__chi)):
            Z = self.__chi[s] * (-gravity_force+self.__g[s]) + Z

        self.__pde.setValue(Y=Data(), X=Z)
        ZT=self.__pde.getSolution()
        return ZT*(-4*PI*self.__G)

1
2	########################################################
3	#
4	# Copyright (c) 2003-2012 by University of Queensland
5	# Earth Systems Science Computational Center (ESSCC)
6	# http://www.uq.edu.au/esscc
7	#
8	# Primary Business: Queensland, Australia
9	# Licensed under the Open Software License version 3.0
10	# http://www.opensource.org/licenses/osl-3.0.php
11	#
12	########################################################
13
14	__copyright__="""Copyright (c) 2003-2012 by University of Queensland
15	Earth Systems Science Computational Center (ESSCC)
16	http://www.uq.edu.au/esscc
17	Primary Business: Queensland, Australia"""
18	__license__="""Licensed under the Open Software License version 3.0
19	http://www.opensource.org/licenses/osl-3.0.php"""
20	__url__="https://launchpad.net/escript-finley"
21
22	__all__ = ['ForwardModel','GravityModel']
23
24	from esys.escript import unitsSI as U
25	from esys.escript import *
26	from esys.escript.linearPDEs import LinearSinglePDE
27
28	PI = 3.14159265358979323846
29	G = 6.6742e-11U.m3/(U.kgU.sec**2)
30
31	class ForwardModel(object):
32	"""
33	An abstract forward model that can be plugged into a cost function.
34	Subclasses should implement getValue() and getGradient().
35	"""
36	def __init__(self):
37	pass
38
39	def getArguments(self, x):
40	return None
41
42	def getValue(self, x, *args):
43	pass
44
45	def getGradient(self, x, *args):
46	pass
47
48	class GravityModel(ForwardModel):
49	"""
50	Forward Model for gravity inversion as described in the inversion cookbook.
51	"""
52	def __init__(self, domain, chi, g, constrain_top=True, gravity_constant=G, tol=1e-8):
53	"""
54	Creates a new gravity model on the given domain with one or more
55	surveys (chi, g).
56	"""
57	self.__domain = domain
58	try:
59	n=len(chi)
60	m=len(g)
61	if m != n:
62	raise ValueError("Length of chi and g must be the same.")
63	self.__chi = chi
64	self.__g = g
65	except TypeError:
66	self.__chi = [chi]
67	self.__g = [g]
68
69	A=0
70	for s in xrange(len(self.__chi)):
71	A2 = integrate(inner(self.__chi[s], self.__g[s]**2))
72	if A2 < 0:
73	raise ValueError("Negative weighting factor for survey %s"%s)
74	A=max(A2, A)
75	if not A > 0:
76	raise ValueError("No reference data set.")
77	self.__G = gravity_constant
78	BX = boundingBox(domain)
79	DIM = domain.getDim()
80	x = domain.getX()
81	self.__pde=LinearSinglePDE(domain)
82	self.__pde.getSolverOptions().setTolerance(tol)
83	if constrain_top:
84	constraint=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0])
85	else:
86	constraint=whereZero(x[DIM-1]-BX[DIM-1][0])
87	for i in xrange(DIM-1):
88	constraint=constraint+whereZero(x[i]-BX[i][1])+whereZero(x[i]-BX[i][0])
89	self.__pde.setValue(A=kronecker(DIM), q=constraint)
90
91	def getArguments(self, rho):
92	"""
93	Returns precomputed values shared by getValue() and getGradient().
94	"""
95	phi = self.getPotential(rho)
96	return phi, -grad(phi)
97
98	def getPotential(self, rho):
99	self.__pde.setValue(Y=(-4PIself.__G) * rho, X=Data())
100	phi=self.__pde.getSolution()
101	return phi
102
103	def getValue(self, rho, phi, gravity_force):
104	A=0.
105	for s in xrange(len(self.__chi)):
106	A = inner(self.__chi[s], (gravity_force-self.__g[s])**2) + A
107	return 0.5*integrate(A)
108
109	def getGradient(self, rho, phi, gravity_force):
110	Z=0.
111	for s in xrange(len(self.__chi)):
112	Z = self.__chi[s] * (-gravity_force+self.__g[s]) + Z
113
114	self.__pde.setValue(Y=Data(), X=Z)
115	ZT=self.__pde.getSolution()
116	return ZT(-4PI*self.__G)
117