In [97]:
from SimPEG import *
%pylab inline
In [98]:
import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 100
In [99]:
cs, ncx, ncy, ncz, npad = 20., 30, 20, 30, 12
hx = [(cs,npad,-1.4), (cs,ncx), (cs,npad,1.4)]
hy = [(cs,npad,-1.4), (cs,ncy), (cs,npad,1.4)]
hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
print ("Padding distance x: %10.5f m") % (np.sum(mesh.hx[:npad]))
print ("Padding distance z: %10.5f m") % (np.sum(mesh.hz[:npad]))
print ("Min dx: %10.5f m") % (mesh.hx.min())
print ("Min dz: %10.5f m") % (mesh.hz.min())
In [100]:
mesh1D = Mesh.TensorMesh([5], 'C')
active1D = mesh1D.gridCC < 0.2
actmap1d = Maps.ActiveCellsTopo(mesh1D, active1D)
In [101]:
print mesh1D.gridCC
print active1D
In [102]:
reg = Regularization.Tikhonov(mesh1D, mapping=actmap1d)
In [103]:
print reg.Wx.todense()
Wxactive = reg.Wx[:,active1D]
print Wxactive.todense()
In [104]:
m_temp = np.ones(5)*2
m_temp[active1D] = 2
In [105]:
Wxx_active = Wxactive.T*Wxactive
Wxx = reg.Wx.T*reg.Wx
In [106]:
# print m_temp
print (Wxx*m_temp)[active1D]
print (Wxx_active*m_temp[active1D])
In [106]:
In [107]:
print Wxx.todense()
In [108]:
print Wxx_active.todense()
In [109]:
mesh1D.gridCC.shape
Out[109]:
In [110]:
sig1D = np.random.randn(mesh1D.nC)
In [111]:
mesh1D.plotImage(actmap1d*sig1D[active1D])
Out[111]:
In [112]:
mesh2D = Mesh.TensorMesh([hx,hz], 'CC')
active2D = mesh2D.gridCC[:,1] < 0.
actmap2d = Maps.ActiveCellsTopo(mesh2D, active2D)
In [113]:
sig2D = Utils.mkvc(Utils.ModelBuilder.randomModel((mesh2D.nCx, mesh2D.nCy)))
In [114]:
mesh2D.plotImage(actmap2d*sig2D[active2D])
Out[114]:
In [115]:
cs, ncx, ncy, ncz, npad = 20., 30, 20, 15, 12
hx = [(cs,npad,-1.4), (cs,ncx), (cs,npad,1.4)]
hy = [(cs,npad,-1.4), (cs,ncy), (cs,npad,1.4)]
hz = [(cs,npad,-1.4), (cs,ncz)]
meshact = Mesh.TensorMesh([hx,hy,hz], 'CCN')
In [116]:
mopt = np.load('./inv3D_realistic_ref2D/model_5.npy')
active = mesh.gridCC[:,2] < 0.
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nC)
actMapreg = Maps.ActiveCellsTopo(mesh, active, nC=mesh.nC)
mapping = Maps.ExpMap(mesh) * actMap
regmapping = Maps.ExpMap(mesh) * actMapreg
regmap = Maps.IdentityMap(meshact)
In [117]:
mesh.dim
Out[117]:
In [118]:
dat1 = mesh.plotSlice(np.log10(regmapping*mopt), ind = 18, normal='Y')
plt.colorbar(dat1[0])
plt.ylim(-1000, 100)
# plt.xlim(-500, 500)
Out[118]:
In [119]:
import simpegEM as EM
In [120]:
for maps in mapping.maps:
print maps
In [121]:
temp = map(lambda m: isinstance(m,Maps.ActiveCells),mapping.maps)
In [122]:
print len(map(bool, temp))
In [123]:
reduce
Out[123]:
In [124]:
print mapping.maps[int(np.where(temp)[0])]
In [125]:
print mapping
In [126]:
x1 = np.arange(30)*10 - 300.
y1 = np.arange(30)*10 - 150.
xyz1 = Utils.ndgrid(x1, y1, np.r_[0.])
x2 = np.arange(30)*10 + 10.
y2 = np.arange(30)*10 - 150.
xyz2 = Utils.ndgrid(x2, y2, np.r_[0.])
In [127]:
ntx = 2
nrx1 = xyz1.shape[0]
time = np.logspace(-4, -2, 31)
In [128]:
rx1 = EM.TDEM.RxTDEM(xyz1, time, 'bz')
tx1 = EM.TDEM.TxTDEM(np.array([0., -150., 0.]), 'CircularLoop_MVP', [rx1])
tx1.radius = 250.
rx2 = EM.TDEM.RxTDEM(xyz2, time, 'bz')
tx2 = EM.TDEM.TxTDEM(np.array([0., 150., 0.]), 'CircularLoop_MVP', [rx2])
tx2.radius = 250.
In [129]:
# survey = EM.TDEM.SurveyTDEM([tx1, tx2])
survey = EM.TDEM.SurveyTDEM([tx1])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping, verbose=True)
# prb.solver = MumpsSolver
# prb.solverOpts = {"symmetric":True}
# prb.timeSteps = [(1e-4/15, 10), (1e-3/15, 10), (1e-2/15, 5)]
prb.timeSteps = [(1e-4/15, 10)]
if prb.ispaired:
prb.unpair()
if survey.ispaired:
survey.unpair()
prb.pair(survey)
In [130]:
dpred = np.load('bz_realistic.npy')
In [131]:
noise = abs(dpred)*np.random.randn(dpred.size)*0.05
dobs = dpred+noise
survey.dobs = dpred
std = 0.05
In [132]:
Imap = Maps.IdentityMap(meshact)
# print Imap.deriv(0.).shape
In [133]:
sigma = mapping*mopt
In [134]:
sigtest = np.load('invTEM2D.npy')
sighalf = np.ones_like(sigtest)*1e-3
m0 = np.log(sigtest[active])
mref = np.log(sighalf[active])
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1/(abs(survey.dobs)*std)
opt = Optimization.InexactGaussNewton(maxIter = 15)
reg = Regularization.Tikhonov(mesh, mapping=mapping)
# reg.active = active
#betaest = Directives.BetaEstimate_ByEig(beta0_ratio = 1e2)
beta = Directives.BetaSchedule(coolingFactor = 8, coolingRate = 3)
invprb = InvProblem.BaseInvProblem(dmis, reg, opt)
# inv = Inversion.BaseInversion(invprb, directiveList = [beta,RememberXC()])
#inv = Inversion.BaseInversion(invprb, directiveList = [beta, ])
reg.alpha_s = 1e-9
reg.alpha_x = 1.
reg.alpha_y = 1.
reg.alpha_z = 1.
reg.active = False
reg.mref = mref
C = Utils.Counter()
prb.counter = C
opt.counter = C
opt.LSshorten = 0.5
In [135]:
print reg.W.shape
print reg.Wx.shape
In [135]:
In [136]:
Wactive = reg.W[:,active]
In [137]:
Wactive.shape
Out[137]:
In [138]:
active.shape
Out[138]:
In [139]:
test = reg.evalDeriv(mref)
test_1 = reg.evalDeriv(m0)
In [140]:
print reg.eval(m0)
print reg.eval(mref)
In [141]:
def ExtendTopo(mesh, sigma, active):
act_temp = active.reshape((mesh.nCx*mesh.nCy, mesh.nCz), order = 'F')
val_temp = sigma.reshape((mesh.nCx*mesh.nCy, mesh.nCz), order = 'F')
out = np.zeros(mesh.nC)
z_temp = mesh.gridCC[:,2].reshape((mesh.nCx*mesh.nCy, mesh.nCz), order = 'F')
for i in range(mesh.nCx*mesh.nCy):
act_tempxy = act_temp[i,:] == 1
val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])]
out[~active] = Utils.mkvc(val_temp)[~active]
return out
In [142]:
plt.plot(actMapreg*test)
Out[142]:
In [143]:
dat1 = mesh.plotSlice(actMapreg*test, ind = 4, normal='X')
plt.colorbar(dat1[0])
# plt.ylim(-500, 0)
# plt.xlim(-500, 500)
Out[143]:
In [144]:
print test.min()
print test.max()
print test_1.min()
print test_1.max()
In [145]:
dum = actMap*test
print dum[62000]
print mesh.gridCC[62000, :]
print mesh.vectorCCz[26]
print mesh.vectorCCy[4]
In [146]:
dat1 = meshact.plotSlice(np.log10(np.exp(m0)), ind = 4, normal='X')
# dat1 = mesh.plotSlice(actMap*m0, ind = 18)
plt.colorbar(dat1[0])
plt.ylim(-500, 0)
Out[146]:
In [146]:
In [147]:
dat1 = mesh.plotSlice(actMapreg*mopt, ind = 14)
plt.colorbar(dat1[0])
Out[147]:
In [174]:
dest = np.load('inv3D_realistic_ref2D/dpred_5.npy')
In [175]:
Dpred = dobs.reshape((900, 31, 2), order='F')
Dpred1 = Dpred[:,:,0]
Dpred2 = Dpred[:,:,1]
Dest = dest.reshape((900, 31, 2), order='F')
Dest1 = Dest[:,:,0]
Dest2 = Dest[:,:,1]
In [50]:
itime = 28
fig, ax = plt.subplots(1,2, figsize = (12, 7))
vmin = Utils.mkvc(Dpred[:,itime,0]).min()
vmax = Utils.mkvc(Dpred[:,itime,0]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'), Dpred[:,itime,0].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'), Dest[:,itime,0].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
for i in range(2):
ax[i].set_xlabel('Easting (m)', fontsize = 16)
ax[i].set_ylabel('Northing (m)', fontsize = 16)
cb = plt.colorbar(dat1, ax=ax[i], orientation='horizontal')
cb.set_label('Magnetic flux density (Wb/m$^2$)', fontsize = 14)
ax[0].set_title('Observed data from Tx#1', fontsize = 16)
ax[1].set_title('Predicted data from Tx#1', fontsize = 16)
fig.savefig('./figures/obspredTD_7_3mstx1.png')
In [51]:
itime = 9
fig, ax = plt.subplots(1,2, figsize = (12, 7))
vmin = Utils.mkvc(Dpred[:,itime,1]).min()
vmax = Utils.mkvc(Dpred[:,itime,1]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'), Dpred[:,itime,1].reshape((30, 30), order='F'), 30)
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'), Dest[:,itime,1].reshape((30, 30), order='F'), 30)
for i in range(2):
ax[i].set_xlabel('Easting (m)', fontsize = 16)
ax[i].set_ylabel('Northing (m)', fontsize = 16)
cb = plt.colorbar(dat1, ax=ax[i], orientation='horizontal')
cb.set_label('Magnetic flux density (Wb/m$^2$)', fontsize = 14)
ax[0].set_title('Observed data from Tx#2', fontsize = 16)
ax[1].set_title('Predicted data from Tx#2', fontsize = 16)
fig.savefig('./figures/obspredTD_7_3mstx2.png')
In [32]:
rxind = 20
plt.loglog(time, Dpred[rxind,:,0])
plt.loglog(time, Dest[rxind,:,0])
Out[32]:
In [33]:
itime = 0
fig, ax = plt.subplots(1,2, figsize = (10, 4))
ax[0].hist(Dpred1[:,itime])
ax[1].hist(Dpred2[:,itime])
Out[33]:
In [35]:
# from JSAnimation import IPython_display
# from matplotlib import animation
# fig, ax = plt.subplots(1,2, figsize = (12, 5))
# for i in range(2):
# ax[i].set_xlabel('x (m)', fontsize = 16)
# ax[i].set_ylabel('y (m)', fontsize = 16)
# def animate(itime):
# frame1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'), Dpred[:,itime,0].reshape((30, 30), order='F'), 30)
# frame2 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'), Dpred[:,itime,1].reshape((30, 30), order='F'), 30)
# # cb1 = plt.colorbar(frame1, ax = ax[0])
# # cb2 = plt.colorbar(frame2, ax = ax[1])
# return frame1, frame2
# animation.FuncAnimation(fig, animate, frames=31, interval=40, blit=True)
In [ ]:
# np.save('bzobs_realistic', dobs)
# dmis = DataMisfit.l2_DataMisfit(survey)
In [ ]: