From 23e6ece0cd9ae76d1ced84e68b0c4ac79187582e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 17 May 2014 21:05:31 -0700 Subject: [PATCH] Maps documentation. --- SimPEG/Maps.py | 65 ++++++++++----- SimPEG/Tests/TestUtils.py | 24 +++--- docs/api_Maps.rst | 168 ++++++++++++++++++++++++++++++++++---- 3 files changed, 213 insertions(+), 44 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 5061191e..ec481145 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -43,59 +43,82 @@ class IdentityMap(object): __metaclass__ = Utils.SimPEGMetaClass - counter = None #: A SimPEG.Utils.Counter object mesh = None #: A SimPEG Mesh def __init__(self, mesh): self.mesh = mesh + @property + def nP(self): + """ + :rtype: int + :return: number of parameters in the model + """ + return self.mesh.nC + + @property + def shape(self): + """ + The default shape is (mesh.nC, nP). + + :rtype: (int,int) + :return: shape of the operator as a tuple + """ + return (self.mesh.nC, self.nP) + def transform(self, m): """ + Changes the model into the physical property. + + .. note:: + + This can be called by the __mul__ property against a numpy.ndarray. + :param numpy.array m: model :rtype: numpy.array :return: transformed model - The *transform* changes the model into the physical property. - """ return m def transformInverse(self, D): """ + Changes the physical property into the model. + + .. note:: + + The *transformInverse* may not be easy to create in general. + :param numpy.array D: physical property :rtype: numpy.array :return: model - The *transformInverse* changes the physical property into the model. - - .. note:: The *transformInverse* may not be easy to create in general. - """ raise NotImplementedError('The transformInverse is not implemented.') def transformDeriv(self, m): """ + The derivative of the transformation. + :param numpy.array m: model :rtype: scipy.csr_matrix :return: derivative of transformed model - The *transform* changes the model into the physical property. - The *transformDeriv* provides the derivative of the *transform*. """ return sp.identity(m.size) - @property - def nP(self): - """Number of parameters in the model.""" - return self.mesh.nC - - def example(self): - return np.random.rand(self.nP) - def test(self, m=None, **kwargs): + """Test the derivative of the mapping. + + :param numpy.array m: model + :param kwargs: key word arguments of :meth:`SimPEG.Tests.checkDerivative` + :rtype: bool + :return: passed the test? + + """ print 'Testing the %s Class!' % self.__class__.__name__ if m is None: - m = self.example() + m = np.random.rand(self.nP) if 'plotIt' not in kwargs: kwargs['plotIt'] = False return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, **kwargs) @@ -419,7 +442,7 @@ class ComboMap(IdentityMap): return self.transform(val) class ComplexMap(IdentityMap): - """docstring for ComplexMap + """ComplexMap default nP is nC in the mesh times 2 [real, imag] @@ -434,6 +457,10 @@ class ComplexMap(IdentityMap): def nP(self): return self._nP + @property + def shape(self): + return (self.nP/2,self.nP) + def transform(self, m): nC = self.mesh.nC return m[:nC] + m[nC:]*1j diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index 0fb1d200..3cd40991 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -202,7 +202,7 @@ def Rosenbrock(x, return_g=True, return_H=True): out += (H,) return out if len(out) > 1 else out[0] -def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10): +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10, ax=None): """ Basic derivative check @@ -231,7 +231,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole """ print "%s checkDerivative %s" % ('='*20, '='*20) - print "iter h |f0-ft| |f0-ft-h*J0*dx| Order\n%s" % ('-'*57) + print "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n%s" % ('-'*57) f0, J0 = fctn(x0) @@ -282,14 +282,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole if plotIt: - plt.figure() - plt.clf() - plt.loglog(h, E0, 'b') - plt.loglog(h, E1, 'g--') - plt.title('checkDerivative') - plt.xlabel('h') - plt.ylabel('error of Taylor approximation') - plt.legend(['0th order', '1st order'], loc='upper left') + ax = ax or plt.subplot(111) + ax.loglog(h, E0, 'b') + ax.loglog(h, E1, 'g--') + ax.set_title('Check Derivative - %s' % ('PASSED :)' if passTest else 'FAILED :(')) + ax.set_xlabel('h') + ax.set_ylabel('Error') + leg = ax.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$'], loc='best', + title="$f(x + h\Delta x) - f(x) - h g(x) \Delta x - \mathcal{O}(h^2) = 0$", + frameon=False) + plt.setp(leg.get_title(),fontsize=15) plt.show() return passTest @@ -328,6 +330,6 @@ if __name__ == '__main__': def simpleFail(x): return np.sin(x), -sdiag(np.cos(x)) - checkDerivative(simplePass, np.random.randn(5), plotIt=False) + checkDerivative(simplePass, np.random.randn(5), plotIt=True) checkDerivative(simpleFunction, np.random.randn(5), plotIt=False) checkDerivative(simpleFail, np.random.randn(5), plotIt=False) diff --git a/docs/api_Maps.rst b/docs/api_Maps.rst index 7974c272..b590203c 100644 --- a/docs/api_Maps.rst +++ b/docs/api_Maps.rst @@ -1,18 +1,74 @@ .. _api_Maps: -Maps -**** +SimPEG Maps +*********** A SimPEG Map operates on a vector and transforms it to another space. We will use an example commonly applied in electromagnetics (EM) of the log-conductivity model. -Electrical conductivity varies over many orders of magnitude, so it is a common -technique when solving the inverse problem to parameterize and optimize in terms -of log conductivity. This makes sense not only because it ensures all conductivities -will be positive, but because this is fundamentally the space where conductivity -lives (i.e. it varies logarithmically). In SimPEG, we use a (:class:`SimPEG.Maps.ExpMap`) to -describe how to map back to conductivity. + +.. math:: + + m = \log(\sigma) + +Here we require a *mapping* to get from \\\(m\\\) to \\\(\\sigma\\\), +we will call this map \\\(\\mathcal{M}\\\). + +.. math:: + + \sigma = \mathcal{M}(m) = \exp(m) + +In SimPEG, we use a (:class:`SimPEG.Maps.ExpMap`) to describe how to map +back to conductivity. This is a relatively trivial example (we are just taking +the exponential!) but by defining maps we can start to combine and manipulate +exactly what we think about as our model, \\\(m\\\). In code, this looks like + +.. plot:: + :include-source: + + from SimPEG import * + import matplotlib.pyplot as plt + M = Mesh.TensorMesh([100]) + expMap = Maps.ExpMap(M) + m = np.zeros(M.nC) + m[M.vectorCCx>0.5] = 1.0 + sig = expMap.transform(m) # or just: expMap * m + plt.plot(M.vectorCCx, np.c_[m,sig]) + plt.legend(['$m=\log(\sigma)$','$\sigma=\exp(m)$'],'best') + plt.ylim([-0.5,3]) + +Taking Derivatives +================== + +Now that we have wrapped up the mapping, we can ensure that it is easy to take +derivatives (or at least have access to them!). In the :class:`SimPEG.Maps.ExpMap` +there are no dependencies between model parameters, so it will be a diagonal matrix: + +.. math:: + + \left(\frac{\partial \mathcal{M}(m)}{\partial m}\right)_{ii} = \frac{\partial \exp(m_i)}{\partial m} = \exp(m_i) + +Or equivalently: + +.. math:: + + \frac{\partial \mathcal{M}(m)}{\partial m} = \text{diag}( \exp(m) ) + +The mapping API makes this really easy to test that you have got the derivative correct. +When these are used in the inverse problem, this is extremely important!! + +.. plot:: + :include-source: + + from SimPEG import * + import matplotlib.pyplot as plt + M = Mesh.TensorMesh([100]) + expMap = Maps.ExpMap(M) + m = np.zeros(M.nC) + m[M.vectorCCx>0.5] = 1.0 + expMap.test(m, plotIt=True) + The API ======= @@ -21,25 +77,109 @@ The API :members: :undoc-members: -.. autoclass:: SimPEG.Maps.NonLinearMap - :members: - :undoc-members: -.. autoclass:: SimPEG.Maps.ComboMap - :members: - :undoc-members: +Combining Maps +============== + +We will use an example where we want a 1D layered earth as +our model, but we want to map this to a 2D discretization to do our forward +modeling. We will also assume that we are working in log conductivity still, +so after the transformation we want to map to conductivity space. +To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.Vertical1DMap`), +which does the first part of what we just described. The second part will be +done by the :class:`SimPEG.Maps.ExpMap` described above. + +:: + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.Vertical1DMap(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + +.. plot:: + + from SimPEG import * + import matplotlib.pyplot as plt + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.Vertical1DMap(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + figs, axs = plt.subplots(1,2) + axs[0].plot(m, M.vectorCCy, 'b-o') + axs[0].set_title('Model') + axs[0].set_ylabel('Depth, y') + axs[0].set_xlabel('Value, $m_i$') + axs[0].set_xlim(0,3) + axs[0].set_ylim(0,1) + clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0]) + axs[1].set_title('Physical Property') + axs[1].set_ylabel('Depth, y') + clbar.set_label('$\sigma = \exp(\mathbf{P}m)$') + plt.tight_layout() + +If you noticed, it was pretty easy to combine maps. What is even cooler is +that the derivatives also are made for you (if everything goes right). +Just to be sure that the derivative is correct, you should always run the test +on the mapping that you create. + Common Maps =========== + +Exponential Map +--------------- + +Electrical conductivity varies over many orders of magnitude, so it is a common +technique when solving the inverse problem to parameterize and optimize in terms +of log conductivity. This makes sense not only because it ensures all conductivities +will be positive, but because this is fundamentally the space where conductivity +lives (i.e. it varies logarithmically). + .. autoclass:: SimPEG.Maps.ExpMap :members: :undoc-members: + +Vertical 1D Map +--------------- + .. autoclass:: SimPEG.Maps.Vertical1DMap :members: :undoc-members: + +Mesh to Mesh Map +---------------- + .. autoclass:: SimPEG.Maps.Mesh2Mesh :members: :undoc-members: + + +Some Extras +=========== + +Combo Map +--------- + +The ComboMap holds the information for multiplying and combining +maps. It also uses the chain rule create the derivative. +Remember, any time that you make your own combination of mappings +be sure to test that the derivative is correct. + +.. autoclass:: SimPEG.Maps.ComboMap + :members: + :undoc-members: + + +Non Linear Map +-------------- + +.. autoclass:: SimPEG.Maps.NonLinearMap + :members: + :undoc-members: