Skip to content

Commit

Permalink
added tests on mapping inverses
Browse files Browse the repository at this point in the history
  • Loading branch information
@lheagy
lheagy committed Jul 24, 2016
1 parent a3c0c10 commit 718742b
Show file tree
Hide file tree
Showing 2 changed files with 137 additions and 69 deletions.
151 changes: 94 additions & 57 deletions SimPEG/Maps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ def test(self, m=None, **kwargs):
m = abs(np.random.rand(self.nP))
if 'plotIt' not in kwargs:
kwargs['plotIt'] = False
return checkDerivative(lambda m : [self * m, self.deriv(m)], m, num=4,
return checkDerivative(lambda m: [self * m, self.deriv(m)], m, num=4,
**kwargs)

def testVec(self, m=None, **kwargs):
Expand Down Expand Up @@ -217,7 +217,9 @@ def deriv(self, m, v=None):
return deriv

def __str__(self):
return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in
self.maps]), self.shape[0],
self.shape[1])


class ExpMap(IdentityMap):
Expand Down Expand Up @@ -247,7 +249,8 @@ def inverse(self, D):
:rtype: numpy.array
:return: model
The *transformInverse* changes the physical property into the model.
The *transformInverse* changes the physical property into the
model.
.. math::
Expand Down Expand Up @@ -287,7 +290,8 @@ def deriv(self, m, v=None):

class ReciprocalMap(IdentityMap):
"""
Reciprocal mapping. For example, electrical resistivity and conductivity.
Reciprocal mapping. For example, electrical resistivity and
conductivity.
.. math::
Expand All @@ -298,11 +302,11 @@ def _transform(self, m):
return 1.0 / Utils.mkvc(m)

def inverse(self, D):
return 1.0 / Utils.mkvc(m)
return 1.0 / Utils.mkvc(D)

def deriv(self, m, v=None):
# TODO: if this is a tensor, you might have a problem.
deriv = Utils.sdiag( - Utils.mkvc(m)**(-2) )
deriv = Utils.sdiag(- Utils.mkvc(m)**(-2))
if v is not None:
return deriv * v
return deriv
Expand All @@ -324,7 +328,8 @@ class LogMap(IdentityMap):
m = \\exp(p)
NOTE: If you have a model which is log conductivity (ie. \\(m = \\log(\\sigma)\\)),
NOTE: If you have a model which is log conductivity
(ie. \\(m = \\log(\\sigma)\\)),
you should be using an ExpMap
"""
Expand Down Expand Up @@ -378,17 +383,19 @@ def deriv(self, m, v=None):
:rtype: numpy.array
:return: derivative of transformed model
"""
deriv = np.ones([self.mesh.nC,1])
deriv = np.ones([self.mesh.nC, 1])
if v is not None:
return deriv * v
return deriv


class FullMap(SurjectFull):
"""FullMap is depreciated. Use SurjectVertical1DMap instead.
"""
def __init__(self, mesh, **kwargs):
warnings.warn(
"`FullMap` is deprecated and will be removed in future versions. Use `SurjectFull` instead",
"`FullMap` is deprecated and will be removed in future versions."
" Use `SurjectFull` instead",
FutureWarning)
SurjectFull.__init__(self, mesh, **kwargs)

Expand Down Expand Up @@ -429,9 +436,9 @@ def deriv(self, m, v=None):
"""
repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
repVec = sp.csr_matrix(
(np.ones(repNum),
(range(repNum), np.zeros(repNum))
), shape=(repNum, 1))
(np.ones(repNum),
(range(repNum), np.zeros(repNum))
), shape=(repNum, 1))
deriv = sp.kron(sp.identity(self.nP), repVec)
if v is not None:
return deriv * v
Expand All @@ -444,7 +451,8 @@ class Vertical1DMap(SurjectVertical1D):
"""
def __init__(self, mesh, **kwargs):
warnings.warn(
"`Vertical1DMap` is deprecated and will be removed in future versions. Use `SurjectVertical1D` instead",
"`Vertical1DMap` is deprecated and will be removed in future"
"versions. Use `SurjectVertical1D` instead",
FutureWarning)
SurjectVertical1D.__init__(self, mesh, **kwargs)

Expand All @@ -459,9 +467,10 @@ class Surject2Dto3D(IdentityMap):
normal = 'Y' #: The normal

def __init__(self, mesh, **kwargs):
assert mesh.dim == 3, 'Only works for a 3D Mesh'
assert mesh.dim == 3, 'Surject2Dto3D Only works for a 3D Mesh'
IdentityMap.__init__(self, mesh, **kwargs)
assert self.normal in ['X', 'Y', 'Z'], 'For now, only "Y" normal is supported'
assert self.normal in ['X', 'Y', 'Z'], ('For now, only "Y"'
' normal is supported')

@property
def nP(self):
Expand All @@ -484,11 +493,17 @@ def _transform(self, m):
"""
m = Utils.mkvc(m)
if self.normal == 'Z':
return Utils.mkvc(m.reshape(self.mesh.vnC[[0,1]], order='F')[:,:,np.newaxis].repeat(self.mesh.nCz,axis=2))
return Utils.mkvc(m.reshape(self.mesh.vnC[[0,1]], order='F'
)[:, :, np.newaxis].repeat(self.mesh.nCz,
axis=2))
elif self.normal == 'Y':
return Utils.mkvc(m.reshape(self.mesh.vnC[[0,2]], order='F')[:,np.newaxis,:].repeat(self.mesh.nCy,axis=1))
return Utils.mkvc(m.reshape(self.mesh.vnC[[0,2]], order='F'
)[:, np.newaxis, :].repeat(self.mesh.nCy,
axis=1))
elif self.normal == 'X':
return Utils.mkvc(m.reshape(self.mesh.vnC[[1,2]], order='F')[np.newaxis,:,:].repeat(self.mesh.nCx,axis=0))
return Utils.mkvc(m.reshape(self.mesh.vnC[[1,2]], order='F'
)[np.newaxis, :, :].repeat(self.mesh.nCx,
axis=0))

