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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4373 - (hide annotations)
Mon Apr 22 03:16:24 2013 UTC (6 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 15991 byte(s)
more on geodetic coordinates
1 caltinay 3946
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4 jfenwick 4154 # Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3946
16 caltinay 4007 """Collection of forward models that define the inversion problem"""
17    
18 jfenwick 4154 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 jfenwick 3981 http://www.uq.edu.au
20 caltinay 3946 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 caltinay 4058 __all__ = ['ForwardModel','ForwardModelWithPotential','GravityModel','MagneticModel']
26 caltinay 3947
27 caltinay 3946 from esys.escript import unitsSI as U
28 caltinay 4095 from esys.escript import Data, Vector, Scalar, Function
29 gross 4037 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
30 gross 4373 from .coordinates import makeTranformation
31 caltinay 4007 from esys.escript.util import *
32 gross 4037 from math import pi as PI
33 gross 4102 import numpy as np
34 caltinay 3946
35 jfenwick 4238
36 caltinay 3946 class ForwardModel(object):
37     """
38     An abstract forward model that can be plugged into a cost function.
39 caltinay 4007 Subclasses need to implement `getValue()`, `getGradient()`, and possibly
40 gross 4373 `getArguments()` and 'getCoordinateTransformation'.
41 caltinay 3946 """
42     def __init__(self):
43     pass
44    
45     def getArguments(self, x):
46 caltinay 4044 return ()
47 caltinay 3946
48     def getValue(self, x, *args):
49 caltinay 4044 raise NotImplementedError
50 caltinay 3946
51     def getGradient(self, x, *args):
52 caltinay 4044 raise NotImplementedError
53 gross 4373
54     def getCoordinateTransformation(self):
55     return None
56    
57 gross 4076
58 caltinay 4044 class ForwardModelWithPotential(ForwardModel):
59 caltinay 3946 """
60 caltinay 4044 Base class for a forward model using a potential such as magnetic or
61 caltinay 4095 gravity. It defines a cost function::
62 caltinay 4044
63 gross 4099 defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) ) ** 2 )
64 caltinay 4044
65     where s runs over the survey, weight_i are weighting factors, data_i are
66     the data, and r_i are the results produced by the forward model.
67     It is assumed that the forward model is produced through postprocessing
68     of the solution of a potential PDE.
69 caltinay 3946 """
70 gross 4373 def __init__(self, domain, w, data, coordinates=None,
71     fixPotentialAtBottom=False,
72     tol=1e-8):
73 caltinay 3946 """
74 caltinay 4095 initializes a new forward model with potential.
75 caltinay 4044
76 gross 4037 :param domain: domain of the model
77 caltinay 4097 :type domain: `Domain`
78 gross 4093 :param w: data weighting factors
79 caltinay 4097 :type w: ``Vector`` or list of ``Vector``
80 gross 4037 :param data: data
81 caltinay 4097 :type data: ``Vector`` or list of ``Vector``
82 gross 4373 :param coordinates: defines coordinate system to be used
83     :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
84 gross 4037 :param tol: tolerance of underlying PDE
85     :type tol: positive ``float``
86 gross 4344 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
87     in addition to the top.
88     :type fixPotentialAtBottom: ``bool``
89 caltinay 3946 """
90 gross 4037 super(ForwardModelWithPotential, self).__init__()
91 caltinay 3946 self.__domain = domain
92 gross 4373 self.__trafo=makeTranformation(domain, coordinates)
93    
94 caltinay 3946 try:
95 gross 4093 n=len(w)
96 gross 4037 m=len(data)
97 gross 4099 if not m == n:
98 caltinay 4077 raise ValueError("Length of weight and data must be the same.")
99 gross 4093 self.__weight = w
100 gross 4037 self.__data = data
101 caltinay 3946 except TypeError:
102 gross 4093 self.__weight = [w]
103 gross 4037 self.__data = [data]
104 caltinay 3946
105     BX = boundingBox(domain)
106     DIM = domain.getDim()
107     x = domain.getX()
108     self.__pde=LinearSinglePDE(domain)
109     self.__pde.getSolverOptions().setTolerance(tol)
110 caltinay 4058 self.__pde.setSymmetryOn()
111 gross 4121 z=x[DIM-1]
112 gross 4344 q0=whereZero(z-BX[DIM-1][1])
113     if fixPotentialAtBottom: q0+=whereZero(z-BX[DIM-1][0])
114     self.__pde.setValue(q=q0)
115 caltinay 4213
116 gross 4102 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
117     self.diameter=1./sqrt(sum(1./self.edge_lengths**2))
118 gross 4373
119     if not self.__trafo.isCartesian():
120     fd=1./self.__trafo().getScalingFactors()
121     fw=self.__trafo().getScalingFactors()*sqrt(self.__trafo().getVolumeFactor())
122    
123     for s in range(len(self.__weight)):
124     self.__weight[s] = fw * self.__weight[s]
125     self.__data[s] = fd * self.__data[s]
126 caltinay 4213
127 gross 4373
128 gross 4102 def _rescaleWeights(self, scale=1., fetch_factor=1.):
129     """
130 caltinay 4213 rescales the weights such that
131    
132 gross 4102 *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
133     """
134     if not scale > 0:
135     raise ValueError("Value for scale must be positive.")
136     A=0
137 caltinay 4263 for s in range(len(self.__weight)):
138 gross 4102 A += integrate( abs(inner(self.__weight[s], self.__data[s])* inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
139     if A > 0:
140     A=sqrt(scale/A)/self.diameter
141 caltinay 4263 for s in range(len(self.__weight)): self.__weight[s]*=A
142 gross 4102 else:
143     raise ValueError("Rescaling of weights failed.")
144 caltinay 4213
145 gross 4037 def getDomain(self):
146     """
147 caltinay 4044 Returns the domain of the forward model.
148 caltinay 4095
149 caltinay 4097 :rtype: `Domain`
150 gross 4037 """
151     return self.__domain
152 gross 4373
153     def getCoordinateTransformation(self):
154     """
155     returns the coordinate transformation being used
156 caltinay 4044
157 gross 4373 :rtype: ``CoordinateTransformation``
158     """
159     return self.__trafo
160    
161 gross 4037 def getPDE(self):
162     """
163 caltinay 4044 Return the underlying PDE.
164    
165 gross 4037 :rtype: `LinearPDE`
166     """
167     return self.__pde
168 caltinay 4044
169 gross 4037 def getDefect(self, result):
170     """
171 caltinay 4044 Returns the defect value.
172    
173 gross 4037 :param result: a result vector
174     :type result: `Vector`
175     :rtype: ``float``
176     """
177     A=0.
178 caltinay 4263 for s in range(len(self.__weight)):
179 gross 4093 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
180 gross 4099 return A/2
181 caltinay 4044
182 gross 4093 def getDefectGradient(self, result):
183     Y=0.
184 caltinay 4094 for s in range(len(self.__weight)):
185     Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
186 gross 4099 return Y
187 caltinay 4097
188 caltinay 3960 def getSurvey(self, index=None):
189     """
190 caltinay 4063 Returns the pair (data_index, weight_index), where data_i is the data
191     of survey i, weight_i is the weighting factor for survey i.
192 caltinay 3960 If index is None, all surveys will be returned in a pair of lists.
193     """
194     if index is None:
195 gross 4037 return self.__data, self.__weight
196     if index>=len(self.__data):
197     raise IndexError("Forward model only has %d surveys"%len(self.__data))
198     return self.__data[index], self.__weight[index]
199 caltinay 4044
200    
201 gross 4074
202 gross 4037 class GravityModel(ForwardModelWithPotential):
203     """
204 caltinay 4044 Forward Model for gravity inversion as described in the inversion
205     cookbook.
206 gross 4037 """
207 caltinay 4097 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
208 gross 4373 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
209 gross 4037 """
210     Creates a new gravity model on the given domain with one or more
211 gross 4093 surveys (w, g).
212 caltinay 4044
213 gross 4037 :param domain: domain of the model
214 caltinay 4097 :type domain: `Domain`
215 gross 4093 :param w: data weighting factors
216 caltinay 4097 :type w: ``Vector`` or list of ``Vector``
217 caltinay 4044 :param g: gravity anomaly data
218 caltinay 4097 :type g: ``Vector`` or list of ``Vector``
219 gross 4373 :param coordinates: defines coordinate system to be used
220     :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
221 gross 4037 :param tol: tolerance of underlying PDE
222     :type tol: positive ``float``
223 gross 4344 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
224     in addition to the top.
225     :type fixPotentialAtBottom: ``bool``
226 gross 4093
227 gross 4102
228     :note: It is advisable to call rescaleWeights() to rescale weights before start inversion.
229 gross 4037 """
230 gross 4373 super(GravityModel, self).__init__(domain, w, g, coordinates, fixPotentialAtBottom, tol)
231 caltinay 4044
232 gross 4373 if not self.getCoordinateTransformation().isCartesian():
233     self.__G = 4*PI*gravity_constant * self.getCoordinateTransformation().getVolumeFactor()
234    
235     fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
236     A=self.getPDE().createCoeffiecient("A")
237     for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
238     self.getPDE().setValue(A=A)
239    
240     else:
241     self.__G = 4*PI*gravity_constant
242     self.getPDE().setValue(A=kronecker(self.getDomain()))
243 caltinay 4213
244 gross 4102 def rescaleWeights(self, scale=1., rho_scale=1.):
245     """
246 caltinay 4213 rescales the weights such that
247    
248 gross 4102 *sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale*
249 caltinay 4213
250 gross 4102 :param scale: scale of data weighting factors
251     :type scale: positive ``float``
252     :param rho_scale: scale of density.
253 caltinay 4213 :type rho_scale: ``Scalar``
254 gross 4102 """
255 gross 4373 self._rescaleWeights(scale, self.__G*rho_scale)
256 caltinay 4213
257 caltinay 3946 def getArguments(self, rho):
258     """
259 caltinay 4097 Returns precomputed values shared by `getValue()` and `getGradient()`.
260 caltinay 4044
261 gross 4037 :param rho: a suggestion for the density distribution
262 caltinay 4097 :type rho: ``Scalar``
263 gross 4037 :return: gravity potential and corresponding gravity field.
264     :rtype: ``Scalar``, ``Vector``
265 caltinay 3946 """
266     phi = self.getPotential(rho)
267 gross 4037 gravity_force = -grad(phi)
268     return phi, gravity_force
269 caltinay 3946
270     def getPotential(self, rho):
271 gross 4037 """
272 caltinay 4044 Calculates the gravity potential for a given density distribution.
273    
274 gross 4037 :param rho: a suggestion for the density distribution
275 caltinay 4097 :type rho: ``Scalar``
276 gross 4037 :return: gravity potential
277     :rtype: ``Scalar``
278     """
279     pde=self.getPDE()
280 caltinay 4044
281 gross 4037 pde.resetRightHandSideCoefficients()
282 gross 4373 pde.setValue(Y=-self.__G*rho)
283 gross 4037 phi=pde.getSolution()
284 caltinay 4044
285 caltinay 3946 return phi
286 caltinay 4044
287 caltinay 3946 def getValue(self, rho, phi, gravity_force):
288 gross 4037 """
289 caltinay 4044 Returns the value of the defect
290    
291 gross 4037 :param rho: density distribution
292 caltinay 4097 :type rho: ``Scalar``
293 gross 4037 :param phi: corresponding potential
294 caltinay 4097 :type phi: ``Scalar``
295 gross 4037 :param gravity_force: gravity force
296 caltinay 4097 :type gravity_force: ``Vector``
297 gross 4037 :rtype: ``float``
298     """
299     return self.getDefect(gravity_force)
300 caltinay 3946
301     def getGradient(self, rho, phi, gravity_force):
302 gross 4037 """
303 caltinay 4044 Returns the gradient of the defect with respect to density.
304    
305 gross 4037 :param rho: density distribution
306 caltinay 4097 :type rho: ``Scalar``
307 gross 4037 :param phi: corresponding potential
308 caltinay 4097 :type phi: ``Scalar``
309 gross 4037 :param gravity_force: gravity force
310 caltinay 4097 :type gravity_force: ``Vector``
311     :rtype: ``Scalar``
312 caltinay 4044 """
313 gross 4037 pde=self.getPDE()
314     pde.resetRightHandSideCoefficients()
315 gross 4093 pde.setValue(X=self.getDefectGradient(gravity_force))
316 gross 4037 ZT=pde.getSolution()
317 gross 4373 return ZT*(-self.__G)
318 caltinay 3946
319 caltinay 4044
320 gross 4037 class MagneticModel(ForwardModelWithPotential):
321     """
322 caltinay 4044 Forward Model for magnetic inversion as described in the inversion
323     cookbook.
324 gross 4037 """
325 gross 4373 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
326 gross 4037 """
327 caltinay 4044 Creates a new magnetic model on the given domain with one or more
328 gross 4093 surveys (w, B).
329 caltinay 4044
330 gross 4037 :param domain: domain of the model
331     :type domain: `Domain`
332 gross 4093 :param w: data weighting factors
333 caltinay 4097 :type w: ``Vector`` or list of ``Vector``
334 gross 4037 :param B: magnetic field data
335 caltinay 4097 :type B: ``Vector`` or list of ``Vector``
336 gross 4037 :param tol: tolerance of underlying PDE
337     :type tol: positive ``float``
338 gross 4373 :param coordinates: defines coordinate system to be used
339     :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
340 gross 4344 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
341     in addition to the top.
342     :type fixPotentialAtBottom: ``bool``
343 gross 4102 """
344 gross 4373 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
345     self.__background_magnetic_flux_density=interpolate(background_magnetic_flux_density,
346     B[0].getFunctionSpace())
347     if not self.getCoordinateTransformation().isCartesian():
348     self.__F = self.getCoordinateTransformation().getVolumeFactor()
349     self.__B_r=self.__background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
350     self.__B_b=self.__background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
351    
352     A=self.getPDE().createCoeffiecient("A")
353     fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
354     for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
355     self.getPDE().setValue(A=A)
356    
357     else:
358     self.getPDE().setValue(A=kronecker(self.getDomain()))
359     self.__B_r=self.__background_magnetic_flux_density
360     self.__B_b=self.__background_magnetic_flux_density
361 caltinay 4213
362 gross 4102 def rescaleWeights(self, scale=1., k_scale=1.):
363     """
364 caltinay 4213 rescales the weights such that
365    
366 gross 4115 *sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale*
367 caltinay 4213
368 gross 4099 :param scale: scale of data weighting factors
369 gross 4093 :type scale: positive ``float``
370 gross 4102 :param k_scale: scale of susceptibility.
371 caltinay 4213 :type k_scale: ``Scalar``
372 gross 4037 """
373 gross 4115 self._rescaleWeights(scale, inner(self.__background_magnetic_flux_density,1/self.edge_lengths ) * k_scale)
374 caltinay 4213
375 gross 4037 def getArguments(self, k):
376     """
377 caltinay 4097 Returns precomputed values shared by `getValue()` and `getGradient()`.
378 caltinay 4044
379     :param k: susceptibility
380 caltinay 4097 :type k: ``Scalar``
381 caltinay 4044 :return: scalar magnetic potential and corresponding magnetic field
382 gross 4037 :rtype: ``Scalar``, ``Vector``
383     """
384     phi = self.getPotential(k)
385 gross 4373 magnetic_flux_density = k * self.__B_b -grad(phi)
386 gross 4121 return phi, magnetic_flux_density
387 gross 4037
388     def getPotential(self, k):
389     """
390 caltinay 4044 Calculates the magnetic potential for a given susceptibility.
391    
392 gross 4037 :param k: susceptibility
393 caltinay 4097 :type k: ``Scalar``
394 caltinay 4044 :return: magnetic potential
395 gross 4037 :rtype: ``Scalar``
396     """
397     pde=self.getPDE()
398 caltinay 4044
399 gross 4037 pde.resetRightHandSideCoefficients()
400 gross 4373 pde.setValue(X = k* self.__B_r)
401 gross 4037 phi=pde.getSolution()
402 caltinay 4044
403 gross 4037 return phi
404 caltinay 4044
405 gross 4121 def getValue(self, k, phi, magnetic_flux_density):
406 gross 4037 """
407 caltinay 4044 Returns the value of the defect.
408    
409 gross 4037 :param k: susceptibility
410 caltinay 4097 :type k: ``Scalar``
411 gross 4037 :param phi: corresponding potential
412 caltinay 4097 :type phi: ``Scalar``
413 gross 4121 :param magnetic_flux_density: magnetic field
414     :type magnetic_flux_density: ``Vector``
415 gross 4037 :rtype: ``float``
416     """
417 gross 4121 return self.getDefect(magnetic_flux_density)
418 gross 4037
419 gross 4121 def getGradient(self, k, phi, magnetic_flux_density):
420 gross 4037 """
421 caltinay 4044 Returns the gradient of the defect with respect to susceptibility.
422    
423 gross 4037 :param k: susceptibility
424 caltinay 4097 :type k: ``Scalar``
425 gross 4037 :param phi: corresponding potential
426 caltinay 4097 :type phi: ``Scalar``
427 gross 4121 :param magnetic_flux_density: magnetic field
428     :type magnetic_flux_density: ``Vector``
429 caltinay 4097 :rtype: ``Scalar``
430 caltinay 4044 """
431 gross 4121 Y=self.getDefectGradient(magnetic_flux_density)
432 gross 4037 pde=self.getPDE()
433     pde.resetRightHandSideCoefficients()
434     pde.setValue(X=Y)
435     YT=pde.getSolution()
436 gross 4373 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
437 gross 4037

  ViewVC Help
Powered by ViewVC 1.1.26