mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 08:46:27 +08:00
Maps documentation.
This commit is contained in:
rowanc1
+46
-19
@@ -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
|
||||
|
||||
+13
-11
@@ -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)
|
||||
|
||||
+154
-14
@@ -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:
|
||||
|
||||
Reference in New Issue
Block a user