def deriv(self, m, v=None):
"""
Expand All @@ -498,10 +513,9 @@ def deriv(self, m, v=None):
"""
inds = self * np.arange(self.nP)
nC, nP = self.mesh.nC, self.nP
P = sp.csr_matrix(
(np.ones(nC),
(range(nC), inds)
), shape=(nC, nP))
P = sp.csr_matrix((np.ones(nC),
(range(nC), inds)
), shape=(nC, nP))
if v is not None:
return P * v
return P
Expand All @@ -513,7 +527,8 @@ class Map2Dto3D(Surject2Dto3D):

def __init__(self, mesh, **kwargs):
warnings.warn(
"`Map2Dto3D` is deprecated and will be removed in future versions. Use `Surject2Dto3D` instead",
"`Map2Dto3D` is deprecated and will be removed in future versions."
" Use `Surject2Dto3D` instead",
FutureWarning)
Surject2Dto3D.__init__(self, mesh, **kwargs)

Expand All @@ -528,12 +543,14 @@ def __init__(self, meshes, **kwargs):

assert type(meshes) is list, "meshes must be a list of two meshes"
assert len(meshes) == 2, "meshes must be a list of two meshes"
assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
assert meshes[0].dim == meshes[1].dim, ("The two meshes must be the "
"same dimension")

self.mesh = meshes[0]
self.mesh2 = meshes[1]

self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC, 'CC', zerosOutside=True)
self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC, 'CC',
zerosOutside=True)

@property
def shape(self):
Expand All @@ -560,12 +577,12 @@ class InjectActiveCells(IdentityMap):
"""

indActive = None #: Active Cells
indActive = None #: Active Cells
valInactive = None #: Values of inactive Cells
nC = None #: Number of cells in the full model
nC = None #: Number of cells in the full model

def __init__(self, mesh, indActive, valInactive, nC=None):
self.mesh = mesh
self.mesh = mesh

self.nC = nC or mesh.nC

Expand All @@ -585,7 +602,7 @@ def __init__(self, mesh, indActive, valInactive, nC=None):

inds = np.nonzero(self.indActive)[0]
self.P = sp.csr_matrix((np.ones(inds.size), (inds, range(inds.size))),
shape=(self.nC, self.nP)
shape=(self.nC, self.nP)
)

@property
Expand Down Expand Up @@ -630,7 +647,7 @@ class Weighting(IdentityMap):
nC = None #: Number of cells in the full model

def __init__(self, mesh, weights=None, nC=None):
self.mesh = mesh
self.mesh = mesh

self.nC = nC or mesh.nC

Expand Down Expand Up @@ -672,7 +689,7 @@ class ComplexMap(IdentityMap):
def __init__(self, mesh, nP=None):
IdentityMap.__init__(self, mesh)
if nP is not None:
assert nP%2 == 0, 'nP must be even.'
assert nP % 2 == 0, 'nP must be even.'
self._nP = nP or (self.mesh.nC * 2)

@property
Expand All @@ -681,7 +698,7 @@ def nP(self):

@property
def shape(self):
return (self.nP/2,self.nP)
return (self.nP/2, self.nP)

def _transform(self, m):
nC = self.mesh.nC
Expand All @@ -700,7 +717,7 @@ def adj(v):
return LinearOperator(shp, matvec=fwd, rmatvec=adj) * v
return LinearOperator(shp, matvec=fwd, rmatvec=adj)

inverse = deriv
# inverse = deriv


class ParametricCircleMap(IdentityMap):
Expand Down Expand Up @@ -762,9 +779,14 @@ def deriv(self, m, v=None):
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))
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))

if v is not None:
return sp.csr_matrix(np.c_[g1, g2, g3, g4, g5]) * v
Expand Down Expand Up @@ -843,9 +865,9 @@ def _transform(self, m):
if self.mesh.dim == 2:
X = self.mesh.gridCC[self.actInd, 0]
Y = self.mesh.gridCC[self.actInd, 1]
if self.normal =='X':
if self.normal == 'X':
f = polynomial.polyval(Y, c) - X
elif self.normal =='Y':
elif self.normal == 'Y':
f = polynomial.polyval(X, c) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
Expand Down Expand Up @@ -949,12 +971,13 @@ class ParametricSplineMap(IdentityMap):

slope = 1e4

def __init__(self, mesh, pts, ptsv=None, order=3, logSigma=True, normal='X'):
def __init__(self, mesh, pts, ptsv=None, order=3, logSigma=True,
normal='X'):
IdentityMap.__init__(self, mesh)
self.logSigma = logSigma
self.order = order
self.normal = normal
self.pts= pts
self.pts = pts
self.npts = np.size(pts)
self.ptsv = ptsv
self.spl = None
Expand Down Expand Up @@ -1002,13 +1025,16 @@ def _transform(self, m):
if np.mod(c.size, 2):
raise(Exception("Put even points!"))

self.spl = {"splb": UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
"splt": UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
self.spl = {"splb": UnivariateSpline(self.pts, c[:npts],
k=self.order, s=0),
"splt": UnivariateSpline(self.pts, c[npts:],
k=self.order, s=0)}

if self.normal == 'X':
zb = self.ptsv[0]
zt = self.ptsv[1]
flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
flines = ((self.spl["splt"](Y) - self.spl["splb"](Y)) *
(Z - zb) / (zt - zb) + self.spl["splb"](Y))
f = flines - X
# elif self.normal =='Y':
# elif self.normal =='Z':
Expand All @@ -1024,7 +1050,7 @@ def deriv(self, m, v=None):
sig1, sig2, c = m[0], m[1], m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
# 2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
Expand All @@ -1035,15 +1061,16 @@ def deriv(self, m, v=None):
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
#3D
# 3D
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
Z = self.mesh.gridCC[:, 2]
if self.normal =='X':
zb = self.ptsv[0]
zt = self.ptsv[1]
flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
flines = ((self.spl["splt"](Y)-self.spl["splb"](Y)) *
(Z - zb) / (zt - zb) + self.spl["splb"](Y))
f = flines - X
# elif self.normal =='Y':
# elif self.normal =='Z':
Expand Down Expand Up @@ -1074,7 +1101,8 @@ def deriv(self, m, v=None):
spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
fderiv = (spla(X)-splb(X))/(2*dy)
g3[:, i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
g3[:, i] = Utils.sdiag(alpha*(sig2-sig1) /
(1.+(alpha*f)**2) / np.pi)*fderiv

elif self.mesh.dim == 3:
g3 = np.zeros((self.mesh.nC, self.npts*2))
Expand All @@ -1091,20 +1119,29 @@ def deriv(self, m, v=None):

# treat bottom boundary
if i < self.npts:
splba = UnivariateSpline(self.pts, ca[:self.npts], k=self.order, s=0)
splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0)
flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X
flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X
splba = UnivariateSpline(self.pts, ca[:self.npts],
k=self.order, s=0)
splbb = UnivariateSpline(self.pts, cb[:self.npts],
k=self.order, s=0)
flinesa = ((self.spl["splt"](Y) - splba(Y)) * (Z-zb) /
(zt-zb) + splba(Y) - X)
flinesb = ((self.spl["splt"](Y) - splbb(Y)) * (Z-zb) /
(zt-zb) + splbb(Y) - X)

# treat top boundary
else:
splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
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
splta = UnivariateSpline(self.pts, ca[self.npts:],
k=self.order, s=0)
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 = (flinesa-flinesb)/(2*dy)
g3[:, i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
else :
g3[:, i] = Utils.sdiag(alpha*(sig2-sig1) /
(1.+(alpha*f)**2) / np.pi)*fderiv
else:
raise(Exception("Not Implemented for Y and Z, your turn :)"))

if v is not None:
Expand Down

0 comments on commit 718742b

Please sign in to comment.