diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 3ca3f780..31842d79 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.11 +current_version = 0.1.12 files = setup.py SimPEG/__init__.py docs/conf.py 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/README.rst b/README.rst index c9b30a78..85777c4c 100644 --- a/README.rst +++ b/README.rst @@ -1,4 +1,4 @@ -.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png +.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png :alt: SimPEG Logo ====== 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/EM/Analytics/FDEMDipolarfields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py new file mode 100644 index 00000000..e65bbdd1 --- /dev/null +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -0,0 +1,302 @@ +from __future__ import division +import numpy as np +from scipy.constants import mu_0, pi, epsilon_0 +from scipy.special import erf +from SimPEG import Utils + +omega = lambda f: 2.*np.pi*f +# TODO: +# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.) +# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + +def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.): + + """ + Computing Analytic Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omega(f)*epsilon + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ey = front*(dx*dy / r**2)*mid + Ez = front*(dx*dz / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ez = front*(dy*dz / r**2)*mid + Ex = front*(dy*dx / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ex = front*(dz*dx / r**2)*mid + Ey = front*(dz*dy / r**2)*mid + return Ex, Ey, Ez + + +def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omega(f)*epsilon + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.)) + Ey_galvanic = front*(dx*dy / r**2)*mid + Ez_galvanic = front*(dx*dz / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.)) + Ez_galvanic = front*(dy*dz / r**2)*mid + Ex_galvanic = front*(dy*dx / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.)) + Ex_galvanic = front*(dz*dx / r**2)*mid + Ey_galvanic = front*(dz*dy / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + +def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omega(f)*epsilon + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Ex_inductive = front*(k**2 * r**2) + Ey_inductive = np.zeros_like(Ex_inductive) + Ez_inductive = np.zeros_like(Ex_inductive) + return Ex_inductive, Ey_inductive, Ez_inductive + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey_inductive = front*(k**2 * r**2) + Ez_inductive = np.zeros_like(Ey_inductive) + Ex_inductive = np.zeros_like(Ey_inductive) + return Ex_inductive, Ey_inductive, Ez_inductive + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez_inductive = front*(k**2 * r**2) + Ex_inductive = np.zeros_like(Ez_inductive) + Ey_inductive = np.zeros_like(Ez_inductive) + return Ex_inductive, Ey_inductive, Ez_inductive + + +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) + Jx = sig*Ex + Jy = sig*Ey + Jz = sig*Ez + return Jx, Jy, Jz + + +def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) + Jx_galvanic = sig*Ex_galvanic + Jy_galvanic = sig*Ey_galvanic + Jz_galvanic = sig*Ez_galvanic + return Jx_galvanic, Jy_galvanic, Jz_galvanic + + +def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) + Jx_inductive = sig*Ex_inductive + Jy_inductive = sig*Ey_inductive + Jz_inductive = sig*Ez_inductive + return Jx_inductive, Jy_inductive, Jz_inductive + + +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Magnetic fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + + front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Hy = front*(-dz / r) + Hz = front*(dy / r) + Hx = np.zeros_like(Hy) + return Hx, Hy, Hz + + elif orientation.upper() == 'Y': + Hx = front*(dz / r) + Hz = front*(-dx / r) + Hy = np.zeros_like(Hx) + return Hx, Hy, Hz + + elif orientation.upper() == 'Z': + Hx = front*(-dy / r) + Hy = front*(dx / r) + Hz = np.zeros_like(Hx) + return Hx, Hy, Hz + + +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Magnetic flux densites from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) + Bx = mu*Hx + By = mu*Hy + Bz = mu*Hz + return Bx, By, Bz + + +def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Electric vector potentials from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + + front = current * length / (4.*np.pi*r) + + if orientation.upper() == 'X': + Ax = front*np.exp(-1j*k*r) + Ay = np.zeros_like(Ax) + Az = np.zeros_like(Ax) + return Ax, Ay, Az + + elif orientation.upper() == 'Y': + Ay = front*np.exp(-1j*k*r) + Ax = np.zeros_like(Ay) + Az = np.zeros_like(Ay) + return Ax, Ay, Az + + elif orientation.upper() == 'Z': + Az = front*np.exp(-1j*k*r) + Ax = np.zeros_like(Ay) + Ay = np.zeros_like(Ay) + return Ax, Ay, Az + + + + + diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 9df2aef7..8331501d 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -2,3 +2,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * from DC import DCAnalyticHalf, DCAnalyticSphere +from FDEMDipolarfields import * 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 dd2fef10..32ad2aae 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -558,7 +558,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/SimPEG/Survey.py b/SimPEG/Survey.py index fbc88276..65c46972 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -311,7 +311,6 @@ class BaseSurvey(object): if f is None: f = self.prob.fields(m) return Utils.mkvc(self.eval(f)) - @Utils.count def eval(self, f): """eval(f) @@ -322,7 +321,7 @@ class BaseSurvey(object): d_\\text{pred} = \mathbf{P} f(m) """ - raise NotImplemented('eval is not yet implemented.') + raise NotImplementedError('eval is not yet implemented.') @Utils.count def evalDeriv(self, f): @@ -334,7 +333,7 @@ class BaseSurvey(object): \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P} """ - raise NotImplemented('eval is not yet implemented.') + raise NotImplementedError('eval is not yet implemented.') @Utils.count def residual(self, m, f=None): diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 23d4de4a..536a8103 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.11' +__version__ = '0.1.12' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 143ea545..d696464d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013 - 2016, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.11' +version = '0.1.12' # The full version, including alpha/beta/rc tags. -release = '0.1.11' +release = '0.1.12' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. 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: diff --git a/docs/index.rst b/docs/index.rst index fb6fb4f3..fd6959bc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,4 +1,4 @@ -.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png +.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png :alt: SimPEG Logo SimPEG Documentation diff --git a/setup.py b/setup.py index 06cd9b27..a52bfe2c 100644 --- a/setup.py +++ b/setup.py @@ -83,7 +83,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.11", + version = "0.1.12", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13',