In [1]:
## Forward model a data

In [2]:
import SimPEG as simpeg
import simpegMT as simpegmt
from simpegMT.Utils.dataUtils import rec2ndarr
import cPickle as pickle

In [3]:
# Make the mesh.
mTensor = simpeg.Utils.meshTensor
cSize = [50,20]
# Cells constant size mesh
hx = mTensor([(cSize[0],50)])
hy = mTensor([(cSize[0],50)])
hz = mTensor([(cSize[1],48)])
x0 = np.array([-1250,-1250,- 30*20])
mesh3dCons = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)
# With padding
hPad = mTensor([(cSize[0],5,1.5)])
aPad = mTensor([(cSize[1],9,1.5)])
bPad = mTensor([(cSize[1],9,-1.5)])
hxPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))
hyPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))
hzPad = np.hstack((bPad,mTensor([(cSize[1],30)]),aPad))
x0Pad = np.array([-(np.sum(hPad)+1000),-(np.sum(hPad)+1000),-(np.sum(bPad)+(20*30))])
mesh3d = simpeg.Mesh.TensorMesh([hxPad,hyPad,hzPad],x0Pad)

In [4]:
# Read the model
modelname = "simpegTDmodel.con"

# Load the model to the uniform cell mesh
modelUniCell = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3dCons)

# Load the model to the mesh with padding cells
modelT = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)
# Adjust the model to reflect changes in the mesh (fewer aircells)
modMat = mesh3d.r(modelT,'CC','CC','M')
modNewMat = np.ones((50,50,48))*modMat[0,0,0]
modNewMat[:,:,9::] = modMat[:,:,:-9]
modelTD = mesh3d.r(modNewMat,'CC','CC','V')

In [5]:
# Load the model to the uniform cell mesh
modelUniCell = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3dCons)
# Save as a vtk file
simpeg.Utils.meshutils.writeVTRFile('modelTDuniMesh.vtr',mesh3dCons,{'S/m':modelUniCell})

In [6]:
# Load the model to the mesh with padding cells
# modelTD = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)
# Save as a vtk file
simpeg.Utils.meshutils.writeVTRFile('modelTDpaddedMesh.vtr',mesh3d,{'S/m':modelTD})

In [7]:
# Define the data locations
xG,yG = np.meshgrid(np.linspace(-700,700,8),np.linspace(-700,700,8))
zG = np.zeros_like(xG)
locs = np.hstack((simpeg.mkvc(xG.ravel(),2),simpeg.mkvc(yG.ravel(),2),simpeg.mkvc(zG.ravel(),2)))

In [8]:
# Make the receiver list
rxList = []
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
    rxList.append(simpegmt.SurveyMT.RxMT(locs,rxType))    
# Source list
srcList =[]
freqs = np.logspace(4,0,17)
for freq in freqs:
    srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
# Survey MT
survey = simpegmt.SurveyMT.SurveyMT(srcList)


# Forward model the data
if False:
    fields = problem.fields(modelTD)
    mtData = survey.projectFields(fields)

In [9]:
# Load the data
mtSeogiData = simpegmt.DataMT.DataMT(survey,np.load('seogiModel_MTdata.npy'))
mtSeogirecData = mtSeogiData.toRecArray('Complex')


/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2834: FutureWarning: Numpy has detected that you (may be) writing to an array returned
by numpy.diagonal or by selecting multiple fields in a record
array. This code will likely break in a future numpy release --
see numpy.diagonal or arrays.indexing reference docs for details.
The quick fix is to make an explicit copy (e.g., do
arr.diagonal().copy() or arr[['f0','f1']].copy()).
  if (obj.__array_interface__["data"][0]
/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2835: FutureWarning: Numpy has detected that you (may be) writing to an array returned
by numpy.diagonal or by selecting multiple fields in a record
array. This code will likely break in a future numpy release --
see numpy.diagonal or arrays.indexing reference docs for details.
The quick fix is to make an explicit copy (e.g., do
arr.diagonal().copy() or arr[['f0','f1']].copy()).
  != self.__array_interface__["data"][0]):

