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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4267 - (show annotations)
Fri Mar 1 04:50:54 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 22804 byte(s)
Updated inversion example scripts, added ER Mapper data, checked results
and corrected some issues. Cookbook does not reflect the changes yet.

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 """Higher-level classes that allow running inversions with minimal set up"""
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__ = ['InversionDriver', 'GravityInversion','MagneticInversion', 'JointGravityMagneticInversion']
26
27 import logging
28 import numpy as np
29
30 from esys.escript import *
31 from esys.escript import unitsSI as U
32 from esys.weipa import createDataset
33
34 from .inversioncostfunctions import InversionCostFunction
35 from .forwardmodels import GravityModel, MagneticModel
36 from .mappings import *
37 from .minimizers import *
38 from .regularizations import Regularization
39 from .datasources import DataSource
40
41 class InversionDriver(object):
42 """
43 Base class for running an inversion
44 """
45 def __init__(self, solverclass=None, debug=False):
46 """
47 creates an instance of an inversion driver.
48
49 :param solverclass: class of the solver to be used.
50 :type solverclass: 'AbstractMinimizer'.
51 """
52 # This configures basic logging to get some meaningful output
53 if debug:
54 level=logging.DEBUG
55 else:
56 level=logging.INFO
57 logging.getLogger('inv').setLevel(level)
58 self.logger=logging.getLogger('inv.%s'%self.__class__.__name__)
59 if solverclass is None: solverclass=MinimizerLBFGS
60 self.__solver=solverclass()
61 self.initial_value = None
62 self.m=None
63
64
65 def setCostFunction(self, costfunction):
66 """
67 sets the cost function of the inversion. This function needs to be called
68 before the inversion iteration can be started.
69
70 :param costfunction: domain of the inversion
71 :type costfunction: 'InversionCostFunction'
72 """
73 self.__solver.setCostFunction(costfunction)
74
75 def getCostFunction(self):
76 """
77 returns the domain of the inversion
78
79 :rtype: 'InversionCostFunction'
80 """
81 if self.isSetUp():
82 return self.getSolver().getCostFunction()
83 else:
84 raise RuntimeError("Inversion is not set up.")
85
86 def getSolver(self):
87 """
88 The solver to be used in the inversion process. See the minimizers
89 module for available solvers. By default, the L-BFGS minimizer is used.
90
91 :rtype: 'AbstractMinimizer'.
92 """
93 return self.__solver
94
95 def isSetUp(self):
96 """
97 returns True if the inversion is set up and is ready to run.
98
99 :rtype: `bool`
100 """
101 if self.getSolver().getCostFunction():
102 return True
103 else:
104 return False
105
106 def getDomain(self):
107 """
108 returns the domain of the inversion
109
110 :rtype: `Domain`
111 """
112 self.getCostFunction().getDomain()
113
114 def setSolverCallback(self, callback=None):
115 """
116 Sets the callback function which is called after every solver
117 iteration.
118 """
119 self.getSolver().setCallback(callback)
120
121 def setSolverMaxIterations(self, maxiter=None):
122 """
123 Sets the maximum number of solver iterations to run.
124 If `maxiter` is reached the iteration is terminated and
125 `MinimizerMaxIterReached` is thrown.
126
127 :param maxiter: maximum number of iteration steps.
128 :type maxiter: positive ``int``
129 """
130 if maxiter == None: maxiter = 200
131 if maxiter>0:
132 self.getSolver().setMaxIterations(maxiter)
133 else:
134 raise ValueError("maxiter must be positive.")
135
136
137 def setSolverTolerance(self, m_tol=None, J_tol=None):
138 """
139 Sets the error tolerance for the solver. An acceptable solution is
140 considered to be found once the tolerance is reached.
141
142 :param m_tol: tolerance for changes to level set function. If `None` changes to the
143 level set function are not checked for convergence during iteration.
144 :param J_tol: tolerance for changes to cost function. If `None` changes to the
145 cost function are not checked for convergence during iteration.
146 :type m_tol: `float` or `None`
147 :type J_tol: `float` or `None`
148 :note: if both arguments are equal to `None` the default setting m_tol=1e-4, J_tol=None is used.
149
150 """
151 if m_tol == None and J_tol==None:
152 m_tol=1e-4
153 self.getSolver().setTolerance(m_tol=m_tol, J_tol=J_tol)
154
155 def setInitialGuess(self, *args):
156 """
157 Sets the initial guess for the inversion iteration.
158 By default zero is used.
159 """
160 self.initial_value=self.getCostFunction().createLevelSetFunction(*args)
161
162 def run(self):
163 """
164 This function runs the inversion.
165
166 :return: physical parameters as result of the inversion
167 :rtype: ``list`` of physical parameters or a physical parameter
168 """
169 if not self.isSetUp():
170 raise RuntimeError("Inversion is not setup.")
171
172 if self.initial_value == None: self.setInitialGuess()
173 self.logger.info("Starting solver, please stand by...")
174 try:
175 self.getSolver().run(self.initial_value)
176 self.m=self.getSolver().getResult()
177 self.p=self.getCostFunction().getProperties(self.m)
178 except MinimizerException as e:
179 self.m=self.getSolver().getResult()
180 self.p=self.getCostFunction().getProperties(self.m)
181 self.logger.info("Solver iteration failed!")
182 raise e
183 self.logger.info("Solver completed successfully!")
184 self.getSolver().logSummary()
185 return self.p
186
187 def getLevelSetFunction(self):
188 """
189 returns the level set function as solution of the optimization problem
190
191 :rtype: `Data`
192 """
193 return self.m
194
195 def setup(self, *args, **k_args):
196 """
197 returns True if the inversion is set up and ready to run.
198 """
199 pass
200
201 class GravityInversion(InversionDriver):
202 """
203 Driver class to perform an inversion of Gravity (Bouguer) anomaly data.
204 The class uses the standard `Regularization` class for a single level set
205 function, `DensityMapping` mapping, and the gravity forward model
206 `GravityModel`.
207 """
208 def setup(self, domainbuilder,
209 rho0=None, drho=None, z0=None, beta=None,
210 w0=None, w1=None):
211 """
212 Sets up the inversion parameters from a `DomainBuilder`.
213
214 :param domainbuilder: Domain builder object with gravity source(s)
215 :type domainbuilder: `DomainBuilder`
216 :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
217 :type rho0: ``float`` or ``Scalar``
218 :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
219 :type drho: ``float`` or ``Scalar``
220 :param z0: reference depth for depth weighting, see `DensityMapping`. If not specified, zero is used.
221 :type z0: ``float`` or ``Scalar``
222 :param beta: exponent for depth weighting, see `DensityMapping`. If not specified, zero is used.
223 :type beta: ``float`` or ``Scalar``
224 :param w0: weighting factor for level set term regularization. If not set zero is assumed.
225 :type w0: ``Scalar`` or ``float``
226 :param w1: weighting factor for the gradient term in the regularization. If not set zero is assumed
227 :type w1: ``Vector`` or list of ``float``
228 """
229 self.logger.info('Retrieving domain...')
230 dom=domainbuilder.getDomain()
231 DIM=dom.getDim()
232 #========================
233 self.logger.info('Creating mapping...')
234 rho_mapping=DensityMapping(dom, rho0=rho0, drho=drho, z0=z0, beta=beta)
235 scale_mapping=rho_mapping.getTypicalDerivative()
236 #========================
237 self.logger.info("Setting up regularization...")
238 if w1 is None:
239 w1=[1.]*DIM
240 rho_mask = domainbuilder.getSetDensityMask()
241 regularization=Regularization(dom, numLevelSets=1,\
242 w0=w0, w1=w1, location_of_set_m=rho_mask)
243 #====================================================================
244 self.logger.info("Retrieving gravity surveys...")
245 surveys=domainbuilder.getGravitySurveys()
246 g=[]
247 w=[]
248 for g_i,sigma_i in surveys:
249 w_i=safeDiv(1., sigma_i)
250 if g_i.getRank()==0:
251 g_i=g_i*kronecker(DIM)[DIM-1]
252 if w_i.getRank()==0:
253 w_i=w_i*kronecker(DIM)[DIM-1]
254 g.append(g_i)
255 w.append(w_i)
256 self.logger.debug("Added gravity survey:")
257 self.logger.debug("g = %s"%g_i)
258 self.logger.debug("sigma = %s"%sigma_i)
259 self.logger.debug("w = %s"%w_i)
260 #====================================================================
261
262 self.logger.info("Setting up model...")
263 forward_model=GravityModel(dom, w, g)
264 forward_model.rescaleWeights(rho_scale=scale_mapping)
265
266 #====================================================================
267 self.logger.info("Setting cost function...")
268 self.setCostFunction(InversionCostFunction(regularization, rho_mapping, forward_model))
269
270 def setInitialGuess(self, rho=None):
271 """
272 set the initial guess *rho* for density the inversion iteration. If no *rho* present
273 then an appropriate initial guess is chosen.
274
275 :param rho: initial value for the density anomaly.
276 :type rho: `Scalar`
277 """
278 if rho:
279 super(GravityInversion,self).setInitialGuess(rho)
280 else:
281 super(GravityInversion,self).setInitialGuess()
282
283 def siloWriterCallback(self, k, m, Jm, g_Jm):
284 """
285 callback function that can be used to track the solution
286
287 :param k: iteration count
288 :param m: current m approximation
289 :param Jm: value of cost function
290 :param g_Jm: gradient of f at x
291 """
292 fn='inv.%d'%k
293 ds=createDataset(density=self.getCostFunction().getProperties(m))
294 ds.setCycleAndTime(k,k)
295 ds.saveSilo(fn)
296 self.logger.debug("J(m) = %e"%Jm)
297
298 class MagneticInversion(InversionDriver):
299 """
300 Driver class to perform an inversion of magnetic anomaly data. The class
301 uses the standard `Regularization` class for a single level set function,
302 `SusceptibilityMapping` mapping and the linear magnetic forward model
303 `MagneticModel`.
304 """
305 def setup(self, domainbuilder,
306 k0=None, dk=None, z0=None, beta=None,
307 w0=None, w1=None):
308 """
309 Sets up the inversion from a `DomainBuilder`.
310 If magnetic data are given as scalar it is assumed that values are
311 collected in direction of the background magnetic field.
312
313 :param domainbuilder: Domain builder object with gravity source(s)
314 :type domainbuilder: `DomainBuilder`
315 :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
316 :type k0: ``float`` or ``Scalar``
317 :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
318 :type dk: ``float`` or ``Scalar``
319 :param z0: reference depth for depth weighting, see `SusceptibilityMapping`. If not specified, zero is used.
320 :type z0: ``float`` or ``Scalar``
321 :param beta: exponent for depth weighting, see `SusceptibilityMapping`. If not specified, zero is used.
322 :type beta: ``float`` or ``Scalar``
323 :param w0: weighting factor for level set term regularization. If not set zero is assumed.
324 :type w0: ``Scalar`` or ``float``
325 :param w1: weighting factor for the gradient term in the regularization. If not set zero is assumed
326 :type w1: ``Vector`` or list of ``float``
327 """
328 self.logger.info('Retrieving domain...')
329 dom=domainbuilder.getDomain()
330 DIM=dom.getDim()
331
332 #========================
333 self.logger.info('Creating mapping ...')
334 k_mapping=SusceptibilityMapping(dom, k0=k0, dk=dk, z0=z0, beta=beta)
335 scale_mapping=k_mapping.getTypicalDerivative()
336 #========================
337 self.logger.info("Setting up regularization...")
338 if w1 is None:
339 w1=[1.]*DIM
340 k_mask = domainbuilder.getSetSusceptibilityMask()
341 regularization=Regularization(dom, numLevelSets=1,w0=w0, w1=w1, location_of_set_m=k_mask)
342
343 #====================================================================
344 self.logger.info("Retrieving magnetic field surveys...")
345 d_b=normalize(domainbuilder.getBackgroundMagneticFluxDensity())
346 surveys=domainbuilder.getMagneticSurveys()
347 B=[]
348 w=[]
349 for B_i,sigma_i in surveys:
350 w_i=safeDiv(1., sigma_i)
351 if B_i.getRank()==0:
352 B_i=B_i*d_b
353 if w_i.getRank()==0:
354 w_i=w_i*d_b
355 B.append(B_i)
356 w.append(w_i)
357 self.logger.debug("Added magnetic survey:")
358 self.logger.debug("B = %s"%B_i)
359 self.logger.debug("sigma = %s"%sigma_i)
360 self.logger.debug("w = %s"%w_i)
361 #====================================================================
362 self.logger.info("Setting up model...")
363 forward_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity())
364 forward_model.rescaleWeights(k_scale=scale_mapping)
365
366 #====================================================================
367 self.logger.info("Setting cost function...")
368 self.setCostFunction(InversionCostFunction(regularization, k_mapping, forward_model))
369
370 def setInitialGuess(self, k=None):
371 """
372 set the initial guess *k* for susceptibility for the inversion iteration. If no *k* present
373 then an appropriate initial guess is chosen.
374
375 :param k: initial value for the susceptibility anomaly.
376 :type k: `Scalar`
377 """
378 if k:
379 super(MagneticInversion,self).setInitialGuess(k)
380 else:
381 super(MagneticInversion,self).setInitialGuess()
382
383 def siloWriterCallback(self, k, m, Jm, g_Jm):
384 """
385 callback function that can be used to track the solution
386
387 :param k: iteration count
388 :param m: current m approximation
389 :param Jm: value of cost function
390 :param g_Jm: gradient of f at x
391 """
392 fn='inv.%d'%k
393 ds=createDataset(susceptibility=self.getCostFunction().getProperties(m))
394 ds.setCycleAndTime(k,k)
395 ds.saveSilo(fn)
396 self.logger.debug("J(m) = %e"%Jm)
397
398 class JointGravityMagneticInversion(InversionDriver):
399 """
400 Driver class to perform a joint inversion of Gravity (Bouguer) and
401 magnetic anomaly data. The class uses the standard `Regularization` class
402 for two level set functions with cross-gradient correlation,
403 `DensityMapping` and `SusceptibilityMapping` mappings, the gravity forward
404 model `GravityModel` and the linear magnetic forward model `MagneticModel`.
405 """
406 DENSITY=0
407 SUSCEPTIBILITY=1
408
409 def setup(self, domainbuilder,
410 rho0=None, drho=None, rho_z0=None, rho_beta=None,
411 k0=None, dk=None, k_z0=None, k_beta=None, w0=None, w1=None,
412 w_gc=None):
413 """
414 Sets up the inversion from an instance ``domainbuilder`` of a
415 `DomainBuilder`. Gravity and magnetic data attached to the
416 ``domainbuilder`` are considered in the inversion.
417 If magnetic data are given as scalar it is assumed that values are
418 collected in direction of the background magnetic field.
419
420 :param domainbuilder: Domain builder object with gravity source(s)
421 :type domainbuilder: `DomainBuilder`
422 :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
423 :type rho0: ``float`` or `Scalar`
424 :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
425 :type drho: ``float`` or `Scalar`
426 :param rho_z0: reference depth for depth weighting for density, see `DensityMapping`. If not specified, zero is used.
427 :type rho_z0: ``float`` or `Scalar`
428 :param rho_beta: exponent for depth weighting for density, see `DensityMapping`. If not specified, zero is used.
429 :type rho_beta: ``float`` or `Scalar`
430 :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
431 :type k0: ``float`` or `Scalar`
432 :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
433 :type dk: ``float`` or `Scalar`
434 :param k_z0: reference depth for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
435 :type k_z0: ``float`` or `Scalar`
436 :param k_beta: exponent for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
437 :type k_beta: ``float`` or `Scalar`
438 :param w0: weighting factors for level set term regularization, see `Regularization`. If not set zero is assumed.
439 :type w0: `Data` or ``ndarray`` of shape (2,)
440 :param w1: weighting factor for the gradient term in the regularization see `Regularization`. If not set zero is assumed
441 :type w1: `Data` or ``ndarray`` of shape (2,DIM)
442 :param w_gc: weighting factor for the cross gradient term in the regularization, see `Regularization`. If not set one is assumed
443 :type w_gc: `Scalar` or `float`
444 """
445 self.logger.info('Retrieving domain...')
446 dom=domainbuilder.getDomain()
447 DIM=dom.getDim()
448
449 #========================
450 self.logger.info('Creating mappings ...')
451 rho_mapping=DensityMapping(dom, rho0=rho0, drho=drho, z0=rho_z0, beta=rho_beta)
452 rho_scale_mapping=rho_mapping.getTypicalDerivative()
453 self.logger.debug("rho_scale_mapping = %s"%rho_scale_mapping)
454 k_mapping=SusceptibilityMapping(dom, k0=k0, dk=dk, z0=k_z0, beta=k_beta)
455 k_scale_mapping=k_mapping.getTypicalDerivative()
456 self.logger.debug("k_scale_mapping = %s"%k_scale_mapping)
457 #========================
458 self.logger.info("Setting up regularization...")
459
460 if w1 is None:
461 w1=np.ones((2,DIM))
462
463 wc=Data(0.,(2,2), Function(dom))
464 if w_gc is None:
465 wc[0,1]=1
466 else:
467 wc[0,1]=w_gc
468
469 reg_mask=Data(0.,(2,), Solution(dom))
470 reg_mask[self.DENSITY] = domainbuilder.getSetDensityMask()
471 reg_mask[self.SUSCEPTIBILITY] = domainbuilder.getSetSusceptibilityMask()
472 regularization=Regularization(dom, numLevelSets=2,\
473 w0=w0, w1=w1,wc=wc, location_of_set_m=reg_mask)
474 #====================================================================
475 self.logger.info("Retrieving gravity surveys...")
476 surveys=domainbuilder.getGravitySurveys()
477 g=[]
478 w=[]
479 for g_i,sigma_i in surveys:
480 w_i=safeDiv(1., sigma_i)
481 if g_i.getRank()==0:
482 g_i=g_i*kronecker(DIM)[DIM-1]
483 if w_i.getRank()==0:
484 w_i=w_i*kronecker(DIM)[DIM-1]
485 g.append(g_i)
486 w.append(w_i)
487 self.logger.debug("Added gravity survey:")
488 self.logger.debug("g = %s"%g_i)
489 self.logger.debug("sigma = %s"%sigma_i)
490 self.logger.debug("w = %s"%w_i)
491
492 self.logger.info("Setting up gravity model...")
493 gravity_model=GravityModel(dom, w, g)
494 gravity_model.rescaleWeights(rho_scale=rho_scale_mapping)
495 #====================================================================
496 self.logger.info("Retrieving magnetic field surveys...")
497 d_b=normalize(domainbuilder.getBackgroundMagneticFluxDensity())
498 surveys=domainbuilder.getMagneticSurveys()
499 B=[]
500 w=[]
501 for B_i,sigma_i in surveys:
502 w_i=safeDiv(1., sigma_i)
503 if B_i.getRank()==0:
504 B_i=B_i*d_b
505 if w_i.getRank()==0:
506 w_i=w_i*d_b
507 B.append(B_i)
508 w.append(w_i)
509 self.logger.debug("Added magnetic survey:")
510 self.logger.debug("B = %s"%B_i)
511 self.logger.debug("sigma = %s"%sigma_i)
512 self.logger.debug("w = %s"%w_i)
513
514 self.logger.info("Setting up magnetic model...")
515 magnetic_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity())
516 magnetic_model.rescaleWeights(k_scale=k_scale_mapping)
517 #====================================================================
518 self.logger.info("Setting cost function...")
519 self.setCostFunction(InversionCostFunction(regularization,
520 ((rho_mapping,self.DENSITY), (k_mapping, self.SUSCEPTIBILITY)),
521 ((gravity_model,0), (magnetic_model,1)) ))
522
523 def setInitialGuess(self, rho=None, k=None):
524 """
525 set the initial guess *rho* for density and *k* for susceptibility for the inversion iteration.
526
527 :param rho: initial value for the density anomaly.
528 :type rho: `Scalar`
529 :param k: initial value for the susceptibility anomaly.
530 :type k: `Scalar`
531 """
532 super(JointGravityMagneticInversion,self).setInitialGuess(rho, k)
533
534 def siloWriterCallback(self, k, m, Jm, g_Jm):
535 """
536 callback function that can be used to track the solution
537
538 :param k: iteration count
539 :param m: current m approximation
540 :param Jm: value of cost function
541 :param g_Jm: gradient of f at x
542 """
543 fn='inv.%d'%k
544 p=self.getCostFunction().getProperties(m)
545 ds=createDataset(density=p[self.DENSITY], susceptibility=p[self.SUSCEPTIBILITY])
546 ds.setCycleAndTime(k,k)
547 ds.saveSilo(fn)
548 self.logger.debug("J(m) = %e"%Jm)
549

  ViewVC Help
Powered by ViewVC 1.1.26