diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index 18d9cfc3..cae85ba4 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -61,7 +61,7 @@ class BaseObjFunction(object): if self.debug: print 'Calling ObjFunction.startup' if self.reg.mref is None: - print 'Regularization has not set mref. SimPEG will set it to m0.' + print 'Regularization has not set mref. SimPEG.ObjFunction will set it to m0.' self.reg.mref = m0 self.phi_d = np.nan diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 8bb5c95c..7a89fa89 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -800,13 +800,15 @@ class NewtonRoot(object): """ tol = 1.000e-06 - solveTol = 0 # Default direct solve. maxIter = 20 stepDcr = 0.5 maxLS = 30 comments = False doLS = True + Solver = Solver + solverOpts = {} + def __init__(self, **kwargs): Utils.setKwargs(self, **kwargs) @@ -828,13 +830,9 @@ class NewtonRoot(object): while True: r, J = fun(x, return_g=True) - if self.solveTol == 0: - Jinv = Solver(J) - dh = - Jinv.solve(r) - else: - raise NotImplementedError('Iterative solve on NewtonRoot is not yet implemented.') - # M = @(x) tril(J)\(diag(J).*(triu(J)\x)); - # [dh, ~] = bicgstab(J,-r,O.solveTol,500,M); + + Jinv = self.Solver(J, **self.solverOpts) + dh = - Jinv.solve(r) muLS = 1. LScnt = 1 @@ -862,7 +860,7 @@ class NewtonRoot(object): if norm(rt) < self.tol: break if self.iter > self.maxIter: - print 'NewtonRoot stopped by maxIters. norm: %4.4e' % norm(rt) + print 'NewtonRoot stopped by maxIters (%d). norm: %4.4e' % (self.maxIter, norm(rt)) break return x diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py index 417051df..fa15ff67 100644 --- a/SimPEG/Parameters.py +++ b/SimPEG/Parameters.py @@ -56,7 +56,7 @@ class Parameter(object): if (self.current is None or not self.opt.iter == self.currentIter): self.current = self.nextIter() - self.currentIter = self.opt.iter + self.currentIter = getattr(self.opt, 'iter', 0) return self.current def nextIter(self): @@ -162,3 +162,16 @@ class BetaSchedule(BetaEstimate): self.beta /= self.coolingFactor return self.beta + + +class UpdateReferenceModel(Parameter): + + mref0 = None + + def nextIter(self): + mref = getattr(self, 'm_prev', None) + if mref is None: + if self.debug: print 'UpdateReferenceModel is using mref0' + mref = self.mref0 + self.m_prev = self.objFunc.m_current + return mref diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index c360f02c..477ba604 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -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) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 9ed4a157..09a2ca6c 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -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)) @@ -205,7 +209,10 @@ class BaseRx(object): def nD(self): return self.locs.shape[0] +<<<<<<< HEAD +======= +>>>>>>> dbd1334e0bf48dedc12f744841e71725a9d98d50 class BaseTx(object): """SimPEG Transmitter Object""" diff --git a/SimPEG/Utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py index c70c7019..462b9cfd 100644 --- a/SimPEG/Utils/ModelBuilder.py +++ b/SimPEG/Utils/ModelBuilder.py @@ -70,20 +70,37 @@ def getIndecesBlock(p0,p1,ccMesh): # Return a tuple return ind -def defineBlockConductivity(ccMesh,p0,p1,condVals): +def defineBlock(ccMesh,p0,p1,vals=[0,1]): """ Build a block with the conductivity specified by condVal. Returns an array. - condVals[0] conductivity of the block - condVals[1] conductivity of the ground + vals[0] conductivity of the block + vals[1] conductivity of the ground """ - sigma = np.zeros(ccMesh.shape[0]) + condVals[1] + sigma = np.zeros(ccMesh.shape[0]) + vals[1] ind = getIndecesBlock(p0,p1,ccMesh) - sigma[ind] = condVals[0] + sigma[ind] = vals[0] return sigma -def defineTwoLayeredConductivity(ccMesh,depth,condVals): +def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.): + G = ccMesh.copy() + dim = ccMesh.shape[1] + for i in range(dim): + G[:, i] = G[:,i] - center[i] + + theta = -theta*np.pi/180 + M = np.array([[np.cos(theta),-np.sin(theta),0],[np.sin(theta),np.cos(theta),0],[0,0,1.]]) + M = M[:dim,:dim] + G = M.dot(G.T).T + + for i in range(dim): + G[:, i] = G[:,i]/anisotropy[i]*2. + + D = np.sqrt(np.sum(G**2,axis=1)) + return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5 + +def defineTwoLayers(ccMesh,depth,vals=[0,1]): """ Define a two layered model. Depth of the first layer must be specified. CondVals vector with the conductivity values of the layers. Eg: @@ -94,7 +111,7 @@ def defineTwoLayeredConductivity(ccMesh,depth,condVals): 0 depth zf 1st layer 2nd layer """ - sigma = np.zeros(ccMesh.shape[0]) + condVals[1] + sigma = np.zeros(ccMesh.shape[0]) + vals[1] dim = np.size(ccMesh[0,:]) @@ -116,7 +133,7 @@ def defineTwoLayeredConductivity(ccMesh,depth,condVals): ind = getIndecesBlock(p0,p1,ccMesh) - sigma[ind] = condVals[0]; + sigma[ind] = vals[0]; return sigma @@ -230,9 +247,9 @@ if __name__ == '__main__': p0 = np.array([0.5,0.5,0.5])[:testDim] p1 = np.array([1.0,1.0,1.0])[:testDim] - condVals = np.array([100,1e-6]) + vals = np.array([100,1e-6]) - sigma = defineBlockConductivity(ccMesh,p0,p1,condVals) + sigma = defineBlockConductivity(ccMesh,p0,p1,vals) # Plot sigma model print sigma.shape @@ -242,10 +259,10 @@ if __name__ == '__main__': # ----------------------------------------- print('Testing the two layered model') - condVals = np.array([100,1e-5]); + vals = np.array([100,1e-5]); depth = 1.0; - sigma = defineTwoLayeredConductivity(ccMesh,depth,condVals) + sigma = defineTwoLayeredConductivity(ccMesh,depth,vals) M.plotImage(sigma) print sigma diff --git a/docs/api_Forward.rst b/docs/api_Forward.rst deleted file mode 100644 index 66d9d36d..00000000 --- a/docs/api_Forward.rst +++ /dev/null @@ -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: - diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index e5a72491..065dc58e 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -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 diff --git a/docs/api_Model.rst b/docs/api_Model.rst new file mode 100644 index 00000000..05963ff0 --- /dev/null +++ b/docs/api_Model.rst @@ -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: diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst new file mode 100644 index 00000000..0c371b55 --- /dev/null +++ b/docs/api_Problem.rst @@ -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: + diff --git a/docs/api_Survey.rst b/docs/api_Survey.rst new file mode 100644 index 00000000..25fd0e37 --- /dev/null +++ b/docs/api_Survey.rst @@ -0,0 +1,11 @@ +.. _api_Survey: + + +Survey +====== + +.. automodule:: SimPEG.Survey + :show-inheritance: + :members: + :undoc-members: + :inherited-members: diff --git a/docs/index.rst b/docs/index.rst index 71aaa282..98f4efb7 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -39,7 +39,9 @@ Forward Problems .. toctree:: :maxdepth: 2 - api_Forward + api_Model + api_Survey + api_Problem Inversion *********