In [10]:
sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf

In [11]:
% matplotlib qt
# Plot the data
mtSeogiMap = iPf.MTinteractiveMap([mtSeogirecData])

In [12]:
mtSeogirecData.shape


Out[12]:
(1088,)

In [13]:
# Add noise to the data
std = 0.05 # 5% std
if os.path.isfile('seogiModel_MTdata.npy') and os.path.isfile('seogiModel_dobs.npy'):
    d_true = np.load('seogiModel_MTdata.npy')
    d_obs = np.load('seogiModel_dobs.npy')
else:
    d_true = simpeg.mkvc(mtSeogiData)
    d_obs = d_true + std*abs(d_true)*np.random.randn(*d_true.shape)
    np.save('seogiModel_dobs.npy',d_obs)
# Assign the dobs
survey.dtrue = d_true
survey.dobs = d_obs
survey.std = np.abs(survey.dobs*std) #+ 0.01*np.linalg.norm(survey.dobs) #survey.dobs*0 + std
# Assign the data weight
Wd = 1/survey.std #(abs(survey.dobs)*survey.std)

In [14]:
# Make an MTdata object with dobs
mtSeogiDobs = simpegmt.DataMT.DataMT(survey,survey.dobs)

In [15]:
# Run the setup
mt1DdataZyxList = simpegmt.Utils.dataUtils.convert3Dto1Dobject(mtSeogiDobs,'zyx')
mt1DdataZxyList = simpegmt.Utils.dataUtils.convert3Dto1Dobject(mtSeogiDobs,'zxy')
mt1DdataZdetList = simpegmt.Utils.dataUtils.convert3Dto1Dobject(mtSeogiDobs,'zdet')

In [16]:
m1d = simpeg.Mesh.TensorMesh([mesh3d.hz],[mesh3d.x0[-1]])
# Setup the problem object
sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model
# Set the mapping
actMap = simpeg.Maps.ActiveCells(m1d, sigma1d > 1e-7, np.log(1e-8), nC=m1d.nCx)
mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
problem = simpegmt.ProblemMT1D.eForm_psField(m1d,sigmaPrimary = sigma1d,mapping=mappingExpAct)
problem.verbose = False

In [17]:
sigma1d


Out[17]:
array([  5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   1.00000000e-03,   1.00000000e-03,
         1.00000000e-03,   1.00000000e-03,   1.00000000e-03,
         1.00000000e-08,   1.00000000e-08,   1.00000000e-08,
         1.00000000e-08,   1.00000000e-08,   1.00000000e-08,
         1.00000000e-08,   1.00000000e-08,   1.00000000e-08])

In [22]:
# Define a function to run the forward problem
def runInversionModel(data,problem,m1d,nameflag):
    problem.unpair()
    data.survey.unpair()
    problem.pair(data.survey)
    # Define a counter
    C =  simpeg.Utils.Counter()
    # Set the optimization
    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 10)
    opt.counter = C
    opt.LSshorten = 0.5
    opt.remember('xc')
    # Data misfit
    dmis = simpeg.DataMisfit.l2_DataMisfit(data.survey)
    dmis.Wd = 1./data.survey.std
    # Regularization
    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
    reg = simpeg.Regularization.Tikhonov(regMesh)
    reg.smoothModel = False
    reg.alpha_s = 1e-6
    reg.alpha_x = 1.
    # reg.alpha_xx = .001
    # Inversion problem
    invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
    invProb.counter = C
    # Beta cooling
    beta = simpeg.Directives.BetaSchedule()
    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
    targmis = simpeg.Directives.TargetMisfit()
