diff --git a/.gitignore b/.gitignore index e1264b7d..5fd423b1 100644 --- a/.gitignore +++ b/.gitignore @@ -40,3 +40,4 @@ nosetests.xml docs/_build/ Makefile docs/warnings.txt +.DS_Store diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 69e85e84..765711b7 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -253,8 +253,7 @@ class SaveOutputDictEveryIteration(SaveEveryIteration): class Update_IRLS(InversionDirective): eps_min = None - eps_p = None - eps_q = None + eps = None norms = [2.,2.,2.,2.] factor = None gamma = None @@ -263,6 +262,7 @@ class Update_IRLS(InversionDirective): f_old = None f_min_change = 1e-2 beta_tol = 5e-2 + prctile = 95 # Solving parameter for IRLS (mode:2) IRLSiter = 0 @@ -297,9 +297,22 @@ class Update_IRLS(InversionDirective): print "Convergence with smooth l2-norm regularization: Start IRLS steps..." self.mode = 2 - print self.eps_p, self.eps_q, self.norms - self.reg.eps_p = self.eps_p - self.reg.eps_q = self.eps_q + + # Either use the supplied epsilon, or fix base on distribution of + # model values + if getattr(self, 'reg.eps', None) is None: + self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile) + else: + self.reg.eps_p = self.eps[0] + + if getattr(self, 'reg.eps', None) is None: + self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile) + else: + self.reg.eps_q = self.eps[1] + + print "L[p qx qy qz]-norm : " + str(self.reg.norms) + print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q) + self.reg.norms = self.norms self.coolingFactor = 1. self.coolingRate = 1 @@ -343,14 +356,14 @@ class Update_IRLS(InversionDirective): else: self.f_old = phim_new - # Cool the threshold parameter if required - if getattr(self, 'factor', None) is not None: - eps = self.reg.eps / self.factor - - if getattr(self, 'eps_min', None) is not None: - self.reg.eps = np.max([self.eps_min,eps]) - else: - self.reg.eps = eps +# # Cool the threshold parameter if required +# if getattr(self, 'factor', None) is not None: +# eps = self.reg.eps / self.factor +# +# if getattr(self, 'eps_min', None) is not None: +# self.reg.eps = np.max([self.eps_min,eps]) +# else: +# self.reg.eps = eps # Get phi_m at the end of current iteration self.phi_m_last = self.invProb.phi_m_last diff --git a/SimPEG/Examples/DC_Analytic_Dipole.py b/SimPEG/Examples/DC_Analytic_Dipole.py index ba999fa3..347446e3 100644 --- a/SimPEG/Examples/DC_Analytic_Dipole.py +++ b/SimPEG/Examples/DC_Analytic_Dipole.py @@ -1,7 +1,7 @@ from SimPEG import * import SimPEG.EM.Static.DC as DC -def run(plotIt=False): +def run(plotIt=True): cs = 25. hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] @@ -65,4 +65,4 @@ def run(plotIt=False): if __name__ == '__main__': - print run(plotIt=True) + print run() diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index afd90525..0a0757cf 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -42,55 +42,33 @@ def run(N=100, plotIt=True): survey = Survey.LinearSurvey() survey.pair(prob) survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk) - #survey.makeSyntheticData(mtrue, std=std_noise) wd = np.ones(nk) * std_noise - #print survey.std[0] - #M = prob.mesh # Distance weighting wr = np.sum(prob.G**2.,axis=0)**0.5 wr = ( wr/np.max(wr) ) -# reg = Regularization.Simple(mesh) -# reg.mref = mref -# reg.cell_weights = wr -# dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd -# -# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) -# invProb = InvProblem.BaseInvProblem(dmis, reg, opt) -# invProb.curModel = m0 -# -# beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) -# target = Directives.TargetMisfit() -# + betaest = Directives.BetaEstimate_ByEig() -# inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) -# -# -# mrec = inv.run(m0) -# ml2 = mrec -# print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) -# -# # Switch regularization to sparse -# phim = invProb.phi_m_last -# phid = invProb.phi_d reg = Regularization.Sparse(mesh) reg.mref = mref reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - eps_p = 5e-2 - eps_q = 5e-2 - norms = [0., 0., 2., 2.] + opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) update_Jacobi = Directives.Update_lin_PreCond() - IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q) + + # Set the IRLS directive, penalize the lowest 25 percentile of model values + # Start with an l2-l2, then switch to lp-norms + norms = [0., 0., 2., 2.] + IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3) inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi]) diff --git a/SimPEG/Examples/Maps_ComboMaps.py b/SimPEG/Examples/Maps_ComboMaps.py new file mode 100644 index 00000000..3e495023 --- /dev/null +++ b/SimPEG/Examples/Maps_ComboMaps.py @@ -0,0 +1,62 @@ +from SimPEG import Mesh, Maps, np + +def run(plotIt=True): + """ + + Maps: ComboMaps + =============== + + We will use an example where we want a 1D layered earth as + our model, but we want to map this to a 2D discretization to do our forward + modeling. We will also assume that we are working in log conductivity still, + so after the transformation we want to map to conductivity space. + To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`), + which does the first part of what we just described. The second part will be + done by the :class:`SimPEG.Maps.ExpMap` described above. + + .. code-block:: python + :linenos: + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.SurjectVertical1D(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + + If you noticed, it was pretty easy to combine maps. What is even cooler is + that the derivatives also are made for you (if everything goes right). + Just to be sure that the derivative is correct, you should always run the test + on the mapping that you create. + + """ + + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.SurjectVertical1D(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + + if not plotIt: return + + import matplotlib.pyplot as plt + figs, axs = plt.subplots(1,2) + axs[0].plot(m, M.vectorCCy, 'b-o') + axs[0].set_title('Model') + axs[0].set_ylabel('Depth, y') + axs[0].set_xlabel('Value, $m_i$') + axs[0].set_xlim(0,3) + axs[0].set_ylim(0,1) + clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0]) + axs[1].set_title('Physical Property') + axs[1].set_ylabel('Depth, y') + clbar.set_label('$\sigma = \exp(\mathbf{P}m)$') + plt.tight_layout() + plt.show() + + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/Maps_Mesh2Mesh.py b/SimPEG/Examples/Maps_Mesh2Mesh.py new file mode 100644 index 00000000..254dc036 --- /dev/null +++ b/SimPEG/Examples/Maps_Mesh2Mesh.py @@ -0,0 +1,41 @@ +from SimPEG import Mesh, Maps, Utils + +def run(plotIt=True): + """ + + Maps: Mesh2Mesh + =============== + + This mapping allows you to go from one mesh to another. + + """ + + M = Mesh.TensorMesh([100,100]) + h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) + h1 = h1/h1.sum() + M2 = Mesh.TensorMesh([h1,h1]) + V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) + v = Utils.mkvc(V) + modh = Maps.Mesh2Mesh([M,M2]) + modH = Maps.Mesh2Mesh([M2,M]) + H = modH * v + h = modh * H + + if not plotIt: return + + import matplotlib.pyplot as plt + ax = plt.subplot(131) + M.plotImage(v, ax=ax) + ax.set_title('Fine Mesh (Original)') + ax = plt.subplot(132) + M2.plotImage(H,clim=[0,1],ax=ax) + ax.set_title('Course Mesh') + ax = plt.subplot(133) + M.plotImage(h,clim=[0,1],ax=ax) + ax.set_title('Fine Mesh (Interpolated)') + plt.show() + + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py index e4f748cb..07e3e4aa 100644 --- a/SimPEG/Examples/Utils_surface2ind_topo.py +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -2,7 +2,7 @@ from SimPEG import * from SimPEG.Utils import surface2ind_topo -def run(plotIt=False, nx=5, ny=5): +def run(plotIt=True, nx=5, ny=5): """ Utils: surface2ind_topo diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 2ee4eb6d..3c9cbca9 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -10,6 +10,8 @@ import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Inversion_IRLS import Inversion_Linear +import Maps_ComboMaps +import Maps_Mesh2Mesh import Mesh_Basic_ForwardDC import Mesh_Basic_PlotImage import Mesh_Basic_Types @@ -22,7 +24,7 @@ import MT_1D_ForwardAndInversion import MT_3D_Foward import Utils_surface2ind_topo -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "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", "Utils_surface2ind_topo"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Maps_ComboMaps", "Maps_Mesh2Mesh", "Mesh_Basic_ForwardDC", "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", "Utils_surface2ind_topo"] ##### AUTOIMPORTS ##### diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index ecfd51e1..77e7e9a3 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -502,7 +502,9 @@ class InjectActiveCells(IdentityMap): if Utils.isScalar(valInactive): self.valInactive = np.ones(self.nC)*float(valInactive) else: - self.valInactive = valInactive.copy() + self.valInactive = np.ones(self.nC) + self.valInactive[self.indInactive] = valInactive.copy() + self.valInactive[self.indActive] = 0 inds = np.nonzero(self.indActive)[0] diff --git a/docs/content/api_core/api_Maps.rst b/docs/content/api_core/api_Maps.rst index a5de8a2f..e1e744c6 100644 --- a/docs/content/api_core/api_Maps.rst +++ b/docs/content/api_core/api_Maps.rst @@ -63,26 +63,8 @@ done by the :class:`SimPEG.Maps.ExpMap` described above. .. plot:: - from SimPEG import * - import matplotlib.pyplot as plt - M = Mesh.TensorMesh([7,5]) - v1dMap = Maps.SurjectVertical1D(M) - expMap = Maps.ExpMap(M) - myMap = expMap * v1dMap - m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! - sig = myMap * m - figs, axs = plt.subplots(1,2) - axs[0].plot(m, M.vectorCCy, 'b-o') - axs[0].set_title('Model') - axs[0].set_ylabel('Depth, y') - axs[0].set_xlabel('Value, $m_i$') - axs[0].set_xlim(0,3) - axs[0].set_ylim(0,1) - clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0]) - axs[1].set_title('Physical Property') - axs[1].set_ylabel('Depth, y') - clbar.set_label('$\sigma = \exp(\mathbf{P}m)$') - plt.tight_layout() + from SimPEG import Examples + Examples.Maps_ComboMaps.run() If you noticed, it was pretty easy to combine maps. What is even cooler is that the derivatives also are made for you (if everything goes right). @@ -167,31 +149,10 @@ Map 2D Cross-Section to 3D Model Mesh to Mesh Map ---------------- - .. plot:: - from SimPEG import * - import matplotlib.pyplot as plt - M = Mesh.TensorMesh([100,100]) - h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) - h1 = h1/h1.sum() - M2 = Mesh.TensorMesh([h1,h1]) - V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) - v = Utils.mkvc(V) - modh = Maps.Mesh2Mesh([M,M2]) - modH = Maps.Mesh2Mesh([M2,M]) - H = modH * v - h = modh * H - ax = plt.subplot(131) - M.plotImage(v, ax=ax) - ax.set_title('Fine Mesh (Original)') - ax = plt.subplot(132) - M2.plotImage(H,clim=[0,1],ax=ax) - ax.set_title('Course Mesh') - ax = plt.subplot(133) - M.plotImage(h,clim=[0,1],ax=ax) - ax.set_title('Fine Mesh (Interpolated)') - plt.show() + from SimPEG import Examples + Examples.Maps_Mesh2Mesh.run() .. autoclass:: SimPEG.Maps.Mesh2Mesh diff --git a/docs/content/examples/Inversion_IRLS.rst b/docs/content/examples/Inversion_IRLS.rst new file mode 100644 index 00000000..192ada50 --- /dev/null +++ b/docs/content/examples/Inversion_IRLS.rst @@ -0,0 +1,26 @@ +.. _examples_Inversion_IRLS: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Inversion: Linear Problem +========================= + +Here we go over the basics of creating a linear problem and inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.Inversion_IRLS.run() + +.. literalinclude:: ../../../SimPEG/Examples/Inversion_IRLS.py + :language: python + :linenos: diff --git a/docs/content/examples/Maps_ComboMaps.rst b/docs/content/examples/Maps_ComboMaps.rst new file mode 100644 index 00000000..69ebbf39 --- /dev/null +++ b/docs/content/examples/Maps_ComboMaps.rst @@ -0,0 +1,48 @@ +.. _examples_Maps_ComboMaps: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + + +Maps: ComboMaps +=============== + +We will use an example where we want a 1D layered earth as +our model, but we want to map this to a 2D discretization to do our forward +modeling. We will also assume that we are working in log conductivity still, +so after the transformation we want to map to conductivity space. +To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`), +which does the first part of what we just described. The second part will be +done by the :class:`SimPEG.Maps.ExpMap` described above. + +.. code-block:: python + :linenos: + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.SurjectVertical1D(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + +If you noticed, it was pretty easy to combine maps. What is even cooler is +that the derivatives also are made for you (if everything goes right). +Just to be sure that the derivative is correct, you should always run the test +on the mapping that you create. + + + +.. plot:: + + from SimPEG import Examples + Examples.Maps_ComboMaps.run() + +.. literalinclude:: ../../../SimPEG/Examples/Maps_ComboMaps.py + :language: python + :linenos: diff --git a/docs/content/examples/Maps_Mesh2Mesh.rst b/docs/content/examples/Maps_Mesh2Mesh.rst new file mode 100644 index 00000000..60233b57 --- /dev/null +++ b/docs/content/examples/Maps_Mesh2Mesh.rst @@ -0,0 +1,27 @@ +.. _examples_Maps_Mesh2Mesh: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + + +Maps: Mesh2Mesh +=============== + +This mapping allows you to go from one mesh to another. + + + +.. plot:: + + from SimPEG import Examples + Examples.Maps_Mesh2Mesh.run() + +.. literalinclude:: ../../../SimPEG/Examples/Maps_Mesh2Mesh.py + :language: python + :linenos: