wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-11 11:26:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Documentation updates for forward problem.

Browse Source
This commit is contained in:
rowanc1
2014-04-06 21:50:11 -06:00
parent 1506575b01
commit 3a9def5ad8
8 changed files with 177 additions and 112 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Problem.py
+24 -66
View File
@@ -4,42 +4,14 @@ import Model
class BaseProblem(object):
"""
Problem is the base class for all geophysical forward problems in SimPEG.
The problem is a partial differential equation of the form:
.. math::
c(m, u) = 0
Here, m is the model and u is the field (or fields).
Given the model, m, we can calculate the fields u(m),
however, the data we collect is a subset of the fields,
and can be defined by a linear projection, P.
.. math::
d_\\text{pred} = Pu(m)
We are interested in how changing the model transforms the data,
as such we can take write the Taylor expansion:
.. math::
Pu(m + hv) = Pu(m) + hP\\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \\right\| )
We can linearize and define the sensitivity matrix as:
.. math::
J = P\\frac{\partial u}{\partial m}
The sensitivity matrix, and it's transpose will be used in the inverse problem
to (locally) find how model parameters change the data, and optimize!
"""
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
surveyPair = Survey.BaseSurvey
modelPair = Model.BaseModel
surveyPair = Survey.BaseSurvey #: A SimPEG.Survey Class
modelPair = Model.BaseModel #: A SimPEG.Model Class
def __init__(self, model, **kwargs):
Utils.setKwargs(self, **kwargs)
@@ -49,7 +21,9 @@ class BaseProblem(object):
self.model = model
@property
def mesh(self): return self.model.mesh
def mesh(self):
"""SimPEG mesh that is associated with the model provided."""
return self.model.mesh
@property
def survey(self):
@@ -73,48 +47,33 @@ class BaseProblem(object):
self._survey = None
@property
def ispaired(self): return self.survey is not None
def ispaired(self):
"""True if the problem is paired to a survey."""
return self.survey is not None
@Utils.timeIt
def Jvec(self, m, v, u=None):
"""
Effect of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: Jv
Working with the general PDE, c(m, u) = 0, where m is the model and u is the field,
the sensitivity is defined as:
.. math::
J = P\\frac{\partial u}{\partial m}
We can take the derivative of the PDE:
.. math::
\\nabla_m c(m, u) \delta m + \\nabla_u c(m, u) \delta u = 0
If the forward problem is invertible, then we can rearrange for du/dm:
.. math::
J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
This can often be computed given a vector (i.e. J(v)) rather than stored, as J is a large dense matrix.
"""
raise NotImplementedError('J is not yet implemented.')
@Utils.timeIt
def Jtvec(self, m, v, u=None):
"""
Effect of transpose of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: JTv
Effect of transpose of J on a vector v.
"""
raise NotImplementedError('Jt is not yet implemented.')
@@ -122,29 +81,26 @@ class BaseProblem(object):
@Utils.timeIt
def Jvec_approx(self, m, v, u=None):
"""
Approximate effect of J(m) on a vector v
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: Jv
Approximate effect of J on a vector v
:return: approxJv
"""
return self.Jvec(m, v, u)
@Utils.timeIt
def Jtvec_approx(self, m, v, u=None):
"""
Approximate effect of transpose of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: JTv
Approximate transpose of J*v
"""
return self.Jtvec(m, v, u)
@@ -152,26 +108,28 @@ class BaseProblem(object):
"""
The field given the model.
.. math::
u(m)
:param numpy.array m: model
:rtype: numpy.array
:return: u, the fields
"""
pass
raise NotImplementedError('fields is not yet implemented.')
#TODO: Rename and refactor to createSyntheticData
def createSyntheticSurvey(self, m, std=0.05, u=None, **geometry_kwargs):
def createSyntheticSurvey(self, m, std=0.05, u=None, **survey_kwargs):
"""
Create synthetic survey given a model, and a standard deviation.
:param numpy.array m: geophysical model
:param numpy.array std: standard deviation
:param numpy.array u: fields for the given model (if pre-calculated)
:param numpy.array survey_kwargs: Keyword arguments for initiating the survey.
:rtype: SurveyObject
:return: survey
Returns the observed data with random Gaussian noise
and Wd which is the same size as data, and can be used to weight the inversion.
"""
survey = self.surveyPair(mtrue=m, **geometry_kwargs)
survey = self.surveyPair(mtrue=m, **survey_kwargs)
survey.pair(self)
survey.dtrue = survey.dpred(m, u=u)
noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
SimPEG/Survey.py
+16 -12
View File
@@ -60,7 +60,8 @@ class BaseSurvey(object):
@Utils.count
@Utils.requires('prob')
def dpred(self, m, u=None):
"""
"""dpred(m, u=None)
Create the projected data from a model.
The field, u, (if provided) will be used for the predicted data
instead of recalculating the fields (which may be expensive!).
@@ -77,9 +78,9 @@ class BaseSurvey(object):
@Utils.count
def projectFields(self, u):
"""
This function projects the fields onto the data space.
"""projectFields(u)
This function projects the fields onto the data space.
.. math::
@@ -89,22 +90,23 @@ class BaseSurvey(object):
@Utils.count
def projectFieldsDeriv(self, u):
"""
This function projects the fields onto the data space.
"""projectFieldsDeriv(u)
This function s the derivative of projects the fields onto the data space.
.. math::
\\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
"""
return sp.identity(u.size)
raise NotImplemented('projectFields is not yet implemented.')
@Utils.count
def residual(self, m, u=None):
"""
"""residual(m, u=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: float
:rtype: numpy.array
:return: data residual
The data residual:
@@ -129,6 +131,7 @@ class BaseSurvey(object):
"""
if getattr(self,'_Wd',None) is None:
print 'SimPEG is making Survey.Wd to be norm of the data plus a floor.'
eps = np.linalg.norm(Utils.mkvc(self.dobs),2)*1e-5
self._Wd = 1/(abs(self.dobs)*self.std+eps)
return self._Wd
@@ -137,11 +140,12 @@ class BaseSurvey(object):
self._Wd = value
def residualWeighted(self, m, u=None):
"""
"""residualWeighted(m, u=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: float
:return: data residual
:rtype: numpy.array
:return: weighted data residual
The weighted data residual:
@@ -149,7 +153,7 @@ class BaseSurvey(object):
\mu_\\text{data}^{\\text{weighted}} = \mathbf{W}_d(\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
Where W_d is a covariance matrix that weights the data residual.
Where \\\\(W_d\\\\) is a covariance matrix that weights the data residual.
"""
return Utils.mkvc(self.Wd*self.residual(m, u=u))
docs/api_Forward.rst
-30
View File
@@ -1,30 +0,0 @@
.. _api_Forward:
Model
=====
.. automodule:: SimPEG.Model
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
Survey
======
.. automodule:: SimPEG.Survey
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
Problem
=======
.. automodule:: SimPEG.Problem
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
docs/api_Mesh.rst
+7 -3
View File
@@ -39,6 +39,7 @@ the implementations.
axes[1].set_title('TreeMesh')
rM.plotGrid(ax=axes[2], **opts)
axes[2].set_title('LogicallyRectMesh')
plt.show()
Variable Locations and Terminology
@@ -58,10 +59,12 @@ of the TensorMesh.
:include-source:
from SimPEG import Mesh, np
import matplotlib.pyplot as plt
hx = np.r_[3,2,1,1,1,1,2,3]
hy = np.r_[3,1,1,3]
M = Mesh.TensorMesh([hx, hy])
M.plotGrid(centers=True)
plt.show()
In this simple mesh, the hx vector defines the widths of the cell
@@ -85,6 +88,7 @@ plotted above as red circles. Other terminology for this mesh are:
M.plotGrid(faces=True, nodes=True)
plt.title('Cell faces in the x- and y-directions.')
plt.legend(('Nodes', 'X-Faces', 'Y-Faces'))
plt.show()
Generally, the faces are used to discretize fluxes, quantities that
leave or enter the cells. As such, these fluxes have a direction to
@@ -102,7 +106,7 @@ and live on the edges(!) of the cell.
:include-source:
from SimPEG import Mesh
Mesh.TensorMesh([1,1,1]).plotGrid(faces=True, edges=True, centers=True)
Mesh.TensorMesh([1,1,1]).plotGrid(faces=True, edges=True, centers=True, showIt=True)
How many of each?
-----------------
@@ -159,7 +163,7 @@ vector grid size.
:include-source:
from SimPEG import Mesh
Mesh.TensorMesh([4,5]).plotGrid(faces=True)
Mesh.TensorMesh([4,5]).plotGrid(faces=True, showIt=True)
Making Tensors
@@ -183,7 +187,7 @@ notation::
from SimPEG import Mesh, Utils
h1 = (5, 10, 1.5), (20, 5), (3, 10)
M = Mesh.TensorMesh(Utils.meshTensors(h1, h1))
M.plotGrid()
M.plotGrid(showIt=True)
Hopefully, you now know how to create TensorMesh objects in SimPEG,
and by extension you are also familiar with how to create and use
docs/api_Model.rst
+45
View File
@@ -0,0 +1,45 @@
.. _api_Model:
Model
*****
A SimPEG model 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 (:class:`SimPEG.Model.LogModel`).
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 the term Model to
describe how to get between these two spaces.
The API
=======
.. autoclass:: SimPEG.Model.BaseModel
:members:
:undoc-members:
.. autoclass:: SimPEG.Model.BaseNonLinearModel
:members:
:undoc-members:
.. autoclass:: SimPEG.Model.ComboModel
:members:
:undoc-members:
Common Models
=============
.. autoclass:: SimPEG.Model.LogModel
:members:
:undoc-members:
.. autoclass:: SimPEG.Model.Vertical1DModel
:members:
:undoc-members:
.. autoclass:: SimPEG.Model.Mesh2Mesh
:members:
:undoc-members:
docs/api_Problem.rst
+71
View File
@@ -0,0 +1,71 @@
.. _api_Problem:
Problem
*******
The problem is a partial differential equation of the form:
.. math::
c(m, u) = 0
Here, \\(m\\) is the model and u is the field (or fields).
Given the model, \\(m\\), we can calculate the fields \\(u(m)\\),
however, the data we collect is a subset of the fields,
and can be defined by a linear projection, \\(P\\).
.. math::
d_\text{pred} = P u(m)
For the inverse problem, we are interested in how changing the model transforms the data,
as such we can take write the Taylor expansion:
.. math::
Pu(m + hv) = Pu(m) + hP\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \right\| )
We can linearize and define the sensitivity matrix as:
.. math::
J = P\frac{\partial u}{\partial m}
The sensitivity matrix, and it's transpose will be used in the inverse problem
to (locally) find how model parameters change the data, and optimize!
Working with the general PDE, \\(c(m, u) = 0\\), where m is the model and u is the field,
the sensitivity is defined as:
.. math::
J = P\frac{\partial u}{\partial m}
We can take the derivative of the PDE:
.. math::
\nabla_m c(m, u) \partial m + \nabla_u c(m, u) \partial u = 0
If the forward problem is invertible, then we can rearrange for \\(\\frac{\\partial u}{\\partial m}\\):
.. math::
J = - P \left( \nabla_u c(m, u) \right)^{-1} \nabla_m c(m, u)
This can often be computed given a vector (i.e. \\(J(v)\\)) rather than stored, as \\(J\\) is a large dense matrix.
.. math::
u(m)
The API
=======
.. automodule:: SimPEG.Problem
:members:
:undoc-members:
docs/api_Survey.rst
+11
View File
@@ -0,0 +1,11 @@
.. _api_Survey:
Survey
======
.. automodule:: SimPEG.Survey
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
docs/index.rst
+3 -1
View File
@@ -39,7 +39,9 @@ Forward Problems
.. toctree::
:maxdepth: 2
api_Forward
api_Model
api_Survey
api_Problem
Inversion
*********