#     targmis.target = .6 * data.survey.nD
#     print 'Target misfit is {:.0f}'.format(targmis.target)
    locs = np.ones((regMesh.nC,1))*data.survey.srcList[0].rxList[0].locs[:,:-1]
    saveModel = simpeg.Directives.SaveModelEveryIteration()
    saveModel.fileName = 'Inversion_modAt{:.0f}_{:.0f}'.format(locs[0,0],locs[0,1])
    print saveModel.fileName
    # Create an inversion object
    inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
    # Run
    m_0 = np.log(1e-2+0*problem.sigmaPrimary[problem.mapping.sigmaMap.maps[-1].indActive])
    problem.survey.mtrue = m_0
    mopt = inv.run(m_0)
    # Save the model
    
    xyzv = np.hstack((locs,simpeg.mkvc(regMesh.gridCC,2),simpeg.mkvc(mopt,2)))
    np.save('xyzmod_{:s}_{:.0f}_{:.0f}.npy'.format(nameflag,locs[0,0],locs[0,1]),xyzv)
    return mopt

In [23]:
dat = mt1DdataZyxList[0]
mopt = runInversionModel(dat,problem,m1d,'zyx')


Inversion_modAt-700_-700
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  1.37e+03  1.07e+03  0.00e+00  1.07e+03    3.77e+02      0              
   1  1.37e+03  4.40e+01  1.93e-02  7.05e+01    7.43e+01      0              
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.0682e+02
1 : |xc-x_last| = 1.4393e+00 <= tolX*(1+|x0|) = 2.9759e+00
0 : |proj(x-g)-x|    = 7.4277e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 7.4277e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      10    <= iter          =      2
------------------------- DONE! -------------------------

In [25]:
dat.survey


Out[25]:
<simpegMT.SurveyMT.SurveyMT at 0x7f3614079310>

In [31]:
if True:
    for dat in mt1DdataZyxList:
        runInversionModel(dat,problem,m1d,'zyx')

    for dat in mt1DdataZxyList:
        runInversionModel(dat,problem,m1d,'zxy')
    for dat in mt1DdataZdetList:
        runInversionModel(dat,problem,m1d,'zdet')


Inversion_modAt-700_-700
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  1.69e+03  1.07e+03  0.00e+00  1.07e+03    3.77e+02      0              
   1  1.69e+03  3.88e+01  1.72e-02  6.80e+01    6.98e+01      0              
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.0682e+02
1 : |xc-x_last| = 1.3066e+00 <= tolX*(1+|x0|) = 2.9759e+00
0 : |proj(x-g)-x|    = 6.9756e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 6.9756e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      10    <= iter          =      2
------------------------- DONE! -------------------------
Inversion_modAt-700_-500
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  1.12e+03  1.10e+03  0.00e+00  1.10e+03    3.79e+02      0              
   1  1.12e+03  5.89e+01  2.11e-02  8.26e+01    9.39e+01      0              
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.1005e+02
1 : |xc-x_last| = 1.8888e+00 <= tolX*(1+|x0|) = 2.9759e+00
0 : |proj(x-g)-x|    = 9.3881e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 9.3881e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      10    <= iter          =      2
------------------------- DONE! -------------------------
Inversion_modAt-700_-300
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  7.28e+02  1.14e+03  0.00e+00  1.14e+03    3.81e+02      0              
   1  7.28e+02  6.54e+01  2.62e-02  8.44e+01    8.78e+01      0              
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.1411e+02
1 : |xc-x_last| = 1.8407e+00 <= tolX*(1+|x0|) = 2.9759e+00
0 : |proj(x-g)-x|    = 8.7820e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 8.7820e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      10    <= iter          =      2
------------------------- DONE! -------------------------
Inversion_modAt-700_-100
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-31-496b6a805e71> in <module>()
      1 if True:
      2     for dat in mt1DdataZyxList:
----> 3         runInversionModel(dat,problem,m1d,'zyx')
      4 
      5     for dat in mt1DdataZxyList:

<ipython-input-22-5072715ac24c> in runInversionModel(data, problem, m1d, nameflag)
     39     m_0 = np.log(1e-2+0*problem.sigmaPrimary[problem.mapping.sigmaMap.maps[-1].indActive])
     40     problem.survey.mtrue = m_0
