From 7a4aa4ddc3140b5c9293f85e7f589e82f4428132 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 15 Jan 2016 10:07:11 -0800 Subject: [PATCH 01/13] Update and rename EM_FDEM_1D_Inversion.py to EM_TDEM_1D_Inversion.py - this is a TDEM inversion (not FDEM) --- .../{EM_FDEM_1D_Inversion.py => EM_TDEM_1D_Inversion.py} | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) rename SimPEG/Examples/{EM_FDEM_1D_Inversion.py => EM_TDEM_1D_Inversion.py} (97%) diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py similarity index 97% rename from SimPEG/Examples/EM_FDEM_1D_Inversion.py rename to SimPEG/Examples/EM_TDEM_1D_Inversion.py index aba70f4b..d4d80e55 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_TDEM_1D_Inversion.py @@ -5,10 +5,10 @@ from scipy.constants import mu_0 def run(plotIt=True): """ - EM: FDEM: 1D: Inversion + EM: TDEM: 1D: Inversion ======================= - Here we will create and run a FDEM 1D inversion. + Here we will create and run a TDEM 1D inversion. """ From 0b971f4a50315d439c82672cbc6a885a9712e106 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 Jan 2016 10:46:31 -0800 Subject: [PATCH 02/13] updated examples init to also update docs --- SimPEG/Examples/__init__.py | 4 ++-- ..._FDEM_1D_Inversion.rst => EM_TDEM_1D_Inversion.rst} | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) rename docs/examples/{EM_FDEM_1D_Inversion.rst => EM_TDEM_1D_Inversion.rst} (64%) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 8431e4ba..78b7c5a6 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,7 +1,7 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### -import EM_FDEM_1D_Inversion +import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent import Inversion_Linear @@ -13,7 +13,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_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 ##### diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_TDEM_1D_Inversion.rst similarity index 64% rename from docs/examples/EM_FDEM_1D_Inversion.rst rename to docs/examples/EM_TDEM_1D_Inversion.rst index acbc8cdc..53f6f9ef 100644 --- a/docs/examples/EM_FDEM_1D_Inversion.rst +++ b/docs/examples/EM_TDEM_1D_Inversion.rst @@ -1,4 +1,4 @@ -.. _examples_EM_FDEM_1D_Inversion: +.. _examples_EM_TDEM_1D_Inversion: .. --------------------------------- .. .. .. @@ -9,18 +9,18 @@ .. --------------------------------- .. -EM: FDEM: 1D: Inversion +EM: TDEM: 1D: Inversion ======================= -Here we will create and run a FDEM 1D inversion. +Here we will create and run a TDEM 1D inversion. .. plot:: from SimPEG import Examples - Examples.EM_FDEM_1D_Inversion.run() + Examples.EM_TDEM_1D_Inversion.run() -.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py +.. literalinclude:: ../../SimPEG/Examples/EM_TDEM_1D_Inversion.py :language: python :linenos: From 7107cd9d94c392d37b007ddea36a05d5fae3bca3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 19 Jan 2016 23:57:41 -0800 Subject: [PATCH 03/13] plot a harmonic mag dipole --- .../EM_Analytic_MagDipoleWholespace.py | 47 +++++++++++++++++++ SimPEG/Examples/__init__.py | 3 +- .../EM_Analytic_MagDipoleWholespace.rst | 26 ++++++++++ 3 files changed, 75 insertions(+), 1 deletion(-) create mode 100644 SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py create mode 100644 docs/examples/EM_Analytic_MagDipoleWholespace.rst diff --git a/SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py b/SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py new file mode 100644 index 00000000..31408bcc --- /dev/null +++ b/SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py @@ -0,0 +1,47 @@ +from SimPEG import * +import SimPEG.EM as EM + +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm + +x = np.arange(-100.5,100.5,step = 1.) #(avoid putting measurement points where source is located) +y = np.r_[0] +z = x + +sig = 1. +freq = 1. + +XYZ = Utils.ndgrid(x,y,z) + + +def run(XYZ=XYZ,sig=sig,freq=freq,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. + + """ + + Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, np.r_[0.,0.,0.], sig, freq,orientation=orientation) + absB = np.sqrt(Bx*Bx.conj()+By*By.conj()+Bz*Bz.conj()).real + + + if plotIt: + 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() \ No newline at end of file diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 78b7c5a6..055feeac 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,6 +1,7 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### +import EM_Analytic_MagDipoleWholespace import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent @@ -13,7 +14,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -__examples__ = ["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"] +__examples__ = ["EM_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 ##### diff --git a/docs/examples/EM_Analytic_MagDipoleWholespace.rst b/docs/examples/EM_Analytic_MagDipoleWholespace.rst new file mode 100644 index 00000000..bfe8d33b --- /dev/null +++ b/docs/examples/EM_Analytic_MagDipoleWholespace.rst @@ -0,0 +1,26 @@ +.. _examples_EM_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_Analytic_MagDipoleWholespace.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py + :language: python + :linenos: From 98f209d1f0a368912a42442aab7a5b9aa03d993f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 20 Jan 2016 08:04:01 -0800 Subject: [PATCH 04/13] specify that example is FDEM --- ...holespace.py => EM_FDEM_Analytic_MagDipoleWholespace.py} | 0 SimPEG/Examples/__init__.py | 4 ++-- ...lespace.rst => EM_FDEM_Analytic_MagDipoleWholespace.rst} | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) rename SimPEG/Examples/{EM_Analytic_MagDipoleWholespace.py => EM_FDEM_Analytic_MagDipoleWholespace.py} (100%) rename docs/examples/{EM_Analytic_MagDipoleWholespace.rst => EM_FDEM_Analytic_MagDipoleWholespace.rst} (73%) diff --git a/SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py b/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py similarity index 100% rename from SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py rename to SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 055feeac..894d6ea0 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,7 +1,7 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### -import EM_Analytic_MagDipoleWholespace +import EM_FDEM_Analytic_MagDipoleWholespace import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent @@ -14,7 +14,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -__examples__ = ["EM_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"] +__examples__ = ["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 ##### diff --git a/docs/examples/EM_Analytic_MagDipoleWholespace.rst b/docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst similarity index 73% rename from docs/examples/EM_Analytic_MagDipoleWholespace.rst rename to docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst index bfe8d33b..0e1b2d93 100644 --- a/docs/examples/EM_Analytic_MagDipoleWholespace.rst +++ b/docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst @@ -1,4 +1,4 @@ -.. _examples_EM_Analytic_MagDipoleWholespace: +.. _examples_EM_FDEM_Analytic_MagDipoleWholespace: .. --------------------------------- .. .. .. @@ -19,8 +19,8 @@ Here we plot the magnetic flux density from a harmonic dipole in a wholespace. .. plot:: from SimPEG import Examples - Examples.EM_Analytic_MagDipoleWholespace.run() + Examples.EM_FDEM_Analytic_MagDipoleWholespace.run() -.. literalinclude:: ../../SimPEG/Examples/EM_Analytic_MagDipoleWholespace.py +.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py :language: python :linenos: From 91e37c5fa737f4215bfe2b73146f09e3d146b6fa Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 20 Jan 2016 12:52:51 -0700 Subject: [PATCH 05/13] Clean up init code, put inside run, remove mpl from TL import. --- .../EM_FDEM_Analytic_MagDipoleWholespace.py | 32 ++++++++----------- 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py b/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py index 31408bcc..791f8c63 100644 --- a/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py +++ b/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py @@ -1,33 +1,29 @@ from SimPEG import * import SimPEG.EM as EM -import matplotlib.pyplot as plt -from matplotlib.colors import LogNorm - -x = np.arange(-100.5,100.5,step = 1.) #(avoid putting measurement points where source is located) -y = np.r_[0] -z = x - -sig = 1. -freq = 1. - -XYZ = Utils.ndgrid(x,y,z) - - -def run(XYZ=XYZ,sig=sig,freq=freq,orientation='Z', plotIt=True): +def run(XYZ=None, 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. + Here we plot the magnetic flux density from a harmonic dipole in a wholespace. """ - Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, np.r_[0.,0.,0.], sig, freq,orientation=orientation) + 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, np.r_[0.,0.,0.], sig, freq, orientation=orientation) absB = np.sqrt(Bx*Bx.conj()+By*By.conj()+Bz*Bz.conj()).real - if plotIt: + 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) @@ -44,4 +40,4 @@ def run(XYZ=XYZ,sig=sig,freq=freq,orientation='Z', plotIt=True): return fig, ax if __name__ == '__main__': - run() \ No newline at end of file + run() From 84400650fa35240c4f2bc3b8b0ab90a995c741cc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 23 Jan 2016 16:37:11 -0700 Subject: [PATCH 06/13] Clean up init file in examples. --- SimPEG/Examples/__init__.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 894d6ea0..7ef15451 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -29,16 +29,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: @@ -53,7 +54,7 @@ if __name__ == '__main__': out += '\n##### AUTOIMPORTS #####\n' f.close() - f = file(__file__, 'w') + f = file(fName, 'w') f.write(out) f.close() From a8539cc23439795d1bb5b84e648fbd8e556dabe4 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 27 Jan 2016 08:41:28 -0800 Subject: [PATCH 07/13] pass loc as a variable --- SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py b/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py index 791f8c63..15ff9f62 100644 --- a/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py +++ b/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py @@ -1,7 +1,7 @@ from SimPEG import * import SimPEG.EM as EM -def run(XYZ=None, sig=1.0, freq=1.0, orientation='Z', plotIt=True): +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 ==================================== @@ -17,7 +17,7 @@ def run(XYZ=None, sig=1.0, freq=1.0, orientation='Z', plotIt=True): XYZ = Utils.ndgrid(x,y,z) - Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, np.r_[0.,0.,0.], sig, freq, orientation=orientation) + 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 From d50232b385e3f7b3b4ae5c96119d235a64b040c1 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 13:04:13 -0800 Subject: [PATCH 08/13] give Zero a transpose --- SimPEG/Utils/matutils.py | 12 ++++++++++++ tests/utils/test_Zero.py | 7 ++++++- 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index b38bb4a1..ba900c72 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -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): diff --git a/tests/utils/test_Zero.py b/tests/utils/test_Zero.py index 594de6a6..7b3c6e5d 100644 --- a/tests/utils/test_Zero.py +++ b/tests/utils/test_Zero.py @@ -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]) From f5333f35a241c4bba3e6dd92ebc0da6637d1fecc Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 13:07:48 -0800 Subject: [PATCH 09/13] use survey.createSyntheticData to make synthetic data --- SimPEG/Examples/EM_TDEM_1D_Inversion.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/SimPEG/Examples/EM_TDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py index d4d80e55..86198064 100644 --- a/SimPEG/Examples/EM_TDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_TDEM_1D_Inversion.py @@ -50,20 +50,15 @@ def run(plotIt=True): prb.Solver = SolverLU prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)] prb.pair(survey) - dtrue = survey.dpred(mtrue) - - survey.dtrue = dtrue + # create observed data 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.dobs = survey.makeSyntheticData(mtrue,std) 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.dtrue, 'b.-') ax.loglog(rx.times, survey.dobs, 'r.-') ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16) ax.set_xlabel('Time (s)', fontsize = 14) @@ -76,6 +71,7 @@ def run(plotIt=True): 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) From c24416ea126dc046c26fa5b4afbecb207cc6778b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 13:20:48 -0800 Subject: [PATCH 10/13] use the std and eps from survey in definition of data misfit if defined (and demonstrate it with the TDEM example) --- SimPEG/DataMisfit.py | 28 ++++++++++--------------- SimPEG/Examples/EM_TDEM_1D_Inversion.py | 3 +++ SimPEG/Survey.py | 1 + 3 files changed, 15 insertions(+), 17 deletions(-) diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index c1203d47..425fe4ce 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -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 diff --git a/SimPEG/Examples/EM_TDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py index 86198064..f217912d 100644 --- a/SimPEG/Examples/EM_TDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_TDEM_1D_Inversion.py @@ -53,7 +53,10 @@ def run(plotIt=True): # 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 diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 88355df1..47a88ae2 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -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 From 4bcc9e08502c99a732646b215adff5f1f54d61ec Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 14:00:47 -0800 Subject: [PATCH 11/13] - name cleanup in Jvec and Jtvec (use u instead of f to be consistent with the rest of SimPEG) - example 1D inversion for FDEM --- SimPEG/EM/FDEM/FDEM.py | 30 +++--- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 116 ++++++++++++++++++++++++ 2 files changed, 131 insertions(+), 15 deletions(-) create mode 100644 SimPEG/Examples/EM_FDEM_1D_Inversion.py diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index f2167fd8..cb6adc12 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -53,13 +53,13 @@ class BaseFDEMProblem(BaseEMProblem): 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,33 +71,33 @@ 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) 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 @@ -113,12 +113,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 @@ -127,7 +127,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 @@ -140,7 +140,7 @@ class BaseFDEMProblem(BaseEMProblem): else: raise Exception('Must be real or imag') - return Jtv + return Utils.mkvc(Jtv) def getSourceTerm(self, freq): """ diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py new file mode 100644 index 00000000..a6478016 --- /dev/null +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -0,0 +1,116 @@ +from SimPEG import * +import SimPEG.EM as EM +from scipy.constants import mu_0 + + +def run(plotIt=True): + """ + EM: FDEM: 1D: Inversion + ======================= + + Here we will create and run a FDEM 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') + + layerz = -100. + + active = mesh.vectorCCz<0. + 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-2 + sig_air = 1e-8 + 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(-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=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.pair(survey) + + std = 0.05 + 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 = (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) + 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 = 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-3 + 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(-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}$'],loc='best') + plt.show() + + +if __name__ == '__main__': + run() From 37d37369d89958f88fb4d9817517739e41578a27 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 22:36:39 -0800 Subject: [PATCH 12/13] added example to docs --- SimPEG/Examples/__init__.py | 3 ++- docs/examples/EM_FDEM_1D_Inversion.rst | 26 ++++++++++++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) create mode 100644 docs/examples/EM_FDEM_1D_Inversion.rst diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 7ef15451..47721ccf 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,6 +1,7 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### +import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 @@ -14,7 +15,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -__examples__ = ["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"] +__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 ##### diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst new file mode 100644 index 00000000..acbc8cdc --- /dev/null +++ b/docs/examples/EM_FDEM_1D_Inversion.rst @@ -0,0 +1,26 @@ +.. _examples_EM_FDEM_1D_Inversion: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +EM: FDEM: 1D: Inversion +======================= + +Here we will create and run a FDEM 1D inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_FDEM_1D_Inversion.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py + :language: python + :linenos: From 1457e51fc1e8949a3cf6e765a92a2c7c6c09aab3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 23:44:02 -0800 Subject: [PATCH 13/13] from SimPEG.EM import mu_0 --- SimPEG/EM/FDEM/FDEM.py | 2 +- SimPEG/EM/__init__.py | 2 +- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 4 ++-- SimPEG/Examples/EM_TDEM_1D_Inversion.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index ed440144..4c8b0c89 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -142,7 +142,7 @@ class BaseFDEMProblem(BaseEMProblem): raise Exception('Must be real or imag') ATinv.clean() - + return Utils.mkvc(Jtv) def getSourceTerm(self, freq): diff --git a/SimPEG/EM/__init__.py b/SimPEG/EM/__init__.py index 6a1ca774..565f63a8 100644 --- a/SimPEG/EM/__init__.py +++ b/SimPEG/EM/__init__.py @@ -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 diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index a6478016..ff87b6a6 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -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): @@ -59,7 +59,7 @@ def run(plotIt=True): prb.Solver = MumpsSolver except ImportError, e: prb.Solver = SolverLU - + prb.pair(survey) std = 0.05 diff --git a/SimPEG/Examples/EM_TDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py index f217912d..65ae6669 100644 --- a/SimPEG/Examples/EM_TDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_TDEM_1D_Inversion.py @@ -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):