diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 297270ef..8dbee17d 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -32,24 +32,6 @@ class BaseInversion(object): self.opt.printers.insert(2,IterationPrinters.phi_d) self.opt.printers.insert(3,IterationPrinters.phi_m) - #TODO: Move this to the data class? - @property - def phi_d_target(self): - """ - target for phi_d - - By default this is the number of data. - - Note that we do not set the target if it is None, but we return the default value. - """ - if getattr(self, '_phi_d_target', None) is None: - return self.data.dobs.size # - return self._phi_d_target - - @phi_d_target.setter - def phi_d_target(self, value): - self._phi_d_target = value - @Utils.timeIt def run(self, m0): """run(m0) @@ -70,50 +52,3 @@ class BaseInversion(object): **finish** is called at the end of the optimization. """ pass - - - def save(self, group): - group.attrs['phi_d'] = self.phi_d - group.attrs['phi_m'] = self.phi_m - group.setArray('m', self.m) - group.setArray('dpred', self.dpred) - - - -# class Inversion(Cooling, Remember, BaseInversion): - -# maxIter = 10 -# name = "SimPEG Inversion" - -# def __init__(self, prob, reg, opt, data, **kwargs): -# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) - -# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) - -# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: -# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) - -# class TimeSteppingInversion(Remember, BaseInversion): -# """ -# A slightly different view on regularization parameters, -# let Beta be viewed as 1/dt, and timestep by updating the -# reference model every optimization iteration. -# """ -# maxIter = 1 -# name = "Time-Stepping SimPEG Inversion" - -# def __init__(self, prob, reg, opt, data, **kwargs): -# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) - -# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) - -# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: -# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) - -# def _startup_TimeSteppingInversion(self, m0): - -# def _doEndIteration_updateMref(self, xt): -# if self.debug: 'Updating the reference model.' -# self.parent.reg.mref = self.xc - -# self.opt.hook(_doEndIteration_updateMref, overwrite=True) diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index 2b1912b8..18d9cfc3 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -1,7 +1,7 @@ -import Utils, Parameters, numpy as np, scipy.sparse as sp +import Utils, Parameters, Survey, Problem, numpy as np, scipy.sparse as sp class BaseObjFunction(object): - """BaseObjFunction(data, reg, **kwargs)""" + """BaseObjFunction(forward, reg, **kwargs)""" __metaclass__ = Utils.SimPEGMetaClass @@ -10,6 +10,9 @@ class BaseObjFunction(object): debug = False #: Print debugging information counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things + surveyPair = Survey.BaseSurvey + problemPair = Problem.BaseProblem + name = 'Base Objective Function' #: Name of the objective function u_current = None #: The most current evaluated field @@ -31,18 +34,19 @@ class BaseObjFunction(object): def objFunc(self): return self @property def opt(self): return getattr(self.parent,'opt',None) - @property - def prob(self): return self.data.prob - @property - def mesh(self): return self.data.prob.mesh - @property - def model(self): return self.data.prob.model - def __init__(self, data, reg, **kwargs): + def __init__(self, forward, reg, **kwargs): Utils.setKwargs(self, **kwargs) - self.data = data + assert forward.ispaired, 'The forward problem and survey must be paired.' + if isinstance(forward, self.surveyPair): + self.survey = forward + self.prob = forward.prob + elif isinstance(forward, self.problemPair): + self.prob = forward + self.survey = forward.survey + self.reg = reg self.reg.parent = self @@ -73,13 +77,13 @@ class BaseObjFunction(object): self.u_current = None self.m_current = m - u = self.data.prob.fields(m) + u = self.prob.fields(m) self.u_current = u phi_d = self.dataObj(m, u=u) phi_m = self.reg.modelObj(m) - self.dpred = self.data.dpred(m, u=u) # This is a cheap matrix vector calculation. + self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation. self.phi_d, self.phi_d_last = phi_d, self.phi_d self.phi_m, self.phi_m_last = phi_m, self.phi_m @@ -124,7 +128,7 @@ class BaseObjFunction(object): u is the field of interest; d_obs is the observed data; and W is the weighting matrix. """ # TODO: ensure that this is a data is vector and Wd is a matrix. - R = self.data.residualWeighted(m, u=u) + R = self.survey.residualWeighted(m, u=u) return 0.5*np.vdot(R, R) @Utils.timeIt @@ -160,11 +164,11 @@ class BaseObjFunction(object): \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} """ - if u is None: u = self.data.prob.fields(m) + if u is None: u = self.prob.fields(m) - R = self.data.residualWeighted(m, u=u) + R = self.survey.residualWeighted(m, u=u) - dmisfit = self.data.prob.Jtvec(m, self.data.Wd * R, u=u) + dmisfit = self.prob.Jtvec(m, self.survey.Wd * R, u=u) return dmisfit @@ -204,12 +208,12 @@ class BaseObjFunction(object): \\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J} """ - if u is None: u = self.data.prob.fields(m) + if u is None: u = self.prob.fields(m) - R = self.data.residualWeighted(m, u=u) + R = self.survey.residualWeighted(m, u=u) # TODO: abstract to different norms a little cleaner. # \/ it goes here. in l2 it is the identity. - dmisfit = self.data.prob.Jtvec_approx(m, self.data.Wd * self.data.Wd * self.data.prob.Jvec_approx(m, v, u=u), u=u) + dmisfit = self.prob.Jtvec_approx(m, self.survey.Wd * self.survey.Wd * self.prob.Jvec_approx(m, v, u=u), u=u) return dmisfit diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py index 5a6760e6..417051df 100644 --- a/SimPEG/Parameters.py +++ b/SimPEG/Parameters.py @@ -41,7 +41,7 @@ class Parameter(object): @property def reg(self): return self.parent.reg @property - def data(self): return self.parent.data + def survey(self): return self.parent.survey @property def prob(self): return self.parent.prob @property @@ -131,13 +131,13 @@ class BetaEstimate(Parameter): :return: beta0 """ objFunc = self.parent - data = objFunc.data + survey = objFunc.survey m = objFunc.m_current u = objFunc.u_current if u is None: - u = data.prob.fields(m) + u = survey.prob.fields(m) x0 = np.random.rand(*m.shape) t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u)) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index a509f6cd..c360f02c 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -1,4 +1,4 @@ -import Utils, Data, numpy as np, scipy.sparse as sp +import Utils, Survey, numpy as np, scipy.sparse as sp import Model class BaseProblem(object): @@ -38,40 +38,42 @@ class BaseProblem(object): counter = None #: A SimPEG.Utils.Counter object - dataPair = Data.BaseData + surveyPair = Survey.BaseSurvey modelPair = Model.BaseModel - def __init__(self, mesh, model, **kwargs): + def __init__(self, model, **kwargs): Utils.setKwargs(self, **kwargs) - self.mesh = mesh assert (isinstance(model, self.modelPair) or isinstance(model, Model.ComboModel) and isinstance(model.models[0], self.modelPair) ), "Model object must be an instance of a %s class."%(self.modelPair.__name__) self.model = model @property - def data(self): + def mesh(self): return self.model.mesh + + @property + def survey(self): """ - The data object for this problem. + The survey object for this problem. """ - return getattr(self, '_data', None) + return getattr(self, '_survey', None) def pair(self, d): - """Bind a data to this problem instance using pointers.""" - assert isinstance(d, self.dataPair), "Data object must be an instance of a %s class."%(self.dataPair.__name__) + """Bind a survey to this problem instance using pointers.""" + assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__) if d.ispaired: - raise Exception("The data object is already paired to a problem. Use data.unpair()") - self._data = d + raise Exception("The survey object is already paired to a problem. Use survey.unpair()") + self._survey = d d._prob = self def unpair(self): - """Unbind a data from this problem instance.""" + """Unbind a survey from this problem instance.""" if not self.ispaired: return - self.data._prob = None - self._data = None + self.survey._prob = None + self._survey = None @property - def ispaired(self): return self.data is not None + def ispaired(self): return self.survey is not None @Utils.timeIt def Jvec(self, m, v, u=None): @@ -156,25 +158,26 @@ class BaseProblem(object): """ pass - def createSyntheticData(self, m, std=0.05, u=None, **geometry_kwargs): + #TODO: Rename and refactor to createSyntheticData + def createSyntheticSurvey(self, m, std=0.05, u=None, **geometry_kwargs): """ - Create synthetic data given a model, and a standard deviation. + Create synthetic survey given a model, and a standard deviation. :param numpy.array m: geophysical model :param numpy.array std: standard deviation - :rtype: numpy.array, numpy.array - :return: dobs, Wd + :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. """ - data = self.dataPair(mtrue=m, **geometry_kwargs) - data.pair(self) - data.dtrue = data.dpred(m, u=u) - noise = std*abs(data.dtrue)*np.random.randn(*data.dtrue.shape) - data.dobs = data.dtrue+noise - data.std = data.dobs*0 + std - return data + survey = self.surveyPair(mtrue=m, **geometry_kwargs) + survey.pair(self) + survey.dtrue = survey.dpred(m, u=u) + noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape) + survey.dobs = survey.dtrue+noise + survey.std = survey.dobs*0 + std + return survey diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8dc32b79..8f7b9ede 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -46,7 +46,7 @@ class BaseRegularization(object): @property def prob(self): return self.parent.prob @property - def data(self): return self.parent.data + def survey(self): return self.parent.survey @property def mesh(self): return self.model.mesh diff --git a/SimPEG/Data.py b/SimPEG/Survey.py similarity index 55% rename from SimPEG/Data.py rename to SimPEG/Survey.py index 12696406..c640d865 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Survey.py @@ -1,8 +1,8 @@ import Utils, numpy as np -class BaseData(object): - """Data holds the observed data, and the standard deviations.""" +class BaseSurvey(object): + """Survey holds the observed data, and the standard deviations.""" __metaclass__ = Utils.SimPEGMetaClass @@ -19,25 +19,25 @@ class BaseData(object): @property def prob(self): """ - The geophysical problem that explains this data, use:: + The geophysical problem that explains this survey, use:: - data.pair(prob) + survey.pair(prob) """ return getattr(self, '_prob', None) def pair(self, p): - """Bind a problem to this data instance using pointers""" - assert hasattr(p, 'dataPair'), "Problem must have an attribute 'dataPair'." - assert isinstance(self, p.dataPair), "Problem requires data object must be an instance of a %s class."%(p.dataPair.__name__) + """Bind a problem to this survey instance using pointers""" + assert hasattr(p, 'surveyPair'), "Problem must have an attribute 'surveyPair'." + assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a %s class."%(p.surveyPair.__name__) if p.ispaired: - raise Exception("The problem object is already paired to a data. Use prob.unpair()") + raise Exception("The problem object is already paired to a survey. Use prob.unpair()") self._prob = p - p._data = self + p._survey = self def unpair(self): - """Unbind a problem from this data instance""" + """Unbind a problem from this survey instance""" if not self.ispaired: return - self.prob._data = None + self.prob._survey = None self._prob = None @property @@ -108,7 +108,7 @@ class BaseData(object): """ Data weighting matrix. This is a covariance matrix used in:: - def data.residualWeighted(m,u=None): + def residualWeighted(m,u=None): return self.Wd*self.residual(m, u=u) By default, this is based on the norm of the data plus a noise floor. @@ -139,20 +139,82 @@ class BaseData(object): """ return Utils.mkvc(self.Wd*self.residual(m, u=u)) - @property - def RHS(self): - """ - Source matrix. - """ - return getattr(self, '_RHS', None) - @RHS.setter - def RHS(self, value): - self._RHS = value - @property def isSynthetic(self): "Check if the data is synthetic." - return (self.mtrue is not None) + return self.mtrue is not None + + + #TODO: Move this to the data class? + # @property + # def phi_d_target(self): + # """ + # target for phi_d + + # By default this is the number of data. + + # Note that we do not set the target if it is None, but we return the default value. + # """ + # if getattr(self, '_phi_d_target', None) is None: + # return self.data.dobs.size # + # return self._phi_d_target + + # @phi_d_target.setter + # def phi_d_target(self, value): + # self._phi_d_target = value + + +class BaseRxList(object): + """SimPEG Receiver List Object""" + + locs = None #: Locations (nRx x 3) + + knownRxTypes = None #: Set this to a list of strings to ensure that txType is known + + def __init__(self, locs, rxType, **kwargs): + self.locs = locs + self.rxType = rxType + Utils.setKwargs(self, **kwargs) + + @property + def rxType(self): + """Receiver Type""" + return getattr(self, '_rxType', None) + @rxType.setter + def rxType(self, value): + known = self.knownRxTypes + if known is not None: + assert value in known, "rxType must be in ['%s']" % ("', '".join(known)) + self._rxType = value + +class BaseTx(object): + """SimPEG Transmitter Object""" + + loc = None #: Location [x,y,z] + + rxList = None #: SimPEG Receiver List + rxListPair = BaseRxList + + knownTxTypes = None #: Set this to a list of strings to ensure that txType is known + + def __init__(self, loc, txType, rxList, **kwargs): + assert isinstance(rxList, self.rxListPair), 'rxList must be a %s'%self.rxListPair.__name__ + + self.loc = loc + self.txType = txType + self.rxList = rxList + Utils.setKwargs(self, **kwargs) + + @property + def txType(self): + """Transmitter Type""" + return getattr(self, '_txType', None) + @txType.setter + def txType(self, value): + known = self.knownTxTypes + if known is not None: + assert value in known, "txType must be in ['%s']" % ("', '".join(known)) + self._txType = value if __name__ == '__main__': d = BaseData() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index a44cfd74..209d1cff 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -159,11 +159,11 @@ def requires(var): .. note:: - To use data.%s(), SimPEG requires that a problem be bound to the data. + To use survey.%s(), SimPEG requires that a problem be bound to the survey. If a problem has not been bound, an Exception will be raised. To bind a problem to the Data object:: - data.pair(myProblem) + survey.pair(myProblem) """ % f.__name__ else: diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 17f64d54..874f4ffe 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -5,7 +5,7 @@ from Solver import Solver import Mesh import Model import Problem -import Data +import Survey import Regularization import ObjFunction import Optimization diff --git a/Tutorials/Linear.py b/Tutorials/Linear.py index 7d3aa805..90714abc 100644 --- a/Tutorials/Linear.py +++ b/Tutorials/Linear.py @@ -5,8 +5,8 @@ import matplotlib.pyplot as plt class LinearProblem(Problem.BaseProblem): """docstring for LinearProblem""" - def __init__(self, mesh, model, G, **kwargs): - Problem.BaseProblem.__init__(self, mesh, model, **kwargs) + def __init__(self, model, G, **kwargs): + Problem.BaseProblem.__init__(self, model, **kwargs) self.G = G def fields(self, m, u=None): @@ -20,8 +20,7 @@ class LinearProblem(Problem.BaseProblem): def example(N): - h = np.ones(N)/N - M = Mesh.TensorMesh([h]) + M = Mesh.TensorMesh([N]) nk = 20 jk = np.linspace(1.,20.,nk) @@ -43,8 +42,8 @@ def example(N): model = Model.BaseModel(M) - prob = LinearProblem(M, model, G) - data = prob.createSyntheticData(mtrue, std=0.01) + prob = LinearProblem(model, G) + data = prob.createSyntheticSurvey(mtrue, std=0.01) return prob, data, model diff --git a/docs/simpeg-framework.png b/docs/simpeg-framework.png index ad8e5c1c..f90f8260 100644 Binary files a/docs/simpeg-framework.png and b/docs/simpeg-framework.png differ