---> 41     mopt = inv.run(m_0)
     42     # Save the model
     43 

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/CounterUtils.pyc in wrapper(self, *args, **kwargs)
     90         counter = getattr(self,'counter',None)
     91         if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
---> 92         out = f(self,*args,**kwargs)
     93         if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
     94         return out

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Inversion.pyc in run(self, m0)
     60         self.invProb.startup(m0)
     61         self.directiveList.call('initialize')
---> 62         self.m = self.opt.minimize(self.invProb.evalFunction, self.invProb.curModel)
     63         self.directiveList.call('finish')
     64 

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/CounterUtils.pyc in wrapper(self, *args, **kwargs)
     90         counter = getattr(self,'counter',None)
     91         if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
---> 92         out = f(self,*args,**kwargs)
     93         if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
     94         return out

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Optimization.pyc in minimize(self, evalFunction, x0)
    206         while True:
    207             self.doStartIteration()
--> 208             self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True)
    209             self.printIter()
    210             if self.stoppingCriteria(): break

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/CounterUtils.pyc in wrapper(self, *args, **kwargs)
     90         counter = getattr(self,'counter',None)
     91         if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
---> 92         out = f(self,*args,**kwargs)
     93         if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
     94         return out

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/InvProblem.pyc in evalFunction(self, m, return_g, return_H)
    124         out = (f,)
    125         if return_g:
--> 126             phi_dDeriv = self.dmisfit.evalDeriv(m, u=u)
    127             phi_mDeriv = self.reg.evalDeriv(m)
    128 

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/CounterUtils.pyc in wrapper(self, *args, **kwargs)
     90         counter = getattr(self,'counter',None)
     91         if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
---> 92         out = f(self,*args,**kwargs)
     93         if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
     94         return out

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/DataMisfit.pyc in evalDeriv(self, m, u)
    156         survey = self.survey
    157         if u is None: u = prob.fields(m)
--> 158         return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u=u)), u=u)
    159 
    160     @Utils.timeIt

/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.py in Jtvec(self, m, v, u)
    125                     # Get the
    126                     dA_duIT = ATinv * PTv
--> 127                     dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
    128                     dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
    129                     # Make du_dmT

/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.pyc in getADeriv_m(self, freq, u, v, adjoint)
     68         #
     69         u_src = u['e_1dSolution']
---> 70         dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
     71         if adjoint:
     72             return 1j * omega(freq) * (  dMf_dsig.T * v )

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/PropMaps.pyc in fget(self)
    107 
    108             m = getattr(self, '%sModel'%prop.name)
--> 109             return mapping.deriv( m )
    110         return property(fget=fget)
    111 

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Maps.pyc in deriv(self, m)
    181         mi = m
    182         for map_i in reversed(self.maps):
--> 183             deriv = map_i.deriv(mi) * deriv
    184             mi = map_i * mi
    185         return deriv

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Maps.pyc in deriv(self, m)
    249                 \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m})
    250         """
--> 251         return Utils.sdiag(np.exp(Utils.mkvc(m)))
    252 
    253 class ReciprocalMap(IdentityMap):

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/matutils.pyc in sdiag(h)
     38 def sdiag(h):
     39     """Sparse diagonal matrix"""
---> 40     return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
     41 
     42 def sdInv(M):

/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/construct.pyc in spdiags(data, diags, m, n, format)
     58 
     59     """
---> 60     return dia_matrix((data, diags), shape=(m,n)).asformat(format)
     61 
     62 

/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/base.pyc in asformat(self, format)
    211             return self
    212         else:
--> 213             return getattr(self,'to' + format)()
    214 
    215     ###################################################################

/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/dia.pyc in tocsr(self)
    225     def tocsr(self):
    226         #this could be faster
--> 227         return self.tocoo().tocsr()
    228 
    229     def tocsc(self):

