In [1]:
from SimPEG import *
import simpegDCIP as DC
from numpy.polynomial import polynomial
In [2]:
%pylab inline
In [3]:
cs = 0.5
mesh = Mesh.TensorMesh([np.ones(100)*cs, np.ones(50)*cs], "CN")
x = mesh.vectorCCx
In [4]:
actx = (x>-15.)&(x<15.)
In [5]:
order = 3
Vobs = polynomial.polyvander(x, order)
dobs = 0.05*x-4
H = np.dot(Vobs.T, Vobs)
g = np.dot(Vobs.T, dobs)
mest = np.linalg.solve(H, g)
In [6]:
dpred = Vobs.dot(mest)
In [7]:
V = polynomial.polyvander(x, order)
m1D = Mesh.TensorMesh([V.shape[1]+2])
weight = (1./(V**2).sum(axis=0))**0.5
weight = weight / weight.max()
weightmap = Maps.Weighting(m1D, weights=np.r_[1., 1., weight])
m0_poly = np.r_[np.log(1e-3), np.log(1e-3), np.r_[-3., np.zeros(V.shape[1]-1)] / weight]
mtrue_poly = np.r_[np.log(1e-3), np.log(1e0), mest / weight]
In [8]:
polymap = Maps.PolyMap(mesh, V.shape[1]-1, logSigma=True, normal='Y')
mappingfwd = polymap*weightmap
mapping = Maps.ExpMap(mesh)
In [9]:
sigma = mappingfwd*mtrue_poly
In [10]:
xr = np.linspace(-15, 15, 20)
xz_A = Utils.ndgrid(xr, np.r_[-0.25])
xz_B = Utils.ndgrid(np.ones_like(xr)*22, np.r_[-0.25])
xz_M = Utils.ndgrid(xr, np.r_[-0.25])
xz_N = Utils.ndgrid(np.ones_like(xr)*-22, np.r_[-0.25])
ntx = xz_A.shape[0]
txList = []
for i in range(ntx):
offset = abs(xz_A[i,0]-xz_M[:,0])
actrx = offset > 5.
rx = DC.RxDipole(xz_M[actrx,:], xz_N[actrx,:])
src = DC.SrcDipole([rx], xz_A[i,:], xz_B[i,:])
txList.append(src)
In [11]:
dat = mesh.plotImage(np.log10(sigma), clim=(-2, 0))
plt.colorbar(dat[0])
plot(xz_A[:,0], xz_A[:,1], 'w.')
plot(xz_B[:,0], xz_B[:,1], 'y.')
plot(xz_N[:,0], xz_N[:,1], 'r.')
Out[11]:
In [12]:
dat = mesh.plotImage(np.log10(mappingfwd*m0_poly), clim=(-2, 0))
plt.colorbar(dat[0])
plot(xz_A[:,0], xz_A[:,1], 'w.')
plot(xz_B[:,0], xz_B[:,1], 'y.')
plot(xz_N[:,0], xz_N[:,1], 'r.')
Out[12]:
In [13]:
from SimPEG import SolverLU
In [14]:
mtrue = np.log(sigma)
m0 = np.log(np.ones(mesh.nC)*1e-3)
survey = DC.SurveyDC(txList)
problem = DC.ProblemDC_CC(mesh, mapping = mapping)
problem.pair(survey)
problem.Solver = SolverLU
dtrue = survey.dpred(mtrue)
d0 = survey.dpred(m0)
j = problem.fields
plot(dtrue)
plot(d0)
Out[14]:
In [15]:
abs(dtrue).min()
Out[15]:
In [16]:
hist(np.log10(abs(dtrue)), bins = 20)
Out[16]:
In [17]:
noise = 0.05*abs(dtrue)*np.random.randn(dtrue.shape[0])
In [18]:
survey.dobs = dtrue +noise
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./(0.05*abs(dtrue)+ 0.1)
reg = Regularization.Tikhonov(mesh)
opt = Optimization.InexactGaussNewton(maxIter=30,maxIterLS=20)
opt.remember('xc')
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
betaSched = Directives.BetaSchedule(coolingFactor=8, coolingRate=3)
targetmis = Directives.TargetMisfit()
savemodel = Directives.SaveModelEveryIteration()
inv = Inversion.BaseInversion(invProb, directiveList=[betaSched,targetmis])
reg.alpha_s = 1e-5
reg.mref = m0
mopt = inv.run(m0)
In [19]:
XC = opt.recall('xc')
In [20]:
from ipywidgets import interact, IntSlider
def viewInv(iteration):
fig = plt.figure(num=0,figsize = (10,5))
ax = plt.subplot(111)
dat = mesh.plotImage(np.log10(mapping*XC[iteration]), grid=True, ax=ax, clim=(-3, 0), gridOpts={'alpha':0.2}, pcolorOpts={'cmap':cm.RdPu})
# ax.set_xlim(mesh.vectorNx.min(), mesh.vectorNx.max())
# ax.set_ylim(mesh.vectorNy.min(), mesh.vectorNy.max())
ax.plot(mesh.vectorCCx, dpred, 'r--')
ax.plot(xz_A[:,0], xz_A[:,1], 'k.')
ax.plot(xz_B[:,0], xz_B[:,1], 'r.')
ax.plot(xz_N[:,0], xz_N[:,1], 'b.')
# plt.colorbar(dat[0], ax=ax)
# ax.set_ylim (-15, 0.)
# ax.set_xlim (-15, 15.)
plt.show()
return True
In [21]:
interact(viewInv, iteration = IntSlider(min=0, max=opt.iter-1,step=1, value=0))
In [22]:
plt.plot(survey.dobs)
plt.plot(invProb.dpred)
Out[22]:
In [ ]: