From c8868177731b70bda818bc3d10d7858bdb87530d Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 14 Sep 2015 13:06:51 -0700 Subject: [PATCH 1/2] Fix bug in CircMap Add PolyMap for 2D, and test --- SimPEG/Maps.py | 61 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 6dca13dd..0380447e 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -696,4 +696,63 @@ class CircleMap(IdentityMap): 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)) - return np.c_[g1,g2,g3,g4,g5] + return sp.csr_matrix(np.c_[g1,g2,g3,g4,g5]) + + +class PolyMap(IdentityMap): + + """PolyMap + + Parameterize the model space using a polynomials in a wholespace. + + ..math:: + + y = \mathbf{V} c + + Define the model as: + + ..math:: + + m = [\sigma_1, \sigma_2, c] + + """ + def __init__(self, mesh, order, logSigma=True): + assert mesh.dim == 2, "Working for a 2D mesh only right now. But it isn't that hard to change.. :)" + IdentityMap.__init__(self, mesh) + self.logSigma = logSigma + self.order = order + + slope = 1e4 + + @property + def nP(self): + return self.order+3 + + def _transform(self, m): + alpha = self.slope + sig1,sig2 = m[0],m[1] + c = m[2:] + if self.logSigma: + sig1, sig2 = np.exp(sig1), np.exp(sig2) + X = self.mesh.gridCC[:,0] + Y = self.mesh.gridCC[:,1] + f = polynomial.polyval(X, c) - Y + 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) + X = self.mesh.gridCC[:,0] + Y = self.mesh.gridCC[:,1] + f = polynomial.polyval(X, c) - Y + V = polynomial.polyvander(X, len(c)-1) + if self.logSigma: + g1 = -(np.arctan(alpha*f)/np.pi + 0.5)*sig1 + sig1 + g2 = (np.arctan(alpha*f)/np.pi + 0.5)*sig2 + else: + 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 + return sp.csr_matrix(np.c_[g1,g2,g3]) From f726a27269054d3b7d02fb7ce337abc8eaba7b5e Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 14 Sep 2015 21:25:15 -0700 Subject: [PATCH 2/2] Fix target misfit --- SimPEG/Directives.py | 2 +- SimPEG/Maps.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 3ec9e96e..48d7abcf 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -149,7 +149,7 @@ class TargetMisfit(InversionDirective): @property def target(self): if getattr(self, '_target', None) is None: - self._target = self.survey.nD + self._target = self.survey.nD*0.5 return self._target @target.setter def target(self, val): diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 0380447e..ae6411fd 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -1,7 +1,7 @@ import Utils, numpy as np, scipy.sparse as sp from Tests import checkDerivative from PropMaps import PropMap, Property - +from numpy.polynomial import polynomial class IdentityMap(object): """