/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/dia.pyc in tocoo(self)
    250 
    251         from .coo import coo_matrix
--> 252         return coo_matrix((data,(row,col)), shape=self.shape)
    253 
    254     # needed by _data_matrix

/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/coo.pyc in __init__(self, arg1, shape, dtype, copy)
    204             self.data = self.data.astype(dtype)
    205 
--> 206         self._check()
    207 
    208     def getnnz(self, axis=None):

/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/coo.pyc in _check(self)
    257 
    258         if nnz > 0:
--> 259             if self.row.max() >= self.shape[0]:
    260                 raise ValueError('row index exceeds matrix dimensions')
    261             if self.col.max() >= self.shape[1]:

KeyboardInterrupt: 

In [ ]:


In [32]:
# import multiprocessing, os, numpy as np
# from glob import glob
# import itertools as it
# import subprocess


# def runInv_star(args):
#     return runInversionModel(*args)


# if __name__ == '__main__':
#     pool = multiprocessing.Pool(processes=6)
#     pool.map(runInversionZxy,it.izip(mt1DdataZxyList,it.repeat( (problem,m1d,'zxy') )))
#     pool.join()

In [34]:
seogizxy1Dmods = np.vstack(([np.load(f) for f in glob('xyzmod_zxy*.npy')]))
seogizyx1Dmods = np.vstack(([np.load(f) for f in glob('xyzmod_zyx*.npy')]))
seogizdet1Dmods = np.vstack(([np.load(f) for f in glob('xyzmod_zdet*.npy')]))

In [ ]:


In [35]:
airInd = modelTD ==1e-8
# Interpolate the 1D models
from scipy.interpolate import griddata
mod3d_zxy1dMod = griddata(seogizxy1Dmods[:,0:3],seogizxy1Dmods[:,3],mesh3d.gridCC,'linear')
mod3d_zxy1dMod[airInd] = np.log(1e-8)
mod3d_zyx1dMod = griddata(seogizyx1Dmods[:,0:3],seogizyx1Dmods[:,3],mesh3d.gridCC,'linear')
mod3d_zyx1dMod[airInd] = np.log(1e-8)
mod3d_zdet1dMod = griddata(seogizdet1Dmods[:,0:3],seogizdet1Dmods[:,3],mesh3d.gridCC,'linear')
mod3d_zdet1dMod[airInd] = np.log(1e-8)

In [36]:
simpeg.Utils.meshutils.writeVTRFile('model3D_zxy1Dinv_paddedMesh.vtr',mesh3d,{'S/m':np.exp(mod3d_zxy1dMod)})
simpeg.Utils.meshutils.writeVTRFile('model3D_zyx1Dinv_paddedMesh.vtr',mesh3d,{'S/m':np.exp(mod3d_zyx1dMod)})
simpeg.Utils.meshutils.writeVTRFile('model3D_zdet1Dinv_paddedMesh.vtr',mesh3d,{'S/m':np.exp(mod3d_zdet1dMod)})

In [52]:
%matplotlib qt
m_0 = np.log(problem.sigmaPrimary[problem.mapping.sigmaMap.maps[-1].indActive])
simpegmt.Utils.dataUtils.plotMT1DModelData(problem,[m_0,mopt])
plt.suptitle('Target misfit-smooth False')
plt.show()

In [50]:
simpegmt.Utils.dataUtils.getAppRes(problem.dataPair(problem.survey,problem.survey.dpred(mopt)))


Out[50]:
[(794.79531859664246, 56.319860569098076),
 (674.66872912184408, 57.635385492661392),
 (559.52741397235366, 58.423832448812171),
 (459.43398962262597, 58.410853192444343),
 (380.28287396790699, 57.529223784875711),
 (322.8521167724395, 56.001248998922648),
 (283.43537890603488, 54.287024515800717),
 (255.71182522194459, 52.81526241219477),
 (233.73967183794295, 51.557597236493933),
 (216.97743323521141, 50.281238562685871),
 (204.41678397612844, 49.05493425669367),
 (196.11922783991341, 47.732975834266547),
 (192.81620390650565, 46.626437679858043),
 (191.83932042041988, 45.994944912824344),
 (191.13999703648773, 45.673623213889364),
 (190.40201746939533, 45.487392219519187),
 (189.72482164986508, 45.362566052988335)]

