/[escript]/trunk/downunder/py_src/forwardmodels/base.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/forwardmodels/base.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (hide annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 7263 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 caltinay 5548 ##############################################################################
2     #
3 jfenwick 5593 # Copyright (c) 2003-2015 by The University of Queensland
4 caltinay 5548 # http://www.uq.edu.au
5     #
6     # Primary Business: Queensland, Australia
7     # Licensed under the Open Software License version 3.0
8     # http://www.opensource.org/licenses/osl-3.0.php
9     #
10     # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11     # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13     #
14     ##############################################################################
15    
16     """Base classes for forward models"""
17    
18 sshaw 5707 from __future__ import division, print_function
19    
20 jfenwick 5593 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
21 caltinay 5548 http://www.uq.edu.au
22     Primary Business: Queensland, Australia"""
23     __license__="""Licensed under the Open Software License version 3.0
24     http://www.opensource.org/licenses/osl-3.0.php"""
25     __url__="https://launchpad.net/escript-finley"
26    
27     __all__ = ['ForwardModel','ForwardModelWithPotential']
28    
29     from esys.downunder.coordinates import makeTranformation
30     from esys.escript.linearPDEs import LinearSinglePDE
31     from esys.escript.util import *
32     import numpy as np
33    
34     class ForwardModel(object):
35     """
36     An abstract forward model that can be plugged into a cost function.
37     Subclasses need to implement `getDefect()`, `getGradient()`, and possibly
38     `getArguments()` and 'getCoordinateTransformation'.
39     """
40     def __init__(self):
41     pass
42    
43     def getArguments(self, x):
44     return ()
45    
46     def getCoordinateTransformation(self):
47     return None
48    
49     def getDefect(self, x, *args):
50     raise NotImplementedError
51    
52     def getGradient(self, x, *args):
53     raise NotImplementedError
54    
55    
56     class ForwardModelWithPotential(ForwardModel):
57     """
58     Base class for a forward model using a potential such as magnetic or
59     gravity. It defines a cost function:
60    
61     defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) )**2 )
62    
63     where s runs over the survey, weight_i are weighting factors, data_i are
64     the data, and r_i are the results produced by the forward model.
65     It is assumed that the forward model is produced through postprocessing
66     of the solution of a potential PDE.
67     """
68     def __init__(self, domain, w, data, coordinates=None,
69     fixPotentialAtBottom=False,
70     tol=1e-8):
71     """
72     initializes a new forward model with potential.
73    
74     :param domain: domain of the model
75     :type domain: `Domain`
76     :param w: data weighting factors
77     :type w: ``Vector`` or list of ``Vector``
78     :param data: data
79     :type data: ``Vector`` or list of ``Vector``
80     :param coordinates: defines coordinate system to be used
81     :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
82     :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
83     in addition to the top.
84     :type fixPotentialAtBottom: ``bool``
85     :param tol: tolerance of underlying PDE
86     :type tol: positive ``float``
87     """
88     super(ForwardModelWithPotential, self).__init__()
89     self.__domain = domain
90     self.__trafo = makeTranformation(domain, coordinates)
91    
92     try:
93     n=len(w)
94     m=len(data)
95     if not m == n:
96     raise ValueError("Length of weight and data must be the same.")
97     self.__weight = w
98     self.__data = data
99     except TypeError:
100     self.__weight = [w]
101     self.__data = [data]
102    
103     BX = boundingBox(domain)
104     DIM = domain.getDim()
105     x = domain.getX()
106     self.__pde=LinearSinglePDE(domain)
107     self.__pde.getSolverOptions().setTolerance(tol)
108     self.__pde.setSymmetryOn()
109     z=x[DIM-1]
110     q0=whereZero(z-BX[DIM-1][1])
111     if fixPotentialAtBottom: q0+=whereZero(z-BX[DIM-1][0])
112     self.__pde.setValue(q=q0)
113    
114     self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
115     self.diameter=1./sqrt(sum(1./self.edge_lengths**2))
116    
117     self.__origweight=[]
118     for s in range(len(self.__weight)):
119     # save a copy of the original weights in case of rescaling
120     self.__origweight.append(1.*self.__weight[s])
121    
122     if not self.__trafo.isCartesian():
123     fd=1./self.__trafo.getScalingFactors()
124     fw=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
125     for s in range(len(self.__weight)):
126     self.__weight[s] = fw * self.__weight[s]
127     self.__data[s] = fd * self.__data[s]
128    
129     def _rescaleWeights(self, scale=1., fetch_factor=1.):
130     """
131     rescales the weights such that
132    
133     *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
134     """
135     if not scale > 0:
136     raise ValueError("Value for scale must be positive.")
137     A=0
138     # copy back original weights before rescaling
139     self.__weight=[1.*ow for ow in self.__origweight]
140    
141     for s in range(len(self.__weight)):
142     A += integrate(abs(inner(self.__weight[s], self.__data[s]) * inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
143     if A > 0:
144     A=sqrt(scale/A)/self.diameter
145     if not self.__trafo.isCartesian():
146     A*=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
147     for s in range(len(self.__weight)):
148     self.__weight[s]*=A
149     else:
150     raise ValueError("Rescaling of weights failed.")
151    
152     def getDomain(self):
153     """
154     Returns the domain of the forward model.
155    
156     :rtype: `Domain`
157     """
158     return self.__domain
159    
160     def getCoordinateTransformation(self):
161     """
162     returns the coordinate transformation being used
163    
164     :rtype: ``CoordinateTransformation``
165     """
166     return self.__trafo
167    
168     def getPDE(self):
169     """
170     Return the underlying PDE.
171    
172     :rtype: `LinearPDE`
173     """
174     return self.__pde
175    
176     def _getDefect(self, result):
177     """
178     Returns the defect value.
179    
180     :param result: a result vector
181     :type result: `Vector`
182     :rtype: ``float``
183     """
184     A=0.
185     for s in range(len(self.__weight)):
186     A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
187     return A/2
188    
189     def getDefectGradient(self, result):
190     Y=0.
191     for s in range(len(self.__weight)):
192     Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
193     return Y
194    
195     def getSurvey(self, index=None):
196     """
197     Returns the pair (data_index, weight_index), where data_i is the data
198     of survey i, weight_i is the weighting factor for survey i.
199     If index is None, all surveys will be returned in a pair of lists.
200     """
201     if index is None:
202     return self.__data, self.__weight
203     if index>=len(self.__data):
204     raise IndexError("Forward model only has %d surveys"%len(self.__data))
205     return self.__data[index], self.__weight[index]
206    

  ViewVC Help
Powered by ViewVC 1.1.26