/[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 4263 - (show annotations)
Thu Feb 28 05:59:04 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 13319 byte(s)
No more xrange in downunder.

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

  ViewVC Help
Powered by ViewVC 1.1.26