In [53]:
ls


model3D_zdet1Dinv_paddedMesh.vtr  xyzmod_zxy_300_500.npy
model3D_zxy1Dinv_paddedMesh.vtr   xyzmod_zxy_-300_-700.npy
model3D_zyx1Dinv_paddedMesh.vtr   xyzmod_zxy_-300_700.npy
modelTDpaddedMesh.vtr             xyzmod_zxy_300_-700.npy
modelTDuniMesh.vtr                xyzmod_zxy_300_700.npy
MT3DforData1Dinv.                 xyzmod_zxy_-500_-100.npy
MT3DforData1Dinv.ipynb            xyzmod_zxy_-500_100.npy
MT3DforData1Dinv.py               xyzmod_zxy_500_-100.npy
seogiModel_dobs.npy               xyzmod_zxy_500_100.npy
seogiModel_MTdata.npy             xyzmod_zxy_-500_-300.npy
simpegTDmodel.con                 xyzmod_zxy_-500_300.npy
xyzmod_zdet_-100_-100.npy         xyzmod_zxy_500_-300.npy
xyzmod_zdet_-100_100.npy          xyzmod_zxy_500_300.npy
xyzmod_zdet_100_-100.npy          xyzmod_zxy_-500_-500.npy
xyzmod_zdet_100_100.npy           xyzmod_zxy_-500_500.npy
xyzmod_zdet_-100_-300.npy         xyzmod_zxy_500_-500.npy
xyzmod_zdet_-100_300.npy          xyzmod_zxy_500_500.npy
xyzmod_zdet_100_-300.npy          xyzmod_zxy_-500_-700.npy
xyzmod_zdet_100_300.npy           xyzmod_zxy_-500_700.npy
xyzmod_zdet_-100_-500.npy         xyzmod_zxy_500_-700.npy
xyzmod_zdet_-100_500.npy          xyzmod_zxy_500_700.npy
xyzmod_zdet_100_-500.npy          xyzmod_zxy_-700_-100.npy
xyzmod_zdet_100_500.npy           xyzmod_zxy_-700_100.npy
xyzmod_zdet_-100_-700.npy         xyzmod_zxy_700_-100.npy
xyzmod_zdet_-100_700.npy          xyzmod_zxy_700_100.npy
xyzmod_zdet_100_-700.npy          xyzmod_zxy_-700_-300.npy
xyzmod_zdet_100_700.npy           xyzmod_zxy_-700_300.npy
xyzmod_zdet_-300_-100.npy         xyzmod_zxy_700_-300.npy
xyzmod_zdet_-300_100.npy          xyzmod_zxy_700_300.npy
xyzmod_zdet_300_-100.npy          xyzmod_zxy_-700_-500.npy
xyzmod_zdet_300_100.npy           xyzmod_zxy_-700_500.npy
xyzmod_zdet_-300_-300.npy         xyzmod_zxy_700_-500.npy
xyzmod_zdet_-300_300.npy          xyzmod_zxy_700_500.npy
xyzmod_zdet_300_-300.npy          xyzmod_zxy_-700_-700.npy
xyzmod_zdet_300_300.npy           xyzmod_zxy_-700_700.npy
xyzmod_zdet_-300_-500.npy         xyzmod_zxy_700_-700.npy
xyzmod_zdet_-300_500.npy          xyzmod_zxy_700_700.npy
xyzmod_zdet_300_-500.npy          xyzmod_zyx_-100_-100.npy
xyzmod_zdet_300_500.npy           xyzmod_zyx_-100_100.npy
xyzmod_zdet_-300_-700.npy         xyzmod_zyx_100_-100.npy
xyzmod_zdet_-300_700.npy          xyzmod_zyx_100_100.npy
xyzmod_zdet_300_-700.npy          xyzmod_zyx_-100_-300.npy
xyzmod_zdet_300_700.npy           xyzmod_zyx_-100_300.npy
xyzmod_zdet_-500_-100.npy         xyzmod_zyx_100_-300.npy
xyzmod_zdet_-500_100.npy          xyzmod_zyx_100_300.npy
xyzmod_zdet_500_-100.npy          xyzmod_zyx_-100_-500.npy
xyzmod_zdet_500_100.npy           xyzmod_zyx_-100_500.npy
xyzmod_zdet_-500_-300.npy         xyzmod_zyx_100_-500.npy
xyzmod_zdet_-500_300.npy          xyzmod_zyx_100_500.npy
xyzmod_zdet_500_-300.npy          xyzmod_zyx_-100_-700.npy
xyzmod_zdet_500_300.npy           xyzmod_zyx_-100_700.npy
xyzmod_zdet_-500_-500.npy         xyzmod_zyx_100_-700.npy
xyzmod_zdet_-500_500.npy          xyzmod_zyx_100_700.npy
xyzmod_zdet_500_-500.npy          xyzmod_zyx_-300_-100.npy
xyzmod_zdet_500_500.npy           xyzmod_zyx_-300_100.npy
xyzmod_zdet_-500_-700.npy         xyzmod_zyx_300_-100.npy
xyzmod_zdet_-500_700.npy          xyzmod_zyx_300_100.npy
xyzmod_zdet_500_-700.npy          xyzmod_zyx_-300_-300.npy
xyzmod_zdet_500_700.npy           xyzmod_zyx_-300_300.npy
xyzmod_zdet_-700_-100.npy         xyzmod_zyx_300_-300.npy
xyzmod_zdet_-700_100.npy          xyzmod_zyx_300_300.npy
xyzmod_zdet_700_-100.npy          xyzmod_zyx_-300_-500.npy
xyzmod_zdet_700_100.npy           xyzmod_zyx_-300_500.npy
xyzmod_zdet_-700_-300.npy         xyzmod_zyx_300_-500.npy
xyzmod_zdet_-700_300.npy          xyzmod_zyx_300_500.npy
xyzmod_zdet_700_-300.npy          xyzmod_zyx_-300_-700.npy
xyzmod_zdet_700_300.npy           xyzmod_zyx_-300_700.npy
xyzmod_zdet_-700_-500.npy         xyzmod_zyx_300_-700.npy
xyzmod_zdet_-700_500.npy          xyzmod_zyx_300_700.npy
xyzmod_zdet_700_-500.npy          xyzmod_zyx_-500_-100.npy
xyzmod_zdet_700_500.npy           xyzmod_zyx_-500_100.npy
xyzmod_zdet_-700_-700.npy         xyzmod_zyx_500_-100.npy
xyzmod_zdet_-700_700.npy          xyzmod_zyx_500_100.npy
xyzmod_zdet_700_-700.npy          xyzmod_zyx_-500_-300.npy
xyzmod_zdet_700_700.npy           xyzmod_zyx_-500_300.npy
xyzmod_zxy_-100_-100.npy          xyzmod_zyx_500_-300.npy
xyzmod_zxy_-100_100.npy           xyzmod_zyx_500_300.npy
xyzmod_zxy_100_-100.npy           xyzmod_zyx_-500_-500.npy
xyzmod_zxy_100_100.npy            xyzmod_zyx_-500_500.npy
xyzmod_zxy_-100_-300.npy          xyzmod_zyx_500_-500.npy
xyzmod_zxy_-100_300.npy           xyzmod_zyx_500_500.npy
xyzmod_zxy_100_-300.npy           xyzmod_zyx_-500_-700.npy
xyzmod_zxy_100_300.npy            xyzmod_zyx_-500_700.npy
xyzmod_zxy_-100_-500.npy          xyzmod_zyx_500_-700.npy
xyzmod_zxy_-100_500.npy           xyzmod_zyx_500_700.npy
xyzmod_zxy_100_-500.npy           xyzmod_zyx_-700_-100.npy
xyzmod_zxy_100_500.npy            xyzmod_zyx_-700_100.npy
xyzmod_zxy_-100_-700.npy          xyzmod_zyx_700_-100.npy
xyzmod_zxy_-100_700.npy           xyzmod_zyx_700_100.npy
xyzmod_zxy_100_-700.npy           xyzmod_zyx_-700_-300.npy
xyzmod_zxy_100_700.npy            xyzmod_zyx_-700_300.npy
xyzmod_zxy_-300_-100.npy          xyzmod_zyx_700_-300.npy
xyzmod_zxy_-300_100.npy           xyzmod_zyx_700_300.npy
xyzmod_zxy_300_-100.npy           xyzmod_zyx_-700_-500.npy
xyzmod_zxy_300_100.npy            xyzmod_zyx_-700_500.npy
xyzmod_zxy_-300_-300.npy          xyzmod_zyx_700_-500.npy
xyzmod_zxy_-300_300.npy           xyzmod_zyx_700_500.npy
xyzmod_zxy_300_-300.npy           xyzmod_zyx_-700_-700.npy
xyzmod_zxy_300_300.npy            xyzmod_zyx_-700_700.npy
xyzmod_zxy_-300_-500.npy          xyzmod_zyx_700_-700.npy
xyzmod_zxy_-300_500.npy           xyzmod_zyx_700_700.npy
xyzmod_zxy_300_-500.npy

