In [1]:

    

import numpy as np
data_directory = "../testDataIP/"
import matplotlib.pyplot as plt


    

In [2]:

    

time_gates = "30:50,50:70,70:100,100:130,130:170,170:210,210:260,260:320,320:390,390:470,470:560,560:660,660:770,770:890"
data = "0.780       0.626       0.522       0.456       0.401       0.347       0.309       0.273       0.236       0.208       0.180       0.155       0.135       0.124"
Vp = 24.388
time_gates = np.array([np.array(gate.split(":"), dtype=float) for gate in time_gates.split(',')])
data = np.array(data.split(), dtype=float)
times = np.sqrt(time_gates[:,0] * time_gates[:,1]) * 1e-3


    

In [3]:

    

from simpegEM1D import DigFilter
def ColeColeSeigel(f, sigmaInf, eta, tau, c):
    w = 2*np.pi*f
    return eta/(1 + (1j*w*tau)**c)


    

In [4]:

    

wt, tbase, omega_int = DigFilter.setFrequency(times)
frequency = omega_int / (2*np.pi)
siginf = 1.
sigma = ColeColeSeigel(frequency, siginf, 0.06, 0.06, 0.6)
sigTCole, _ = DigFilter.transFilt(sigma, wt, tbase, omega_int, times, tol=1e-12)


    

In [5]:

    

plt.loglog(times*1e3, data/Vp * 1e3, '.-')
plt.loglog(times*1e3, sigTCole * 1e3, '.-')
plt.xlabel("time (ms)")
plt.ylabel("Voltage (mV/V)")


    

    
Out[5]:





Text(0,0.5,'Voltage (mV/V)')






    

In [21]:

    

from simpegEMIP.StretchedExponential import SEInvProblem, SESurvey
from SimPEG import *
def fit_with_se(time, dobs, eta0=0.01, tau0=0.1, c0=0.5):
    siginf = 1.
    wires = Maps.Wires(('eta', 1), ('tau', 1), ('c', 1))
    taumap = Maps.ExpMap(nP=1)*wires.tau
    etamap = Maps.ExpMap(nP=1)*wires.eta
    cmap = Maps.ExpMap(nP=1)*wires.c    
    survey = SESurvey()
    survey.dobs = dobs
    m1D = Mesh.TensorMesh([np.ones(3)])
    prob = SEInvProblem(m1D, etaMap = etamap, tauMap = taumap, cMap=cmap)
    prob.time = time
    prob.pair(survey)
    eta0 = dobs[1]
    m0 = np.log(np.r_[eta0, tau0, c0])
    perc = 0.05
    dmisfitpeta = DataMisfit.l2_DataMisfit(survey)
    dmisfitpeta.W = 1/(abs(survey.dobs)*perc)
    reg = Regularization.Simple(m1D)
    opt = Optimization.ProjectedGNCG(maxIter = 20)
    invProb = InvProblem.BaseInvProblem(dmisfitpeta, reg, opt)
    target = Directives.TargetMisfit()
    betaSch = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
    opt.upper = np.log(np.r_[1., 10., 1.])
    opt.lower = np.log(np.r_[1e-5, 0.001, 1e-2])
    
    invProb.beta=0.
    inv = Inversion.BaseInversion(invProb, directiveList=[target])
#     reg.mref = 0.*m0
    prob.counter = opt.counter = Utils.Counter()
    opt.LSshorten = 0.5
    opt.remember('xc')
    opt.tolX = 1e-20
    opt.tolF = 1e-20
    opt.tolG = 1e-20
    opt.eps = 1e-20
    mopt = inv.run(m0)   
    return np.exp(mopt), survey.dobs, invProb.dpred

mopt, dobs, dpred = fit_with_se(times, data/Vp)


    

    




SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
SimPEG.InvProblem will set Regularization.mref to m0.

    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
    ***Done using same Solver and solverOpts as the problem***
model has any nan: 0
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  0.00e+00  9.90e+02  2.22e+00  9.90e+02    7.07e+00      0              
   1  0.00e+00  7.13e+02  1.36e+00  7.13e+02    1.01e+01      0              
   2  0.00e+00  5.54e+01  1.41e+00  5.54e+01    9.79e+00      0              
   3  0.00e+00  4.46e+01  1.67e+00  4.46e+01    9.72e+00      2              
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 9.9086e-18
0 : |xc-x_last| = 2.7431e-01 <= tolX*(1+|x0|) = 5.3814e-20
0 : |proj(x-g)-x|    = 9.8067e+00 <= tolG          = 1.0000e-20
0 : |proj(x-g)-x|    = 9.8067e+00 <= 1e3*eps       = 1.0000e-17
0 : maxIter   =      20    <= iter          =      4
------------------------- DONE! -------------------------

In [29]:

    

plt.semilogx(times*1e3, dobs * 1e3, '-')
plt.semilogx(times*1e3, dpred * 1e3, 'x')
plt.xlabel("time (ms)")
plt.ylabel("Voltage (mV/V)")
plt.text(
    40, 6, ("eta: %.1e mV/V, tau: %.1e s, c: %.1e") % (mopt[0]*1e3, mopt[1], mopt[2]),
    fontsize=14
)


    

    
Out[29]:





Text(40,6,'eta: 6.0e+01 mV/V, tau: 8.6e-02 s, c: 4.1e-01')






    

In [23]:

    

current = np.loadtxt(data_directory+"P0-R110-In.xyz")
dt = 1./150.
plt.plot(current[:,0]*dt, current[:,1])


    

    
Out[23]:





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






    

In [24]:

    

plt.semilogx((current[150:,0]*dt-1)*1e3, current[150:,1])
plt.xlabel("Time (ms)")
plt.ylabel("Current (mA)")


    

    
Out[24]:





Text(0,0.5,'Current (mA)')






    

In [25]:

    

plt.semilogx(current[:150,0]*dt*1e3, current[:150,1])
plt.xlabel("Time (s)")
plt.ylabel("Current (mA)")


    

    
Out[25]:





Text(0,0.5,'Current (mA)')