simpeg/simpegPF/notebooks/MagInversion.ipynb

In [5]:

from SimPEG import *
from simpegPF import BaseMag
from scipy.constants import mu_0
from simpegPF.MagAnalytics import spheremodel, CongruousMagBC
from simpegPF.Magnetics import MagneticsDiffSecondary, MagneticsDiffSecondaryInv
# import SeogiUtils as SeUtils
# import simpegEM.Utils.Solver.Mumps as Mumps
%pylab inline

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['linalg']
`%matplotlib` prevents importing * from pylab and numpy

Mag Inversion

Step1: Generating mesh

In [5]:

cs = 25.
hxind = [(cs,5,-1.3), (cs, 31),(cs,5,1.3)]
hyind = [(cs,5,-1.3), (cs, 31),(cs,5,1.3)]
hzind = [(cs,5,-1.3), (cs, 30),(cs,5,1.3)]
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')

Step2: Generating Model: Use Combo model

Here we combined \mu model$^1$, Depth model$^2$ and Active model$^3$

In [6]:

chibkg = 1e-5
chiblk = 0.1
chi = np.ones(mesh.nC)*chibkg
sph_ind = spheremodel(mesh, 0., 0., -150., 80)
chi[sph_ind] = chiblk
active = mesh.gridCC[:,2]<0
actMap = Maps.ActiveCells(mesh, active, chibkg)
dweight = np.ones(mesh.nC)
dweight[active] = (1/abs(mesh.gridCC[active, 2]-13.)**1.5)
baseMap = BaseMag.BaseMagMap(mesh)
depthMap = BaseMag.WeightMap(mesh, dweight)
dmap = baseMap*actMap
rmap = depthMap*actMap
model = (chi)[active]

In [12]:

sph_ind_ini = spheremodel(mesh, 0., 0., -200., 150)

In [8]:

chi_ini = np.ones_like(chi)*chibkg
chi_ini[sph_ind_ini] = chiblk*0.1

In [15]:

fig, ax = plt.subplots(1,1, figsize = (5, 5))
dat1 = mesh.plotSlice(chi, ax = ax, normal = 'X')
plt.colorbar(dat1[0], orientation="horizontal", ax = ax)
ax.set_ylim(-500, 0)

Out [15]:

(-500, 0)

In [16]:

print model.shape
print chi.shape

(33620L,)
(67240L,)

Step3: Generating Data

In [17]:

survey = BaseMag.BaseMagSurvey()
const = 20
Inc = 90.
Dec = 0.
Btot = 51000
survey.setBackgroundField(Inc, Dec, Btot)

In [ ]:

xr = np.linspace(-300, 300, 81)
yr = np.linspace(-300, 300, 81)
X, Y = np.meshgrid(xr, yr)
Z = np.ones((xr.size, yr.size))*(0.)
rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
survey.rxLoc = rxLoc
prob = MagneticsDiffSecondary(mesh, mapping = dmap)
prob.pair(survey)
prob.Solver = Utils.SolverUtils.SolverWrapD(Mumps, factorize=True)

In [7]:

dsyn = survey.dpred(model)

In [8]:

survey.dtrue = Utils.mkvc(dsyn)
std = 0.05
noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
survey.dobs = survey.dtrue+noise
survey.std = survey.dobs*0 + std

In [9]:

fig, ax = plt.subplots(1,2, figsize = (8,5) )
dat = ax[0].imshow(np.reshape(noise, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)])
plt.colorbar(dat, ax = ax[0], orientation="horizontal")
dat2 = ax[1].imshow(np.reshape(survey.dobs, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)])
plt.colorbar(dat2, ax = ax[1], orientation="horizontal")
plt.show()

In [36]:

# m0 = (1e-5*np.ones(mesh.nC))[active]
m0 = chi_ini[active]/dweight[active] 
dmisfit = DataMisfit.l2_DataMisfit(survey)
valmin = abs(survey.dobs).max()
dmisfit.Wd = 1/(np.ones(survey.dobs.size)*valmin)

In [37]:

d_ini = survey.dpred(m0)
fig, ax = plt.subplots(1,2, figsize = (8,5) )
dat1 = ax[0].imshow(np.reshape(d_ini, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)])
vmin = d_ini.min()
vmax = d_ini.max()
plt.colorbar(dat1, ax = ax[0], orientation="horizontal", ticks=[np.linspace(vmin, vmax, 3)], format = FormatStrFormatter('$%5.5f$'))
dat2 = ax[1].imshow(np.reshape(survey.dobs, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)])
vmin = survey.dobs.min()
vmax = survey.dobs.max()
plt.colorbar(dat2, ax = ax[1], orientation="horizontal", ticks=[np.linspace(vmin, vmax, 5)])
plt.show()

In [40]:

reg = Regularization.Tikhonov(mesh, mapping = rmap)
opt = Optimization.ProjectedGNCG(maxIter = 30)
opt.lower = 1e-10
opt.maxIterLS = 50
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
beta = Directives.BetaSchedule(coolingFactor=8, coolingRate=2)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=10**0)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
opt.tolG = 1e-20
opt.eps = 1e-20
reg.alpha_s = 1e-9
reg.alpha_x = 1.
reg.alpha_y = 1.
reg.alpha_z = 1.
prob.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.1
opt.remember('xc')

In [41]:

mopt = inv.run(m0)

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***
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0 -8.78e+04  1.08e+06  2.70e-01  1.06e+06    1.38e+05      0              
   1 -8.78e+04  3.89e+05  1.52e+09 -1.34e+14    1.57e+08      0              
   2 -1.10e+04  2.49e+05  1.52e+09 -1.67e+13    1.96e+07     15   Skip BFGS  
   3 -1.10e+04  2.49e+05  1.52e+09 -1.67e+13    1.96e+07     19   Skip BFGS  
   4 -1.37e+03  2.49e+05  1.52e+09 -2.09e+12    2.45e+06     19   Skip BFGS  

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-41-e70b86582598> in <module>()
----> 1 mopt = inv.run(m0)

/home/seogi/Documents/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

/home/seogi/Documents/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 

/home/seogi/Documents/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

/home/seogi/Documents/simpeg/SimPEG/Optimization.pyc in minimize(self, evalFunction, x0)
    179             self.printIter()
    180             if self.stoppingCriteria(): break
--> 181             self.searchDirection = self.findSearchDirection()
    182             del self.H #: Doing this saves memory, as it is not needed in the rest of the computations.
    183             p = self.scaleSearchDirection(self.searchDirection)

/home/seogi/Documents/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

/home/seogi/Documents/simpeg/SimPEG/Optimization.pyc in findSearchDirection(self)
    974 
    975                 #  Form product Hessian*pc.
--> 976                 Hp = self.H*pc
    977                 Hp = (1-Active)*Hp
    978 

/usr/local/lib/python2.7/dist-packages/scipy/sparse/linalg/interface.pyc in __mul__(self, x)
    201 
    202             if x.ndim == 1 or x.ndim == 2 and x.shape[1] == 1:
--> 203                 return self.matvec(x)
    204             elif x.ndim == 2:
    205                 return self.matmat(x)

/usr/local/lib/python2.7/dist-packages/scipy/sparse/linalg/interface.pyc in matvec(self, x)
    132             raise ValueError('dimension mismatch')
    133 
--> 134         y = self._matvec(x)
    135 
    136         if isinstance(x, np.matrix):

/home/seogi/Documents/simpeg/SimPEG/InvProblem.pyc in H_fun(v)
    132         if return_H:
    133             def H_fun(v):
--> 134                 phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u)
    135                 phi_m2Deriv = self.reg.eval2Deriv(m, v=v)
    136 

/home/seogi/Documents/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

/home/seogi/Documents/simpeg/SimPEG/DataMisfit.pyc in eval2Deriv(self, m, v, u)
    135         prob   = self.prob
    136         if u is None: u = prob.fields(m)
--> 137         return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u)

/home/seogi/Documents/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

/home/seogi/Documents/simpeg/SimPEG/Problem.pyc in Jtvec_approx(self, m, v, u)
    399             :return: JTv
    400         """
--> 401         return self.Jtvec(m, v, u)
    402 
    403     def fields(self, m):

/home/seogi/Documents/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

/home/seogi/Documents/simpegpf/simpegPF/Magnetics.pyc in Jtvec(self, m, v, u)
    328         Atemp = Utils.sdiag(self.MfMu0*B0)*(dMfMuI * (dmudm))
    329         Btemp = Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm))
--> 330         Jtv = Atemp.T*(P.T*v) - Btemp.T*(P.T*v) - Ctv
    331 
    332         return Utils.mkvc(Jtv)

/usr/local/lib/python2.7/dist-packages/scipy/sparse/base.pyc in __mul__(self, other)
    294             # Fast path for the most common case
    295             if other.shape == (N,):
--> 296                 return self._mul_vector(other)
    297             elif other.shape == (N, 1):
    298                 return self._mul_vector(other.ravel()).reshape(M, 1)

/usr/local/lib/python2.7/dist-packages/scipy/sparse/compressed.pyc in _mul_vector(self, other)
    445         # csr_matvec or csc_matvec
    446         fn = getattr(_sparsetools,self.format + '_matvec')
--> 447         fn(M, N, self.indptr, self.indices, self.data, other, result)
    448 
    449         return result

KeyboardInterrupt: 

In [20]:

opt.counter.summary()
xc = opt.recall('xc')

Counters:
  ProjectedGNCG.activeSet                 :       30
  ProjectedGNCG.doEndIteration            :       30
  ProjectedGNCG.doStartIteration          :       31
  ProjectedGNCG.projection                :      394
  ProjectedGNCG.scaleSearchDirection      :       30

Times:                                        mean      sum
  MagneticsDiffSecondary.Jtvec            : 3.10e-01, 5.64e+01,  182x
  MagneticsDiffSecondary.Jtvec_approx     : 3.48e-01, 5.26e+01,  151x
  MagneticsDiffSecondary.Jvec             : 3.48e-01, 5.26e+01,  151x
  MagneticsDiffSecondary.Jvec_approx      : 3.48e-01, 5.26e+01,  151x
  ProjectedGNCG.findSearchDirection       : 3.66e+00, 1.10e+02,   30x
  ProjectedGNCG.minimize                  : 2.57e+02, 2.57e+02,    1x
  ProjectedGNCG.modifySearchDirection     : 4.71e+00, 1.41e+02,   30x

In [21]:

from JSAnimation import IPython_display
from matplotlib import animation
from SimPEG import *

fig, ax = subplots(1,2, figsize = (16, 5))
ax[0].set_xlabel('Easting (m)')
ax[0].set_ylabel('Depth (m)')
ax[1].set_xlabel('Easting (m)')
ax[1].set_ylabel('Depth (m)')


def animate(i_id):
    indx = 18
    temp = dmap*(xc[i_id])
    minval = (temp).min()
    maxval = (temp).max()
    
    frame1 = mesh.plotSlice(temp, vType='CC', ind=indx, normal='X',ax = ax[1], grid=False, gridOpts={'color':'b','lw':0.3, 'alpha':0.5}, )
    frame2 = mesh.plotSlice(chi, vType='CC', ind=indx, normal='X',ax = ax[0], grid=False, gridOpts={'color':'b','lw':0.3, 'alpha':0.5}, );
    ax[0].set_title('True model', fontsize = 16)
    ax[1].set_title('Estimated model at iteration = ' + str(i_id+1), fontsize = 16)
    ax[0].set_ylim(-500, 0)
    ax[1].set_ylim(-500, 0)
    return frame1[0]

animation.FuncAnimation(fig, animate, frames=10, interval=40, blit=True)

Out [21]:

– + Once Loop Reflect

In [16]:

import matplotlib
matplotlib.rcParams.update({'font.size': 14, 'text.usetex': True, 'font.family': 'arial'})

In [17]:

indx = 18
iteration = 9
fig, axes = subplots(1,2, figsize = (12, 5))
vmin = chi.min()
vmax = chi.max()
ps1 = mesh.plotSlice(chi, vType='CC', ind=indx, normal='X',ax = axes[0], grid=True, gridOpts={'color':'b','lw':0.3, 'alpha':0.5});
axes[0].set_title('$\chi_{true}$', fontsize = 16)
axes[0].set_ylim(-500, 0.)
cb1 = colorbar(ps1[0], ax =  axes[0], orientation="horizontal", ticks=[np.linspace(vmin, vmax, 5)], format = FormatStrFormatter('$%5.3f$'))
axes[0].set_xlabel('Easting (m)')
axes[0].set_ylabel('Depth (m)')

vmin = (actMap*xc[iteration]).min()
vmax = (actMap*xc[iteration]).max()
ps2 = mesh.plotSlice(actMap*xc[iteration], vType='CC', ind=indx, normal='X', ax = axes[1], grid=True, gridOpts={'color':'b','lw':0.3, 'alpha':0.5});
axes[1].set_title('$\chi_{pred}$', fontsize = 16)
axes[1].set_ylim(-500, 0.)
cb2 = colorbar(ps2[0], ax =  axes[1], orientation="horizontal", ticks=[np.linspace(vmin, vmax, 5)], format = FormatStrFormatter('$%5.3f$'))
cb1.set_label('Susceptibility (dimensionless)')
cb2.set_label('Susceptibility (dimensionless)')
axes[1].set_xlabel('Easting (m)')
axes[1].set_ylabel('Depth (m)')
fig.savefig('model.png', dpi = 200)

/usr/local/lib/python2.7/dist-packages/matplotlib/lines.py:503: RuntimeWarning: invalid value encountered in greater_equal
  return np.alltrue(x[1:] - x[0:-1] >= 0)

In [18]:

dpred_xc = survey.dpred(xc[iteration])
fig, ax = plt.subplots(1,2, figsize = (12,7) )
vmin = survey.dobs.min()
vmax = survey.dobs.max()
dat2 = ax[0].imshow(np.reshape(survey.dobs, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)], vmin = vmin, vmax = vmax)
cb1 = plt.colorbar(dat2, ax = ax[0], orientation="horizontal", ticks=[np.linspace(vmin, vmax, 5)])
dat = ax[1].imshow(np.reshape(dpred_xc, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)], vmin = vmin, vmax = vmax)
cb2 = plt.colorbar(dat, ax = ax[1], orientation="horizontal", ticks=[np.linspace(vmin, vmax, 5)])
ax[0].plot(rxLoc[:,0],rxLoc[:,1],'w.', ms=1)
ax[1].plot(rxLoc[:,0],rxLoc[:,1],'w.', ms=1)
ax[0].set_title('Observed', fontsize = 16)
ax[1].set_title('Predicted', fontsize = 16)
ax[0].set_xlabel('Easting (m)')
ax[0].set_ylabel('Northing (m)')
ax[1].set_xlabel('Easting (m)')
ax[1].set_ylabel('Northing (m)')
cb1.set_label('Total magnetic intensity (nT)')
cb2.set_label('Total magnetic intensity (nT)')
fig.savefig('obspred.png', dpi = 200)

In [18]: