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

  ViewVC Help
Powered by ViewVC 1.1.26