notebook.community

Edit and run 





In [1]:



    


import numpy as np
%pylab inline



    


















    






Populating the interactive namespace from numpy and matplotlib

$A(t,T) = \Sigma_i A_i e^{-t/\tau_i} / (1 + e^{-T/2\tau_i})$





In [2]:



    


def AofT(time,T, ai, taui):
    return ai*np.exp(-time/taui)/(1.+np.exp(-T/(2*taui)))



    











In [3]:



    


from SimPEG import *
import sys
sys.path.append("./DoubleLog/")
from plotting import mapDat



    











In [4]:



    


class LinearSurvey(Survey.BaseSurvey):
    nD = None
    def __init__(self, time, **kwargs): 
        self.time = time
        self.nD = time.size
        
    def projectFields(self, u):
        return u

class LinearProblem(Problem.BaseProblem):

    surveyPair = LinearSurvey

    def __init__(self, mesh, G, **kwargs):
        Problem.BaseProblem.__init__(self, mesh, **kwargs)
        self.G = G

    def fields(self, m, u=None):
        return self.G.dot(m)

    def Jvec(self, m, v, u=None):
        return self.G.dot(v)

    def Jtvec(self, m, v, u=None):
        return self.G.T.dot(v)

Simple exponential basis

$$ \mathbf{A}\mathbf{\alpha} = \mathbf{d}$$ 




In [11]:



    


# time = np.cumsum(np.r_[0., 1e-5*np.ones(10), 5e-5*np.ones(10), 1e-4*np.ones(10), 5e-4*np.ones(10), 1e-3*np.ones(10)])
time = np.cumsum(np.r_[0., 1e-5*np.ones(10), 5e-5*np.ones(10),1e-4*np.ones(10), 5e-4*np.ones(10)])
M = 41
tau = np.logspace(-5.1, -1, M)



    











In [12]:



    


N = time.size
A = np.zeros((N, M))
for j in range(M):
    A[:,j] = np.exp(-time/tau[j])//tau[j]



    











In [13]:



    


mtrue = np.zeros(M)



    











In [14]:



    


np.random.seed(1)
inds = np.random.random_integers(0, 41, size=5)
mtrue[inds] = np.r_[-3, 2, 1, 4, 5]



    











In [15]:



    


out = np.dot(A,mtrue)



    











In [16]:



    


fig = plt.figure(figsize=(6,4.5))
ax = plt.subplot(111)

for i, ind in enumerate(inds):
    temp, dum, dum = mapDat(mtrue[inds][i]*np.exp(-time/tau[ind])/tau[j], 1e-1, stretch=2)
    plt.semilogx(time, temp, 'k', alpha = 0.5)    
outmap, ticks, tickLabels = mapDat(out, 1e-1, stretch=2)    
ax.semilogx(time, outmap, 'k', lw=2)
ax.set_yticks(ticks)
ax.set_yticklabels(tickLabels)
# ax.set_ylim(ticks.min(), ticks.max())
ax.set_ylim(ticks.min(), ticks.max())
ax.set_xlim(time.min(), time.max())
ax.grid(True)



    


















    
























In [159]:



    


# from pymatsolver import MumpsSolver



    











In [168]:



    


abs(survey.dobs).min()



    


















    
Out[168]:








6.0558934897099022

















In [169]:



    


mesh = Mesh.TensorMesh([M])
prob = LinearProblem(mesh, A)
survey = LinearSurvey(time)
survey.pair(prob)
survey.makeSyntheticData(mtrue, std=0.01)
# survey.dobs = out
reg = Regularization.BaseRegularization(mesh)
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./(0.05*abs(survey.dobs)+0.05*100.)
opt = Optimization.ProjectedGNCG(maxIter=20)
# opt = Optimization.InexactGaussNewton(maxIter=20)
opt.lower = -1e-10
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
invProb.beta = 1e-4
beta = Directives.BetaSchedule()
beta.coolingFactor = 2
target = Directives.TargetMisfit()
inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
m0 = np.zeros_like(survey.mtrue)
mrec = 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  1.00e-04  4.62e+03  0.00e+00  4.62e+03    4.12e+03      0              
   1  1.00e-04  4.15e+03  9.34e-02  4.15e+03    1.76e+03      3              
   2  1.00e-04  3.74e+03  4.20e+00  3.74e+03    6.75e+01      0   Skip BFGS  
   3  5.00e-05  2.39e+03  4.12e+00  2.39e+03    9.07e+01      0              
   4  5.00e-05  2.16e+03  5.14e+00  2.16e+03    5.20e+01      2              
   5  5.00e-05  2.11e+03  1.45e+01  2.11e+03    5.38e+00      0   Skip BFGS  
   6  2.50e-05  8.34e+02  9.78e+00  8.34e+02    9.15e+01      0              
   7  2.50e-05  8.16e+02  1.10e+01  8.16e+02    2.36e+01      3              
   8  2.50e-05  7.24e+02  2.30e+01  7.24e+02    2.01e+00      1              
   9  1.25e-05  7.07e+02  3.14e+01  7.07e+02    5.63e+00      2   Skip BFGS  
  10  1.25e-05  7.02e+02  3.00e+01  7.02e+02    1.49e+00      0   Skip BFGS  
  11  1.25e-05  7.01e+02  3.35e+01  7.01e+02    8.18e+00      2              
  12  6.25e-06  7.01e+02  3.29e+01  7.01e+02    2.46e-01      0              
  13  6.25e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      0              
  14  6.25e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      1              
  15  3.13e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      1   Skip BFGS  
  16  3.13e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      0              
  17  3.13e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      2   Skip BFGS  
  18  1.56e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      1   Skip BFGS  
  19  1.56e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      0              
  20  1.56e-06  7.01e+02  3.28e+01  7.01e+02    1.66e-01      2   Skip BFGS  
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 4.6234e+02
1 : |xc-x_last| = 0.0000e+00 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 1.6602e-01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.6602e-01 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      20    <= iter          =     20
------------------------- DONE! -------------------------

















In [170]:



    


plt.semilogx(tau, mtrue, '.')
plt.semilogx(tau, mrec, '.')



    


















    
Out[170]:








[<matplotlib.lines.Line2D at 0x10ada9510>]










    
























In [171]:



    


fig = plt.figure(figsize=(6,4.5))
ax = plt.subplot(111)
obsmap, ticks, tickLabels = mapDat(survey.dobs, 1e0, stretch=2)    
predmap, dum, dum = mapDat(invProb.dpred, 1e0, stretch=2)    
ax.loglog(time, survey.dobs, 'k', lw=2)
ax.loglog(time, invProb.dpred, 'k.', lw=2)
# ax.set_yticks(ticks)
# ax.set_yticklabels(tickLabels)
# ax.set_ylim(ticks.min(), ticks.max())
# ax.set_ylim(ticks.min(), ticks.max())
ax.set_xlim(time.min(), time.max())
ax.grid(True)



    


















    
























In [172]:



    


time = np.cumsum(np.r_[0., 1e-5*np.ones(10), 5e-5*np.ones(10), 1e-4*np.ones(10), 5e-4*np.ones(10), 1e-3*np.ones(10)])
N = time.size
A = np.zeros((N, M))
for j in range(M):
    A[:,j] = np.exp(-time/tau[j]) /tau[j]



    











In [173]:



    


mfund = mtrue.copy()
mfund[mfund<0.] = 0.
obs = np.dot(A, mtrue)
fund = np.dot(A, mfund)
pred = np.dot(A, mrec)



    











In [174]:



    


ip = obs-fund
ipobs = obs-pred



    











In [175]:



    


plt.loglog(time, obs, 'k.-', lw=2)
plt.loglog(time, -obs, 'k--', lw=2)
plt.loglog(time, fund, 'b.', lw=2)
plt.loglog(time, pred, 'b-', lw=2)
plt.loglog(time, -ip, 'r--', lw=2)
plt.loglog(time, abs(ipobs), 'r.', lw=2)
plt.ylim(abs(obs).min(), abs(obs).max())
plt.xlim(time.min(), time.max())



    


















    
Out[175]:








(1.0000000000000001e-05, 0.016600000000000007)

Range of tau is really important to fit ...





In [ ]:

Content source: sgkang/simpegAIP

Similar notebooks:

notebook.community | gallery | about