fixed merge conflict in examples __init__.py

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -0,0 +1,116 @@
+from SimPEG import *
+import SimPEG.EM as EM
+from SimPEG.EM 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()

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):
@@ -50,20 +50,18 @@ 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)
+    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, 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 +74,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)

SimPEG/Examples/__init__.py

@@ -1,22 +1,23 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import Mesh_QuadTree_Creation
-import EM_TDEM_1D_Inversion
-import Mesh_QuadTree_FaceDiv
-import Mesh_Tensor_Creation
-import FLOW_Richards_1D_Celia1990
-import Mesh_Operators_CahnHilliard
-import Mesh_Basic_Types
-import Inversion_Linear
-import MT_3D_Foward
-import MT_1D_ForwardAndInversion
-import Forward_BasicDirectCurrent
+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
 import Mesh_Basic_PlotImage
+import Mesh_Basic_Types
+import Mesh_Operators_CahnHilliard
+import Mesh_QuadTree_Creation
+import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
+import Mesh_Tensor_Creation
+import MT_1D_ForwardAndInversion
+import MT_3D_Foward
 
-__examples__ = ["Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_Linear", "MT_3D_Foward", "MT_1D_ForwardAndInversion", "Forward_BasicDirectCurrent", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "Mesh_QuadTree_HangingNodes"]
+__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", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
 
 ##### AUTOIMPORTS #####