In [105]:
import numpy as np
%pylab inline
In [106]:
def AofT(time,T, ai, taui):
return ai*np.exp(-time/taui)/(1.+np.exp(-T/(2*taui)))
In [107]:
from SimPEG import *
from simpegem1d.Waveform import CausalConv
import sys
sys.path.append("./DoubleLog/")
from plotting import mapDat
import cPickle as pickle
In [108]:
milligan = pickle.load(open("milligan.p"))
In [109]:
DAT = milligan['DAT']
Header = milligan['header']
tindst = 0
time_vtem_mus = np.r_[99,120,141,167,198,234,281,339,406,484,573,682,818,974,1151,1370,1641,1953,2307,2745,3286,3911,4620,5495,6578,7828,9245]
time = time_vtem_mus[tindst:]*1e-6
obsdata = DAT[:,13:13+27]
xyz = np.c_[DAT[:,0], DAT[:,1], DAT[:,3]]
line = milligan['Line']
In [110]:
linenum = 6
Line = np.unique(line)
indline = line == Line[linenum]
print Line[linenum]
# plt.plot(xyz[:,0], xyz[:,1], '.')
# plt.plot(xyz[indline,0], xyz[indline,1], 'r.')
In [111]:
obs_line = obsdata[indline,tindst:][::10,:]
xyz_line = xyz[indline,:][::10,:]
In [112]:
class LinearSurvey(Survey.BaseSurvey):
nD = None
xyz = None
uncert = None
def __init__(self, time, xyz,**kwargs):
self.time = time
self.xyz = xyz
self.ntime = time.size
self.ntx = self.xyz.shape[0]
self.nD = self.ntime*self.ntx
def projectFields(self, u):
return u
def setUncertainty(self, dobs, perc=0.1, floor=0.):
self.uncert = np.zeros((self.ntime, self.ntx))
self.dobs = dobs
dobs = dobs.reshape((self.ntime, self.ntx), order='F')
for itx in range(self.ntx):
ipind = dobs[:,itx]<0.
if (ipind).sum() > 3:
ip = dobs[ipind,itx]
self.uncert[:,itx] = perc*abs(dobs[:,itx]) + abs(ip).max()
else:
self.uncert[:,itx] = perc*abs(dobs[:,itx])+floor
self.uncert = Utils.mkvc(self.uncert)
return self.uncert
class LinearProblem(Problem.BaseProblem):
surveyPair = LinearSurvey
tau = None
def __init__(self, mesh,**kwargs):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
self.setTau()
#TODO: make general (not it works for db/dt)
def getG(self, currentderiv, timeconv, t0=0.):
dt = timeconv[1]-timeconv[0]
ntime = time_conv.size
meshtime = Mesh.TensorMesh([dt*np.ones(ntime)], x0=[-dt/2.])
P = meshtime.getInterpolationMat(time+t0, 'CC')
self.G = np.zeros((self.survey.ntime, self.ntau))
for j in range(self.ntau):
self.G[:,j] = P*(CausalConv(1./self.tau[j]*np.exp(-time_conv/self.tau[j]), -currentderiv, time_conv))
def setTau(self, minlogtau=-5, maxlogtau=-2, ntau = 81):
self.tau = np.logspace(minlogtau, maxlogtau, ntau)
self.ntau = ntau
def fields(self, m, u=None):
m = m.reshape((self.ntau, self.survey.ntx), order='F')
pred = self.G.dot(m)
return Utils.mkvc(pred)
def Jvec(self, m, v, u=None):
v = v.reshape((self.ntau, self.survey.ntx), order='F')
jvec = self.G.dot(v)
return Utils.mkvc(jvec)
def Jtvec(self, m, v, u=None):
v = v.reshape((self.survey.ntime, self.survey.ntx), order='F')
jtvec = self.G.T.dot(v)
return Utils.mkvc(jtvec)
In [113]:
Regularization.BaseRegularization
Out[113]:
In [ ]:
In [114]:
meshline = Mesh.TensorMesh([2])
Gx = meshline.cellGradx
In [115]:
Gx.todense()
Out[115]:
In [116]:
tind = 12
# time = np.load('./exampledata/timevtem.npy')[4:]
# obs = np.load('./exampledata/souding.npy')[4:]
wave = np.loadtxt('/Users/sgkang/Dropbox/Shared/SeogiDikun/Milligan/Data/7042_106_wform.xyz', skiprows=7)
M = 81
tau = np.logspace(-5, -2, M)
In [117]:
twave = (np.arange(10000)+1)*1e-5
indstart = 4439
indend = 6000
t0_wave = twave[indstart:indend].min()
time_conv= twave[indstart:indend]-t0_wave
In [118]:
currentderiv = wave[indstart:indend]
currentderiv[time_conv>4.4e-3] = 0.
current = CausalConv(wave[indstart:indend], np.ones_like(wave[indstart:indend]), time_conv)
In [119]:
# figsize(6, 4)
plt.plot(time_conv, wave[indstart:indend], 'b.', lw=2, ms=4)
plt.plot(time_conv, wave[indstart:indend]*0., 'k--')
plt.plot(np.r_[4.4000000e-03, 4.4000000e-03], np.r_[-4.5, 4.5], 'k:')
plt.ylim(-4.5, 4.5)
# plt.xlim(-2e-4, 5.5e-3)
Out[119]:
In [120]:
figsize(6, 4)
plt.plot(time_conv, current/current.max(), 'k-')
# plt.plot(time_conv, wave[indstart:indend]*0., 'k--')
plt.plot(np.r_[4.4000000e-03, 4.4000000e-03], np.r_[1e-6, 1e-2], 'k:')
plt.plot(time, np.zeros_like(time), 'r.')
plt.ylim(0, 1.)
Out[120]:
In [121]:
from SimPEG import Mesh
from simpegem1d.Waveform import SineFun, SineFunDeriv, CausalConv
dt = 1e-5
t0 = 4.4e-03
ntime = time_conv.size
meshtime = Mesh.TensorMesh([dt*np.ones(ntime)], x0=[-dt/2.])
P = meshtime.getInterpolationMat(time+t0, 'CC')
# time_conv = meshtime.gridN
# currentderiv = SineFunDeriv(time_conv, t0)
# current = SineFun(time_conv, t0)
In [122]:
temp = np.exp(-time_conv/1e-2)/1e-2
out = CausalConv(temp, currentderiv, time_conv)
# plt.plot(time_conv, currentderiv)
In [123]:
plt.plot(time_conv, out)
Out[123]:
In [124]:
time_conv.min(), time_conv.max()
Out[124]:
In [125]:
P.shape
Out[125]:
In [126]:
# N = time.size
# A = np.zeros((N, M))
# for j in range(M):
# A[:,j] = P*(CausalConv(1./tau[j]*np.exp(-time_conv/tau[j]), -currentderiv, time_conv))
In [210]:
meshline = Mesh.TensorMesh([xyz_line.shape[0]])
prob = LinearProblem(meshline)
survey = LinearSurvey(time, xyz_line)
prob.pair(survey)
prob.getG(currentderiv, time_conv, t0=4.4e-3)
uncert = survey.setUncertainty(obs_line.flatten(), perc=0.1, floor=1e-4)
In [211]:
mtrue = np.zeros(prob.ntau)
np.random.seed(1)
inds = np.random.random_integers(0, 41, size=5)
mtrue[inds] = np.r_[0.1, 2, 1, 4, 5]
m0 = (np.repeat(mtrue.reshape([1,-1]), survey.ntx, axis=0)).flatten()
In [212]:
from SimPEG import Maps
In [213]:
weightone = np.sqrt(np.diag(np.dot(prob.G.T, prob.G)))
weight = (np.repeat(weightone.reshape([1,-1]), survey.ntx, axis=0)).flatten()
meshreg = Mesh.TensorMesh([weight.size])
wmap = Maps.Weighting(meshreg, weights=weight)
In [214]:
plt.semilogx(prob.tau, weightone)
Out[214]:
In [215]:
print uncert.shape
print survey.dobs.shape
In [216]:
from SimPEG.Regularization import BaseRegularization as BaseRegularization
class SmoothAEMRegularization(BaseRegularization):
xyz_line = None
ntau = 81
alpha_s = 1e0
alpha_x = 1e-2
def __init__(self, mesh, mapping=None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
@property
def W(self):
"""Regularization matrix W"""
if getattr(self, '_W', None) is None:
ntx = self.xyz_line.shape[0]
meshline = Mesh.TensorMesh([ntx])
Gx = meshline.cellGradx
Wx = np.sqrt(self.alpha_x)*sp.kron(Gx, Utils.speye(self.ntau))
Ws = np.sqrt(self.alpha_s)*Utils.speye(self.ntau*ntx)
wlist = (Wx, Ws)
self._W = sp.vstack(wlist)
return self._W
In [217]:
meshline = Mesh.TensorMesh([survey.ntx])
In [218]:
meshline.nF*prob.ntau
Out[218]:
In [219]:
# test = reg.W*Utils.mkvc(mrec)
# test = test.reshape((prob.ntau, meshline.nF), order='F')
# itau = 11
# plt.plot(meshline.vectorCCx, mrec[itau,:])
# plt.plot(meshline.vectorNx, test[itau,:])
In [232]:
# reg = Regularization.BaseRegularization(meshreg, mapping=wmap)
# reg = Regularization.BaseRegularization(meshreg)
reg = SmoothAEMRegularization(meshreg, mapping=wmap)
reg.xyz_line = xyz_line
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./(uncert)
opt = Optimization.ProjectedGNCG(maxIter=200, maxIterLS = 40)
# opt = Optimization.InexactGaussNewton(maxIter=20)
opt.lower = 1e-20
# opt.upper = 1e-10
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e-12)
beta = Directives.BetaSchedule()
beta.coolingFactor = 1
target = Directives.TargetMisfit()
inv = Inversion.BaseInversion(invProb, directiveList=[beta, target, betaest])
opt.tolG = 1e-10
opt.tolG = 1e-10
m0 = np.ones_like(weight)*1e-10
reg.mref = np.zeros_like(m0)
mrec = inv.run(m0)
In [233]:
mrec = mrec.reshape((prob.ntau, survey.ntx), order='F')
In [ ]:
In [234]:
target.target
Out[234]:
In [235]:
14*234
Out[235]:
In [236]:
pred = prob.fields(mrec)
In [237]:
pred = pred.reshape((survey.ntime, survey.ntx), order='F')
obs = survey.dobs.reshape((survey.ntime, survey.ntx), order='F')
UNCERT = survey.uncert.reshape((survey.ntime, survey.ntx), order='F')
In [238]:
predmap, ticks, ticklabels = mapDat(pred, 1e-3, stretch=3)
obsmap, ticks, ticklabels = mapDat(obs, 1e-3, stretch=3)
In [239]:
itx = 80
fig = plt.figure(figsize=(10,3))
ax = plt.subplot(111)
for i in range(20):
ax.plot(xyz_line[:,0], obsmap[i,:], 'k')
ax.plot(xyz_line[:,0], predmap[i,:], 'b.')
ax.set_xlim(xyz_line[:,0].min(), xyz_line[:,0].max())
ax.set_ylim(ticks.min(), ticks.max())
ax.set_yticks(ticks)
# ax.set_ylim(0., ticks.max())
plt.plot(xyz_line[itx,0]*np.ones(2), np.r_[ticks.min(), ticks.max()],'k-', lw=2)
ax.set_xlabel("Easting (m)")
ax.set_yticklabels(ticklabels)
Out[239]:
In [240]:
figsize(10, 3)
# for itau in range(9,11,1):
itau = 31
plt.plot(xyz_line[:,0], mrec[itau,:], 'k.')
xlim(xyz_line[:,0].min(), xyz_line[:,0].max())
# ylim(1e-1, 1e5)
Out[240]:
In [229]:
# for itime in range(time.size):
# plt.semilogy(xyz_line[:,0], obs[itime,:], 'k-')
# plt.semilogy(xyz_line[:,0], pred[itime,:], 'b-')
# plt.plot(xyz_line[itx,0]*np.ones(2), np.r_[1e-8, 1e1],'k-', lw=2)
# ylim(1e-6, 2e0)
# xlim(xyz_line[:,0].min(), xyz_line[:,0].max())
In [231]:
figsize(6,5)
plt.loglog(time, pred[:,itx], 'b-')
plt.loglog(time, obs[:,itx], 'k.-')
plt.loglog(time, -obs[:,itx], 'k--')
plt.loglog(time, -(obs[:,itx]-pred[:,itx]), 'r--')
plt.loglog(time, UNCERT[:,itx], 'g-')
ylim(1e-5, 1e1)
Out[231]:
In [54]:
# for itx in range(survey.ntx):
# plt.semilogx(prob.tau, mrec[:,itx], '.')
In [55]:
hist(np.log10(tau), bins=30)
Out[55]:
In [ ]:
In [ ]:
In [56]:
# mfund = mrec.copy()
# mip = mrec.copy()
# mfund[mfund<0.] = 0.
# mip[mip>0.] = 0.
# fund = np.dot(A, mfund)
In [58]:
tind = 7
print obs[tind], invProb.dpred[tind], ip[tind]
In [59]:
fig = plt.figure(figsize=(7,4.5))
ax = plt.subplot(111)
ax.plot(time, obs, 'k.-', lw=2)
# ax.set_ylim(1e-4, 1e0)
ax.set_xscale('linear')
ax.set_yscale('linear')
ax.set_xlim(time.min(), time.max())
ax.grid(True)
In [60]:
fig = plt.figure(figsize=(7,4.5))
ax = plt.subplot(111)
ax.plot(time, obs, 'k.-', lw=2)
ax.plot(time, -obs, 'k--', lw=2)
ax.plot(time, invProb.dpred, 'b-', lw=2)
ax.plot(time, -ip, 'r.-', lw=2, ms=10)
ax.plot(time, uncert, 'g.-', lw=2, ms=10)
ax.set_ylim(1e-4, 1e0)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(time.min(), time.max())
ax.grid(True)
In [61]:
predmap, ticks, tickLabels = mapDat(invProb.dpred,1e-3, stretch=3)
In [62]:
fig = plt.figure(figsize=(5,4))
ax = plt.subplot(111)
# for i, ind in enumerate(inds):
# temp, dum, dum = mapDat(mtrue[inds][i]*np.exp(-time/tau[ind]), 1e-5, stretch=2)
# plt.semilogx(time, temp, 'k', alpha = 0.5)
outmap, ticks, tickLabels = mapDat(obs,1e-3, stretch=3)
ax.plot(time*1e6, outmap, 'k', lw=2)
ax.plot(time*1e6, predmap, 'b', lw=2)
ax.plot(time[7:]*1e6, outmap[7:]-predmap[7:], 'r', lw=2)
ax.plot(np.r_[-0.002, 0.]*1e6, np.zeros(2), 'k', lw=2)
ax.plot(np.zeros(2), np.r_[ticks.min(), ticks.max()], 'k--', lw=1)
ax.set_yticks(ticks)
ax.set_yticklabels(tickLabels)
ax.set_ylim(ticks.min(), ticks.max())
ax.plot(np.r_[-0.002, time.max()]*1e6, np.zeros(2), 'k--')
ax.set_xlim(-0.002*1e6, time.max()*1e6)
ax.set_xlabel("Time (micro-s)", fontsize = 16)
ax.set_ylabel("$db_z/dt (pV/A$-$m^4)$ ", fontsize = 16)
# ax.grid(True)
In [63]:
weight_d = np.sqrt(np.diag(np.dot(np.dot(np.diag(1./uncert), A), (np.dot(np.diag(1./uncert), A)).T)))
# weight_d = np.sqrt(np.diag( np.dot(A, A.T)))
plt.semilogx(time, weight_d)
In [ ]:
In [ ]: