/[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 4059 by caltinay, Wed Oct 31 07:27:35 2012 UTC revision 4060 by caltinay, Fri Nov 9 03:50:19 2012 UTC
# Line 22  __license__="""Licensed under the Open S Line 22  __license__="""Licensed under the Open S
22  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
23  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
24    
25  import logging  __all__ = ['InversionBase','SingleParameterInversionBase','GravityInversion','MagneticInversion']
26    
27    import logging
28  from esys.escript import *  from esys.escript import *
29  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
30  from esys.weipa import createDataset  from esys.weipa import createDataset
31    
32  from costfunctions import SimpleCostFunction  from .costfunctions import SimpleCostFunction
33  from forwardmodels import GravityModel, ForwardModel, MagneticModel  from .forwardmodels import GravityModel, MagneticModel
34  from mappings import *  from .mappings import *
35  from minimizers import *  from .minimizers import *
36  from regularizations import Regularization  from .regularizations import Regularization
37  from datasources import DataSource  from .datasources import DataSource
38    
39  class InversionBase(object):  class InversionBase(object):
40      """      """
# Line 257  class GravityInversion(SingleParameterIn Line 258  class GravityInversion(SingleParameterIn
258      """      """
259      Inversion of Gravity (Bouguer) anomaly data.      Inversion of Gravity (Bouguer) anomaly data.
260      """      """
261      def setup(self, source, rho_ref=None, w0=None, w=None):      def setup(self, domainbuilder, rho_ref=None, w0=None, w=None):
262          """          """
263          Sets up the inversion parameters from a downunder.`DataSource`.          Sets up the inversion parameters from a downunder.`DomainBuilder`.
264    
265          :param source: suitable data source          :param domainbuilder: Domain builder object with gravity source(s)
266          :type source: `DataSource`          :type domainbuilder: `DomainBuilder`
267          :param rho_ref: reference density. If not specified, zero is used.          :param rho_ref: reference density. If not specified, zero is used.
268          :type rho_ref: ``float`` or `Scalar`          :type rho_ref: ``float`` or `Scalar`
269          :param w0: weighting factor for L2 term in the regularization          :param w0: weighting factor for L2 term in the regularization
# Line 271  class GravityInversion(SingleParameterIn Line 272  class GravityInversion(SingleParameterIn
272          :type w: ``list`` of float or `Vector`          :type w: ``list`` of float or `Vector`
273          """          """
274          if rho_ref is None:          if rho_ref is None:
275            rho_ref=0.              rho_ref=0.
276    
277          self.logger.info('Retrieving domain...')          self.logger.info('Retrieving domain...')
278          self.setDomain(source.getDomain())          self.setDomain(domainbuilder.getDomain())
279          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
280    
281          self.logger.info("Retrieving density mask...")          self.logger.info("Retrieving density mask...")
282          rho_mask = source.getSetDensityMask()          rho_mask = domainbuilder.getSetDensityMask()
283          if w is None:          if w is None:
284              w=[1.]*DIM              w=[1.]*DIM
285          m_ref=self.getMapping().getInverse(rho_ref)          m_ref=self.getMapping().getInverse(rho_ref)
# Line 287  class GravityInversion(SingleParameterIn Line 288  class GravityInversion(SingleParameterIn
288          self.setRegularization(Regularization(self.getDomain(),\          self.setRegularization(Regularization(self.getDomain(),\
289                  m_ref=m_ref, w0=w0, w=w, location_of_set_m=rho_mask))                  m_ref=m_ref, w0=w0, w=w, location_of_set_m=rho_mask))
290    
291          self.logger.info("Retrieving gravity and standard deviation data...")          self.logger.info("Retrieving gravity surveys...")
292          g, sigma=source.getGravityAndStdDev()          surveys=domainbuilder.getGravitySurveys()
293          chi=safeDiv(1., sigma*sigma)          g=[]
294          self.logger.debug("g = %s"%g)          chi=[]
295          self.logger.debug("sigma = %s"%sigma)          for g_i,sigma_i in surveys:
296          self.logger.debug("chi = %s"%chi)              chi_i=safeDiv(1., sigma_i*sigma_i)
297          chi=chi*kronecker(DIM)[DIM-1]              if g_i.getRank()==0:
298                    g_i=g_i*kronecker(DIM)[DIM-1]
299                if chi_i.getRank()==0:
300                    chi_i=chi_i*kronecker(DIM)[DIM-1]
301                g.append(g_i)
302                chi.append(chi_i)
303                self.logger.debug("Added gravity survey:")
304                self.logger.debug("g = %s"%g_i)
305                self.logger.debug("sigma = %s"%sigma_i)
306                self.logger.debug("chi = %s"%chi_i)
307    
308          self.logger.info("Setting up model...")          self.logger.info("Setting up model...")
309          self.setForwardModel(GravityModel(self.getDomain(), chi, g))          self.setForwardModel(GravityModel(self.getDomain(), chi, g))
# Line 312  class MagneticInversion(SingleParameterI Line 322  class MagneticInversion(SingleParameterI
322      """      """
323      Inversion of magnetic data.      Inversion of magnetic data.
324      """      """
325      def setup(self, source, k_ref=None, w0=None, w=None):      def setup(self, domainbuilder, k_ref=None, w0=None, w=None):
326          """          """
327          Sets up the inversion parameters from a downunder.`DataSource`.          Sets up the inversion parameters from a downunder.`DomainBuilder`.
328    
329          :param source: suitable data source          :param domainbuilder: Domain builder object with magnetic source(s)
330          :type source: `DataSource`          :type domainbuilder: `DomainBuilder`
331          :param k_ref: reference susceptibility. If not specified, zero is used.          :param k_ref: reference susceptibility. If not specified, zero is used.
332          :type k_ref: ``float`` or `Scalar`          :type k_ref: ``float`` or `Scalar`
333          :param w0: weighting factor for the L2 term in the regularization          :param w0: weighting factor for the L2 term in the regularization
# Line 329  class MagneticInversion(SingleParameterI Line 339  class MagneticInversion(SingleParameterI
339              k_ref=0.              k_ref=0.
340    
341          self.logger.info('Retrieving domain...')          self.logger.info('Retrieving domain...')
342          self.setDomain(source.getDomain())          self.setDomain(domainbuilder.getDomain())
343          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
344    
345          self.logger.info("Retrieving susceptibility mask...")          self.logger.info("Retrieving susceptibility mask...")
346          k_mask = source.getSetSusceptibilityMask()          k_mask = domainbuilder.getSetSusceptibilityMask()
347    
348          if w is None:          if w is None:
349              w=[1.]*DIM              w=[1.]*DIM
# Line 343  class MagneticInversion(SingleParameterI Line 353  class MagneticInversion(SingleParameterI
353          self.setRegularization(Regularization(self.getDomain(),\          self.setRegularization(Regularization(self.getDomain(),\
354                  m_ref=m_ref, w0=w0, w=w, location_of_set_m=k_mask))                  m_ref=m_ref, w0=w0, w=w, location_of_set_m=k_mask))
355    
356          self.logger.info("Retrieving magnetic field and standard deviation data...")          self.logger.info("Retrieving magnetic field surveys...")
357          B, sigma=source.getMagneticFieldAndStdDev()          surveys=domainbuilder.getMagneticSurveys()
358          chi=safeDiv(1., sigma*sigma)          B=[]
359          self.logger.debug("B = %s"%B)          chi=[]
360          self.logger.debug("sigma = %s"%sigma)          for B_i,sigma_i in surveys:
361          self.logger.debug("chi = %s"%chi)              chi_i=safeDiv(1., sigma_i*sigma_i)
362          chi=chi*kronecker(DIM)[DIM-1]              if B_i.getRank()==0:
363                    B_i=B_i*kronecker(DIM)[DIM-1]
364                if chi_i.getRank()==0:
365                    chi_i=chi_i*kronecker(DIM)[DIM-1]
366                B.append(B_i)
367                chi.append(chi_i)
368                self.logger.debug("Added magnetic survey:")
369                self.logger.debug("B = %s"%B_i)
370                self.logger.debug("sigma = %s"%sigma_i)
371                self.logger.debug("chi = %s"%chi_i)
372    
373          self.logger.info("Setting up model...")          self.logger.info("Setting up model...")
374          self.setForwardModel(MagneticModel(self.getDomain(), chi, B, source.getBackgroundMagneticField()))          self.setForwardModel(MagneticModel(self.getDomain(), chi, B, domainbuilder.getBackgroundMagneticField()))
375    
376          if self._mu_reg is None:          if self._mu_reg is None:
377              x=self.getDomain().getX()              x=self.getDomain().getX()

Legend:
Removed from v.4059  
changed lines
  Added in v.4060

  ViewVC Help
Powered by ViewVC 1.1.26