Skip to content

Commit

Permalink
Remove the Meshless Identity Map.
Browse files Browse the repository at this point in the history 
 - This is now default functionality in the IdentityMap.
  • Loading branch information
Rowan Cockett committed Dec 4, 2015
1 parent cfc921b commit c298ebe
Show file tree
Hide file tree
Showing 2 changed files with 57 additions and 78 deletions.
127 changes: 53 additions & 74 deletions SimPEG/Maps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,33 +10,41 @@ class IdentityMap(object):
SimPEG Map
"""

__metaclass__ = Utils.SimPEGMetaClass

mesh = None #: A SimPEG Mesh

def __init__(self, mesh, **kwargs):
def __init__(self, mesh=None, nP=None, **kwargs):
Utils.setKwargs(self, **kwargs)

if nP is not None:
assert type(nP) in [int, long], ' Number of parameters must be an integer.'

self.mesh = mesh
self._nP = nP

@property
def nP(self):
"""
:rtype: int
:return: number of parameters in the model
:return: number of parameters that the mapping accepts
"""
if self._nP is not None:
return self._nP
if self.mesh is None:
return '*'
return self.mesh.nC

@property
def shape(self):
"""
The default shape is (mesh.nC, nP).
The default shape is (mesh.nC, nP) if the mesh is defined.
If this is a meshless mapping (i.e. nP is defined independently)
the shape will be the the shape (nP,nP).
:rtype: (int,int)
:return: shape of the operator as a tuple
"""
if self._nP is not None:
return (self.nP, self.nP)
if self.mesh is None:
return ('*', self.nP)
return (self.mesh.nC, self.nP)
Expand Down Expand Up @@ -119,35 +127,6 @@ def __str__(self):
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])


class IdentityMap_Meshless(IdentityMap):

def __init__(self, nP=None, **kwargs):
IdentityMap.__init__(self, None, **kwargs)
self._nP = nP

@property
def nP(self):
"""
:rtype: int
:return: number of parameters in the model
"""
if self._nP is None:
return '*'
return self._nP

@property
def shape(self):
"""
The default shape is (mesh.nC, nP).
:rtype: (int,int)
:return: shape of the operator as a tuple
"""
if self._nP is None:
return ('*', '*')
return (self.nP, self.nP)


class ComboMap(IdentityMap):
"""Combination of various maps."""

Expand Down Expand Up @@ -505,7 +484,7 @@ def __init__(self, mesh, indActive, valInactive, nC=None):
else:
self.valInactive = valInactive.copy()
self.valInactive[self.indActive] = 0

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

Expand Down Expand Up @@ -738,7 +717,7 @@ class PolyMap(IdentityMap):
Parameterize the model space using a polynomials in a wholespace.
..math::
y = \mathbf{V} c
Define the model as:
Expand Down Expand Up @@ -782,10 +761,10 @@ def _transform(self, m):
else:
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
elif self.normal =='Y':
Expand All @@ -796,43 +775,43 @@ def _transform(self, m):
raise(Exception("Input for normal = X or Y or Z"))
else:
raise(Exception("Only supports 2D"))


return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)

def deriv(self, m):
alpha = self.slope
sig1,sig2, c = m[0],m[1],m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
if self.mesh.dim == 2:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]

if self.normal =='X':
f = polynomial.polyval(Y, c) - X
V = polynomial.polyvander(Y, len(c)-1)
V = polynomial.polyvander(Y, len(c)-1)
elif self.normal =='Y':
f = polynomial.polyval(X, c) - Y
V = polynomial.polyvander(X, len(c)-1)
V = polynomial.polyvander(X, len(c)-1)
else:
raise(Exception("Input for normal = X or Y or Z"))
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]

if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
V = polynomial.polyvander2d(Y, Z, self.order)
V = polynomial.polyvander2d(Y, Z, self.order)
elif self.normal =='Y':
f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
V = polynomial.polyvander2d(X, Z, self.order)
V = polynomial.polyvander2d(X, Z, self.order)
elif self.normal =='Z':
f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
V = polynomial.polyvander2d(X, Y, self.order)
V = polynomial.polyvander2d(X, Y, self.order)
else:
raise(Exception("Input for normal = X or Y or Z"))

Expand All @@ -845,16 +824,16 @@ def deriv(self, m):

g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V

return sp.csr_matrix(np.c_[g1,g2,g3])
return sp.csr_matrix(np.c_[g1,g2,g3])

class SplineMap(IdentityMap):

"""SplineMap
Parameterize the boundary of two geological units using a spline interpolation
Parameterize the boundary of two geological units using a spline interpolation
..math::
g = f(x)-y
Define the model as:
Expand All @@ -879,7 +858,7 @@ def __init__(self, mesh, pts, ptsv=None,order=3, logSigma=True, normal='X'):
def nP(self):
if self.mesh.dim == 2:
return np.size(self.pts)+2
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
return np.size(self.pts)*2+2
else:
raise(Exception("Only supports 2D and 3D"))
Expand All @@ -896,28 +875,28 @@ def _transform(self, m):
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0)
if self.normal =='X':
if self.normal =='X':
f = self.spl(Y) - X
elif self.normal =='Y':
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))

