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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4212 by jfenwick, Tue Jan 22 09:30:23 2013 UTC revision 4213 by caltinay, Tue Feb 19 01:16:29 2013 UTC
# Line 25  __url__="https://launchpad.net/escript-f Line 25  __url__="https://launchpad.net/escript-f
25  __all__ = ['InversionDriver', 'GravityInversion','MagneticInversion', 'JointGravityMagneticInversion']  __all__ = ['InversionDriver', 'GravityInversion','MagneticInversion', 'JointGravityMagneticInversion']
26    
27  import logging  import logging
28    import numpy as np
29    
30  from esys.escript import *  from esys.escript import *
31  import esys.escript.unitsSI as U  from esys.escript import unitsSI as U
32  from esys.weipa import createDataset  from esys.weipa import createDataset
 import numpy as np  
33    
34  from .inversioncostfunctions import InversionCostFunction  from .inversioncostfunctions import InversionCostFunction
35  from .forwardmodels import GravityModel, MagneticModel  from .forwardmodels import GravityModel, MagneticModel
# Line 41  class InversionDriver(object): Line 42  class InversionDriver(object):
42      """      """
43      Base class for running an inversion      Base class for running an inversion
44      """      """
45      def __init__(self, solverclass=None):      def __init__(self, solverclass=None, debug=False):
46          """          """
47          creates an instance of an inversion driver.          creates an instance of an inversion driver.
48            
49          :param solverclass: class of the solver to be used.          :param solverclass: class of the solver to be used.
50          :type solverclass: 'AbstractMinimizer'.          :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__)          self.logger=logging.getLogger('inv.%s'%self.__class__.__name__)
59          if solverclass is None: solverclass=MinimizerLBFGS          if solverclass is None: solverclass=MinimizerLBFGS
60          self.__solver=solverclass()          self.__solver=solverclass()
# Line 57  class InversionDriver(object): Line 64  class InversionDriver(object):
64    
65      def setCostFunction(self, costfunction):      def setCostFunction(self, costfunction):
66          """          """
67          sets the cost function of the inversion. This function needs to be called          sets the cost function of the inversion. This function needs to be called
68          before the inversion iteration can be started.          before the inversion iteration can be started.
69    
70          :param costfunction: domain of the inversion          :param costfunction: domain of the inversion
# Line 75  class InversionDriver(object): Line 82  class InversionDriver(object):
82              return self.getSolver().getCostFunction()              return self.getSolver().getCostFunction()
83          else:          else:
84             raise RuntimeError("Inversion is not set up.")             raise RuntimeError("Inversion is not set up.")
85        
86      def getSolver(self):      def getSolver(self):
87          """          """
88          The solver to be used in the inversion process. See the minimizers          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.          module for available solvers. By default, the L-BFGS minimizer is used.
90            
91          :rtype: 'AbstractMinimizer'.          :rtype: 'AbstractMinimizer'.
92          """          """
93          return self.__solver          return self.__solver
94                
95      def isSetUp(self):      def isSetUp(self):
96          """          """
97          returns True if the inversion is set up and is ready to run.          returns True if the inversion is set up and is ready to run.
98            
99          :rtype: `bool`          :rtype: `bool`
100          """          """
101          if  self.getSolver().getCostFunction():          if  self.getSolver().getCostFunction():
102              return True              return True
103          else:          else:
104              return False              return False
105                
106      def getDomain(self):      def getDomain(self):
107          """          """
108          returns the domain of the inversion          returns the domain of the inversion
# Line 103  class InversionDriver(object): Line 110  class InversionDriver(object):
110          :rtype: `Domain`          :rtype: `Domain`
111          """          """
112          self.getCostFunction().getDomain()          self.getCostFunction().getDomain()
113            
114      def setSolverCallback(self, callback=None):      def setSolverCallback(self, callback=None):
115          """          """
116          Sets the callback function which is called after every solver          Sets the callback function which is called after every solver
117          iteration.          iteration.
118          """          """
119          self.getSolver().setCallback(callback)          self.getSolver().setCallback(callback)
120            
           
121      def setSolverMaxIterations(self, maxiter=None):      def setSolverMaxIterations(self, maxiter=None):
122          """          """
123          Sets the maximum number of solver iterations to run. If ``maxiter`` is reached the iteration          Sets the maximum number of solver iterations to run.
124          is terminated and `MinimizerMaxIterReached` is thrown.          If `maxiter` is reached the iteration is terminated and
125                    `MinimizerMaxIterReached` is thrown.
126    
127          :param maxiter: maximum number of iteration steps.          :param maxiter: maximum number of iteration steps.
128          :type maxiter: positive ``int``          :type maxiter: positive ``int``
129          """          """
# Line 129  class InversionDriver(object): Line 136  class InversionDriver(object):
136    
137      def setSolverTolerance(self, m_tol=None, J_tol=None):      def setSolverTolerance(self, m_tol=None, J_tol=None):
138          """          """
139          Sets the error tolerance for the self.getSolver(). An acceptable solution is          Sets the error tolerance for the solver. An acceptable solution is
140          considered to be found once the tolerance is reached. is None          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          :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.                      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          :param J_tol: tolerance for changes to cost function. If `None` changes to the
145                      cost function are not checked for convergence during iteration.                      cost function are not checked for convergence during iteration.
146          :type m_tol: `float` or `None`          :type m_tol: `float` or `None`
147          :type J_tol: `float` or `None`          :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.          :note: if both arguments are equal to `None` the default setting m_tol=1e-4, J_tol=None is used.
149    
150          """          """
# Line 145  class InversionDriver(object): Line 152  class InversionDriver(object):
152              m_tol=1e-4              m_tol=1e-4
153          self.getSolver().setTolerance(m_tol=m_tol, J_tol=J_tol)          self.getSolver().setTolerance(m_tol=m_tol, J_tol=J_tol)
154    
   
155      def setInitialGuess(self, *args):      def setInitialGuess(self, *args):
156          """          """
157          set the initial guess for the inversion iteration. By default zero is used.          Sets the initial guess for the inversion iteration.
158                    By default zero is used.
159          """          """
160          self.initial_value=self.getCostFunction().createLevelSetFunction(*args)          self.initial_value=self.getCostFunction().createLevelSetFunction(*args)
161            
162      def run(self):      def run(self):
163          """          """
164          This function runs the inversion.          This function runs the inversion.
165            
166          :return: physical parameters as result of the inversion          :return: physical parameters as result of the inversion
167          :rtype: ``list`` of physical parameters or a physical parameter          :rtype: ``list`` of physical parameters or a physical parameter
168          """          """
169          if not self.isSetUp():          if not self.isSetUp():
170              raise RuntimeError("Inversion is not setup.")                        raise RuntimeError("Inversion is not setup.")
171    
172          if self.initial_value == None: self.setInitialGuess()          if self.initial_value == None: self.setInitialGuess()
173          self.logger.info("Starting self.getSolver()...")          self.logger.info("Starting solver, please stand by...")
174          try:          try:
175              self.getSolver().run(self.initial_value)              self.getSolver().run(self.initial_value)
176              self.m=self.getSolver().getResult()              self.m=self.getSolver().getResult()
177              self.p=self.getCostFunction().getProperties(self.m)              self.p=self.getCostFunction().getProperties(self.m)
178          except MinimizerException as e:          except MinimizerException as e:
179              self.m=self.getSolver().getResult()              self.m=self.getSolver().getResult()
180              self.p=self.getCostFunction().getProperties(self.m)              self.p=self.getCostFunction().getProperties(self.m)
181              self.logger.info("iteration failed.")              self.logger.info("Solver iteration failed!")
182              raise e              raise e
183          self.logger.info("iteration completed.")          self.logger.info("Solver completed successfully!")
184          self.getSolver().logSummary()          self.getSolver().logSummary()
185          return self.p          return self.p
186            
187      def getLevelSetFunction(self):      def getLevelSetFunction(self):
188          """          """
189          returns the level set function as solution of the optimization problem          returns the level set function as solution of the optimization problem
190            
191          :rtype: `Data`          :rtype: `Data`
192          """          """
193          return self.m          return self.m
194            
195      def setup(self, *args, **k_args):      def setup(self, *args, **k_args):
196          """          """
197          returns True if the inversion is set up and ready to run.          returns True if the inversion is set up and ready to run.
198          """          """
# Line 256  class GravityInversion(InversionDriver): Line 262  class GravityInversion(InversionDriver):
262          forward_model=GravityModel(dom, w, g)          forward_model=GravityModel(dom, w, g)
263          forward_model.rescaleWeights(rho_scale=scale_mapping)          forward_model.rescaleWeights(rho_scale=scale_mapping)
264    
   
