wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-29 23:25:49 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #217 from simpeg/em/inversioncleanup

Browse Source 
Em/inversioncleanup
This commit is contained in:
Lindsey
2016-02-01 19:23:39 -08:00
parent 8da4d78bde 1457e51fc1
commit 28a8a2ea80
12 changed files with 294 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DataMisfit.py
+11 -17
View File
@@ -59,20 +59,6 @@ class BaseDataMisfit(object):
"""
raise NotImplementedError('This method should be overwritten.')
# TODO: implement target misfit as a property, or possibly as an inversion directive.
# def target(self, forward):
# """target(forward)
# Target for data misfit. By default this is the number of data,
# which satisfies the Discrepancy Principle.
# :rtype: float
# :return: data misfit target
# """
# prob, survey = self.splitForward(forward)
# return survey.nD
class l2_DataMisfit(BaseDataMisfit):
@@ -103,10 +89,18 @@ class l2_DataMisfit(BaseDataMisfit):
"""
if getattr(self, '_Wd', None) is None:
print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.'
survey = self.survey
eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
if getattr(survey,'std', None) is None:
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
survey.std = 0.05
if getattr(survey, 'eps', None) is None:
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
return self._Wd
@Wd.setter
SimPEG/EM/FDEM/FDEM.py
+16 -15
View File
@@ -53,13 +53,13 @@ class BaseFDEMProblem(BaseEMProblem):
Ainv.clean()
return F
def Jvec(self, m, v, f=None):
def Jvec(self, m, v, u=None):
"""
Sensitivity times a vector
"""
if f is None:
f = self.fields(m)
if u is None:
u = self.fields(m)
self.curModel = m
@@ -71,34 +71,34 @@ class BaseFDEMProblem(BaseEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = f[src, ftype]
u_src = u[src, ftype]
dA_dm = self.getADeriv_m(freq, u_src, v)
dRHS_dm = self.getRHSDeriv_m(freq, src, v)
du_dm = Ainv * ( - dA_dm + dRHS_dm )
for rx in src.rxList:
df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
Jv[src, rx] = P(Df_Dm)
Ainv.clean()
return Utils.mkvc(Jv)
def Jtvec(self, m, v, f=None):
def Jtvec(self, m, v, u=None):
"""
Sensitivity transpose times a vector
"""
if f is None:
f = self.fields(m)
if u is None:
u = self.fields(m)
self.curModel = m
@@ -114,12 +114,12 @@ class BaseFDEMProblem(BaseEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = f[src, ftype]
u_src = u[src, ftype]
for rx in src.rxList:
PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_duT = df_duTFun(src, PTv, adjoint=True)
ATinvdf_duT = ATinv * df_duT
@@ -128,7 +128,7 @@ class BaseFDEMProblem(BaseEMProblem):
dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
du_dmT = -dA_dmT + dRHS_dmT
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
dfT_dm = df_dmFun(src, PTv, adjoint=True)
du_dmT += dfT_dm
@@ -142,7 +142,8 @@ class BaseFDEMProblem(BaseEMProblem):
raise Exception('Must be real or imag')
ATinv.clean()
return Jtv
return Utils.mkvc(Jtv)
def getSourceTerm(self, freq):
"""
SimPEG/EM/__init__.py
+1 -1
View File
@@ -1,6 +1,6 @@
# from EM import *
import TDEM
import FDEM
import Base
import Analytics
import Utils
from scipy.constants import mu_0, epsilon_0
SimPEG/Examples/EM_FDEM_1D_Inversion.py
+38 -29
View File
@@ -1,6 +1,6 @@
from SimPEG import *
import SimPEG.EM as EM
from scipy.constants import mu_0
from SimPEG.EM import mu_0
def run(plotIt=True):
@@ -17,54 +17,62 @@ def run(plotIt=True):
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')
layerz = -100.
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-100.)
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
sig_half = 2e-3
sig_half = 2e-2
sig_air = 1e-8
sig_layer = 1e-3
sig_layer = 1e-2
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
sigma[layer] = sig_layer
mtrue = np.log(sigma[active])
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
ax.set_ylim(-600, 0)
ax.set_xlim(1e-4, 1e-2)
ax.set_ylim(-500, 0)
ax.set_xlim(1e-3, 1e-1)
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
ax.set_ylabel('Depth (m)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
rxOffset=1e-3
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
survey = EM.TDEM.SurveyTDEM([src])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
rxOffset=10.
bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
freqs = np.logspace(1,3,10)
srcLoc = np.array([0., 0., 10.])
srcList = []
[srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
survey = EM.FDEM.Survey(srcList)
prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
prb.Solver = SolverLU
prb.Solver = SolverLU
prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
prb.pair(survey)
dtrue = survey.dpred(mtrue)
survey.dtrue = dtrue
std = 0.05
noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
survey.dobs = survey.dtrue+noise
survey.std = survey.dobs*0 + std
survey.Wd = 1/(abs(survey.dobs)*std)
survey.makeSyntheticData(mtrue, std)
survey.std = std
survey.eps = np.linalg.norm(survey.dtrue)*1e-5
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (10, 6))
ax.loglog(rx.times, dtrue, 'b.-')
ax.loglog(rx.times, survey.dobs, 'r.-')
fig, ax = plt.subplots(1,1, figsize = (6, 6))
ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-')
ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-')
ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
ax.set_xlabel('Time (s)', fontsize = 14)
ax.set_ylabel('$B_z$ (T)', fontsize = 16)
@@ -74,14 +82,15 @@ def run(plotIt=True):
dmisfit = DataMisfit.l2_DataMisfit(survey)
regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
reg = Regularization.Tikhonov(regMesh)
opt = Optimization.InexactGaussNewton(maxIter = 5)
opt = Optimization.InexactGaussNewton(maxIter = 6)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
m0 = np.log(np.ones(mtrue.size)*sig_half)
reg.alpha_s = 1e-2
reg.alpha_s = 1e-3
reg.alpha_x = 1.
prb.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.5
@@ -94,12 +103,12 @@ def run(plotIt=True):
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
ax.set_ylim(-600, 0)
ax.set_xlim(1e-4, 1e-2)
ax.set_ylim(-500, 0)
ax.set_xlim(1e-3, 1e-1)
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
ax.set_ylabel('Depth (m)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'])
plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best')
plt.show()
SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
+43
View File
@@ -0,0 +1,43 @@
from SimPEG import *
import SimPEG.EM as EM
def run(XYZ=None, loc=np.r_[0.,0.,0.], sig=1.0, freq=1.0, orientation='Z', plotIt=True):
"""
EM: Magnetic Dipole in a Whole-Space
====================================
Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
"""
if XYZ is None:
x = np.arange(-100.5,100.5,step = 1.) #(avoid putting measurement points where source is located)
y = np.r_[0]
z = x
XYZ = Utils.ndgrid(x,y,z)
Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, loc, sig, freq, orientation=orientation)
absB = np.sqrt(Bx*Bx.conj()+By*By.conj()+Bz*Bz.conj()).real
if plotIt:
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
fig, ax = plt.subplots(1,1,figsize=(6,5))
bxplt = Bx.reshape(x.size,z.size)
bzplt = Bz.reshape(x.size,z.size)
pc = ax.pcolor(x,z,absB.reshape(x.size,z.size),norm=LogNorm())
ax.streamplot(x,z,bxplt.real,bzplt.real,color='k',density=1)
ax.set_xlim([x.min(),x.max()])
ax.set_ylim([z.min(),z.max()])
ax.set_xlabel('x')
ax.set_ylabel('z')
cb = plt.colorbar(pc,ax = ax)
cb.set_label('|B| (T)')
plt.show()
return fig, ax
if __name__ == '__main__':
run()
SimPEG/Examples/EM_TDEM_1D_Inversion.py
+106
View File
@@ -0,0 +1,106 @@
from SimPEG import *
import SimPEG.EM as EM
from SimPEG.EM import mu_0
def run(plotIt=True):
"""
EM: TDEM: 1D: Inversion
=======================
Here we will create and run a TDEM 1D inversion.
"""
cs, ncx, ncz, npad = 5., 25, 15, 15
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-100.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
sig_half = 2e-3
sig_air = 1e-8
sig_layer = 1e-3
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
sigma[layer] = sig_layer
mtrue = np.log(sigma[active])
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
ax.set_ylim(-600, 0)
ax.set_xlim(1e-4, 1e-2)
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
ax.set_ylabel('Depth (m)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
rxOffset=1e-3
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
survey = EM.TDEM.SurveyTDEM([src])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
prb.Solver = SolverLU
prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
prb.pair(survey)
# create observed data
std = 0.05
survey.dobs = survey.makeSyntheticData(mtrue,std)
survey.std = std
survey.eps = 1e-5*np.linalg.norm(survey.dobs)
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (10, 6))
ax.loglog(rx.times, survey.dtrue, 'b.-')
ax.loglog(rx.times, survey.dobs, 'r.-')
ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
ax.set_xlabel('Time (s)', fontsize = 14)
ax.set_ylabel('$B_z$ (T)', fontsize = 16)
ax.set_xlabel('Time (s)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
dmisfit = DataMisfit.l2_DataMisfit(survey)
regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
reg = Regularization.Tikhonov(regMesh)
opt = Optimization.InexactGaussNewton(maxIter = 5)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
m0 = np.log(np.ones(mtrue.size)*sig_half)
reg.alpha_s = 1e-2
reg.alpha_x = 1.
prb.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.5
opt.remember('xc')
mopt = inv.run(m0)
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
ax.set_ylim(-600, 0)
ax.set_xlim(1e-4, 1e-2)
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
ax.set_ylabel('Depth (m)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'])
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/__init__.py
+8 -5
View File
@@ -2,6 +2,8 @@
##### AUTOIMPORTS #####
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Forward_BasicDirectCurrent
import Inversion_Linear
@@ -13,7 +15,7 @@ import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
__examples__ = ["EM_FDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
##### AUTOIMPORTS #####
@@ -28,16 +30,17 @@ if __name__ == '__main__':
from SimPEG import Examples
# Create the examples dir in the docs folder.
docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples'])
fName = os.path.realpath(__file__)
docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
shutil.rmtree(docExamplesDir)
os.makedirs(docExamplesDir)
# Get all the python examples in this folder
thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1])
thispath = os.path.sep.join(fName.split(os.path.sep)[:-1])
exfiles = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')]
# Add the imports to the top in the AUTOIMPORTS section
f = file(__file__, 'r')
f = file(fName, 'r')
inimports = False
out = ''
for line in f:
@@ -52,7 +55,7 @@ if __name__ == '__main__':
out += '\n##### AUTOIMPORTS #####\n'
f.close()
f = file(__file__, 'w')
f = file(fName, 'w')
f.write(out)
f.close()
SimPEG/Survey.py
+1
View File
@@ -205,6 +205,7 @@ class BaseSurvey(object):
__metaclass__ = Utils.SimPEGMetaClass
std = None #: Estimated Standard Deviations
eps = None #: Estimated Noise Floor
dobs = None #: Observed data
dtrue = None #: True data, if data is synthetic
mtrue = None #: True model, if data is synthetic
SimPEG/Utils/matutils.py
+12
View File
@@ -26,6 +26,9 @@ def mkvc(x, numDims=1):
if hasattr(x, 'tovec'):
x = x.tovec()
if isinstance(x, Zero):
return x
assert isinstance(x, np.ndarray), "Vector must be a numpy array"
if numDims == 1:
@@ -37,6 +40,9 @@ def mkvc(x, numDims=1):
def sdiag(h):
"""Sparse diagonal matrix"""
if isinstance(h, Zero):
return h
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
def sdInv(M):
@@ -417,6 +423,12 @@ class Zero(object):
def __ge__(self, v):return 0 >= v
def __gt__(self, v):return 0 > v
@property
def transpose(self): return Zero()
@property
def T(self): return Zero()
class Identity(object):
_positive = True
def __init__(self, positive=True):
docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst
+26
View File
@@ -0,0 +1,26 @@
.. _examples_EM_FDEM_Analytic_MagDipoleWholespace:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
EM: Magnetic Dipole in a Whole-Space
====================================
Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
.. plot::
from SimPEG import Examples
Examples.EM_FDEM_Analytic_MagDipoleWholespace.run()
.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
:language: python
:linenos:
docs/examples/EM_TDEM_1D_Inversion.rst
+26
View File
@@ -0,0 +1,26 @@
.. _examples_EM_TDEM_1D_Inversion:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
EM: TDEM: 1D: Inversion
=======================
Here we will create and run a TDEM 1D inversion.
.. plot::
from SimPEG import Examples
Examples.EM_TDEM_1D_Inversion.run()
.. literalinclude:: ../../SimPEG/Examples/EM_TDEM_1D_Inversion.py
:language: python
:linenos:
tests/utils/test_Zero.py
+6 -1
View File
@@ -1,5 +1,5 @@
import unittest
from SimPEG.Utils import Zero, Identity, sdiag
from SimPEG.Utils import Zero, Identity, sdiag, mkvc
from SimPEG import np, sp
class Tests(unittest.TestCase):
@@ -29,6 +29,11 @@ class Tests(unittest.TestCase):
assert a == 1
self.assertRaises(ZeroDivisionError, lambda:3/z)
assert mkvc(z) == 0
assert sdiag(z)*a == 0
assert z.T == 0
assert z.transpose == 0
def test_mat_zero(self):
z = Zero()
S = sdiag(np.r_[2,3])