# 3D:
# Comments:
# 3D:
# Comments:
# Make two spline functions and link them using linear interpolation.
# This is not quite direct extension of 2D to 3D case
# Using 2D interpolation is possible

elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]

npts = np.size(self.pts)
npts = np.size(self.pts)
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)}

Expand All @@ -932,7 +911,7 @@ def _transform(self, m):
raise(Exception("Input for normal = X or Y or Z"))
else:
raise(Exception("Only supports 2D and 3D"))


return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)

Expand All @@ -942,7 +921,7 @@ def deriv(self, m):
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
if self.mesh.dim == 2:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]

Expand All @@ -951,17 +930,17 @@ def deriv(self, m):
elif self.normal =='Y':
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
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)
f = flines - X
f = flines - X
# elif self.normal =='Y':
# elif self.normal =='Z':
else:
Expand All @@ -974,7 +953,7 @@ def deriv(self, m):
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:
g3 = np.zeros((self.mesh.nC, self.npts))
if self.normal =='Y':
Expand All @@ -988,7 +967,7 @@ def deriv(self, m):
cb = c.copy()
dy = self.mesh.hy[ind]*1.5
ca[i] = ctemp+dy
cb[i] = ctemp-dy
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 = (spla(X)-splb(X))/(2*dy)
Expand All @@ -998,7 +977,7 @@ def deriv(self, m):
g3 = np.zeros((self.mesh.nC, self.npts*2))
if self.normal =='X':
# Here we use perturbation to compute sensitivity
for i in range(self.npts*2):
for i in range(self.npts*2):
ctemp = c[i]
ind = np.argmin(abs(self.mesh.vectorCCy-ctemp))
ca = c.copy()
Expand All @@ -1012,20 +991,20 @@ def deriv(self, m):
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
#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
fderiv = (flinesa-flinesb)/(2*dy)
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 :
raise(Exception("Not Implemented for Y and Z, your turn :)"))
return sp.csr_matrix(np.c_[g1,g2,g3])
return sp.csr_matrix(np.c_[g1,g2,g3])







8 changes: 4 additions & 4 deletions tests/base/test_regularization.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ def test_regularization(self):
if not inspect.isclass(r): continue
if not issubclass(r, Regularization.BaseRegularization):
continue

for i, mesh in enumerate(self.meshlist):

print 'Testing %iD'%mesh.dim
Expand All @@ -32,7 +32,7 @@ def test_regularization(self):
m = np.random.rand(mapping.nP)
reg.mref = np.ones_like(m)*np.mean(m)

print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
passed = reg.eval(reg.mref) < TOL
self.assertTrue(passed)

Expand All @@ -50,7 +50,7 @@ def test_regularization_ActiveCells(self):
if not inspect.isclass(r): continue
if not issubclass(r, Regularization.BaseRegularization):
continue

for i, mesh in enumerate(self.meshlist):

print 'Testing Active Cells %iD'%(mesh.dim)
Expand All @@ -62,7 +62,7 @@ def test_regularization_ActiveCells(self):
elif mesh.dim == 3:
indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)

mapping = Maps.IdentityMap_Meshless(nP=indAct.nonzero()[0].size)
mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)

reg = r(mesh, mapping=mapping, indActive=indAct)
m = np.random.rand(mesh.nC)[indAct]
Expand Down

0 comments on commit c298ebe

Please sign in to comment.