In [ ]:
import multiprocessing, os, numpy as np
from glob import glob
import subprocess

def runMT3Dfwd(path):
    """ Worker function """
    orgDir = os.getcwd()
    print 'Starting in {:s}\n'.format(orgDir)
    
    os.chdir(path)
    print '########### \nRunning in {:s}\n###########'.format(os.getcwd())
    cmd = '/tera_raid/gudni/Codes/ZTEM_MT3Dinv_06082014/MT3Dfwd'
    # os.system(cmd)
    subprocess.call(cmd)
    os.chdir(orgDir)


if __name__ == '__main__':
    pool = multiprocessing.Pool(processes=12)
    foldList = glob('MT3Dfwd/model100w*/Freq*')
    
    pool.map(runMT3Dfwd,foldList)
    pool.cself.mesh.getEdgeInnerProductlose()
    pool.join()

In [ ]:
locs=rec2ndarr(uniLocs)

In [ ]:
%debug

In [ ]:
simpeg.mkvc(locs[0,:],2).T

In [ ]:
simpegmt.DataMT.DataMT

In [ ]:
%matplotlib qt
sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf
dataRec = mt1DdataList[0].toRecArray('Complex')
iPf.MTinteractiveMap([dataRec])

In [42]:
simpegmt.Utils.dataUtils.getAppRes(mt1DdataZdetList[0])


Out[42]:
[(979.21054686559228, 60.117584762903384),
 (742.46348350358323, 60.178574585305192),
 (575.52692078307155, 59.504482832217768),
 (448.01468147227973, 57.739339607718549),
 (380.47444523645635, 58.103945329921899),
 (329.83191047669334, 55.337196869237211),
 (289.04291928642527, 54.655448967096305),
 (266.75073534275322, 53.22914208405799),
 (233.25010687113831, 51.415310785084358),
 (235.86124553266205, 49.932520459157615),
 (204.17536410665218, 49.937485696638511),
 (216.58363015060209, 48.71243562166476),
 (213.79479472770714, 45.446839687730105),
 (194.30255638458473, 46.234641994868568),
 (198.73561291019692, 47.04819962721875),
 (206.89598365667399, 44.89666240461478),
 (213.17403683446534, 43.140057028057178)]

In [ ]: