simpeg/simpegPF/notebooks/MagInversion.ipynb

In [2]:

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
from pymatsolver import MumpsSolver
%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 [13]:

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 [14]:

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 [15]:

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

In [16]:

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

In [17]:

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

Out [17]:

(-500, 0)

In [18]:

print model.shape
print chi.shape

(33620,)
(67240,)

Step3: Generating Data

In [19]:

survey = BaseMag.BaseMagSurvey()
const = 20
Inc = 90.
Dec = 0.
Btot = 51000
survey.setBackgroundField(Inc, Dec, Btot)
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 = MumpsSolver

In [12]:

dsyn = survey.dpred(model)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-12-63f359d726f9> in <module>()
----> 1 dsyn = survey.dpred(model)

/Users/sgkang/Projects/simpeg/SimPEG/Utils/CounterUtils.pyc in wrapper(self, *args, **kwargs)
     81         counter = getattr(self,'counter',None)
     82         if type(counter) is Counter: counter.count(self.__class__.__name__+'.'+f.__name__)
---> 83         out = f(self,*args,**kwargs)
     84         return out
     85     return wrapper

/Users/sgkang/Projects/simpeg/SimPEG/Utils/codeutils.pyc in requiresVarWrapper(self, *args, **kwargs)
    221             if getattr(self, var, None) is None:
    222                 raise Exception(extra)
--> 223             return f(self,*args,**kwargs)
    224 
    225         doc = requiresVarWrapper.__doc__

/Users/sgkang/Projects/simpeg/SimPEG/Survey.pyc in dpred(self, m, u)
    307         """
    308         if u is None: u = self.prob.fields(m)
--> 309         return Utils.mkvc(self.eval(u))
    310 
    311 

/Users/sgkang/Projects/simpeg/SimPEG/Utils/CounterUtils.pyc in wrapper(self, *args, **kwargs)
     81         counter = getattr(self,'counter',None)
     82         if type(counter) is Counter: counter.count(self.__class__.__name__+'.'+f.__name__)
---> 83         out = f(self,*args,**kwargs)
     84         return out
     85     return wrapper

/Users/sgkang/Projects/simpeg/SimPEG/Survey.pyc in eval(self, u)
    320                 d_\\text{pred} = \mathbf{P} u(m)
    321         """
--> 322         raise NotImplemented('eval is not yet implemented.')
    323 
    324     @Utils.count

TypeError: 'NotImplementedType' object is not callable

In [ ]:

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 [ ]:

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 [15]:

# 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 [16]:

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 [22]:

reg = Regularization.Tikhonov(mesh, mapping = rmap)
opt = Optimization.ProjectedGNCG(maxIter = 2)
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**2)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
opt.tolG = 1e-20
opt.eps = 1e-20
reg.alpha_s = 1e-5
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 [23]:

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 -4.89e+06  1.16e+06  2.70e-01 -1.59e+05    1.45e+05      0              
   1 -4.89e+06  8.02e+05  2.15e+05 -1.05e+12    2.87e+07      0              
   2 -6.12e+05  8.00e+05  2.15e+05 -1.32e+11    3.59e+06     10   Skip BFGS  
------------------------- STOP! -------------------------
0 : |fc-fOld| = 9.2136e+11 <= tolF*(1+|f0|) = 1.5946e+04
1 : |xc-x_last| = 4.6116e-05 <= tolX*(1+|x0|) = 1.0566e+02
0 : |proj(x-g)-x|    = 3.5921e+06 <= tolG          = 1.0000e-20
0 : |proj(x-g)-x|    = 3.5921e+06 <= 1e3*eps       = 1.0000e-17
1 : maxIter   =       2    <= iter          =      2
------------------------- DONE! -------------------------

In [24]:

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

Counters:
  ProjectedGNCG.activeSet                 :        2
  ProjectedGNCG.doEndIteration            :        2
  ProjectedGNCG.doStartIteration          :        3
  ProjectedGNCG.projection                :       24
  ProjectedGNCG.scaleSearchDirection      :        2

Times:                                        mean      sum
  MagneticsDiffSecondary.Jtvec            : 5.36e-01, 5.89e+00,   11x
  MagneticsDiffSecondary.Jtvec_approx     : 6.45e-01, 5.16e+00,    8x
  MagneticsDiffSecondary.Jvec             : 6.31e-01, 5.05e+00,    8x
  MagneticsDiffSecondary.Jvec_approx      : 6.31e-01, 5.05e+00,    8x
  ProjectedGNCG.findSearchDirection       : 4.77e+00, 9.54e+00,    2x
  ProjectedGNCG.minimize                  : 2.02e+01, 2.02e+01,    1x
  ProjectedGNCG.modifySearchDirection     : 4.53e+00, 9.06e+00,    2x

In [26]:

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=2, interval=40, blit=True)

Out [26]:

– + Once Loop Reflect

In [16]:

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

In [28]:

indx = 18
iteration = 1
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)

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]: