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

Possibly deal with a good chunk of the old_divs

Browse Source
This commit is contained in:
Brendan Smithyman
2016-07-17 16:02:43 -05:00
parent 7189ec5b2f
commit f7a70aa6a7
85 changed files with 369 additions and 454 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/BaseDC.py
+3 -4
View File
@@ -6,7 +6,6 @@ from builtins import super
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import *
class FieldsDC_CC(Problem.Fields):
@@ -70,7 +69,7 @@ class SrcDipole(Survey.BaseSrc):
pts = [self.loc[0], self.loc[1]]
inds = Utils.closestPoints(prob.mesh, pts)
q = np.zeros(prob.mesh.nC)
q[inds] = - self.current * ( old_div(np.r_[1., -1.], prob.mesh.vol[inds]) )
q[inds] = - self.current * (np.r_[1., -1.] / prob.mesh.vol[inds])
# self._rhsDict[mesh] = q
# return self._rhsDict[mesh]
return q
@@ -154,7 +153,7 @@ class ProblemDC_CC(Problem.BaseProblem):
if getattr(self, '_Msig', None) is None:
sigma = self.curModel.transform
Av = self.mesh.aveF2CC
self._Msig = Utils.sdiag(old_div(1,(self.mesh.dim * Av.T * (old_div(1,sigma)))))
self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
return self._Msig
@property
@@ -162,7 +161,7 @@ class ProblemDC_CC(Problem.BaseProblem):
if getattr(self, '_dMdsig', None) is None:
sigma = self.curModel.transform
Av = self.mesh.aveF2CC
dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(old_div(1.,sigma**2))
dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
return self._dMdsig
SimPEG/DCIP/BaseIP.py
+2 -3
View File
@@ -5,7 +5,6 @@ from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import *
from .BaseDC import SurveyDC, FieldsDC_CC
@@ -60,7 +59,7 @@ class ProblemIP(Problem.BaseProblem):
# sigma = self.curModel.transform
sigma = self.sigma
Av = self.mesh.aveF2CC
self._Msig = Utils.sdiag(old_div(1,(self.mesh.dim * Av.T * (old_div(1,sigma)))))
self._Msig = Utils.sdiag(1//(self.mesh.dim * Av.T * (1/sigma)))
return self._Msig
@property
@@ -69,7 +68,7 @@ class ProblemIP(Problem.BaseProblem):
# sigma = self.curModel.transform
sigma = self.sigma
Av = self.mesh.aveF2CC
dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(old_div(1.,sigma**2))
dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
return self._dMdsig
SimPEG/DCIP/DCIPUtils.py
+15 -16
View File
@@ -8,7 +8,6 @@ from future import standard_library
standard_library.install_aliases()
from builtins import map
from builtins import range
from past.utils import old_div
from SimPEG import np, Utils
from . import BaseDC as DC
from . import BaseDC as IP
@@ -222,8 +221,8 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
MN = np.abs(Rx[1][:,0] - Rx[0][:,0])
# Create mid-point location
Cmid = old_div((Tx[0][0] + Tx[1][0]),2)
Pmid = old_div((Rx[0][:,0] + Rx[1][:,0]),2)
Cmid = (Tx[0][0] + Tx[1][0])/2
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
# Change output for unitType
if unitType == 'volt':
@@ -239,7 +238,7 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
elif surveyType == 'dipole-dipole':
leg = data * 2*np.pi / ( old_div(1,MA) - old_div(1,MB) - old_div(1,NB) + old_div(1,NA) )
leg = data * 2*np.pi / (1/MA - 1/MB - 1/NB + 1/NA)
else:
print("""unitType must be 'pole-dipole' | 'dipole-dipole' """)
@@ -248,7 +247,7 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
if unitType == 'appConductivity':
leg = np.log10(abs(old_div(1.,leg)))
leg = np.log10(abs(1./leg))
rho = np.hstack([rho,leg])
elif unitType == 'appResistivity':
@@ -260,8 +259,8 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
print("""unitType must be 'appResistivity' | 'appConductivity' | 'volt' """)
break
midx = np.hstack([midx, old_div(( Cmid + Pmid ),2) ])
midz = np.hstack([midz, old_div(-np.abs(Cmid-Pmid),2) + old_div((Tx[0][2] + Tx[1][2]),2) ])
midx = np.hstack([midx, (Cmid + Pmid)/2])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2])
# Grid points
grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
@@ -348,10 +347,10 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
# Mesure survey length and direction
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
dl_x = old_div(( endl[1,0] - endl[0,0] ), dl_len)
dl_y = old_div(( endl[1,1] - endl[0,1] ), dl_len)
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
nstn = np.floor( old_div(dl_len, AM_sep) )
nstn = np.floor(dl_len / AM_sep)
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*AM_sep
@@ -387,7 +386,7 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([np.floor( old_div((AB - MN_sep), AM_sep) ) , nrx])
nstn = np.min([(AB - MN_sep) // AM_sep, nrx])
# Check if there is enough space, else break the loop
if nstn <= 0:
@@ -428,16 +427,16 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
max_y = endl[1,1] - dl_y * MN_sep
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = old_div(box_l,2.)
box_w = box_l/2.
nstn = np.floor( old_div(box_l, AM_sep) )
nstn = box_l // AM_sep
# Compute discrete pole location along line
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*AM_sep
# Define number of cross lines
nlin = int(np.floor( old_div(box_w, AM_sep) ))
nlin = int(box_w // AM_sep)
lind = list(range(-nlin,nlin+1))
ngrad = nstn * len(lind)
@@ -994,7 +993,7 @@ def xy_2_lineID(DCsurvey):
ang2 = np.abs(vec3.dot(vec4))
# If the angles are smaller then 45d, than next point is on a new line
if ((ang1 < np.cos(old_div(np.pi,4.))) | (ang2 < np.cos(old_div(np.pi,4.)))) & (np.all(np.r_[r1,r2,r3,r4] > 0)):
if ((ang1 < np.cos(np.pi/4.)) | (ang2 < np.cos(np.pi/4.))) & (np.all(np.r_[r1,r2,r3,r4] > 0)):
# Re-initiate start and mid-point location
xy0 = A[:2]
@@ -1034,7 +1033,7 @@ def r_unit(p1,p2):
if r!=0:
vec = old_div(dx,r)
vec = dx/r
else:
vec = np.zeros(len(p1))
SimPEG/DCIP/Utils.py
+1 -2
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
@@ -13,7 +12,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
import SimPEG.DCIP as DC
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= old_div((nElecs*aSpacing - aSpacing),2)
elocs -= (nElecs*aSpacing - aSpacing) / 2
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
SimPEG/DataMisfit.py
+1 -2
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from builtins import object
from . import Utils, Survey, Problem
import numpy as np, scipy.sparse as sp, gc
@@ -108,7 +107,7 @@ class l2_DataMisfit(BaseDataMisfit):
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(old_div(1,(abs(survey.dobs)*survey.std+survey.eps)))
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
return self._Wd
@Wd.setter
SimPEG/Directives.py
+6 -7
View File
@@ -7,7 +7,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import str
from past.utils import old_div
from builtins import object
from . import Utils
import numpy as np
@@ -140,7 +139,7 @@ class BetaEstimate_ByEig(InversionDirective):
x0 = np.random.rand(*m.shape)
t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f))
b = x0.dot(self.reg.eval2Deriv(m, v=x0))
self.beta0 = self.beta0_ratio*(old_div(t,b))
self.beta0 = self.beta0_ratio*(t/b)
self.invProb.beta = self.beta0
@@ -350,7 +349,7 @@ class Update_IRLS(InversionDirective):
self.IRLSiter += 1
phim_new = self.reg.eval(self.invProb.curModel)
self.f_change = old_div(np.abs(self.f_old - phim_new), self.f_old)
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
print("Regularization decrease: %6.3e" % (self.f_change))
@@ -397,7 +396,7 @@ class Update_IRLS(InversionDirective):
phim_new = self.reg.eval(self.invProb.curModel)
# Update gamma to scale the regularization between IRLS iterations
self.reg.gamma = old_div(self.phi_m_last, phim_new)
self.reg.gamma = self.phi_m_last / phim_new
# Reset the regularization matrices again for new gamma
self.reg._Wsmall = None
@@ -406,7 +405,7 @@ class Update_IRLS(InversionDirective):
self.reg._Wz = None
# Check if misfit is within the tolerance, otherwise scale beta
val = old_div(self.invProb.phi_d, (self.survey.nD*0.5))
val = self.invProb.phi_d / (self.survey.nD*0.5)
if np.abs(1.-val) > self.beta_tol:
self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
@@ -450,7 +449,7 @@ class Update_Wj(InversionDirective):
m = self.invProb.curModel
if self.k is None:
self.k = int(old_div(self.survey.nD,10))
self.k = int(self.survey.nD/10)
def JtJv(v):
@@ -459,6 +458,6 @@ class Update_Wj(InversionDirective):
return self.prob.Jtvec(m,Jv)
JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k)
JtJdiag = old_div(JtJdiag, max(JtJdiag))
JtJdiag = JtJdiag / max(JtJdiag)
self.reg.wght = JtJdiag
SimPEG/EM/Analytics/DC.py
+4 -5
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
from scipy.constants import mu_0, pi
from scipy import special
@@ -32,8 +31,8 @@ def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"):
rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 )
rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 )
phiM = old_div(1.,(4*np.pi*rM*sigma))
phiN = old_div(1.,(4*np.pi*rN*sigma))
phiM = 1./(4*np.pi*rM*sigma)
phiN = 1./(4*np.pi*rN*sigma)
phi = phiM - phiN
if earth_type == "halfspace":
@@ -77,8 +76,8 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
Pleg.append(special.legendre(i, monic=0))
rho = old_div(1.,sigma)
rho1 = old_div(1.,sigma1)
rho = 1./sigma
rho1 = 1./sigma1
# Center of the sphere should be aligned in txloc in y-direction
yc = txloc[1]
SimPEG/EM/Analytics/FDEMcasing.py
+5 -6
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Utils, np
from scipy.constants import mu_0, epsilon_0
from SimPEG.EM.Utils.EMUtils import k
@@ -13,7 +12,7 @@ def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0):
a = float(a)
b = float(b)
# return 1./(2*np.pi) * np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
return np.sqrt(old_div(b, a)) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
return np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
def _r2(xyz):
return np.sum(xyz**2,1)
@@ -41,7 +40,7 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
k2 = k(freq,sigma[2],mu[2],eps)
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + old_div(1., sqrtr2z2))
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1. / sqrtr2z2)
def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
@@ -54,7 +53,7 @@ def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
sqrtr2z2 = np.sqrt(r2z2)
k2 = k(freq,sigma[2],mu[2],eps)
return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + old_div(1.,sqrtr2z2))
return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2)
def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
@@ -69,7 +68,7 @@ def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.o
sqrtr2z2 = np.sqrt(r2 + z**2)
k2 = k(freq,sigma[2],mu[2],eps)
return dHertzZdr*(old_div(-z,sqrtr2z2))*(1j*k2+old_div(1.,sqrtr2z2)) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + old_div(2.,sqrtr2z2))
return dHertzZdr*(-z/sqrtr2z2)*(1j*k2+1./sqrtr2z2) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + 2./sqrtr2z2)
def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
@@ -84,7 +83,7 @@ def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.o
sqrtr2z2 = np.sqrt(r2 + z**2)
k2 = k(freq,sigma[2],mu[2],eps)
return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - old_div(1.,sqrtr2z2)) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2)
return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - 1./sqrtr2z2) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2)
def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
return 1j * omega(freq) * mu * _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
SimPEG/EM/Analytics/TDEM.py
+4 -5
View File
@@ -4,16 +4,15 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
def hzAnalyticDipoleT(r, t, sigma):
theta = np.sqrt(old_div((sigma*mu_0),(4*t)))
theta = np.sqrt((sigma*mu_0)/(4*t))
tr = theta*r
etr = erf(tr)
t1 = (old_div(9,(2*tr**2)) - 1)*etr
t2 = (old_div(1,np.sqrt(pi)))*(old_div(9,tr) + 4*tr)*np.exp(-tr**2)
hz = old_div((t1 - t2),(4*pi*r**3))
t1 = (9/(2*tr**2) - 1)*etr
t2 = (1/np.sqrt(pi))*(9/tr + 4*tr)*np.exp(-tr**2)
hz = (t1 - t2)/(4*pi*r**3)
return hz
SimPEG/EM/Base.py
+1 -2
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
@@ -18,7 +17,7 @@ class EMPropMap(Maps.PropMap):
mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = old_div(1.,mu_0), propertyLink=('mu', Maps.ReciprocalMap))
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
class BaseEMProblem(Problem.BaseProblem):
SimPEG/EM/FDEM/FieldsFDEM.py
+44 -45
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import scipy.sparse as sp
import SimPEG
@@ -293,7 +292,7 @@ class Fields3D_e(FieldsFDEM):
C = self._edgeCurl
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= old_div(- 1.,(1j*omega(src.freq)))
b[:,i] *= -1./(1j*omega(src.freq))
s_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
return b
@@ -339,8 +338,8 @@ class Fields3D_e(FieldsFDEM):
:return: current density
"""
aveE2CCV = self._aveE2CCV
n = int(old_div(aveE2CCV.shape[0], self._nC)) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (aveE2CCV * (self._MeSigma * self._e(eSolution,srcList) ) )
def _jDeriv_u(self, src, du_dm_v, adjoint = False):
@@ -353,8 +352,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint)
@@ -372,8 +371,8 @@ class Fields3D_e(FieldsFDEM):
:return: product of the current density derivative with respect to the inversion model with a vector
"""
e = self[src, 'e']
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
@@ -390,8 +389,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # Number of Components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1. // self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
@@ -405,8 +404,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # Number of Components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
return self._bDeriv_u(src, v, adjoint=adjoint)
@@ -422,8 +421,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # Number of Components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
return self._bDeriv_m(src, v, adjoint=adjoint)
@@ -615,8 +614,8 @@ class Fields3D_b(FieldsFDEM):
:return: primary current density
"""
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
@@ -632,8 +631,8 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
return VI * (self._aveE2CCV * (self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) )
@@ -660,8 +659,8 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
@@ -675,8 +674,8 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfMui.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v) )
@@ -835,7 +834,7 @@ class Fields3D_j(FieldsFDEM):
h = (self._edgeCurl.T * (self._MfRho * jSolution) )
for i, src in enumerate(srcList):
h[:,i] *= old_div(-1.,(1j*omega(src.freq)))
h[:,i] *= -1./(1j*omega(src.freq))
s_m,_ = src.eval(self.prob)
h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m)
return self._MeMuI * h
@@ -898,8 +897,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
@@ -912,8 +911,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
return VI * (self._aveF2CCV * (self._MfRho * du_dm_v))
@@ -929,8 +928,8 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the electric field with respect to the model with a vector
"""
jSolution = Utils.mkvc(self[src,'jSolution'])
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
return VI * (self._aveF2CCV * (self._MfRhoDeriv(jSolution) * v))
@@ -944,8 +943,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
@@ -959,8 +958,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return -1./(1j*omega(src.freq)) * self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
@@ -977,8 +976,8 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the magnetic flux density with respect to the model with a vector
"""
jSolution = self[src,'jSolution']
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
if adjoint:
@@ -1159,8 +1158,8 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
@@ -1173,8 +1172,8 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v ))
@@ -1190,8 +1189,8 @@ class Fields3D_h(FieldsFDEM):
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
hSolution = Utils.mkvc(self[src,'hSolution'])
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v ))
@@ -1206,8 +1205,8 @@ class Fields3D_h(FieldsFDEM):
:return: magnetic flux density
"""
h = self._h(hSolution, srcList)
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1221,8 +1220,8 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
return VI * (self._aveE2CCV * (self._MeMu * du_dm_v))
SimPEG/EM/FDEM/SrcFDEM.py
+3 -4
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
@@ -382,7 +381,7 @@ class MagDipole(BaseSrc):
formulation = prob._formulation
if formulation is 'EB':
mui_s = prob.curModel.mui - old_div(1.,self.mu)
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif formulation is 'HJ':
@@ -496,7 +495,7 @@ class MagDipole_Bfield(BaseSrc):
formulation = prob._formulation
if formulation is 'EB':
mui_s = prob.curModel.mui - old_div(1.,self.mu)
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif formulation is 'HJ':
@@ -608,7 +607,7 @@ class CircularLoop(BaseSrc):
formulation = prob._formulation
if formulation is 'EB':
mui_s = prob.curModel.mui - old_div(1.,self.mu)
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
SimPEG/EM/Static/DC/BoundaryUtils.py
+30 -31
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
def getxBCyBC_CC(mesh, alpha, beta, gamma):
@@ -38,15 +37,15 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
# h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
h_xm, h_xp = mesh.hx[0], mesh.hx[-1]
a_xm = old_div(gamma_xm,(0.5*alpha_xm-old_div(beta_xm,h_xm)))
b_xm = old_div((0.5*alpha_xm+old_div(beta_xm,h_xm)),(0.5*alpha_xm-old_div(beta_xm,h_xm)))
a_xp = old_div(gamma_xp,(0.5*alpha_xp-old_div(beta_xp,h_xp)))
b_xp = old_div((0.5*alpha_xp+old_div(beta_xp,h_xp)),(0.5*alpha_xp-old_div(beta_xp,h_xp)))
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
xBC_xm = 0.5*a_xm
xBC_xp = 0.5*a_xp/b_xp
yBC_xm = 0.5*(1.-b_xm)
yBC_xp = 0.5*(1.-old_div(1.,b_xp))
yBC_xp = 0.5*(1.-1./b_xp)
xBC = np.r_[xBC_xm, xBC_xp]
yBC = np.r_[yBC_xm, yBC_xp]
@@ -69,24 +68,24 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
a_xm = old_div(gamma_xm,(0.5*alpha_xm-old_div(beta_xm,h_xm)))
b_xm = old_div((0.5*alpha_xm+old_div(beta_xm,h_xm)),(0.5*alpha_xm-old_div(beta_xm,h_xm)))
a_xp = old_div(gamma_xp,(0.5*alpha_xp-old_div(beta_xp,h_xp)))
b_xp = old_div((0.5*alpha_xp+old_div(beta_xp,h_xp)),(0.5*alpha_xp-old_div(beta_xp,h_xp)))
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
a_ym = old_div(gamma_ym,(0.5*alpha_ym-old_div(beta_ym,h_ym)))
b_ym = old_div((0.5*alpha_ym+old_div(beta_ym,h_ym)),(0.5*alpha_ym-old_div(beta_ym,h_ym)))
a_yp = old_div(gamma_yp,(0.5*alpha_yp-old_div(beta_yp,h_yp)))
b_yp = old_div((0.5*alpha_yp+old_div(beta_yp,h_yp)),(0.5*alpha_yp-old_div(beta_yp,h_yp)))
a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
xBC_xm = 0.5*a_xm
xBC_xp = 0.5*a_xp/b_xp
yBC_xm = 0.5*(1.-b_xm)
yBC_xp = 0.5*(1.-old_div(1.,b_xp))
yBC_xp = 0.5*(1.-1./b_xp)
xBC_ym = 0.5*a_ym
xBC_yp = 0.5*a_yp/b_yp
yBC_ym = 0.5*(1.-b_ym)
yBC_yp = 0.5*(1.-old_div(1.,b_yp))
yBC_yp = 0.5*(1.-1./b_yp)
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
@@ -121,33 +120,33 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
a_xm = old_div(gamma_xm,(0.5*alpha_xm-old_div(beta_xm,h_xm)))
b_xm = old_div((0.5*alpha_xm+old_div(beta_xm,h_xm)),(0.5*alpha_xm-old_div(beta_xm,h_xm)))
a_xp = old_div(gamma_xp,(0.5*alpha_xp-old_div(beta_xp,h_xp)))
b_xp = old_div((0.5*alpha_xp+old_div(beta_xp,h_xp)),(0.5*alpha_xp-old_div(beta_xp,h_xp)))
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
a_ym = old_div(gamma_ym,(0.5*alpha_ym-old_div(beta_ym,h_ym)))
b_ym = old_div((0.5*alpha_ym+old_div(beta_ym,h_ym)),(0.5*alpha_ym-old_div(beta_ym,h_ym)))
a_yp = old_div(gamma_yp,(0.5*alpha_yp-old_div(beta_yp,h_yp)))
b_yp = old_div((0.5*alpha_yp+old_div(beta_yp,h_yp)),(0.5*alpha_yp-old_div(beta_yp,h_yp)))
a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
a_zm = old_div(gamma_zm,(0.5*alpha_zm-old_div(beta_zm,h_zm)))
b_zm = old_div((0.5*alpha_zm+old_div(beta_zm,h_zm)),(0.5*alpha_zm-old_div(beta_zm,h_zm)))
a_zp = old_div(gamma_zp,(0.5*alpha_zp-old_div(beta_zp,h_zp)))
b_zp = old_div((0.5*alpha_zp+old_div(beta_zp,h_zp)),(0.5*alpha_zp-old_div(beta_zp,h_zp)))
a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
xBC_xm = 0.5*a_xm
xBC_xp = 0.5*a_xp/b_xp
yBC_xm = 0.5*(1.-b_xm)
yBC_xp = 0.5*(1.-old_div(1.,b_xp))
yBC_xp = 0.5*(1.-1./b_xp)
xBC_ym = 0.5*a_ym
xBC_yp = 0.5*a_yp/b_yp
yBC_ym = 0.5*(1.-b_ym)
yBC_yp = 0.5*(1.-old_div(1.,b_yp))
yBC_yp = 0.5*(1.-1./b_yp)
xBC_zm = 0.5*a_zm
xBC_zp = 0.5*a_zp/b_zp
yBC_zm = 0.5*(1.-b_zm)
yBC_zp = 0.5*(1.-old_div(1.,b_zp))
yBC_zp = 0.5*(1.-1./b_zp)
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
SimPEG/EM/Static/DC/ProblemDC_2D.py
+2 -3
View File
@@ -5,7 +5,6 @@ from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from .SurveyDC import Survey, Survey_ky
@@ -190,8 +189,8 @@ class Problem2D_CC(BaseDCProblem_2D):
MfRhoIDeriv = self.MfRhoIDeriv
rho = self.curModel.rho
if adjoint:
return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(old_div(-1.,rho**2)))*v
return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(old_div(-1.,rho**2)))*v
return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
def getRHS(self, ky):
"""
SimPEG/EM/Static/DC/Utils.py
+1 -2
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
@@ -13,7 +12,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
import SimPEG.EM.Static.DC as DC
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= old_div((nElecs*aSpacing - aSpacing),2)
elocs -= (nElecs*aSpacing - aSpacing)/2
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
SimPEG/EM/Static/SIP/ProblemSIP.py
+1 -2
View File
@@ -6,7 +6,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -114,7 +113,7 @@ class BaseSIPProblem(BaseEMProblem):
JvAll = []
#Assume only eta and tau (eta first then tau)
# v = [2*Mx1]
v = v.reshape((int(old_div(v.size,2)), 2), order='F')
v = v.reshape((v.size//2), 2), order='F')
for tind in range(len(self.survey.times)):
t = self.survey.times[tind]
SimPEG/EM/Static/Utils/StaticUtils.py
+18 -19
View File
@@ -6,7 +6,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import np
from SimPEG.EM.Static import DC, IP
@@ -60,7 +59,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
# Create mid-point location
Cmid = Tx[0]
Pmid = old_div((Rx[0][:,0] + Rx[1][:,0]),2)
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
if DCsurvey.mesh.dim == 2:
zsrc = Tx[1]
elif DCsurvey.mesh.dim ==3:
@@ -73,12 +72,12 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
NB = np.abs(Tx[1][0] - Rx[1][:,0])
# Create mid-point location
Cmid = old_div((Tx[0][0] + Tx[1][0]),2)
Pmid = old_div((Rx[0][:,0] + Rx[1][:,0]),2)
Cmid = (Tx[0][0] + Tx[1][0])/2
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
if DCsurvey.mesh.dim == 2:
zsrc = old_div((Tx[0][1] + Tx[1][1]),2)
zsrc = (Tx[0][1] + Tx[1][1])/2
elif DCsurvey.mesh.dim ==3:
zsrc = old_div((Tx[0][2] + Tx[1][2]),2)
zsrc = (Tx[0][2] + Tx[1][2])/2
# Change output for dtype
if dtype == 'volt':
@@ -94,8 +93,8 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
elif stype == 'dpdp':
leg = data * 2*np.pi / ( old_div(1,MA) - old_div(1,MB) + old_div(1,NB) - old_div(1,NA) )
LEG.append(1./(2*np.pi) *( old_div(1,MA) - old_div(1,MB) + old_div(1,NB) - old_div(1,NA) ))
leg = data * 2*np.pi / (1/MA - 1/MB + 1/NB - 1/NA)
LEG.append(1./(2*np.pi) * (1/MA - 1/MB + 1/NB - 1/NA))
else:
print("""dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """)
break
@@ -103,7 +102,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
if dtype == 'appc':
leg = np.log10(abs(old_div(1.,leg)))
leg = np.log10(abs(1./leg))
rho = np.hstack([rho,leg])
elif dtype == 'appr':
@@ -116,11 +115,11 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
break
midx = np.hstack([midx, old_div(( Cmid + Pmid ),2) ])
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
if DCsurvey.mesh.dim==3:
midz = np.hstack([midz, old_div(-np.abs(Cmid-Pmid),2) + zsrc ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
elif DCsurvey.mesh.dim==2:
midz = np.hstack([midz, old_div(-np.abs(Cmid-Pmid),2) + zsrc ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
ax = axs
# Grid points
@@ -193,10 +192,10 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
# Mesure survey length and direction
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
dl_x = old_div(( endl[1,0] - endl[0,0] ), dl_len)
dl_y = old_div(( endl[1,1] - endl[0,1] ), dl_len)
dl_x = (endl[1,0] - endl[0,0]) / dl_len
dl_y = (endl[1,1] - endl[0,1]) / dl_len
nstn = np.floor( old_div(dl_len, a) )
nstn = np.floor(dl_len / a)
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*a
@@ -239,7 +238,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([np.floor( old_div((AB - b), a) ) , n])
nstn = np.min([(AB - b) // a, n])
# Check if there is enough space, else break the loop
if nstn <= 0:
@@ -284,16 +283,16 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
max_y = endl[1,1] - dl_y * b
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = old_div(box_l,2.)
box_w = box_l / 2.
nstn = np.floor( old_div(box_l, a) )
nstn = np.floor(box_l / a)
# Compute discrete pole location along line
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*a
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
# Define number of cross lines
nlin = int(np.floor( old_div(box_w, a) ))
nlin = int(box_w // a)
lind = list(range(-nlin,nlin+1))
ngrad = nstn * len(lind)
SimPEG/EM/TDEM/TDEM_b.py
+2 -3
View File
@@ -5,7 +5,6 @@ from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from SimPEG.Utils import mkvc, sdiag
import numpy as np
@@ -77,14 +76,14 @@ class ProblemTDEM_b(BaseTDEMProblem):
:return: A
"""
dt = self.timeSteps[tInd]
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (old_div(1.0,dt))*self.MfMui
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1.0/dt)*self.MfMui
def getRHS(self, tInd, F):
dt = self.timeSteps[tInd]
B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T
if B_n.shape[0] is not 1:
raise NotImplementedError('getRHS not implemented for this shape of B_n')
RHS = (old_div(1.0,dt))*self.MfMui*B_n[0,:,:] #TODO: This is a hack
RHS = (1.0/dt)*self.MfMui*B_n[0,:,:] #TODO: This is a hack
return RHS
####################################################
SimPEG/EM/Utils/AnalyticUtils.py
+8 -9
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import *
from scipy.special import ellipk, ellipe
from scipy.constants import mu_0, pi
@@ -57,7 +56,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
dR = obsLoc - srcLoc[i, np.newaxis].repeat(nEdges, axis=0)
mCr = np.cross(m, dR)
r = np.sqrt((dR**2).sum(axis=1))
A[:, i] = +(old_div(mu,(4*pi))) * mCr[:,dimInd]/(r**3)
A[:, i] = +(mu/(4*pi)) * mCr[:,dimInd]/(r**3)
if nSrc == 1:
return A.flatten()
return A
@@ -98,11 +97,11 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0):
dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0)
r = np.sqrt((dR**2).sum(axis=1))
if dimInd == 0:
B[:, i] = +(old_div(mu,(4*pi))) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
elif dimInd == 1:
B[:, i] = +(old_div(mu,(4*pi))) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
elif dimInd == 2:
B[:, i] = +(old_div(mu,(4*pi))) /(r**3) * (3*dR[:,2]**2/r**2-1)
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
else:
raise Exception("Not Implemented")
if nSrc == 1:
@@ -156,7 +155,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
y = obsLoc[:, 1] - srcLoc[i, 1]
z = obsLoc[:, 2] - srcLoc[i, 2]
r = np.sqrt(x**2 + y**2)
m = old_div((4 * radius * r), ((radius + r)**2 + z**2))
m = (4 * radius * r) / ((radius + r)**2 + z**2)
m[m > 1.] = 1.
# m might be slightly larger than 1 due to rounding errors
# but ellipke requires 0 <= m <= 1
@@ -166,11 +165,11 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
# % 1/r singular at r = 0 and K(m) singular at m = 1
Aphi = np.zeros(n)
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(old_div(radius, r[ind])) *((1. - old_div(m[ind], 2.)) * K[ind] - E[ind])
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius/ r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
if component == 'x':
A[ind, i] = Aphi[ind] * (old_div(-y[ind], r[ind]) )
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind])
elif component == 'y':
A[ind, i] = Aphi[ind] * ( old_div(x[ind], r[ind]) )
A[ind, i] = Aphi[ind] * (x[ind] / r[ind])
else:
raise ValueError('Invalid component')
SimPEG/EM/Utils/EMUtils.py
+2 -3
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from scipy.constants import mu_0, epsilon_0
@@ -16,8 +15,8 @@ def omega(freq):
def k(freq, sigma, mu=mu_0, eps=epsilon_0):
""" Eq 1.47 - 1.49 in Ward and Hohmann """
w = omega(freq)
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (old_div(sigma, (eps*w)))**2 ) + 1) )
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (old_div(sigma, (eps*w)))**2 ) - 1) )
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) + 1)
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) - 1)
return alp - 1j*beta
SimPEG/EM/Utils/testingUtils.py
+1 -2
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
@@ -32,7 +31,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
else:
mapping = Maps.ExpMap(mesh)
x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + old_div(cs,4.) #don't sample right by the source, slightly off alignment from either staggered grid
x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0])
if comp[2] == 'r':
SimPEG/Examples/DC_Analytic_Dipole.py
+2 -3
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import *
import SimPEG.EM.Static.DC as DC
@@ -45,7 +44,7 @@ def run(plotIt=True):
rn = (srclocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
return I/(sigma*2.*np.pi)*(old_div(1,rP)-old_div(1,rN))
return I/(sigma*2.*np.pi)*(1/rP-1/rN)
data_anaP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
data_anaN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
@@ -68,7 +67,7 @@ def run(plotIt=True):
ax[0].set_title('Computed')
plt.show()
return old_div(np.linalg.norm(data-data_ana),np.linalg.norm(data_ana))
return np.linalg.norm(data-data_ana) / np.linalg.norm(data_ana)
if __name__ == '__main__':
SimPEG/Examples/DC_Forward_PseudoSection.py
+7 -8
View File
@@ -6,7 +6,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import Mesh, Utils, np, sp
import SimPEG.DCIP as DC
import time
@@ -66,7 +65,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
model[ind] = sig[2]
# Get index of the center
indy = int(old_div(mesh.nCy,2))
indy = int(mesh.nCy // 2)
# Plot the model for reference
# Define core mesh extent
@@ -87,8 +86,8 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
# Define some global geometry
dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
dl_x = old_div(( Tx[-1][0,1] - Tx[0][0,0] ), dl_len)
dl_y = old_div(( Tx[-1][1,1] - Tx[0][1,0] ), dl_len)
dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
dl_y = ( Tx[-1][1,1] - Tx[0][1,0] ) / dl_len
#azm = np.arctan(dl_y/dl_x)
#Set boundary conditions
@@ -98,7 +97,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
# line source for simplicity.
Div = mesh.faceDiv
Grad = mesh.cellGrad
Msig = Utils.sdiag(old_div(1.,(mesh.aveF2CC.T*(old_div(1.,model)))))
Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
A = Div*Msig*Grad
@@ -109,7 +108,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
# We will solve the system iteratively, so a pre-conditioner is helpful
# This is simply a Jacobi preconditioner (inverse of the main diagonal)
dA = A.diagonal()
P = sp.spdiags(old_div(1,dA),0,A.shape[0],A.shape[0])
P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
# Now we can solve the system for all the transmitters
# We want to store the data
@@ -133,10 +132,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
tx = np.squeeze(Tx[ii][:,0:1])
tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( old_div([-1], mesh.vol[inds]) )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1] / mesh.vol[inds])
else:
inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( old_div([-1,1], mesh.vol[inds]) )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1,1] / mesh.vol[inds])
# Iterative Solve
Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
+15 -16
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import *
from SimPEG.EM import FDEM, Analytics, mu_0
import time
@@ -75,8 +74,8 @@ def run(plotIt=True):
casing_l = 300 # length of the casing
casing_r = 0.1
casing_a = casing_r - old_div(casing_t,2.) # inner radius
casing_b = casing_r + old_div(casing_t,2.) # outer radius
casing_a = casing_r - casing_t / 2. # inner radius
casing_b = casing_r + casing_t / 2. # outer radius
casing_z = np.r_[-casing_l,0.]
@@ -86,25 +85,25 @@ def run(plotIt=True):
src_loc = np.r_[0.,0.,dsz]
inf_loc = np.r_[0.,0.,1e4]
print('Skin Depth: ', [(old_div(500.,np.sqrt(sigmaback*_))) for _ in freqs])
print('Skin Depth: ', [(500. / np.sqrt(sigmaback*_)) for _ in freqs])
# ------------------ MESH ------------------
# fine cells near well bore
csx1, csx2 = 2e-3, 60.
pfx1, pfx2 = 1.3, 1.3
ncx1 = np.ceil(old_div(casing_b,csx1)+2)
ncx1 = np.ceil(casing_b/csx1)+2
# pad nicely to second cell size
npadx1 = np.floor(old_div(np.log(old_div(csx2,csx1)), np.log(pfx1)))
npadx1 = np.log(csx2/csx1) // np.log(pfx1)
hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
dx1 = sum(hx1a)+sum(hx1b)
dx1 = np.floor(old_div(dx1,csx2))
hx1b *= old_div((dx1*csx2 - sum(hx1a)),sum(hx1b))
dx1 = dx1 // csx2
hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
# second chunk of mesh
dx2 = 300. # uniform mesh out to here
ncx2 = np.ceil(old_div((dx2 - dx1),csx2))
ncx2 = np.ceil((dx2 - dx1) / csx2)
npadx2 = 45
hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
@@ -112,7 +111,7 @@ def run(plotIt=True):
# z-direction
csz = 0.05
nza = 10
ncz, npadzu, npadzd = np.int(np.ceil(old_div(np.diff(casing_z)[0],csz)))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction
# Mesh
@@ -190,7 +189,7 @@ def run(plotIt=True):
# assemble the source
sg = np.hstack([sg_x,sg_y,sg_z])
sg_p = [FDEM.Src.RawVec_e([],_,old_div(sg,mesh.area)) for _ in freqs]
sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs]
# downhole source
dg_x = np.zeros(mesh.vnF[0],dtype=complex)
@@ -199,7 +198,7 @@ def run(plotIt=True):
# vertically directed wire
dgv_indx = (mesh.gridFz[:,0] < csx1) # go through the center of the well
dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + old_div(csz,2.))
dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.)
dgv_ind = dgv_indx & dgv_indz
dg_z[dgv_ind] = -1.
@@ -221,7 +220,7 @@ def run(plotIt=True):
# assemble the source
dg = np.hstack([dg_x,dg_y,dg_z])
dg_p = [FDEM.Src.RawVec_e([],_,old_div(dg,mesh.area)) for _ in freqs]
dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs]
# ------------ Problem and Survey ---------------
survey = FDEM.Survey(sg_p + dg_p)
@@ -259,9 +258,9 @@ def run(plotIt=True):
in1_in = in1[np.r_[inds]]
z_in = mesh.gridFz[inds_fz,2]
in0_in = in0_in.reshape([old_div(in0_in.shape[0],3),3])
in1_in = in1_in.reshape([old_div(in1_in.shape[0],3),3])
z_in = z_in.reshape([old_div(z_in.shape[0],3),3])
in0_in = in0_in.reshape([in0_in.shape[0]//3,3])
in1_in = in1_in.reshape([in1_in.shape[0]//3,3])
z_in = z_in.reshape([z_in.shape[0]//3,3])
I0 = in0_in.sum(1).real
I1 = in1_in.sum(1).real
SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+1 -2
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import *
from SimPEG.FLOW import Richards
@@ -53,7 +52,7 @@ def run(plotIt=True):
def getFields(timeStep,method):
timeSteps = np.ones(old_div(360,timeStep))*timeStep
timeSteps = np.ones(360/timeStep)*timeStep
prob = Richards.RichardsProblem(M, mapping=E, timeSteps=timeSteps,
boundaryConditions=bc, initialConditions=h,
doNewton=False, method=method)
SimPEG/Examples/Inversion_IRLS.py
+5 -6
View File
@@ -6,7 +6,6 @@ from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
from past.utils import old_div
from SimPEG import *
@@ -56,10 +55,10 @@ def run(N=100, plotIt=True):
# Distance weighting
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( old_div(wr,np.max(wr)) )
wr = ( wr/np.max(wr))
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = old_div(1.,wd)
dmis.Wd = 1./wd
betaest = Directives.BetaEstimate_ByEig()
@@ -68,15 +67,15 @@ def run(N=100, plotIt=True):
reg.cell_weights = wr
reg.mref = np.zeros(mesh.nC)
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()
# 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.]
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])
SimPEG/Examples/MT_1D_ForwardAndInversion.py
+1 -2
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import SimPEG as simpeg
import numpy as np
import SimPEG.MT as MT
@@ -91,7 +90,7 @@ def run(plotIt=True):
std = 0.05 # 5% std
survey.std = np.abs(survey.dobs*std)
# Assign the data weight
Wd = old_div(1.,survey.std)
Wd = 1./survey.std
## Setup the inversion proceedure
# Define a counter
SimPEG/Examples/Maps_Mesh2Mesh.py
+1 -2
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Mesh, Maps, Utils
def run(plotIt=True):
@@ -19,7 +18,7 @@ def run(plotIt=True):
M = Mesh.TensorMesh([100,100])
h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
h1 = old_div(h1,h1.sum())
h1 = h1/h1.sum()
M2 = Mesh.TensorMesh([h1,h1])
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
v = Utils.mkvc(V)
SimPEG/FLOW/Richards/Empirical.py
+12 -13
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from builtins import object
from SimPEG import Mesh, Maps, Utils, np
from future.utils import with_metaclass
@@ -193,8 +192,8 @@ class _haverkamp_theta(NonLinearMap):
def transformDerivU(self, u, m):
self.setModel(m)
g = (self.alpha*(old_div((self.theta_s - self.theta_r),
(self.alpha + abs(u)**self.beta)**2))
g = (self.alpha*((self.theta_s - self.theta_r) /
(self.alpha + abs(u)**self.beta)**2)
*(-self.beta*abs(u)**(self.beta-1)*np.sign(u)))
g[u >= 0] = 0
g = Utils.sdiag(g)
@@ -237,7 +236,7 @@ class _haverkamp_k(NonLinearMap):
def transformDerivU(self, u, m):
self.setModel(m)
g = old_div(-(np.exp(self.Ks)*self.A*self.gamma*abs(u)**(self.gamma-1)*np.sign(u)),((self.A+abs(u)**self.gamma)**2))
g = -(np.exp(self.Ks)*self.A*self.gamma*abs(u)**(self.gamma-1)*np.sign(u))/((self.A+abs(u)**self.gamma)**2)
g[u >= 0] = 0
g = Utils.sdiag(g)
return g
@@ -279,9 +278,9 @@ class _vangenuchten_theta(NonLinearMap):
def transform(self, u, m):
self.setModel(m)
m = 1 - old_div(1.0,self.n)
f = (old_div(( self.theta_s - self.theta_r ),
((1+abs(self.alpha*u)**self.n)**m)) + self.theta_r)
m = 1 - 1.0/self.n
f = (( self.theta_s - self.theta_r ) /
((1+abs(self.alpha*u)**self.n)**m) + self.theta_r)
if Utils.isScalar(self.theta_s):
f[u >= 0] = self.theta_s
else:
@@ -293,7 +292,7 @@ class _vangenuchten_theta(NonLinearMap):
self.setModel(m)
def transformDerivU(self, u, m):
g = -self.alpha*self.n*abs(self.alpha*u)**(self.n - 1)*np.sign(self.alpha*u)*(old_div(1.,self.n) - 1)*(self.theta_r - self.theta_s)*(abs(self.alpha*u)**self.n + 1)**(old_div(1.,self.n) - 2)
g = -self.alpha*self.n*abs(self.alpha*u)**(self.n - 1)*np.sign(self.alpha*u)*(1./self.n - 1)*(self.theta_r - self.theta_s)*(abs(self.alpha*u)**self.n + 1)**(1./self.n - 2)
g[u >= 0] = 0
g = Utils.sdiag(g)
return g
@@ -322,10 +321,10 @@ class _vangenuchten_k(NonLinearMap):
I = self.I
n = self.n
Ks = self.Ks
m = 1.0 - old_div(1.0,n)
m = 1.0 - 1.0/n
theta_e = old_div(1.0,((1.0+abs(alpha*u)**n)**m))
f = np.exp(Ks)*theta_e**I* ( ( 1.0 - ( 1.0 - theta_e**(old_div(1.0,m)) )**m )**2 )
theta_e = 1.0/((1.0+abs(alpha*u)**n)**m)
f = np.exp(Ks)*theta_e**I* ( ( 1.0 - ( 1.0 - theta_e**(1.0/m) )**m )**2 )
if Utils.isScalar(self.Ks):
f[u >= 0] = np.exp(self.Ks)
else:
@@ -348,9 +347,9 @@ class _vangenuchten_k(NonLinearMap):
I = self.I
n = self.n
Ks = self.Ks
m = 1.0 - old_div(1.0,n)
m = 1.0 - 1.0/n
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(old_div(1.0,n) - 1.0)*((abs(alpha*u)**n + 1)**(old_div(1.0,n) - 1))**(I - 1)*((1 - old_div(1.0,((abs(alpha*u)**n + 1)**(old_div(1.0,n) - 1))**(old_div(1.0,(old_div(1.0,n) - 1)))))**(1 - old_div(1.0,n)) - 1)**2*(abs(alpha*u)**n + 1)**(old_div(1.0,n) - 2) - old_div((2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(old_div(1.0,n) - 1)*((abs(alpha*u)**n + 1)**(old_div(1.0,n) - 1))**I*((1 - old_div(1.0,((abs(alpha*u)**n + 1)**(old_div(1.0,n) - 1))**(old_div(1.0,(old_div(1.0,n) - 1)))))**(1 - old_div(1.0,n)) - 1)*(abs(alpha*u)**n + 1)**(old_div(1.0,n) - 2)),(((abs(alpha*u)**n + 1)**(old_div(1.0,n) - 1))**(old_div(1.0,(old_div(1.0,n) - 1)) + 1)*(1 - old_div(1.0,((abs(alpha*u)**n + 1)**(old_div(1.0,n) - 1))**(old_div(1.0,(old_div(1.0,n) - 1)))))**(old_div(1.0,n))))
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1))+ 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)
g[u >= 0] = 0
g = Utils.sdiag(g)
return g
SimPEG/FLOW/Richards/RichardsProblem.py
+7 -8
View File
@@ -5,7 +5,6 @@ from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import *
from .Empirical import RichardsMap
import time
@@ -176,7 +175,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
# DIV*diag(GRAD*hn1+BC*bc)*(AV*(1.0/K))^-1
DdiagGh1 = DIV*Utils.sdiag(GRAD*hn1+BC*bc)
diagAVk2_AVdiagK2 = Utils.sdiag((AV*(old_div(1.,K1)))**(-2)) * AV*Utils.sdiag(K1**(-2))
diagAVk2_AVdiagK2 = Utils.sdiag((AV*(1./K1))**(-2)) * AV*Utils.sdiag(K1**(-2))
# The matrix that we are computing has the form:
#
@@ -188,12 +187,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
# | Asub Adiag | | hn | | bn |
# - - - - - -
Asub = (old_div(-1.0,dt))*dT
Asub = (-1.0/dt)*dT
Adiag = (
(old_div(1.0,dt))*dT1
(1.0/dt)*dT1
-DdiagGh1*diagAVk2_AVdiagK2*dK1
-DIV*Utils.sdiag(old_div(1.,(AV*(old_div(1.,K1)))))*GRAD
-DIV*Utils.sdiag(1./(AV*(1./K1)))*GRAD
-Dz*diagAVk2_AVdiagK2*dK1
)
@@ -223,17 +222,17 @@ class RichardsProblem(Problem.BaseTimeProblem):
K = self.mapping.k(h, m)
dK = self.mapping.kDerivU(h, m)
aveK = old_div(1.,(AV*(old_div(1.,K))))
aveK = 1./(AV*(1./K))
RHS = DIV*Utils.sdiag(aveK)*(GRAD*h+BC*bc) + Dz*aveK
if self.method == 'mixed':
r = old_div((T-Tn),dt) - RHS
r = (T-Tn)/dt - RHS
elif self.method == 'head':
r = dT*(h - hn)/dt - RHS
if not return_g: return r
J = old_div(dT,dt) - DIV*Utils.sdiag(aveK)*GRAD
J = dT/dt - DIV*Utils.sdiag(aveK)*GRAD
if self.doNewton:
DDharmAve = Utils.sdiag(aveK**2)*AV*Utils.sdiag(K**(-2)) * dK
J = J - DIV*Utils.sdiag(GRAD*h + BC*bc)*DDharmAve - Dz*DDharmAve
SimPEG/MT/FieldsMT.py
+6 -7
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
from scipy.constants import mu_0
import sys
@@ -70,7 +69,7 @@ class Fields1D_e(BaseMTFields):
C = self.mesh.nodalGrad
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= old_div(- 1.,(1j*omega(src.freq)))
b[:,i] *= -1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -195,7 +194,7 @@ class Fields3D_e(BaseMTFields):
# adjoint: returns a 2*nE long vector with zero's for py
return np.vstack((v,np.zeros_like(v)))
# Not adjoint: return only the px part of the vector
return v[:old_div(len(v),2)]
return v[:len(v)//2]
def _e_pyDeriv_u(self, src, v, adjoint = False):
'''
@@ -205,7 +204,7 @@ class Fields3D_e(BaseMTFields):
# adjoint: returns a 2*nE long vector with zero's for px
return np.vstack((np.zeros_like(v),v))
# Not adjoint: return only the px part of the vector
return v[old_div(len(v),2)::]
return v[len(v)//2::]
def _e_pxDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
@@ -234,7 +233,7 @@ class Fields3D_e(BaseMTFields):
C = self.mesh.edgeCurl
b = (C * e_pxSolution)
for i, src in enumerate(srcList):
b[:,i] *= old_div(- 1.,(1j*omega(src.freq)))
b[:,i] *= -1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -245,7 +244,7 @@ class Fields3D_e(BaseMTFields):
C = self.mesh.edgeCurl
b = (C * e_pySolution)
for i, src in enumerate(srcList):
b[:,i] *= old_div(- 1.,(1j*omega(src.freq)))
b[:,i] *= -1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -355,4 +354,4 @@ class Fields3D_e(BaseMTFields):
This function stacks the fields derivatives appropriately
"""
# The fields have no dependance to the model.
return None
return None
SimPEG/MT/Problem1D/Probs.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG.EM.Utils import omega
from SimPEG import mkvc
from scipy.constants import mu_0
SimPEG/MT/SrcMT.py
+2 -3
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
@@ -85,7 +84,7 @@ class polxy_1Dprimary(BaseMTSrc):
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(old_div(1,( 1j*omega(self.freq) )))
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
return bBG_bp
def S_e(self,problem):
@@ -162,7 +161,7 @@ class polxy_3Dprimary(BaseMTSrc):
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(old_div(1,( 1j*omega(self.freq) )))
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
return bBG_bp
def S_e(self,problem):
SimPEG/MT/SurveyMT.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Survey as SimPEGsurvey, Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
SimPEG/MT/Utils/MT1Danalytic.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import zip
from past.utils import old_div
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
SimPEG/MT/Utils/MT1Dsolutions.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np, SimPEG as simpeg
from .MT1Danalytic import getEHfields
from scipy.constants import mu_0
SimPEG/MT/Utils/dataUtils.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import division
from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
# Utils used for the data,
import numpy as np, matplotlib.pyplot as plt, sys
import SimPEG as simpeg
SimPEG/MT/Utils/ediFilesUtils.py
-1
View File
@@ -7,7 +7,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import object
from past.utils import old_div
# Functions to import and export MT EDI files.
from SimPEG import mkvc
from scipy.constants import mu_0
SimPEG/MT/Utils/plotDataTypes.py
+20 -21
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from matplotlib import pyplot as plt, colors, numpy as np
@@ -42,7 +41,7 @@ def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm=
cmap = plt.get_cmap('RdYlBu')
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
@@ -80,7 +79,7 @@ def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
@@ -95,7 +94,7 @@ def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
@@ -110,7 +109,7 @@ def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
@@ -181,13 +180,13 @@ def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='Sy
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = old_div(1,(8*np.pi**2*10**(-7)))
treshold = 1.0 # 1 meter
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
freq = array['freq'][indUniSta]
if par == 'abs':
zPlot = np.abs(array[flag][indUniSta])
elif par == 'real':
@@ -205,7 +204,7 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
pColor = 'g'
elif 'xy' in flag:
lab = 'XY'
pColor = 'r'
pColor = 'r'
elif 'yx' in flag:
lab = 'YX'
pColor = 'b'
@@ -214,7 +213,7 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
pColor = 'y'
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
@@ -236,8 +235,8 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
else:
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
elif par == 'ares':
elif par == 'ares':
zPlot = old_div(np.abs(array[flag][indSect])**2,(8*np.pi**2*10**(-7)*array['freq'][indSect]))
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
@@ -278,7 +277,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
elif par == 'imag':
zPlot = np.imag(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
@@ -288,7 +287,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
@@ -329,7 +328,7 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
x1 = arr1['x']
y0 = arr0['freq']
y1 = arr1['freq']
if par == 'abs':
if useLog:
@@ -342,7 +341,7 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zPlot[maskInd] = mask
cmap = plt.get_cmap('RdYlBu')#seismic)
elif par == 'ares':
arF = old_div(1,(8*np.pi**2*10**(-7)))
arF = old_div(1,(8*np.pi**2*10**(-7)))
if useLog:
zPlot = old_div((np.log10((old_div(arF,arr0['freq']))*np.abs(arr0[flag])**2) - np.log10((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2)),np.log10((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2))
else:
@@ -352,7 +351,7 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'aphs':
if useLog:
zPlot = old_div((np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(old_div(180,np.pi))) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))) ),np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))))
@@ -373,7 +372,7 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('Spectral')
elif par == 'imag':
if useLog:
zPlot = old_div((np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag)),np.log10(arr1[flag].imag))
@@ -384,7 +383,7 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
@@ -392,14 +391,14 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
if colorNorm=='SymLog':
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
if cLevel:
level = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),50.))
level = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),50.))
clevel = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),10.))
else:
level = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),50.))
SimPEG/Maps.py
+22 -23
View File
@@ -6,7 +6,6 @@ from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
from past.utils import old_div
from builtins import object
from . import Utils
import numpy as np, scipy.sparse as sp
@@ -259,10 +258,10 @@ class ReciprocalMap(IdentityMap):
"""
def _transform(self, m):
return old_div(1.0, Utils.mkvc(m))
return 1.0 / Utils.mkvc(m)
def inverse(self, D):
return old_div(1.0, Utils.mkvc(m))
return 1.0 / Utils.mkvc(m)
def deriv(self, m):
# TODO: if this is a tensor, you might have a problem.
@@ -302,7 +301,7 @@ class LogMap(IdentityMap):
deriv = np.zeros(mod.shape)
tol = 1e-16 # zero
ind = np.greater_equal(np.abs(mod),tol)
deriv[ind] = old_div(1.0,mod[ind])
deriv[ind] = 1.0 / mod[ind]
return Utils.sdiag(deriv)
def inverse(self, m):
@@ -606,14 +605,14 @@ class ComplexMap(IdentityMap):
@property
def shape(self):
return (old_div(self.nP,2),self.nP)
return (self.nP // 2,self.nP)
def _transform(self, m):
nC = self.mesh.nC
return m[:nC] + m[nC:]*1j
def deriv(self, m):
nC = old_div(self.nP,2)
nC = self.nP//2
shp = (nC, nC*2)
def fwd(v):
return v[:nC] + v[nC:]*1j
@@ -658,7 +657,7 @@ class CircleMap(IdentityMap):
sig1, sig2 = np.exp(sig1), np.exp(sig2)
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
return sig1 + (sig2 - sig1)*(old_div(np.arctan(a*(np.sqrt((X-x)**2 + (Y-y)**2) - r)),np.pi) + 0.5)
return sig1 + (sig2 - sig1)*(np.arctan(a*(np.sqrt((X-x)**2 + (Y-y)**2) - r))/np.pi + 0.5)
def deriv(self, m):
a = self.slope
@@ -668,11 +667,11 @@ class CircleMap(IdentityMap):
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
if self.logSigma:
g1 = -(old_div(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2))),np.pi) + 0.5)*sig1 + sig1
g2 = (old_div(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2))),np.pi) + 0.5)*sig2
g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)*sig1 + sig1
g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)*sig2
else:
g1 = -(old_div(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2))),np.pi) + 0.5) + 1.0
g2 = (old_div(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2))),np.pi) + 0.5)
g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5) + 1.0
g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)
g3 = a*(-X + x)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2))
g4 = a*(-Y + y)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2))
g5 = -a*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1))
@@ -761,7 +760,7 @@ class PolyMap(IdentityMap):
raise Exception
return sig1+(sig2-sig1)*(old_div(np.arctan(alpha*f),np.pi)+0.5)
return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
def deriv(self, m):
alpha = self.slope
@@ -800,11 +799,11 @@ class PolyMap(IdentityMap):
raise Exception
if self.logSigma:
g1 = -(old_div(np.arctan(alpha*f),np.pi) + 0.5)*sig1 + sig1
g2 = (old_div(np.arctan(alpha*f),np.pi) + 0.5)*sig2
g1 = -(np.arctan(alpha*f)/np.pi + 0.5)*sig1 + sig1
g2 = (np.arctan(alpha*f)/np.pi + 0.5)*sig2
else:
g1 = -(old_div(np.arctan(alpha*f),np.pi) + 0.5) + 1.0
g2 = (old_div(np.arctan(alpha*f),np.pi) + 0.5)
g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
g2 = (np.arctan(alpha*f)/np.pi + 0.5)
g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V
@@ -897,7 +896,7 @@ class SplineMap(IdentityMap):
raise Exception
return sig1+(sig2-sig1)*(old_div(np.arctan(alpha*f),np.pi)+0.5)
return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
def deriv(self, m):
alpha = self.slope
@@ -931,11 +930,11 @@ class SplineMap(IdentityMap):
raise Exception
if self.logSigma:
g1 = -(old_div(np.arctan(alpha*f),np.pi) + 0.5)*sig1 + sig1
g2 = (old_div(np.arctan(alpha*f),np.pi) + 0.5)*sig2
g1 = -(np.arctan(alpha*f)/np.pi + 0.5)*sig1 + sig1
g2 = (np.arctan(alpha*f)/np.pi + 0.5)*sig2
else:
g1 = -(old_div(np.arctan(alpha*f),np.pi) + 0.5) + 1.0
g2 = (old_div(np.arctan(alpha*f),np.pi) + 0.5)
g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
g2 = (np.arctan(alpha*f)/np.pi + 0.5)
if self.mesh.dim ==2:
@@ -954,7 +953,7 @@ class SplineMap(IdentityMap):
cb[i] = ctemp-dy
spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
fderiv = old_div((spla(X)-splb(X)),(2*dy))
fderiv = (spla(X)-splb(X))/(2*dy)
g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
elif self.mesh.dim==3:
@@ -981,7 +980,7 @@ class SplineMap(IdentityMap):
spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X
flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
fderiv = old_div((flinesa-flinesb),(2*dy))
fderiv = (flinesa-flinesb)/(2*dy)
g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
else :
raise Exception
SimPEG/Mesh/CurvilinearMesh.py
+12 -13
View File
@@ -5,7 +5,6 @@ from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import Utils, np
from .BaseMesh import BaseRectangularMesh
from .DiffOperators import DiffOperators
@@ -16,8 +15,8 @@ from future.utils import with_metaclass
# Some helper functions.
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
normalize2D = lambda x: old_div(x,np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)))
normalize3D = lambda x: old_div(x,np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)))
normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (BaseRectangularMesh, DiffOperators, InnerProducts, CurvView), {}))):
@@ -227,7 +226,7 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
Utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
Utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
self._vol = old_div((vol1 + vol2),2)
self._vol = (vol1 + vol2)/2
return self._vol
return locals()
_vol = None
@@ -289,16 +288,16 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
if self.dim == 2:
return normalize2D(np.r_[self._normals[0], self._normals[1]])
elif self.dim == 3:
normal1 = old_div((self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3]),4)
normal2 = old_div((self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3]),4)
normal3 = old_div((self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3]),4)
normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4
normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
return normalize3D(np.r_[normal1, normal2, normal3])
return locals()
_normals = None
normals = property(**normals())
def edge():
doc = "Edge legnths."
doc = "Edge lengths."
def fget(self):
if(self._edge is None or self._tangents is None):
@@ -309,7 +308,7 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
edge2 = xy[B, :] - xy[A, :]
self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
self._tangents = old_div(np.r_[edge1, edge2],np.c_[self._edge, self._edge])
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
elif(self.dim == 3):
xyz = self.gridN
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz]))
@@ -319,7 +318,7 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz]))
edge3 = xyz[E, :] - xyz[A, :]
self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
self._tangents = old_div(np.r_[edge1, edge2, edge3],np.c_[self._edge, self._edge, self._edge])
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
return self._edge
return locals()
_edge = None
@@ -340,10 +339,10 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
if __name__ == '__main__':
nc = 5
h1 = np.cumsum(np.r_[0, old_div(np.ones(nc),(nc))])
h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
nc = 7
h2 = np.cumsum(np.r_[0, old_div(np.ones(nc),(nc))])
h3 = np.cumsum(np.r_[0, old_div(np.ones(nc),(nc))])
h2 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
dee3 = True
if dee3:
X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False)
SimPEG/Mesh/CylMesh.py
+8 -9
View File
@@ -4,7 +4,6 @@ from __future__ import division
from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import scipy.sparse as sp
from scipy.constants import pi
@@ -75,7 +74,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
Number of x-faces in each direction
:rtype: numpy.array
:rtype: numpy.array
:return: vnFx, (dim, )
"""
return self.vnC
@@ -85,7 +84,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
Number of y-edges in each direction
:rtype: numpy.array
:rtype: numpy.array
:return: vnEy or None if dim < 2, (dim, )
"""
nNx = self.nNx if self.isSymmetric else self.nNx - 1
@@ -96,7 +95,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
Number of z-edges in each direction
:rtype: numpy.array
:rtype: numpy.array
:return: vnEz or None if nCy > 1, (dim, )
"""
if self.isSymmetric:
@@ -200,7 +199,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
D1 = kron3(speye(self.nCz), speye(self.nCy), ddx(self.nCx)[:,1:])
S = self.r(self.area, 'F', 'Fx', 'V')
V = self.vol
self._faceDivx = sdiag(old_div(1,V))*D1*sdiag(S)
self._faceDivx = sdiag(1/V)*D1*sdiag(S)
return self._faceDivx
@property
@@ -212,7 +211,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
D2 = kron3(speye(self.nCz), ddx(self.nCy), speye(self.nCx))
S = self.r(self.area, 'F', 'Fy', 'V')
V = self.vol
self._faceDivy = sdiag(old_div(1,V))*D2*sdiag(S)
self._faceDivy = sdiag(1/V)*D2*sdiag(S)
return self._faceDivy
@property
@@ -222,7 +221,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
D3 = kron3(ddx(self.nCz), speye(self.nCy), speye(self.nCx))
S = self.r(self.area, 'F', 'Fz', 'V')
V = self.vol
self._faceDivz = sdiag(old_div(1,V))*D3*sdiag(S)
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
return self._faceDivz
@@ -261,7 +260,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
A = self.area
E = self.edge
#Edge curl operator
self._edgeCurl = sdiag(old_div(1,A))*sp.vstack((Dz, Dr))*sdiag(E)
self._edgeCurl = sdiag(1/A)*sp.vstack((Dz, Dr))*sdiag(E)
return self._edgeCurl
# @property
@@ -362,7 +361,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
grid = getattr(Mrect, 'grid' + locTypeTo)
# This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + old_div(np.pi,2)
theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
theta[theta < 0] += np.pi*2.0
r = ((grid[:,0] - self.cartesianOrigin[0])**2 + (grid[:,1] - self.cartesianOrigin[1])**2)**0.5
SimPEG/Mesh/DiffOperators.py
+20 -21
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import object
from past.utils import old_div
import numpy as np
from scipy import sparse as sp
from SimPEG.Utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap
@@ -153,7 +152,7 @@ class DiffOperators(object):
# Compute areas of cell faces & volumes
S = self.area
V = self.vol
self._faceDiv = sdiag(old_div(1,V))*D*sdiag(S)
self._faceDiv = sdiag(1/V)*D*sdiag(S)
return self._faceDiv
return locals()
@@ -177,7 +176,7 @@ class DiffOperators(object):
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F', 'Fx', 'V')
V = self.vol
self._faceDivx = sdiag(old_div(1,V))*D1*sdiag(S)
self._faceDivx = sdiag(1/V)*D1*sdiag(S)
return self._faceDivx
return locals()
@@ -200,7 +199,7 @@ class DiffOperators(object):
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F', 'Fy', 'V')
V = self.vol
self._faceDivy = sdiag(old_div(1,V))*D2*sdiag(S)
self._faceDivy = sdiag(1/V)*D2*sdiag(S)
return self._faceDivy
return locals()
@@ -220,7 +219,7 @@ class DiffOperators(object):
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F', 'Fz', 'V')
V = self.vol
self._faceDivz = sdiag(old_div(1,V))*D3*sdiag(S)
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
return self._faceDivz
return locals()
@@ -248,7 +247,7 @@ class DiffOperators(object):
G = sp.vstack((D1, D2, D3), format="csr")
# Compute lengths of cell edges
L = self.edge
self._nodalGrad = sdiag(old_div(1,L))*G
self._nodalGrad = sdiag(1/L)*G
return self._nodalGrad
return locals()
_nodalGrad = None
@@ -264,18 +263,18 @@ class DiffOperators(object):
n = self.vnC
# Compute divergence operator on faces
if(self.dim == 1):
D1 = sdiag(old_div(1.,self.hx)) * ddx(mesh.nCx)
D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
L = - D1.T*D1
elif(self.dim == 2):
D1 = sdiag(old_div(1.,self.hx)) * ddx(n[0])
D2 = sdiag(old_div(1.,self.hy)) * ddx(n[1])
D1 = sdiag(1./self.hx) * ddx(n[0])
D2 = sdiag(1./self.hy) * ddx(n[1])
L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
L = L1 + L2
elif(self.dim == 3):
D1 = sdiag(old_div(1.,self.hx)) * ddx(n[0])
D2 = sdiag(old_div(1.,self.hy)) * ddx(n[1])
D3 = sdiag(old_div(1.,self.hz)) * ddx(n[2])
D1 = sdiag(1./self.hx) * ddx(n[0])
D2 = sdiag(1./self.hy) * ddx(n[1])
D3 = sdiag(1./self.hz) * ddx(n[2])
L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
@@ -340,7 +339,7 @@ class DiffOperators(object):
# Compute areas of cell faces & volumes
S = self.area
V = self.aveCC2F*self.vol # Average volume between adjacent cells
self._cellGrad = sdiag(old_div(S,V))*G
self._cellGrad = sdiag(S/V)*G
return self._cellGrad
return locals()
_cellGrad = None
@@ -367,7 +366,7 @@ class DiffOperators(object):
# Compute areas of cell faces & volumes
S = self.area
V = self.aveCC2F*self.vol # Average volume between adjacent cells
self._cellGradBC = sdiag(old_div(S,V))*G
self._cellGradBC = sdiag(S/V)*G
return self._cellGradBC
return locals()
_cellGradBC = None
@@ -393,7 +392,7 @@ class DiffOperators(object):
G1 = self._cellGradxStencil()
# Compute areas of cell faces & volumes
V = self.aveCC2F*self.vol
L = self.r(old_div(self.area,V), 'F','Fx', 'V')
L = self.r(self.area/V, 'F','Fx', 'V')
self._cellGradx = sdiag(L)*G1
return self._cellGradx
return locals()
@@ -417,7 +416,7 @@ class DiffOperators(object):
G2 = self._cellGradyStencil()
# Compute areas of cell faces & volumes
V = self.aveCC2F*self.vol
L = self.r(old_div(self.area,V), 'F','Fy', 'V')
L = self.r(self.area/V, 'F','Fy', 'V')
self._cellGrady = sdiag(L)*G2
return self._cellGrady
return locals()
@@ -438,7 +437,7 @@ class DiffOperators(object):
G3 = self._cellGradzStencil()
# Compute areas of cell faces & volumes
V = self.aveCC2F*self.vol
L = self.r(old_div(self.area,V), 'F','Fz', 'V')
L = self.r(self.area/V, 'F','Fz', 'V')
self._cellGradz = sdiag(L)*G3
return self._cellGradz
return locals()
@@ -465,7 +464,7 @@ class DiffOperators(object):
D21 = sp.kron(ddx(n[1]), speye(n[0]))
D12 = sp.kron(speye(n[1]), ddx(n[0]))
C = sp.hstack((-D21, D12), format="csr")
self._edgeCurl = C*sdiag(old_div(1,S))
self._edgeCurl = C*sdiag(1/S)
elif self.dim == 3:
@@ -484,7 +483,7 @@ class DiffOperators(object):
sp.hstack((D31, O2, -D13)),
sp.hstack((-D21, D12, O3))), format="csr")
self._edgeCurl = sdiag(old_div(1,S))*(C*sdiag(L))
self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
return self._edgeCurl
return locals()
@@ -663,7 +662,7 @@ class DiffOperators(object):
elif(self.dim == 2):
return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
elif(self.dim == 3):
return (old_div(1.,3.))*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
@property
def aveF2CCV(self):
@@ -735,7 +734,7 @@ class DiffOperators(object):
elif(self.dim == 2):
return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr")
elif(self.dim == 3):
return (old_div(1.,3))*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
@property
def aveE2CCV(self):
SimPEG/Mesh/TensorMesh.py
+4 -5
View File
@@ -7,7 +7,6 @@ from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
from past.utils import old_div
from SimPEG import Utils, np, sp
from .BaseMesh import BaseMesh, BaseRectangularMesh
from .View import TensorView
@@ -307,7 +306,7 @@ class BaseTensorMesh(with_metaclass(Utils.SimPEGMetaClass, BaseMesh)):
prop = np.ones(self.nC)
if invProp:
prop = old_div(1.,prop)
prop = 1./prop
if Utils.isScalar(prop):
prop = prop*np.ones(self.nC)
@@ -352,7 +351,7 @@ class BaseTensorMesh(with_metaclass(Utils.SimPEGMetaClass, BaseMesh)):
if not invMat and not invProp:
dMdprop = self.dim * Av.T * V * ones
elif invMat and invProp:
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(old_div(1.,prop**2))
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2)
if tensorType == 1:
Av = getattr(self, 'ave'+projType+'2CC')
@@ -360,7 +359,7 @@ class BaseTensorMesh(with_metaclass(Utils.SimPEGMetaClass, BaseMesh)):
if not invMat and not invProp:
dMdprop = self.dim * Av.T * V
elif invMat and invProp:
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(old_div(1.,prop**2))
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
if tensorType == 2: # anisotropic
Av = getattr(self, 'ave'+projType+'2CCV')
@@ -368,7 +367,7 @@ class BaseTensorMesh(with_metaclass(Utils.SimPEGMetaClass, BaseMesh)):
if not invMat and not invProp:
dMdprop = Av.T * V
elif invMat and invProp:
dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(old_div(1.,prop**2))
dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
if dMdprop is not None:
def innerProductDeriv(v=None):
SimPEG/Mesh/TreeMesh.py
+28 -29
View File
@@ -10,7 +10,6 @@ from builtins import str
from builtins import zip
from builtins import range
from builtins import object
from past.utils import old_div
# ___ ___ ___ ___ ___
# /\ \ ___ /\__\ /\ \ /\ \ /\ \
# /::\ \ /\ \ /::| | /::\ \ /::\ \ /::\ \
@@ -178,7 +177,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
@property
def fill(self):
"""How filled is the mesh compared to a TensorMesh? As a fraction: [0,1]."""
return old_div(float(self.nC),((2**self.maxLevel)**self.dim))
return float(self.nC)/((2**self.maxLevel)**self.dim)
@property
def maxLevel(self):
@@ -590,7 +589,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
def _cellC(self, p):
"""Cell center of a single cell (without origin correction), given a pointer."""
return (old_div(np.array(self._cellH(p)),2.0) + self._cellN(p)).tolist()
return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist()
def _levelWidth(self, level):
return 2**(self.levels - level)
@@ -878,10 +877,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
p = self._pointer(fx)
n, h = self._cellN(p), self._cellH(p)
if self.dim == 2:
gridFx.append( [n[0], n[1] + old_div(h[1],2.0)] )
gridFx.append( [n[0], n[1] + h[1]/2.0] )
areaFx.append( h[1] )
elif self.dim == 3:
gridFx.append( [n[0], n[1] + old_div(h[1],2.0), n[2] + old_div(h[2],2.0)] )
gridFx.append( [n[0], n[1] + h[1]/2.0, n[2] + h[2]/2.0] )
areaFx.append( h[1]*h[2] )
self._gridFx = np.array(gridFx)
self._areaFxFull = np.array(areaFx)
@@ -894,10 +893,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
p = self._pointer(fy)
n, h = self._cellN(p), self._cellH(p)
if self.dim == 2:
gridFy.append( [n[0] + old_div(h[0],2.0), n[1]] )
gridFy.append( [n[0] + h[0]/2.0, n[1]] )
areaFy.append( h[0] )
elif self.dim == 3:
gridFy.append( [n[0] + old_div(h[0],2.0), n[1], n[2] + old_div(h[2],2.0)] )
gridFy.append( [n[0] + h[0]/2.0, n[1], n[2] + h[2]/2.0] )
areaFy.append( h[0]*h[2] )
self._gridFy = np.array(gridFy)
self._areaFyFull = np.array(areaFy)
@@ -913,7 +912,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
self._fz2i[fz] = ii
p = self._pointer(fz)
n, h = self._cellN(p), self._cellH(p)
gridFz.append( [n[0] + old_div(h[0],2.0), n[1] + old_div(h[1],2.0), n[2]] )
gridFz.append( [n[0] + h[0]/2.0, n[1] + h[1]/2.0, n[2]] )
areaFz.append(h[0]*h[1])
self._gridFz = np.array(gridFz)
self._areaFzFull = np.array(areaFz)
@@ -934,7 +933,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
self._ex2i[ex] = ii
p = self._pointer(ex)
n, h = self._cellN(p), self._cellH(p)
gridEx.append( [n[0] + old_div(h[0],2.0), n[1], n[2]] )
gridEx.append( [n[0] + h[0]/2.0, n[1], n[2]] )
edgeEx.append( h[0] )
self._gridEx = np.array(gridEx)
self._edgeExFull = np.array(edgeEx)
@@ -946,7 +945,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
self._ey2i[ey] = ii
p = self._pointer(ey)
n, h = self._cellN(p), self._cellH(p)
gridEy.append( [n[0], n[1] + old_div(h[1],2.0), n[2]] )
gridEy.append( [n[0], n[1] + h[1]/2.0, n[2]] )
edgeEy.append( h[1] )
self._gridEy = np.array(gridEy)
self._edgeEyFull = np.array(edgeEy)
@@ -958,7 +957,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
self._ez2i[ez] = ii
p = self._pointer(ez)
n, h = self._cellN(p), self._cellH(p)
gridEz.append( [n[0], n[1], n[2] + old_div(h[2],2.0)] )
gridEz.append( [n[0], n[1], n[2] + h[2]/2.0] )
edgeEz.append( h[2] )
self._gridEz = np.array(gridEz)
self._edgeEzFull = np.array(edgeEz)
@@ -998,12 +997,12 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
chy1 = self._cellH([p[0] , p[1] + w, sl])[1]
A = (chy0 + chy1)
self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], old_div(chy0, A)], )
self._hangingFx[self._fx2i[self._index([p[0] , p[1] + w, sl])]] = ([self._fx2i[fx], old_div(chy1, A)], )
self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], chy0 / A], )
self._hangingFx[self._fx2i[self._index([p[0] , p[1] + w, sl])]] = ([self._fx2i[fx], chy1 / A], )
n0, n1 = fx, self._index([p[0], p[1] + 2*w, p[-1]])
self._hangingN[self._n2i[test ]] = ([self._n2i[n0], 1.0], )
self._hangingN[self._n2i[self._index([p[0] , p[1] + w, sl])]] = ([self._n2i[n0], 1.0 - old_div(chy0, A)], [self._n2i[n1], 1.0 - old_div(chy1, A)])
self._hangingN[self._n2i[self._index([p[0] , p[1] + w, sl])]] = ([self._n2i[n0], 1.0 - chy0 / A], [self._n2i[n1], 1.0 - chy1 / A])
self._hangingN[self._n2i[self._index([p[0] , p[1] + 2*w, sl])]] = ([self._n2i[n1], 1.0], )
elif self.dim == 3:
@@ -1094,8 +1093,8 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
chx0 = self._cellH([p[0] , p[1] , sl])[0]
chx1 = self._cellH([p[0] + w, p[1] , sl])[0]
self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], old_div(chx0, (chx0 + chx1))], )
self._hangingFy[self._fy2i[self._index([p[0] + w, p[1] , sl])]] = ([self._fy2i[fy], old_div(chx1, (chx0 + chx1))], )
self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], chx0 / (chx0 + chx1)], )
self._hangingFy[self._fy2i[self._index([p[0] + w, p[1] , sl])]] = ([self._fy2i[fy], chx1 / (chx0 + chx1)], )
n0, n1 = fy, self._index([p[0] + 2*w, p[1], p[-1]])
self._hangingN[self._n2i[test ]] = ([self._n2i[n0], 1.0], )
@@ -1356,7 +1355,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
S = np.r_[self._areaFxFull, self._areaFyFull]
elif self.dim == 3:
S = np.r_[self._areaFxFull, self._areaFyFull, self._areaFzFull]
self._faceDiv = Utils.sdiag(old_div(1.0,VOL))*D*Utils.sdiag(S)*R
self._faceDiv = Utils.sdiag(1.0/VOL)*D*Utils.sdiag(S)*R
return self._faceDiv
@property
@@ -1431,12 +1430,12 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
Rf = self._deflationMatrix('F', withHanging=True, asOnes=False)
Re = self._deflationMatrix('E')
Rf_ave = Utils.sdiag(old_div(1.,Rf.sum(axis=0))) * Rf.T
Rf_ave = Utils.sdiag(1./Rf.sum(axis=0)) * Rf.T
C = sp.csr_matrix((V,(I,J)), shape=(self.ntF, self.ntE))
S = np.r_[self._areaFxFull, self._areaFyFull, self._areaFzFull]
L = np.r_[self._edgeExFull, self._edgeEyFull, self._edgeEzFull]
self._edgeCurl = Rf_ave*Utils.sdiag(old_div(1.0,S))*C*Utils.sdiag(L)*Re
self._edgeCurl = Rf_ave*Utils.sdiag(1.0/S)*C*Utils.sdiag(L)*Re
return self._edgeCurl
@property
@@ -1495,9 +1494,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
Rn = self._deflationMatrix('N')
Re = self._deflationMatrix('E', withHanging=True, asOnes=False)
Re_ave = Utils.sdiag(old_div(1.,Re.sum(axis=0))) * Re.T
Re_ave = Utils.sdiag(1./Re.sum(axis=0)) * Re.T
self._nodalGrad = Re_ave*Utils.sdiag(old_div(1,L))*G*Rn
self._nodalGrad = Re_ave*Utils.sdiag(1/L)*G*Rn
return self._nodalGrad
@property
@@ -1509,7 +1508,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
raise Exception('aveEx2CC not implemented in 2D')
if self.dim == 3:
PM = [old_div(1.,4.)]*4
PM = [0.25]*4
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -1543,7 +1542,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
raise NotImplementedError('aveEy2CC not implemented in 2D')
if self.dim == 3:
PM = [old_div(1.,4.)]*4 # plus / plus
PM = [0.25]*4 # plus / plus
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -1578,7 +1577,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
raise Exception('There are no z edges in 2D')
if self.dim == 3:
PM = [old_div(1.,4.)]*4 # plus / plus
PM = [0.25]*4 # plus / plus
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -1629,7 +1628,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
def aveFx2CC(self):
if getattr(self, '_aveFx2CC', None) is None:
I, J, V = [], [], []
PM = [old_div(1.,2.)]*self.dim # 0.5, 0.5
PM = [0.5]*self.dim # 0.5, 0.5
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -1662,7 +1661,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
def aveFy2CC(self):
if getattr(self, '_aveFy2CC', None) is None:
I, J, V = [], [], []
PM = [old_div(1.,2.)]*2 # 0.5, 0.5
PM = [0.5]*2 # 0.5, 0.5
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -1694,7 +1693,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
def aveFz2CC(self):
if getattr(self, '_aveFz2CC', None) is None:
I, J, V = [], [], []
PM = [old_div(1.,2.)]*2 # 0.5, 0.5
PM = [0.5]*2 # 0.5, 0.5
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -1742,7 +1741,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
def aveN2CC(self):
if getattr(self, '_aveN2CC', None) is None:
I, J, V = [], [], []
PM = [old_div(1.,2.**self.dim)] * 2**self.dim
PM = [1./2.**self.dim] * 2**self.dim
for ii, ind in enumerate(self._sortedCells):
p = self._pointer(ind)
@@ -2165,7 +2164,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
import matplotlib.cm as cmx
szSliceDim = len(getattr(self, 'h'+normal.lower())) #: Size of the sliced dimension
if ind is None: ind = int(old_div(szSliceDim,2))
if ind is None: ind = int(szSliceDim // 2)
assert type(ind) in [int, int], 'ind must be an integer'
indLoc = getattr(self,'vectorCC'+normal.lower())[ind]
normalInd = {'X':0,'Y':1,'Z':2}[normal]
SimPEG/Mesh/View.py
+4 -5
View File
@@ -8,7 +8,6 @@ standard_library.install_aliases()
from builtins import zip
from builtins import range
from builtins import object
from past.utils import old_div
import numpy as np
from SimPEG.Utils import mkvc
try:
@@ -139,7 +138,7 @@ class TensorView(object):
# determine number oE slices in x and y dimension
nX = np.ceil(np.sqrt(self.nCz))
nY = np.ceil(old_div(self.nCz,nX))
nY = np.ceil(self.nCz/nX)
# allocate space for montage
nCx = self.nCx
@@ -239,7 +238,7 @@ class TensorView(object):
assert type(grid) is bool, 'grid must be a boolean'
szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension
if ind is None: ind = int(old_div(szSliceDim,2))
if ind is None: ind = int(szSliceDim//2)
assert type(ind) in [int, int], 'ind must be an integer'
assert not (v.dtype == complex and view == 'vec'), 'Can not plot a complex vector.'
@@ -373,8 +372,8 @@ class TensorView(object):
# spaced vectors at the moment. So we will
# Interpolate down to a regular mesh at the
# smallest mesh size in this 2D slice.
nxi = int(old_div(self.hx.sum(),self.hx.min()))
nyi = int(old_div(self.hy.sum(),self.hy.min()))
nxi = int(self.hx.sum()/self.hx.min())
nyi = int(self.hy.sum()/self.hy.min())
tMi = self.__class__([np.ones(nxi)*self.hx.sum()/nxi,
np.ones(nyi)*self.hy.sum()/nyi], self.x0)
P = self.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True)
SimPEG/Optimization.py
+5 -6
View File
@@ -5,7 +5,6 @@ from __future__ import unicode_literals
from builtins import super
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from builtins import object
from . import Utils
import numpy as np, scipy.sparse as sp
@@ -703,10 +702,10 @@ class BFGS(Minimize, Remember):
d = self.bfgsH0 * d #Assume that bfgsH0 is a SimPEG.Solver
else:
khat = 0 if nn is 0 else np.mod(n-nn+k,nn)
gamma = old_div(np.vdot(S[:,khat],d),np.vdot(Y[:,khat],S[:,khat]))
gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])
d = d - gamma*Y[:,khat]
d = self.bfgsrec(k-1,n,nn,S,Y,d)
d = d + (gamma - old_div(np.vdot(Y[:,khat],d),np.vdot(Y[:,khat],S[:,khat])))*S[:,khat]
d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat]
return d
def findSearchDirection(self):
@@ -984,7 +983,7 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
if cgiter == 1:
pc = dc
else:
betak = old_div(rd, rdlast)
betak = rd / rdlast
pc = dc + betak * pc
# Form product Hessian*pc.
@@ -992,12 +991,12 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
Hp = (1-Active)*Hp
# Update delx and residual.
alphak = old_div(rd, np.dot(pc, Hp))
alphak = rd / np.dot(pc, Hp)
delx = delx + alphak*pc
resid = resid - alphak*Hp
rdlast = rd
if np.logical_or(old_div(norm(resid),normResid0) <= self.tolCG, cgiter == self.maxIterCG):
if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG):
cgFlag = 1
# End CG Iterations
SimPEG/Regularization.py
+15 -16
View File
@@ -5,7 +5,6 @@ from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import object
from past.utils import old_div
from . import Utils, Maps, Mesh
import numpy as np
import scipy.sparse as sp
@@ -143,7 +142,7 @@ class RegularizationMesh(object):
:return: averaging matrix from active x-faces to active cell centers
"""
if getattr(self, '_aveCC2Fx', None) is None:
self._aveCC2Fx = Utils.sdiag(old_div(1.,(self.aveFx2CC.T).sum(1))) * self.aveFx2CC.T
self._aveCC2Fx = Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T
return self._aveCC2Fx
@property
@@ -165,7 +164,7 @@ class RegularizationMesh(object):
:return: averaging matrix from active y-faces to active cell centers
"""
if getattr(self, '_aveCC2Fy', None) is None:
self._aveCC2Fy = Utils.sdiag(old_div(1.,(self.aveFy2CC.T).sum(1))) * self.aveFy2CC.T
self._aveCC2Fy = Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
return self._aveCC2Fy
@property
@@ -187,7 +186,7 @@ class RegularizationMesh(object):
:return: averaging matrix from active z-faces to active cell centers
"""
if getattr(self, '_aveCC2Fz', None) is None:
self._aveCC2Fz = Utils.sdiag(old_div(1.,(self.aveFz2CC.T).sum(1))) * self.aveFz2CC.T
self._aveCC2Fz = Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T
return self._aveCC2Fz
@property
@@ -899,28 +898,28 @@ class Tikhonov(Simple):
class Sparse(Simple):
"""
The regularization is:
.. math::
R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
where the IRLS weight
.. math::
R = \eta TO FINISH LATER!!!
So the derivative is straight forward:
.. math::
R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
The IRLS weights are recomputed after each beta solves.
It is strongly recommended to do a few Gauss-Newton iterations
before updating.
"""
# set default values
eps_p = 1e-1 # Threshold value for the model norm
eps_q = 1e-1 # Threshold value for the model gradient norm
@@ -1003,7 +1002,7 @@ class Sparse(Simple):
def R(self, f_m , eps, exponent):
# Eta scaling is important for mix-norms...do not mess with it
eta = (eps**(1.-old_div(exponent,2.)))**0.5
r = old_div(eta, (f_m**2.+ eps**2.)**(old_div((1.-old_div(exponent,2.)),2.)))
eta = (eps**(1.-exponent/2.))**0.5
r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
return r
SimPEG/Tests.py
+8 -9
View File
@@ -6,7 +6,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
from numpy.linalg import norm
from SimPEG.Utils import mkvc, sdiag, diagEst
@@ -102,7 +101,7 @@ class OrderTest(unittest.TestCase):
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
h = [old_div(hi,np.sum(hi)) for hi in [h1, h2, h3]] # normalize
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
else:
raise Exception('Unexpected meshType')
@@ -139,7 +138,7 @@ class OrderTest(unittest.TestCase):
elif self.meshDimension == 3:
X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd)
self.M = CurvilinearMesh([X, Y, Z])
return old_div(1.,nc)
return 1./nc
elif 'Tree' in self._meshType:
nc *= 2
@@ -149,11 +148,11 @@ class OrderTest(unittest.TestCase):
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
h = [old_div(hi,np.sum(hi)) for hi in [h1, h2, h3]] # normalize
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
else:
raise Exception('Unexpected meshType')
levels = int(old_div(np.log(nc),np.log(2)))
levels = int(np.log(nc) / np.log(2))
self.M = Tree(h[:self.meshDimension], levels=levels)
def function(cell):
if 'notatree' in self._meshType:
@@ -211,8 +210,8 @@ class OrderTest(unittest.TestCase):
print('~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~')
print('%4i | %8.2e |' % (nc, err))
else:
order.append(old_div(np.log(old_div(err,err_old)),np.log(old_div(max_h,max_h_old))))
print('%4i | %8.2e | %6.4f | %6.4f' % (nc, err, old_div(err_old,err), order[-1]))
order.append(np.log(err/err_old)/np.log(max_h/max_h_old))
print('%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]))
err_old = err
max_h_old = max_h
print('---------------------------------------------')
@@ -297,8 +296,8 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
# We assume it is a numpy.ndarray
E1[i] = l2norm( ft - f0 - h[i]*J0.dot(dx) )
order0 = np.log10(old_div(E0[:-1],E0[1:]))
order1 = np.log10(old_div(E1[:-1],E1[1:]))
order0 = np.log10(E0[:-1]/E0[1:])
order1 = np.log10(E1[:-1]/E1[1:])
print(" %d %1.2e %1.3e %1.3e %1.3f" % (i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]))
# Ensure we are about precision
SimPEG/Utils/ModelBuilder.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
import scipy.ndimage as ndi
import scipy.sparse as sp
SimPEG/Utils/SolverUtils.py
+1 -2
View File
@@ -6,7 +6,6 @@ from future import standard_library
standard_library.install_aliases()
from builtins import range
from builtins import object
from past.utils import old_div
import numpy as np, scipy.sparse as sp
from .matutils import mkvc
import warnings
@@ -81,7 +80,7 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6, name=
if factorize and hasattr(self.solver, 'clean'):
return self.solver.clean()
return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
return type(str(name if name is not None else fun.__name__), (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
SimPEG/Utils/codeutils.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import types
import time
import numpy as np
SimPEG/Utils/coordutils.py
+3 -4
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from SimPEG.Utils import mkvc
@@ -30,7 +29,7 @@ def rotationMatrixFromNormals(v0,v1,tol=1e-20):
n0 = v0*1./np.linalg.norm(v0)
n1 = v1*1./np.linalg.norm(v1)
n0dotn1 = n0.dot(n1)
n0dotn1 = n0.dot(n1)
# define the rotation axis, which is the cross product of the two vectors
rotAx = np.cross(n0,n1)
@@ -54,7 +53,7 @@ def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]):
:param numpy.array n0: vector of length 3, should have norm 1
:param numpy.array n1: vector of length 3, should have norm 1
:param numpy.array x0: vector of length 3, point about which we perform the rotation
:param numpy.array x0: vector of length 3, point about which we perform the rotation
:rtype: numpy.array, 3x3
:return: rotation matrix which rotates the frame so that n0 is aligned with n1
"""
@@ -66,4 +65,4 @@ def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]):
X0 = np.ones([XYZ.shape[0],1])*mkvc(x0)
return (XYZ - X0).dot(R.T) + X0 # equivalent to (R*(XYZ - X0)).T + X0
return (XYZ - X0).dot(R.T) + X0 # equivalent to (R*(XYZ - X0)).T + X0
SimPEG/Utils/curvutils.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from scipy import sparse as sp
from .matutils import mkvc, ndgrid, sub2ind, sdiag
SimPEG/Utils/interputils.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
import scipy.sparse as sp
from .matutils import mkvc, sub2ind, spzeros
SimPEG/Utils/matutils.py
-1
View File
@@ -6,7 +6,6 @@ from future import standard_library
standard_library.install_aliases()
from builtins import range
from builtins import object
from past.utils import old_div
import numpy as np
import scipy.sparse as sp
from .codeutils import isScalar
SimPEG/Utils/meshutils.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import print_function
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from scipy import sparse as sp
from .matutils import mkvc, ndgrid, sub2ind, sdiag
tests/base/test_PropMaps.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from scipy.constants import mu_0
tests/base/test_Solver.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG.Mesh import TensorMesh
tests/base/test_regularization.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import unittest
from SimPEG import *
tests/dcip/test_forward_DCproblem.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
import SimPEG.DCIP as DC
tests/dcip/test_forward_IPproblem.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
import SimPEG.DCIP as DC
from SimPEG import *
tests/dcip/test_sens_IPproblem.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
import SimPEG.DCIP as DC
tests/em/fdem/forward/test_FDEM_analytics.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
tests/em/static/test_DC_2D_analytic.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import Mesh, Utils, EM, Maps, np
import SimPEG.EM.Static.DC as DC
tests/em/static/test_DC_analytic.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import Mesh, Utils, EM, Maps, np
import SimPEG.EM.Static.DC as DC
tests/em/static/test_DC_jvecjtvecadj.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
import SimPEG.EM.Static.DC as DC
tests/em/static/test_IP_fwd.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import Mesh, Utils, EM, Maps, np
import SimPEG.EM.Static.DC as DC
tests/em/static/test_IP_jvecjtvecadj.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
import SimPEG.EM.Static.DC as DC
tests/em/static/test_SIP_jvecjtvecadj.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
import SimPEG
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
tests/em/tdem/test_TDEM_combos.py
-1
View File
@@ -5,7 +5,6 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
tests/em/tdem/test_TDEM_forward_Analytic.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
tests/mesh/test_Mixed_boundaryPoisson.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import scipy.sparse as sp
import unittest
tests/mesh/test_TreeMesh.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Mesh, Tests
from SimPEG.Mesh.TreeMesh import CellLookUpException
import numpy as np
tests/mesh/test_TreeOperators.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import unittest
from SimPEG import Utils, Tests
tests/mesh/test_cylMesh.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
import sys
from SimPEG import *
tests/mesh/test_innerProduct.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import unittest
from SimPEG import Utils, Tests
tests/mt/test_ApparentResistivityAnalytic.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import MT
tests/mt/test_Problem1D_againstAnalyticHalfspace.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
import SimPEG as simpeg
from SimPEG import MT
tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
import SimPEG as simpeg
from SimPEG import MT
tests/mt/test_Problem3D_againstAnalytic.py
-1
View File
@@ -6,7 +6,6 @@ from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import str
from past.utils import old_div
# Test functions
from glob import glob
import numpy as np, sys, os, time, scipy, subprocess
tests/utils/test_Zero.py
-1
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG.Utils import Zero, Identity, sdiag, mkvc
from SimPEG import np, sp
tests/utils/test_coordutils.py
+6 -7
View File
@@ -4,7 +4,6 @@ from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest, os
import numpy as np
from SimPEG import Utils
@@ -34,11 +33,11 @@ class coorUtilsTest(unittest.TestCase):
v0*= old_div(1.,np.linalg.norm(v0))
np.random.seed(15)
v1 = np.random.rand(3)
v1*= old_div(1.,np.linalg.norm(v1))
v1 = np.random.rand(3)
v1*= old_div(1.,np.linalg.norm(v1))
v2 = Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1))
v2 = Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1))
self.assertTrue(np.linalg.norm(v2-v1) < tol)
def test_rotateMatrixFromNormals(self):
@@ -47,8 +46,8 @@ class coorUtilsTest(unittest.TestCase):
n0 *= old_div(1.,np.linalg.norm(n0))
np.random.seed(25)
n1 = np.random.rand(3)
n1 *= old_div(1.,np.linalg.norm(n1))
n1 = np.random.rand(3)
n1 *= old_div(1.,np.linalg.norm(n1))
np.random.seed(30)
scale = np.random.rand(100,1)