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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4477 - (show annotations)
Thu Jun 20 02:41:28 2013 UTC (6 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 16567 byte(s)
escript: use scoped_array instead of scoped_ptr
ripley: netcdf returns a copy of an array so use scoped_array
downunder: s/Sperical/Spherical/
scripts: added a few more suppressions

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

  ViewVC Help
Powered by ViewVC 1.1.26