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

  ViewVC Help
Powered by ViewVC 1.1.26