265          #====================================================================          #====================================================================
266          self.logger.info("Setting cost function...")          self.logger.info("Setting cost function...")
267          self.setCostFunction(InversionCostFunction(regularization, rho_mapping, forward_model))          self.setCostFunction(InversionCostFunction(regularization, rho_mapping, forward_model))
268            
269      def setInitialGuess(self, rho=None):      def setInitialGuess(self, rho=None):
270          """          """
271          set the initial guess *rho* for density the inversion iteration. If no *rho* present          set the initial guess *rho* for density the inversion iteration. If no *rho* present
272          then an appropriate initial guess is chosen.          then an appropriate initial guess is chosen.
273            
274          :param rho: initial value for the density anomaly.          :param rho: initial value for the density anomaly.
275          :type rho: `Scalar`          :type rho: `Scalar`
276          """          """
# Line 273  class GravityInversion(InversionDriver): Line 278  class GravityInversion(InversionDriver):
278              super(GravityInversion,self).setInitialGuess(rho)              super(GravityInversion,self).setInitialGuess(rho)
279          else:          else:
280              super(GravityInversion,self).setInitialGuess()              super(GravityInversion,self).setInitialGuess()
281            
282      def siloWriterCallback(self, k, m, Jm, g_Jm):      def siloWriterCallback(self, k, m, Jm, g_Jm):
283          """          """
284          callback function that can be used to track the solution          callback function that can be used to track the solution
# Line 291  class GravityInversion(InversionDriver): Line 296  class GravityInversion(InversionDriver):
296    
297  class MagneticInversion(InversionDriver):  class MagneticInversion(InversionDriver):
298      """      """
299      Driver class to perform an inversion of magnetic anomaly data. The class uses the standard      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      `Regularization` class for a single level set function. 'SusceptibilityMapping` mapping
301      and the linear magnetic forward model `MagneticModel`      and the linear magnetic forward model `MagneticModel`
302      """      """
# Line 359  class MagneticInversion(InversionDriver) Line 364  class MagneticInversion(InversionDriver)
364          #====================================================================          #====================================================================
365          self.logger.info("Setting cost function...")          self.logger.info("Setting cost function...")
366          self.setCostFunction(InversionCostFunction(regularization, k_mapping, forward_model))          self.setCostFunction(InversionCostFunction(regularization, k_mapping, forward_model))
367            
368      def setInitialGuess(self, k=None):      def setInitialGuess(self, k=None):
369          """          """
370          set the initial guess *k* for susceptibility for the inversion iteration. If no *k* present          set the initial guess *k* for susceptibility for the inversion iteration. If no *k* present
371          then an appropriate initial guess is chosen.          then an appropriate initial guess is chosen.
372            
373          :param k: initial value for the susceptibility anomaly.          :param k: initial value for the susceptibility anomaly.
374          :type k: `Scalar`          :type k: `Scalar`
375          """          """
# Line 372  class MagneticInversion(InversionDriver) Line 377  class MagneticInversion(InversionDriver)
377              super(MagneticInversion,self).setInitialGuess(k)              super(MagneticInversion,self).setInitialGuess(k)
378          else:          else:
379              super(MagneticInversion,self).setInitialGuess()              super(MagneticInversion,self).setInitialGuess()
380                
381      def siloWriterCallback(self, k, m, Jm, g_Jm):      def siloWriterCallback(self, k, m, Jm, g_Jm):
382          """          """
383          callback function that can be used to track the solution          callback function that can be used to track the solution
# Line 387  class MagneticInversion(InversionDriver) Line 392  class MagneticInversion(InversionDriver)
392          ds.setCycleAndTime(k,k)          ds.setCycleAndTime(k,k)
393          ds.saveSilo(fn)          ds.saveSilo(fn)
394          self.logger.debug("J(m) = %e"%Jm)          self.logger.debug("J(m) = %e"%Jm)
395            
396  class JointGravityMagneticInversion(InversionDriver):  class JointGravityMagneticInversion(InversionDriver):
397      """      """
398       Driver class to perform a joint inversion of  Gravity (Bouguer) and magnetic anomaly data.  The class uses       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. '       the standard `Regularization` class for two level set functions with cross-gradient correlation. '
400       DensityMapping' and 'SusceptibilityMapping' mappings, the gravity forward model 'GravityModel'       DensityMapping' and 'SusceptibilityMapping' mappings, the gravity forward model 'GravityModel'
401       and the linear magnetic forward model 'MagneticModel.       and the linear magnetic forward model 'MagneticModel.
402      """      """
403      DENSITY=0      DENSITY=0
404      SUSCEPTIBILITY=1      SUSCEPTIBILITY=1
405        
406      def setup(self, domainbuilder,      def setup(self, domainbuilder,
407                      rho0=None, drho=None, rho_z0=None, rho_beta=None,                      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,                      k0=None, dk=None, k_z0=None, k_beta=None, w0=None, w1=None,
# Line 407  class JointGravityMagneticInversion(Inve Line 412  class JointGravityMagneticInversion(Inve
412          Gravity and magnetic data attached to the \member{domainbuilder} are considered in the inversion.          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          If magnetic data are given as scalar it is assumed that values are collected in direction of
414          the background magnetic field.          the background magnetic field.
415            
416          :param domainbuilder: Domain builder object with gravity source(s)          :param domainbuilder: Domain builder object with gravity source(s)
417          :type domainbuilder: `DomainBuilder`          :type domainbuilder: `DomainBuilder`
418            
419          :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.          :param rho0: reference density, see `DensityMapping`. If not specified, zero is used.
420          :type rho0: ``float`` or `Scalar`          :type rho0: ``float`` or `Scalar`
421          :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.          :param drho: density scale, see `DensityMapping`. If not specified, 2750kg/m^3 is used.
422          :type drho: ``float`` or `Scalar`            :type drho: ``float`` or `Scalar`
423          :param rho_z0: reference depth for depth weighting for density, see `DensityMapping`. If not specified, zero is used.          :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`          :type rho_z0: ``float`` or `Scalar`
425          :param rho_beta: exponent for  depth weighting  for density, see `DensityMapping`. If not specified, zero is used.          :param rho_beta: exponent for  depth weighting  for density, see `DensityMapping`. If not specified, zero is used.
426          :type rho_beta: ``float`` or `Scalar`                  :type rho_beta: ``float`` or `Scalar`
427          :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.          :param k0: reference susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
428          :type k0: ``float`` or `Scalar`          :type k0: ``float`` or `Scalar`
429          :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.          :param dk: susceptibility scale, see `SusceptibilityMapping`. If not specified, 2750kg/m^3 is used.
430          :type dk: ``float`` or `Scalar`            :type dk: ``float`` or `Scalar`
431          :param k_z0: reference depth for depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.          :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`          :type k_z0: ``float`` or `Scalar`
433          :param k_beta: exponent for  depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.          :param k_beta: exponent for  depth weighting for susceptibility, see `SusceptibilityMapping`. If not specified, zero is used.
434          :type k_beta: ``float`` or `Scalar`            :type k_beta: ``float`` or `Scalar`
435          :param w0: weighting factors for level set term regularization, see `Regularization`. If not set zero is assumed.          :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,)          :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          :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)            :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          :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`          :type w_gc: `Scalar` or `float`
441          """          """
442          self.logger.info('Retrieving domain...')          self.logger.info('Retrieving domain...')
443          dom=domainbuilder.getDomain()          dom=domainbuilder.getDomain()
# Line 442  class JointGravityMagneticInversion(Inve Line 447  class JointGravityMagneticInversion(Inve
447          self.logger.info('Creating mappings ...')          self.logger.info('Creating mappings ...')
448          rho_mapping=DensityMapping(dom, rho0=rho0, drho=drho, z0=rho_z0, beta=rho_beta)          rho_mapping=DensityMapping(dom, rho0=rho0, drho=drho, z0=rho_z0, beta=rho_beta)
449          rho_scale_mapping=rho_mapping.getTypicalDerivative()          rho_scale_mapping=rho_mapping.getTypicalDerivative()
450          print " rho_scale_mapping = ",rho_scale_mapping          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)          k_mapping=SusceptibilityMapping(dom, k0=k0, dk=dk, z0=k_z0, beta=k_beta)
452          k_scale_mapping=k_mapping.getTypicalDerivative()          k_scale_mapping=k_mapping.getTypicalDerivative()
453          print " rho_scale_mapping = ",k_scale_mapping          self.logger.debug("k_scale_mapping = %s"%k_scale_mapping)
454          #========================          #========================
455          self.logger.info("Setting up regularization...")          self.logger.info("Setting up regularization...")
456            
457          if w1 is None:          if w1 is None:
458              w1=np.ones((2,DIM))              w1=np.ones((2,DIM))
459            
460          wc=Data(0.,(2,2), Function(dom))          wc=Data(0.,(2,2), Function(dom))
461          if w_gc is  None:          if w_gc is  None:
462              wc[0,1]=1              wc[0,1]=1
463          else:          else:
464              wc[0,1]=w_gc              wc[0,1]=w_gc
465            
466                        reg_mask=Data(0.,(2,), Solution(dom))
         reg_mask=Data(0.,(2,), Solution(dom))    
467          reg_mask[self.DENSITY] = domainbuilder.getSetDensityMask()          reg_mask[self.DENSITY] = domainbuilder.getSetDensityMask()
468          reg_mask[self.SUSCEPTIBILITY] = domainbuilder.getSetSusceptibilityMask()          reg_mask[self.SUSCEPTIBILITY] = domainbuilder.getSetSusceptibilityMask()
469          regularization=Regularization(dom, numLevelSets=2,\          regularization=Regularization(dom, numLevelSets=2,\
# Line 503  class JointGravityMagneticInversion(Inve Line 507  class JointGravityMagneticInversion(Inve
507              self.logger.debug("B = %s"%B_i)              self.logger.debug("B = %s"%B_i)
508              self.logger.debug("sigma = %s"%sigma_i)              self.logger.debug("sigma = %s"%sigma_i)
509              self.logger.debug("w = %s"%w_i)              self.logger.debug("w = %s"%w_i)
510                
511          self.logger.info("Setting up magnetic model...")          self.logger.info("Setting up magnetic model...")
512          magnetic_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity())          magnetic_model=MagneticModel(dom, w, B, domainbuilder.getBackgroundMagneticFluxDensity())
513          magnetic_model.rescaleWeights(k_scale=k_scale_mapping)          magnetic_model.rescaleWeights(k_scale=k_scale_mapping)
514          #====================================================================          #====================================================================
515          self.logger.info("Setting cost function...")          self.logger.info("Setting cost function...")
516          self.setCostFunction(InversionCostFunction(regularization,          self.setCostFunction(InversionCostFunction(regularization,
517               ((rho_mapping,self.DENSITY), (k_mapping, self.SUSCEPTIBILITY)),               ((rho_mapping,self.DENSITY), (k_mapping, self.SUSCEPTIBILITY)),
518                 ((gravity_model,0), (magnetic_model,1)) ))                 ((gravity_model,0), (magnetic_model,1)) ))
519            
520      def setInitialGuess(self, rho=None, k=None):      def setInitialGuess(self, rho=None, k=None):
521          """          """
522          set the initial guess *rho* for density and *k* for susceptibility for the inversion iteration.          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.          :param rho: initial value for the density anomaly.
525          :type rho: `Scalar`          :type rho: `Scalar`
526          :param k: initial value for the susceptibility anomaly.          :param k: initial value for the susceptibility anomaly.
527          :type k: `Scalar`          :type k: `Scalar`
528          """          """
529          super(JointGravityMagneticInversion,self).setInitialGuess(rho, k)          super(JointGravityMagneticInversion,self).setInitialGuess(rho, k)
530            
531      def siloWriterCallback(self, k, m, Jm, g_Jm):      def siloWriterCallback(self, k, m, Jm, g_Jm):
532          """          """
533          callback function that can be used to track the solution          callback function that can be used to track the solution
# Line 539  class JointGravityMagneticInversion(Inve Line 543  class JointGravityMagneticInversion(Inve
543          ds.setCycleAndTime(k,k)          ds.setCycleAndTime(k,k)
544          ds.saveSilo(fn)          ds.saveSilo(fn)
545          self.logger.debug("J(m) = %e"%Jm)          self.logger.debug("J(m) = %e"%Jm)
546            

Legend:
Removed from v.4212  
changed lines
  Added in v.4213

  ViewVC Help
Powered by ViewVC 1.1.26