/[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 4213 - (show annotations)
Tue Feb 19 01:16:29 2013 UTC (7 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 22788 byte(s)
Some cleanup and more consistent logging.

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. The class uses the standard
204 'Regularization` class for a single level set function, `DensityMapping` mapping,
205 and the gravity forward model `GravityModel`.
206 """
207 def setup(self, domainbuilder,
208 rho0=None, drho=None, z0=None, beta=None,
209 w0=None, w1=None):
210 """
211 Sets up the inversion parameters from a `DomainBuilder`.
212
213 :param domainbuilder: Domain builder object with gravity source(s)
214 :type domainbuilder: `DomainBuilder`
215 :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
216 :type rho0: ``float`` or ``Scalar``
217 :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
218 :type drho: ``float`` or ``Scalar``
219 :param z0: reference depth for depth weighting, see `DensityMapping`. If not specified, zero is used.
220 :type z0: ``float`` or ``Scalar``
221 :param beta: exponent for depth weighting, see `DensityMapping`. If not specified, zero is used.
222 :type beta: ``float`` or ``Scalar``
223 :param w0: weighting factor for level set term regularization. If not set zero is assumed.
224 :type w0: ``Scalar`` or ``float``
225 :param w1: weighting factor for the gradient term in the regularization. If not set zero is assumed
226 :type w1: ``Vector`` or list of ``float``
227 """
228 self.logger.info('Retrieving domain...')
229 dom=domainbuilder.getDomain()
230 DIM=dom.getDim()
231 #========================
232 self.logger.info('Creating mapping...')
233 rho_mapping=DensityMapping(dom, rho0=rho0, drho=drho, z0=z0, beta=beta)
234 scale_mapping=rho_mapping.getTypicalDerivative()
235 #========================
236 self.logger.info("Setting up regularization...")
237 if w1 is None:
238 w1=[1.]*DIM
239 rho_mask = domainbuilder.getSetDensityMask()
240 regularization=Regularization(dom, numLevelSets=1,\
241 w0=w0, w1=w1, location_of_set_m=rho_mask)
242 #====================================================================
243 self.logger.info("Retrieving gravity surveys...")
244 surveys=domainbuilder.getGravitySurveys()
245 g=[]
246 w=[]
247 for g_i,sigma_i in surveys:
248 w_i=safeDiv(1., sigma_i)
249 if g_i.getRank()==0:
250 g_i=g_i*kronecker(DIM)[DIM-1]
251 if w_i.getRank()==0:
252 w_i=w_i*kronecker(DIM)[DIM-1]
253 g.append(g_i)
254 w.append(w_i)
255 self.logger.debug("Added gravity survey:")
256 self.logger.debug("g = %s"%g_i)
257 self.logger.debug("sigma = %s"%sigma_i)
258 self.logger.debug("w = %s"%w_i)
259 #====================================================================
260
261 self.logger.info("Setting up model...")
262 forward_model=GravityModel(dom, w, g)
263 forward_model.rescaleWeights(rho_scale=scale_mapping)
264
265 #====================================================================
266 self.logger.info("Setting cost function...")
267 self.setCostFunction(InversionCostFunction(regularization, rho_mapping, forward_model))
268
269 def setInitialGuess(self, rho=None):
270 """
271 set the initial guess *rho* for density the inversion iteration. If no *rho* present
272 then an appropriate initial guess is chosen.
273
274 :param rho: initial value for the density anomaly.
275 :type rho: `Scalar`
276 """
277 if rho:
278 super(GravityInversion,self).setInitialGuess(rho)
279 else:
280 super(GravityInversion,self).setInitialGuess()
281
282 def siloWriterCallback(self, k, m, Jm, g_Jm):
283 """
284 callback function that can be used to track the solution
285
286 :param k: iteration count
287 :param m: current m approximation
288 :param Jm: value of cost function
289 :param g_Jm: gradient of f at x
290 """
291 fn='inv.%d'%k
292 ds=createDataset(density=self.getCostFunction().getProperties(m))
293 ds.setCycleAndTime(k,k)
294 ds.saveSilo(fn)
295 self.logger.debug("J(m) = %e"%Jm)
296
297 class MagneticInversion(InversionDriver):
298 """
299 Driver class to perform an inversion of magnetic anomaly data. The class uses the standard
300 `Regularization` class for a single level set function. 'SusceptibilityMapping` mapping
301 and the linear magnetic forward model `MagneticModel`
302 """
303 def setup(self, domainbuilder,
304 k0=None, dk=None, z0=None, beta=None,
305 w0=None, w1=None):
306 """
307 Sets up the inversion from a `DomainBuilder`.
308 If magnetic data are given as scalar it is assumed that values are collected in direction of
309 the background magnetic field.
310
311 :param domainbuilder: Domain builder object with gravity source(s)
312 :type domainbuilder: `DomainBuilder`
313 :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
314 :type k0: ``float`` or ``Scalar``
315 :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
316 :type dk: ``float`` or ``Scalar``
317 :param z0: reference depth for depth weighting, see `SusceptibilityMapping`. If not specified, zero is used.
318 :type z0: ``float`` or ``Scalar``
319 :param beta: exponent for depth weighting, see `SusceptibilityMapping`. If not specified, zero is used.
320 :type beta: ``float`` or ``Scalar``
321 :param w0: weighting factor for level set term regularization. If not set zero is assumed.
322 :type w0: ``Scalar`` or ``float``
323 :param w1: weighting factor for the gradient term in the regularization. If not set zero is assumed
324 :type w1: ``Vector`` or list of ``float``
325 """
326 self.logger.info('Retrieving domain...')
327 dom=domainbuilder.getDomain()
328 DIM=dom.getDim()
329
330 #========================
331 self.logger.info('Creating mapping ...')
332 k_mapping=SusceptibilityMapping(dom, k0=k0, dk=dk, z0=z0, beta=beta)
333 scale_mapping=k_mapping.getTypicalDerivative()
334 #========================
335 self.logger.info("Setting up regularization...")
336 if w1 is None:
337 w1=[1.]*DIM
338 k_mask = domainbuilder.getSetSusceptibilityMask()
339 regularization=Regularization(dom, numLevelSets=1,w0=w0, w1=w1, location_of_set_m=k_mask)
340
341 #====================================================================
342 self.logger.info("Retrieving magnetic field surveys...")
343 d_b=normalize(domainbuilder.getBackgroundMagneticFluxDensity())
344 surveys=domainbuilder.getMagneticSurveys()
345 B=[]
346 w=[]
347 for B_i,sigma_i in surveys:
348 w_i=safeDiv(1., sigma_i)
349 if B_i.getRank()==0:
350 B_i=B_i*d_b
351 if w_i.getRank()==0:
352 w_i=w_i*d_b
353 B.append(B_i)
354 w.append(w_i)
355 self.logger.debug("Added magnetic survey:")
356 self.logger.debug("B = %s"%B_i)
357 self.logger.debug("sigma = %s"%sigma_i)
358 self.logger.debug("w = %s"%w_i)
359 #====================================================================
360 self.logger.info("Setting up model...")
361 forward_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity())
362 forward_model.rescaleWeights(k_scale=scale_mapping)
363
364 #====================================================================
365 self.logger.info("Setting cost function...")
366 self.setCostFunction(InversionCostFunction(regularization, k_mapping, forward_model))
367
368 def setInitialGuess(self, k=None):
369 """
370 set the initial guess *k* for susceptibility for the inversion iteration. If no *k* present
371 then an appropriate initial guess is chosen.
372
373 :param k: initial value for the susceptibility anomaly.
374 :type k: `Scalar`
375 """
376 if k:
377 super(MagneticInversion,self).setInitialGuess(k)
378 else:
379 super(MagneticInversion,self).setInitialGuess()
380
381 def siloWriterCallback(self, k, m, Jm, g_Jm):
382 """
383 callback function that can be used to track the solution
384
385 :param k: iteration count
386 :param m: current m approximation
387 :param Jm: value of cost function
388 :param g_Jm: gradient of f at x
389 """
390 fn='inv.%d'%k
391 ds=createDataset(susceptibility=self.getCostFunction().getProperties(m))
392 ds.setCycleAndTime(k,k)
393 ds.saveSilo(fn)
394 self.logger.debug("J(m) = %e"%Jm)
395
396 class JointGravityMagneticInversion(InversionDriver):
397 """
398 Driver class to perform a joint inversion of Gravity (Bouguer) and magnetic anomaly data. The class uses
399 the standard `Regularization` class for two level set functions with cross-gradient correlation. '
400 DensityMapping' and 'SusceptibilityMapping' mappings, the gravity forward model 'GravityModel'
401 and the linear magnetic forward model 'MagneticModel.
402 """
403 DENSITY=0
404 SUSCEPTIBILITY=1
405
406 def setup(self, domainbuilder,
407 rho0=None, drho=None, rho_z0=None, rho_beta=None,
408 k0=None, dk=None, k_z0=None, k_beta=None, w0=None, w1=None,
409 w_gc=None):
410 """
411 Sets up the inversion from an instance ``domainbuilder`` of a `DomainBuilder`.
412 Gravity and magnetic data attached to the \member{domainbuilder} are considered in the inversion.
413 If magnetic data are given as scalar it is assumed that values are collected in direction of
414 the background magnetic field.
415
416 :param domainbuilder: Domain builder object with gravity source(s)
417 :type domainbuilder: `DomainBuilder`
418
419 :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
420 :type rho0: ``float`` or `Scalar`
421 :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
422 :type drho: ``float`` or `Scalar`
423 :param rho_z0: reference depth for depth weighting for density, see `DensityMapping`. If not specified, zero is used.
424 :type rho_z0: ``float`` or `Scalar`
425 :param rho_beta: exponent for depth weighting for density, see `DensityMapping`. If not specified, zero is used.
426 :type rho_beta: ``float`` or `Scalar`
427 :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
428 :type k0: ``float`` or `Scalar`
429 :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
430 :type dk: ``float`` or `Scalar`
431 :param k_z0: reference depth for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
432 :type k_z0: ``float`` or `Scalar`
433 :param k_beta: exponent for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
434 :type k_beta: ``float`` or `Scalar`
435 :param w0: weighting factors for level set term regularization, see `Regularization`. If not set zero is assumed.
436 :type w0: `Data` or ``ndarray`` of shape (2,)
437 :param w1: weighting factor for the gradient term in the regularization see `Regularization`. If not set zero is assumed
438 :type w1: `Data` or ``ndarray`` of shape (2,DIM)
439 :param w_gc: weighting factor for the cross gradient term in the regularization, see `Regularization`. If not set one is assumed
440 :type w_gc: `Scalar` or `float`
441 """
442 self.logger.info('Retrieving domain...')
443 dom=domainbuilder.getDomain()
444 DIM=dom.getDim()
445
446 #========================
447 self.logger.info('Creating mappings ...')
448 rho_mapping=DensityMapping(dom, rho0=rho0, drho=drho, z0=rho_z0, beta=rho_beta)
449 rho_scale_mapping=rho_mapping.getTypicalDerivative()
450 self.logger.debug("rho_scale_mapping = %s"%rho_scale_mapping)
451 k_mapping=SusceptibilityMapping(dom, k0=k0, dk=dk, z0=k_z0, beta=k_beta)
452 k_scale_mapping=k_mapping.getTypicalDerivative()
453 self.logger.debug("k_scale_mapping = %s"%k_scale_mapping)
454 #========================
455 self.logger.info("Setting up regularization...")
456
457 if w1 is None:
458 w1=np.ones((2,DIM))
459
460 wc=Data(0.,(2,2), Function(dom))
461 if w_gc is None:
462 wc[0,1]=1
463 else:
464 wc[0,1]=w_gc
465
466 reg_mask=Data(0.,(2,), Solution(dom))
467 reg_mask[self.DENSITY] = domainbuilder.getSetDensityMask()
468 reg_mask[self.SUSCEPTIBILITY] = domainbuilder.getSetSusceptibilityMask()
469 regularization=Regularization(dom, numLevelSets=2,\
470 w0=w0, w1=w1,wc=wc, location_of_set_m=reg_mask)
471 #====================================================================
472 self.logger.info("Retrieving gravity surveys...")
473 surveys=domainbuilder.getGravitySurveys()
474 g=[]
475 w=[]
476 for g_i,sigma_i in surveys:
477 w_i=safeDiv(1., sigma_i)
478 if g_i.getRank()==0:
479 g_i=g_i*kronecker(DIM)[DIM-1]
480 if w_i.getRank()==0:
481 w_i=w_i*kronecker(DIM)[DIM-1]
482 g.append(g_i)
483 w.append(w_i)
484 self.logger.debug("Added gravity survey:")
485 self.logger.debug("g = %s"%g_i)
486 self.logger.debug("sigma = %s"%sigma_i)
487 self.logger.debug("w = %s"%w_i)
488
489 self.logger.info("Setting up gravity model...")
490 gravity_model=GravityModel(dom, w, g)
491 gravity_model.rescaleWeights(rho_scale=rho_scale_mapping)
492 #====================================================================
493 self.logger.info("Retrieving magnetic field surveys...")
494 d_b=normalize(domainbuilder.getBackgroundMagneticFluxDensity())
495 surveys=domainbuilder.getMagneticSurveys()
496 B=[]
497 w=[]
498 for B_i,sigma_i in surveys:
499 w_i=safeDiv(1., sigma_i)
500 if B_i.getRank()==0:
501 B_i=B_i*d_b
502 if w_i.getRank()==0:
503 w_i=w_i*d_b
504 B.append(B_i)
505 w.append(w_i)
506 self.logger.debug("Added magnetic survey:")
507 self.logger.debug("B = %s"%B_i)
508 self.logger.debug("sigma = %s"%sigma_i)
509 self.logger.debug("w = %s"%w_i)
510
511 self.logger.info("Setting up magnetic model...")
512 magnetic_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity())
513 magnetic_model.rescaleWeights(k_scale=k_scale_mapping)
514 #====================================================================
515 self.logger.info("Setting cost function...")
516 self.setCostFunction(InversionCostFunction(regularization,
517 ((rho_mapping,self.DENSITY), (k_mapping, self.SUSCEPTIBILITY)),
518 ((gravity_model,0), (magnetic_model,1)) ))
519
520 def setInitialGuess(self, rho=None, k=None):
521 """
522 set the initial guess *rho* for density and *k* for susceptibility for the inversion iteration.
523
524 :param rho: initial value for the density anomaly.
525 :type rho: `Scalar`
526 :param k: initial value for the susceptibility anomaly.
527 :type k: `Scalar`
528 """
529 super(JointGravityMagneticInversion,self).setInitialGuess(rho, k)
530
531 def siloWriterCallback(self, k, m, Jm, g_Jm):
532 """
533 callback function that can be used to track the solution
534
535 :param k: iteration count
536 :param m: current m approximation
537 :param Jm: value of cost function
538 :param g_Jm: gradient of f at x
539 """
540 fn='inv.%d'%k
541 p=self.getCostFunction().getProperties(m)
542 ds=createDataset(density=p[self.DENSITY], susceptibility=p[self.SUSCEPTIBILITY])
543 ds.setCycleAndTime(k,k)
544 ds.saveSilo(fn)
545 self.logger.debug("J(m) = %e"%Jm)
546

  ViewVC Help
Powered by ViewVC 1.1.26