From ca4932fd19d3f26b300e28392d05a8028de27d88 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 16 Jan 2014 13:19:29 -0800 Subject: [PATCH] Moved things around! Packages should now all be capitalized. may need to to tweak git to ensure this... --- SimPEG/{forward/SimPEGData.py => Data.py} | 121 +++++++++------ SimPEG/Model.py | 94 ++++++++++++ SimPEG/{forward => }/Problem.py | 94 ++---------- SimPEG/__init__.py | 17 ++- SimPEG/examples/DC.py | 82 +++++++---- SimPEG/examples/Linear.py | 10 +- SimPEG/forward/ModelTransforms.py | 49 ------- SimPEG/forward/__init__.py | 3 - SimPEG/inverse/Inversion.py | 34 ++--- SimPEG/inverse/Optimize.py | 66 ++++----- SimPEG/inverse/Regularization.py | 40 ++--- SimPEG/mesh/BaseMesh.py | 8 +- SimPEG/mesh/Cyl1DMesh.py | 4 +- SimPEG/mesh/DiffOperators.py | 2 +- SimPEG/mesh/InnerProducts.py | 2 +- SimPEG/mesh/LogicallyOrthogonalMesh.py | 94 ++++++------ SimPEG/mesh/LomView.py | 2 +- SimPEG/mesh/TensorMesh.py | 60 ++++---- SimPEG/mesh/TensorView.py | 2 +- SimPEG/setup.py | 4 +- SimPEG/tests/TestUtils.py | 12 +- SimPEG/tests/test_LogicallyOrthogonalMesh.py | 4 +- SimPEG/tests/test_Solver.py | 4 +- SimPEG/tests/test_basemesh.py | 2 +- SimPEG/tests/test_forward_DCproblem.py | 138 +++++++++--------- SimPEG/tests/test_interpolation.py | 2 +- SimPEG/tests/test_model.py | 27 ++++ SimPEG/tests/test_optimizers.py | 18 +-- ...est_forward_problem.py => test_problem.py} | 16 +- SimPEG/tests/test_tensorMesh.py | 2 +- SimPEG/tests/test_utils.py | 4 +- SimPEG/utils/ModelBuilder.py | 2 +- SimPEG/utils/Save.py | 4 +- SimPEG/utils/__init__.py | 2 +- SimPEG/utils/interputils.py | 4 +- SimPEG/utils/meshutils.py | 4 +- SimPEG/visualize/vtk/vtkTools.py | 2 +- 37 files changed, 546 insertions(+), 489 deletions(-) rename SimPEG/{forward/SimPEGData.py => Data.py} (51%) create mode 100644 SimPEG/Model.py rename SimPEG/{forward => }/Problem.py (63%) delete mode 100644 SimPEG/forward/ModelTransforms.py delete mode 100644 SimPEG/forward/__init__.py create mode 100644 SimPEG/tests/test_model.py rename SimPEG/tests/{test_forward_problem.py => test_problem.py} (51%) diff --git a/SimPEG/forward/SimPEGData.py b/SimPEG/Data.py similarity index 51% rename from SimPEG/forward/SimPEGData.py rename to SimPEG/Data.py index 934fd47a..bc874de2 100644 --- a/SimPEG/forward/SimPEGData.py +++ b/SimPEG/Data.py @@ -1,4 +1,4 @@ -from SimPEG import utils +import Utils def requiresProblem(f): @@ -32,58 +32,28 @@ def requiresProblem(f): return requiresProblemWrapper -class Data(object): +class BaseData(object): """Data holds the observed data, and the standard deviations.""" - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable - std = None #: Estimated Standard Deviations - dobs = None #: Observed data - dtrue = None #: True data, if data is synthetic - mtrue = None #: True model, if data is synthetic - prob = None #: The geophysical problem that explains this data + std = None #: Estimated Standard Deviations + dobs = None #: Observed data + dtrue = None #: True data, if data is synthetic + mtrue = None #: True model, if data is synthetic + prob = None #: The geophysical problem that explains this data + counter = None #: A SimPEG.Utils.Counter object def __init__(self, **kwargs): - utils.setKwargs(self, **kwargs) - - def isSynthetic(self): - "Check if the data is synthetic." - return (self.mtrue is not None) + Utils.setKwargs(self, **kwargs) def setProblem(self, prob): self.prob = prob - @property - def Wd(self): - """ - Standard deviation weighting matrix. - - By default, this is based on the norm of the data plus a noise floor. - - """ - if getattr(self,'_Wd',None) is None: - eps = np.linalg.norm(utils.mkvc(self.dobs),2)*1e-5 - self._Wd = 1/(abs(self.dobs)*self.std+eps) - return self._Wd - @Wd.setter - def Wd(self, value): - self._Wd = value - + @Utils.count @requiresProblem def dpred(self, m, u=None): - if u is None: u = self.prob.field(m) - - @requiresProblem - def residual(self, m, u=None): - if u is None: u = self.prob.field(m) - - @requiresProblem - def residualWeighted(self, m, u=None): - if u is None: u = self.prob.field(m) - - @requiresProblem - def projectField(self, m, u=None): """ Projection matrix. @@ -93,7 +63,74 @@ class Data(object): if u is None: u = self.prob.field(m) return self.P*u + @Utils.count + def residual(self, m, u=None): + """ + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: float + :return: data residual + + The data residual: + + .. math:: + + \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} + + """ + return self.dpred(m, u=u) - self.dobs + + + @property + def Wd(self): + """ + Data weighting matrix. This is a covariance matrix used in:: + + def data.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. + + """ + if getattr(self,'_Wd',None) is None: + eps = np.linalg.norm(Utils.mkvc(self.dobs),2)*1e-5 + self._Wd = 1/(abs(self.dobs)*self.std+eps) + return self._Wd + @Wd.setter + def Wd(self, value): + self._Wd = value + + def residualWeighted(self, m, u=None): + """ + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: float + :return: data residual + + The weighted data residual: + + .. math:: + + \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. + """ + return self.Wd*self.residual(m, u=u) + + @property + def RHS(self): + """ + Source matrix. + """ + return self._RHS + @RHS.setter + def RHS(self, value): + self._RHS = value + + def isSynthetic(self): + "Check if the data is synthetic." + return (self.mtrue is not None) if __name__ == '__main__': - d = SimPEGData() + d = BaseData() d.dpred() diff --git a/SimPEG/Model.py b/SimPEG/Model.py new file mode 100644 index 00000000..75912deb --- /dev/null +++ b/SimPEG/Model.py @@ -0,0 +1,94 @@ +from SimPEG import Utils, np, sp + + +class BaseModel(object): + """SimPEG Model""" + + __metaclass__ = Utils.Save.Savable + + counter = None #: A SimPEG.Utils.Counter object + + def __init__(self): + pass + + def transform(self, m): + """ + :param numpy.array m: model + :rtype: numpy.array + :return: transformed model + + The *transform* changes the model into the physical property. + + A common example of this is to invert for electrical conductivity + in log space. In this case, your model will be log(sigma) and to + get back to sigma, you can take the exponential: + + """ + return m + + def transformDeriv(self, m): + """ + :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) + + def example(self, mesh, type=None): + return np.random.rand(mesh.nC) + + + +class LogModel(BaseModel): + """SimPEG LogModel""" + + def __init__(self, **kwargs): + BaseModel.__init__(self, **kwargs) + + def transform(self, m): + """ + :param numpy.array m: model + :rtype: numpy.array + :return: transformed model + + The *transform* changes the model into the physical property. + + A common example of this is to invert for electrical conductivity + in log space. In this case, your model will be log(sigma) and to + get back to sigma, you can take the exponential: + + .. math:: + + m = \log{\sigma} + + \exp{m} = \exp{\log{\sigma}} = \sigma + """ + return np.exp(Utils.mkvc(m)) + + def transformDeriv(self, m): + """ + :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*. + + If the model *transform* is: + + .. math:: + + m = \log{\sigma} + + \exp{m} = \exp{\log{\sigma}} = \sigma + + Then the derivative is: + + .. math:: + + \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m}) + """ + return Utils.sdiag(np.exp(Utils.mkvc(m))) diff --git a/SimPEG/forward/Problem.py b/SimPEG/Problem.py similarity index 63% rename from SimPEG/forward/Problem.py rename to SimPEG/Problem.py index 1fb5de1f..e3e5fed7 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/Problem.py @@ -1,9 +1,8 @@ -from SimPEG import utils, np, sp -import SimPEGData +from SimPEG import Utils, np, sp, Data norm = np.linalg.norm -class Problem(object): +class BaseProblem(object): """ Problem is the base class for all geophysical forward problems in SimPEG. @@ -36,58 +35,19 @@ class Problem(object): to (locally) find how model parameters change the data, and optimize! """ - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable - counter = None #: A SimPEG.utils.Counter object + counter = None #: A SimPEG.Utils.Counter object + + dataPair = Data.BaseData - def __init__(self, mesh, *args, **kwargs): - utils.setKwargs(self, **kwargs) + def __init__(self, mesh, model, *args, **kwargs): + Utils.setKwargs(self, **kwargs) self.mesh = mesh + self.model = model - @property - def RHS(self): - """ - Source matrix. - """ - return self._RHS - @RHS.setter - def RHS(self, value): - self._RHS = value - - @utils.count - def dpred(self, m, u=None): - """ - Predicted data. - - .. math:: - d_\\text{pred} = Pu(m) - """ - if u is None: - u = self.field(m) - return self.P*u - - @utils.count - def dataResidual(self, m, data, u=None): - """ - :param numpy.array m: geophysical model - :param numpy.array u: fields - :rtype: float - :return: data misfit - - The data misfit: - - .. math:: - - \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data. - """ - - return self.dpred(m, u=u) - data.dobs - - @utils.timeIt + @Utils.timeIt def J(self, m, v, u=None): """ :param numpy.array m: model @@ -117,7 +77,7 @@ class Problem(object): """ raise NotImplementedError('J is not yet implemented.') - @utils.timeIt + @Utils.timeIt def Jt(self, m, v, u=None): """ :param numpy.array m: model @@ -131,7 +91,7 @@ class Problem(object): raise NotImplementedError('Jt is not yet implemented.') - @utils.timeIt + @Utils.timeIt def J_approx(self, m, v, u=None): """ @@ -146,7 +106,7 @@ class Problem(object): """ return self.J(m, v, u) - @utils.timeIt + @Utils.timeIt def Jt_approx(self, m, v, u=None): """ :param numpy.array m: model @@ -170,32 +130,6 @@ class Problem(object): """ pass - def modelTransform(self, m): - """ - :param numpy.array m: model - :rtype: numpy.array - :return: transformed model - - The modelTransform changes the model into the physical property. - - A common example of this is to invert for electrical conductivity - in log space. In this case, your model will be log(sigma) and to - get back to sigma, you can take the exponential: - - """ - return m - - def modelTransformDeriv(self, m): - """ - :param numpy.array m: model - :rtype: scipy.csr_matrix - :return: derivative of transformed model - - The modelTransform changes the model into the physical property. - The modelTransformDeriv provides the derivative of the modelTransform. - """ - return sp.identity(m.size) - def createSyntheticData(self, m, std=0.05, u=None): """ Create synthetic data given a model, and a standard deviation. @@ -212,7 +146,7 @@ class Problem(object): noise = std*abs(dtrue)*np.random.randn(*dtrue.shape) dobs = dtrue+noise stdev = dobs*0 + std - return SimPEGData.Data(dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m) + return self.dataPair(dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 07397233..904766d2 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,14 +1,15 @@ import numpy as np import scipy.sparse as sp -import utils -Solver = utils.Solver -import mesh -import forward -import inverse -import examples -import tests +import Utils +Solver = Utils.Solver +import Mesh +import Model +import Problem +import Data +import Inverse +import Examples +import Tests -Data = forward.Data import scipy.version as _v if _v.version < '0.13.0': diff --git a/SimPEG/examples/DC.py b/SimPEG/examples/DC.py index f975b108..fa50cd1c 100644 --- a/SimPEG/examples/DC.py +++ b/SimPEG/examples/DC.py @@ -1,20 +1,56 @@ from SimPEG import * -class DCProblem(forward.ModelTransforms.LogModel, forward.Problem): + + +class DCData(Data.BaseData): + """ + **DCData** + + Geophysical DC resistivity data. + + """ + + def __init__(self, mesh, model, **kwargs): + problem.BaseProblem.__init__(self, mesh, model) + self.mesh.setCellGradBC('neumann') + Utils.setKwargs(self, **kwargs) + + def reshapeFields(self, u): + if len(u.shape) == 1: + u = u.reshape([-1, self.RHS.shape[1]], order='F') + return u + + def dpred(self, m, u=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pu(m) + """ + if u is None: + u = self.field(m) + + u = self.reshapeFields(u) + + return Utils.mkvc(self.P*u) + + + +class DCProblem(Problem.BaseProblem): """ **DCProblem** Geophysical DC resistivity problem. """ - def __init__(self, mesh): - forward.Problem.__init__(self, mesh) - self.mesh.setCellGradBC('neumann') - def reshapeFields(self, u): - if len(u.shape) == 1: - u = u.reshape([-1, self.RHS.shape[1]], order='F') - return u + dataPair = DCData + + def __init__(self, mesh, model, **kwargs): + problem.BaseProblem.__init__(self, mesh, model) + self.mesh.setCellGradBC('neumann') + Utils.setKwargs(self, **kwargs) + def createMatrix(self, m): """ @@ -31,30 +67,16 @@ class DCProblem(forward.ModelTransforms.LogModel, forward.Problem): """ D = self.mesh.faceDiv G = self.mesh.cellGrad - sigma = self.modelTransform(m) + sigma = self.model.transform(m) Msig = self.mesh.getFaceMass(sigma) A = D*Msig*G return A.tocsc() - def dpred(self, m, u=None): - """ - Predicted data. - - .. math:: - d_\\text{pred} = Pu(m) - """ - if u is None: - u = self.field(m) - - u = self.reshapeFields(u) - - return utils.mkvc(self.P*u) - def field(self, m): A = self.createMatrix(m) solve = Solver(A) phi = solve.solve(self.RHS) - return utils.mkvc(phi) + return Utils.mkvc(phi) def J(self, m, v, u=None): """ @@ -88,17 +110,17 @@ class DCProblem(forward.ModelTransforms.LogModel, forward.Problem): G = self.mesh.cellGrad A = self.createMatrix(m) Av_dm = self.mesh.getFaceMassDeriv() - mT_dm = self.modelTransformDeriv(m) + mT_dm = self.model.transformDeriv(m) dCdu = A dCdm = np.empty_like(u) for i, ui in enumerate(u.T): # loop over each column - dCdm[:, i] = D * ( utils.sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) ) + dCdm[:, i] = D * ( Utils.sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) ) solve = Solver(dCdu) Jv = - P * solve.solve(dCdm) - return utils.mkvc(Jv) + return Utils.mkvc(Jv) def Jt(self, m, v, u=None): """Takes data, turns it into a model..ish""" @@ -114,7 +136,7 @@ class DCProblem(forward.ModelTransforms.LogModel, forward.Problem): G = self.mesh.cellGrad A = self.createMatrix(m) Av_dm = self.mesh.getFaceMassDeriv() - mT_dm = self.modelTransformDeriv(m) + mT_dm = self.model.transformDeriv(m) dCdu = A.T solve = Solver(dCdu) @@ -123,7 +145,7 @@ class DCProblem(forward.ModelTransforms.LogModel, forward.Problem): Jtv = 0 for i, ui in enumerate(u.T): # loop over each column - Jtv += utils.sdiag( G * ui ) * ( D.T * w[:,i] ) + Jtv += Utils.sdiag( G * ui ) * ( D.T * w[:,i] ) Jtv = - mT_dm.T * ( Av_dm.T * Jtv ) return Jtv @@ -174,7 +196,7 @@ if __name__ == '__main__': p0 = [5, 10] p1 = [15, 50] condVals = [sig1, sig2] - mSynth = utils.ModelBuilder.defineBlockConductivity(p0,p1,M.gridCC,condVals) + mSynth = Utils.ModelBuilder.defineBlockConductivity(p0,p1,M.gridCC,condVals) plt.colorbar(M.plotImage(mSynth)) plt.show() diff --git a/SimPEG/examples/Linear.py b/SimPEG/examples/Linear.py index dbee7dff..618ff260 100644 --- a/SimPEG/examples/Linear.py +++ b/SimPEG/examples/Linear.py @@ -1,12 +1,12 @@ -from SimPEG import mesh, forward, inverse, np +from SimPEG import Mesh, Model, Problem, Data, Inverse, np import matplotlib.pyplot as plt -class LinearProblem(forward.Problem): +class LinearProblem(Problem.BaseProblem): """docstring for LinearProblem""" def __init__(self, *args, **kwargs): - forward.Problem.__init__(self, *args, **kwargs) + problem.BaseProblem.__init__(self, *args, **kwargs) def dpred(self, m, u=None): return self.G.dot(m) @@ -39,7 +39,9 @@ def example(N): mtrue[M.vectorCCx > 0.45] = -0.5 mtrue[M.vectorCCx > 0.6] = 0 - prob = LinearProblem(M) + + + prob = LinearProblem(M, None) prob.G = G data = prob.createSyntheticData(mtrue, std=0.01) diff --git a/SimPEG/forward/ModelTransforms.py b/SimPEG/forward/ModelTransforms.py deleted file mode 100644 index ea89b974..00000000 --- a/SimPEG/forward/ModelTransforms.py +++ /dev/null @@ -1,49 +0,0 @@ -import numpy as np -from SimPEG.utils import mkvc, sdiag - -class LogModel(object): - """docstring for LogModel""" - def modelTransform(self, m): - """ - :param numpy.array m: model - :rtype: numpy.array - :return: transformed model - - The modelTransform changes the model into the physical property. - - A common example of this is to invert for electrical conductivity - in log space. In this case, your model will be log(sigma) and to - get back to sigma, you can take the exponential: - - .. math:: - - m = \log{\sigma} - - \exp{m} = \exp{\log{\sigma}} = \sigma - """ - return np.exp(mkvc(m)) - - def modelTransformDeriv(self, m): - """ - :param numpy.array m: model - :rtype: scipy.csr_matrix - :return: derivative of transformed model - - The modelTransform changes the model into the physical property. - The modelTransformDeriv provides the derivative of the modelTransform. - - If the model transform is: - - .. math:: - - m = \log{\sigma} - - \exp{m} = \exp{\log{\sigma}} = \sigma - - Then the derivative is: - - .. math:: - - \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m}) - """ - return sdiag(np.exp(mkvc(m))) diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py deleted file mode 100644 index 47153f5c..00000000 --- a/SimPEG/forward/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from Problem import * -import ModelTransforms -from SimPEGData import * diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 44b3e57b..1679c789 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,14 +1,14 @@ import SimPEG -from SimPEG import utils, sp, np +from SimPEG import Utils, sp, np from Optimize import Remember from BetaSchedule import Cooling -from SimPEG.inverse import IterationPrinters, StoppingCriteria +from SimPEG.Inverse import IterationPrinters, StoppingCriteria class BaseInversion(object): """BaseInversion(prob, reg, opt, data, **kwargs) """ - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable maxIter = 1 #: Maximum number of iterations name = 'BaseInversion' @@ -16,13 +16,13 @@ class BaseInversion(object): debug = False #: Print debugging information comment = '' #: Used by some functions to indicate what is going on in the algorithm - counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things + counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things beta0 = None #: The initial Beta (regularization parameter) beta0_ratio = 0.1 #: When beta0 is set to None, estimateBeta0 is used with this ratio def __init__(self, prob, reg, opt, data, **kwargs): - utils.setKwargs(self, **kwargs) + Utils.setKwargs(self, **kwargs) self.prob = prob self.reg = reg self.opt = opt @@ -59,7 +59,7 @@ class BaseInversion(object): def phi_d_target(self, value): self._phi_d_target = value - @utils.timeIt + @Utils.timeIt def run(self, m0): """run(m0) @@ -78,7 +78,7 @@ class BaseInversion(object): return self.m - @utils.callHooks('startup') + @Utils.callHooks('startup') def startup(self, m0): """ **startup** is called at the start of any new run call. @@ -98,7 +98,7 @@ class BaseInversion(object): self.phi_d_last = np.nan self.phi_m_last = np.nan - @utils.callHooks('doStartIteration') + @Utils.callHooks('doStartIteration') def doStartIteration(self): """ **doStartIteration** is called at the end of each run iteration. @@ -109,7 +109,7 @@ class BaseInversion(object): self._beta = self.getBeta() - @utils.callHooks('doEndIteration') + @Utils.callHooks('doEndIteration') def doEndIteration(self): """ **doEndIteration** is called at the end of each run iteration. @@ -168,7 +168,7 @@ class BaseInversion(object): def stoppingCriteria(self): if self.debug: print 'checking stoppingCriteria' - return utils.checkStoppers(self, self.stoppers) + return Utils.checkStoppers(self, self.stoppers) def printDone(self): @@ -176,9 +176,9 @@ class BaseInversion(object): **printDone** is called at the end of the inversion routine. """ - utils.printStoppers(self, self.stoppers) + Utils.printStoppers(self, self.stoppers) - @utils.callHooks('finish') + @Utils.callHooks('finish') def finish(self): """finish() @@ -186,7 +186,7 @@ class BaseInversion(object): """ pass - @utils.timeIt + @Utils.timeIt def evalFunction(self, m, return_g=True, return_H=True): """evalFunction(m, return_g=True, return_H=True) @@ -226,7 +226,7 @@ class BaseInversion(object): out += (operator,) return out if len(out) > 1 else out[0] - @utils.timeIt + @Utils.timeIt def dataObj(self, m, u=None): """dataObj(m, u=None) @@ -246,10 +246,10 @@ class BaseInversion(object): """ # TODO: ensure that this is a data is vector and Wd is a matrix. R = self.Wd*self.prob.dataResidual(m, self.data, u=u) - R = utils.mkvc(R) + R = Utils.mkvc(R) return 0.5*np.vdot(R, R) - @utils.timeIt + @Utils.timeIt def dataObjDeriv(self, m, u=None): """dataObjDeriv(m, u=None) @@ -291,7 +291,7 @@ class BaseInversion(object): return dmisfit - @utils.timeIt + @Utils.timeIt def dataObj2Deriv(self, m, v, u=None): """dataObj2Deriv(m, v, u=None) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 9d5edcc2..f11cd7b1 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -1,4 +1,4 @@ -from SimPEG import Solver, utils, sp, np +from SimPEG import Solver, Utils, sp, np import matplotlib.pyplot as plt norm = np.linalg.norm @@ -82,7 +82,7 @@ class Minimize(object): Minimize is a general class for derivative based optimization. """ - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable name = "General Optimization Algorithm" #: The name of the optimization algorithm @@ -100,7 +100,7 @@ class Minimize(object): debugLS = False #: Print debugging information for the line-search comment = '' #: Used by some functions to indicate what is going on in the algorithm - counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things + counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things def __init__(self, **kwargs): self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration] @@ -109,9 +109,9 @@ class Minimize(object): self.printers = [IterationPrinters.iteration, IterationPrinters.f, IterationPrinters.norm_g, IterationPrinters.totalLS] self.printersLS = [IterationPrinters.iterationLS, IterationPrinters.LS_ft, IterationPrinters.LS_t, IterationPrinters.LS_armijoGoldstein] - utils.setKwargs(self, **kwargs) + Utils.setKwargs(self, **kwargs) - @utils.timeIt + @Utils.timeIt def minimize(self, evalFunction, x0): """minimize(evalFunction, x0) @@ -189,7 +189,7 @@ class Minimize(object): def parent(self, value): self._parent = value - @utils.callHooks('startup') + @Utils.callHooks('startup') def startup(self, x0): """ **startup** is called at the start of any new minimize call. @@ -214,8 +214,8 @@ class Minimize(object): self.f_last = np.nan self.x_last = x0 - @utils.count - @utils.callHooks('doStartIteration') + @Utils.count + @Utils.callHooks('doStartIteration') def doStartIteration(self): """doStartIteration() @@ -237,9 +237,9 @@ class Minimize(object): """ pad = ' '*10 if inLS else '' name = self.name if not inLS else self.nameLS - utils.printTitles(self, self.printers if not inLS else self.printersLS, name, pad) + Utils.printTitles(self, self.printers if not inLS else self.printersLS, name, pad) - @utils.callHooks('printIter') + @Utils.callHooks('printIter') def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -249,7 +249,7 @@ class Minimize(object): """ pad = ' '*10 if inLS else '' - utils.printLine(self, self.printers if not inLS else self.printersLS, pad=pad) + Utils.printLine(self, self.printers if not inLS else self.printersLS, pad=pad) def printDone(self, inLS=False): """ @@ -262,10 +262,10 @@ class Minimize(object): pad = ' '*10 if inLS else '' stop, done = (' STOP! ', ' DONE! ') if not inLS else ('----------------', ' End Linesearch ') stoppers = self.stoppers if not inLS else self.stoppersLS - utils.printStoppers(self, stoppers, pad='', stop=stop, done=done) + Utils.printStoppers(self, stoppers, pad='', stop=stop, done=done) - @utils.callHooks('finish') + @Utils.callHooks('finish') def finish(self): """finish() @@ -281,10 +281,10 @@ class Minimize(object): if self._iter == 0: self.f0 = self.f self.g0 = self.g - return utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) + return Utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) - @utils.timeIt - @utils.callHooks('projection') + @Utils.timeIt + @Utils.callHooks('projection') def projection(self, p): """projection(p) @@ -298,7 +298,7 @@ class Minimize(object): """ return p - @utils.timeIt + @Utils.timeIt def findSearchDirection(self): """findSearchDirection() @@ -329,7 +329,7 @@ class Minimize(object): """ return -self.g - @utils.count + @Utils.count def scaleSearchDirection(self, p): """scaleSearchDirection(p) @@ -348,7 +348,7 @@ class Minimize(object): nameLS = "Armijo linesearch" #: The line-search name - @utils.timeIt + @Utils.timeIt def modifySearchDirection(self, p): """modifySearchDirection(p) @@ -386,7 +386,7 @@ class Minimize(object): return self._LS_xt, self._iterLS < self.maxIterLS - @utils.count + @Utils.count def modifySearchDirectionBreak(self, p): """modifySearchDirectionBreak(p) @@ -408,8 +408,8 @@ class Minimize(object): print 'The linesearch got broken. Boo.' return p, False - @utils.count - @utils.callHooks('doEndIteration') + @Utils.count + @Utils.callHooks('doEndIteration') def doEndIteration(self, xt): """doEndIteration(xt) @@ -527,7 +527,7 @@ class ProjectedGradient(Minimize, Remember): self.aSet_prev = self.activeSet(x0) - @utils.count + @Utils.count def projection(self, x): """projection(x) @@ -536,7 +536,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.median(np.c_[self.lower,x,self.upper],axis=1) - @utils.count + @Utils.count def activeSet(self, x): """activeSet(x) @@ -545,7 +545,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.logical_or(x == self.lower, x == self.upper) - @utils.count + @Utils.count def inactiveSet(self, x): """inactiveSet(x) @@ -554,7 +554,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.logical_not(self.activeSet(x)) - @utils.count + @Utils.count def bindingSet(self, x): """bindingSet(x) @@ -567,7 +567,7 @@ class ProjectedGradient(Minimize, Remember): bind_low = np.logical_and(x == self.upper, self.g <= 0) return np.logical_or(bind_up, bind_low) - @utils.timeIt + @Utils.timeIt def findSearchDirection(self): """findSearchDirection() @@ -612,7 +612,7 @@ class ProjectedGradient(Minimize, Remember): # aSet_after = self.activeSet(self.xc+p) return p - @utils.timeIt + @Utils.timeIt def _doEndIteration_ProjectedGradient(self, xt): """_doEndIteration_ProjectedGradient(xt)""" aSet = self.activeSet(xt) @@ -718,7 +718,7 @@ class GaussNewton(Minimize, Remember): def __init__(self, **kwargs): Minimize.__init__(self, **kwargs) - @utils.timeIt + @Utils.timeIt def findSearchDirection(self): return Solver(self.H).solve(-self.g) @@ -765,7 +765,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value - @utils.timeIt + @Utils.timeIt def findSearchDirection(self): Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG}) p = Hinv.solve(-self.g) @@ -778,7 +778,7 @@ class SteepestDescent(Minimize, Remember): def __init__(self, **kwargs): Minimize.__init__(self, **kwargs) - @utils.timeIt + @Utils.timeIt def findSearchDirection(self): return -self.g @@ -811,7 +811,7 @@ class NewtonRoot(object): doLS = True def __init__(self, **kwargs): - utils.setKwargs(self, **kwargs) + Utils.setKwargs(self, **kwargs) def root(self, fun, x): """root(fun, x) @@ -885,7 +885,7 @@ if __name__ == '__main__': print 'test the newtonRoot finding.' - fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), utils.sdiag( np.cos(x) ) ) + fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), Utils.sdiag( np.cos(x) ) ) x = np.array([np.pi-0.3, np.pi+0.1, 0]) pnt = NewtonRoot(comments=True).root(fun,x) print pnt diff --git a/SimPEG/inverse/Regularization.py b/SimPEG/inverse/Regularization.py index 8c6f19ef..1569d767 100644 --- a/SimPEG/inverse/Regularization.py +++ b/SimPEG/inverse/Regularization.py @@ -1,4 +1,4 @@ -from SimPEG import utils, np, sp +from SimPEG import Utils, np, sp class Regularization(object): """**Regularization** @@ -83,20 +83,20 @@ class Regularization(object): """ - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable - alpha_s = utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight") - alpha_x = utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") - alpha_y = utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") - alpha_z = utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") - alpha_xx = utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") - alpha_yy = utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") - alpha_zz = utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") + alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight") + alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") + alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") + alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") counter = None def __init__(self, mesh, **kwargs): - utils.setKwargs(self, **kwargs) + Utils.setKwargs(self, **kwargs) self.mesh = mesh @property @@ -112,7 +112,7 @@ class Regularization(object): def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) return self._Ws @property @@ -120,7 +120,7 @@ class Regularization(object): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFv[0]].T*self.mesh.vol - self._Wx = utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx + self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx return self._Wx @property @@ -128,7 +128,7 @@ class Regularization(object): """Regularization matrix Wy""" if getattr(self, '_Wy', None) is None: Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2])].T*self.mesh.vol - self._Wy = utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady + self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady return self._Wy @property @@ -136,28 +136,28 @@ class Regularization(object): """Regularization matrix Wz""" if getattr(self, '_Wz', None) is None: Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.nFv[:2]):].T*self.mesh.vol - self._Wz = utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz + self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz return self._Wz @property def Wxx(self): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: - self._Wxx = utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx return self._Wxx @property def Wyy(self): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: - self._Wyy = utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady return self._Wyy @property def Wzz(self): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: - self._Wzz = utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz return self._Wzz @@ -174,12 +174,12 @@ class Regularization(object): return self._W - @utils.timeIt + @Utils.timeIt def modelObj(self, m): r = self.W * (m - self.mref) return 0.5*r.dot(r) - @utils.timeIt + @Utils.timeIt def modelObjDeriv(self, m): """ @@ -198,7 +198,7 @@ class Regularization(object): """ return self.W.T * ( self.W * (m - self.mref) ) - @utils.timeIt + @Utils.timeIt def modelObj2Deriv(self): """ diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 6151d20b..f034aa5f 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -1,5 +1,5 @@ import numpy as np -from SimPEG import utils +from SimPEG import Utils class BaseMesh(object): @@ -78,7 +78,7 @@ class BaseMesh(object): x_array = np.ones((x.size, len(x))) # Unwrap it and put it in a np array for i, xi in enumerate(x): - x_array[:, i] = utils.mkvc(xi) + x_array[:, i] = Utils.mkvc(xi) x = x_array assert type(x) == np.ndarray, "x must be a numpy array" @@ -91,7 +91,7 @@ class BaseMesh(object): if format == 'M': return xx.reshape(nn, order='F') elif format == 'V': - return utils.mkvc(xx) + return Utils.mkvc(xx) def switchKernal(xx): """Switches over the different options.""" @@ -101,7 +101,7 @@ class BaseMesh(object): return outKernal(xx, nn) elif xType in ['F', 'E']: # This will only deal with components of fields, not full 'F' or 'E' - xx = utils.mkvc(xx) # unwrap it in case it is a matrix + xx = Utils.mkvc(xx) # unwrap it in case it is a matrix nn = self.nFv if xType == 'F' else self.nEv nn = np.r_[0, nn] diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index e22e12b9..2b291c54 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -1,7 +1,7 @@ import numpy as np import scipy.sparse as sp from scipy.constants import pi -from SimPEG.utils import mkvc, ndgrid, sdiag +from SimPEG.Utils import mkvc, ndgrid, sdiag class Cyl1DMesh(object): """ @@ -84,7 +84,7 @@ class Cyl1DMesh(object): doc = "Total number of cells in each direction" fget = lambda self: np.array([self.nCx, self.nCz]) return locals() - nCv = property(**nCv()) + nCv = property(**nCv()) def nNr(): doc = "Number of nodes in the radial direction" diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index d384ca17..1e8c1ccb 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -1,6 +1,6 @@ import numpy as np from scipy import sparse as sp -from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap +from SimPEG.Utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap def checkBC(bc): diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index 45853071..9e7c1d93 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal import numpy as np diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index b3dcd095..a0495630 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -1,4 +1,4 @@ -from SimPEG import utils, np +from SimPEG import Utils, np from BaseMesh import BaseMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts @@ -7,8 +7,8 @@ from LomView import LomView # Some helper functions. length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 -normalize2D = lambda x: x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2)) -normalize3D = lambda x: x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) +normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)) +normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): @@ -21,7 +21,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): """ - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable _meshType = 'LOM' @@ -40,7 +40,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): # Save nodes to private variable _gridN as vectors self._gridN = np.ones((nodes[0].size, self.dim)) for i, node_i in enumerate(nodes): - self._gridN[:, i] = utils.mkvc(node_i.astype(float)) + self._gridN[:, i] = Utils.mkvc(node_i.astype(float)) def gridCC(): doc = "Cell-centered grid." @@ -71,10 +71,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridFx is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] + XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] self._gridFx = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] + XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridFx return locals() @@ -88,10 +88,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridFy is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] + XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] self._gridFy = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] + XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridFy return locals() @@ -104,7 +104,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridFz is None and self.dim == 3: N = self.r(self.gridN, 'N', 'N', 'M') - XYZ = [utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] + XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridFz return locals() @@ -118,10 +118,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridEx is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] + XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] self._gridEx = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] + XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridEx return locals() @@ -135,10 +135,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridEy is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] + XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] self._gridEy = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] + XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridEy return locals() @@ -151,7 +151,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridEz is None and self.dim == 3: N = self.r(self.gridN, 'N', 'N', 'M') - XYZ = [utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] + XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridEz return locals() @@ -194,25 +194,25 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if(self._vol is None): if self.dim == 2: - A, B, C, D = utils.indexCube('ABCD', self.n+1) - normal, area = utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D) + A, B, C, D = Utils.indexCube('ABCD', self.n+1) + normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D) self._vol = area elif self.dim == 3: # Each polyhedron can be decomposed into 5 tetrahedrons # However, this presents a choice so we may as well divide in two ways and average. - A, B, C, D, E, F, G, H = utils.indexCube('ABCDEFGH', self.n+1) + A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.n+1) - vol1 = (utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top - utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top - utils.volTetra(self.gridN, B, D, E, G) + # middle - utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom - utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom + vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top + Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top + Utils.volTetra(self.gridN, B, D, E, G) + # middle + Utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom + Utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom - vol2 = (utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top - utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top - utils.volTetra(self.gridN, A, H, F, C) + # middle - utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom - utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom + vol2 = (Utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top + Utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top + Utils.volTetra(self.gridN, A, H, F, C) + # middle + Utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom + Utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom self._vol = (vol1 + vol2)/2 return self._vol @@ -228,30 +228,30 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): # Compute areas of cell faces if(self.dim == 2): xy = self.gridN - A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) + A, B = Utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) edge1 = xy[B, :] - xy[A, :] normal1 = np.c_[edge1[:, 1], -edge1[:, 0]] area1 = length2D(edge1) - A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) + A, D = Utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) # Note that we are doing A-D to make sure the normal points the right way. # Think about it. Look at the picture. Normal points towards C iff you do this. edge2 = xy[A, :] - xy[D, :] normal2 = np.c_[edge2[:, 1], -edge2[:, 0]] area2 = length2D(edge2) - self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)] + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] self._normals = [normalize2D(normal1), normalize2D(normal2)] elif(self.dim == 3): - A, E, F, B = utils.indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz])) - normal1, area1 = utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False) + A, E, F, B = Utils.indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz])) + normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False) - A, D, H, E = utils.indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz])) - normal2, area2 = utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False) + A, D, H, E = Utils.indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz])) + normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False) - A, B, C, D = utils.indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz])) - normal3, area3 = utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False) + A, B, C, D = Utils.indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz])) + normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False) - self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)] + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] self._normals = [normal1, normal2, normal3] return self._area return locals() @@ -291,21 +291,21 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if(self._edge is None or self._tangents is None): if(self.dim == 2): xy = self.gridN - A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) + A, D = Utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) edge1 = xy[D, :] - xy[A, :] - A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) + A, B = Utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) edge2 = xy[B, :] - xy[A, :] - self._edge = np.r_[utils.mkvc(length2D(edge1)), utils.mkvc(length2D(edge2))] + self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))] self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge] elif(self.dim == 3): xyz = self.gridN - A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz])) + A, D = Utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz])) edge1 = xyz[D, :] - xyz[A, :] - A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz])) + A, B = Utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz])) edge2 = xyz[B, :] - xyz[A, :] - A, E = utils.indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz])) + A, E = Utils.indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz])) edge3 = xyz[E, :] - xyz[A, :] - self._edge = np.r_[utils.mkvc(length3D(edge1)), utils.mkvc(length3D(edge2)), utils.mkvc(length3D(edge3))] + self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))] self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge] return self._edge return locals() @@ -331,10 +331,10 @@ if __name__ == '__main__': h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)]) dee3 = True if dee3: - X, Y, Z = utils.ndgrid(h1, h2, h3, vector=False) + X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False) M = LogicallyOrthogonalMesh([X, Y, Z]) else: - X, Y = utils.ndgrid(h1, h2, vector=False) + X, Y = Utils.ndgrid(h1, h2, vector=False) M = LogicallyOrthogonalMesh([X, Y]) print M.r(M.normals, 'F', 'Fx', 'V') diff --git a/SimPEG/mesh/LomView.py b/SimPEG/mesh/LomView.py index 2a9b242b..8243a248 100644 --- a/SimPEG/mesh/LomView.py +++ b/SimPEG/mesh/LomView.py @@ -2,7 +2,7 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib from mpl_toolkits.mplot3d import Axes3D -from SimPEG.utils import mkvc +from SimPEG.Utils import mkvc class LomView(object): diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index d4ab86b0..544ec716 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -1,4 +1,4 @@ -from SimPEG import utils, np, sp +from SimPEG import Utils, np, sp from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators @@ -23,8 +23,8 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): .. plot:: - from SimPEG import mesh, utils - M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + from SimPEG import mesh, Utils + M = mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) M.plotGrid() For a quick tensor mesh on a (10x12x15) unit cube:: @@ -33,7 +33,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): """ - __metaclass__ = utils.Save.Savable + __metaclass__ = Utils.Save.Savable _meshType = 'TENSOR' @@ -52,7 +52,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" # Ensure h contains 1D vectors - self._h = [utils.mkvc(x.astype(float)) for x in h] + self._h = [Utils.mkvc(x.astype(float)) for x in h] def __str__(self): outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) @@ -170,7 +170,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridCC is None: - self._gridCC = utils.ndgrid(self.getTensor('CC')) + self._gridCC = Utils.ndgrid(self.getTensor('CC')) return self._gridCC return locals() _gridCC = None # Store grid by default @@ -181,7 +181,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridN is None: - self._gridN = utils.ndgrid(self.getTensor('N')) + self._gridN = Utils.ndgrid(self.getTensor('N')) return self._gridN return locals() _gridN = None # Store grid by default @@ -192,7 +192,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFx is None: - self._gridFx = utils.ndgrid(self.getTensor('Fx')) + self._gridFx = Utils.ndgrid(self.getTensor('Fx')) return self._gridFx return locals() _gridFx = None # Store grid by default @@ -203,7 +203,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFy is None and self.dim > 1: - self._gridFy = utils.ndgrid(self.getTensor('Fy')) + self._gridFy = Utils.ndgrid(self.getTensor('Fy')) return self._gridFy return locals() _gridFy = None # Store grid by default @@ -214,7 +214,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFz is None and self.dim > 2: - self._gridFz = utils.ndgrid(self.getTensor('Fz')) + self._gridFz = Utils.ndgrid(self.getTensor('Fz')) return self._gridFz return locals() _gridFz = None # Store grid by default @@ -225,7 +225,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEx is None: - self._gridEx = utils.ndgrid(self.getTensor('Ex')) + self._gridEx = Utils.ndgrid(self.getTensor('Ex')) return self._gridEx return locals() _gridEx = None # Store grid by default @@ -236,7 +236,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEy is None and self.dim > 1: - self._gridEy = utils.ndgrid(self.getTensor('Ey')) + self._gridEy = Utils.ndgrid(self.getTensor('Ey')) return self._gridEy return locals() _gridEy = None # Store grid by default @@ -247,7 +247,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEz is None and self.dim > 2: - self._gridEz = utils.ndgrid(self.getTensor('Ez')) + self._gridEz = Utils.ndgrid(self.getTensor('Ez')) return self._gridEz return locals() _gridEz = None # Store grid by default @@ -262,13 +262,13 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): vh = self.h # Compute cell volumes if(self.dim == 1): - self._vol = utils.mkvc(vh[0]) + self._vol = Utils.mkvc(vh[0]) elif(self.dim == 2): # Cell sizes in each direction - self._vol = utils.mkvc(np.outer(vh[0], vh[1])) + self._vol = Utils.mkvc(np.outer(vh[0], vh[1])) elif(self.dim == 3): # Cell sizes in each direction - self._vol = utils.mkvc(np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) + self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) return self._vol return locals() _vol = None @@ -289,12 +289,12 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): elif(self.dim == 2): area1 = np.outer(np.ones(n[0]+1), vh[1]) area2 = np.outer(vh[0], np.ones(n[1]+1)) - self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)] + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] elif(self.dim == 3): - area1 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], vh[2]))) - area2 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - area3 = np.outer(vh[0], utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)] + area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] return self._area return locals() _area = None @@ -311,16 +311,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): n = self.n # Compute edge lengths if(self.dim == 1): - self._edge = utils.mkvc(vh[0]) + self._edge = Utils.mkvc(vh[0]) elif(self.dim == 2): l1 = np.outer(vh[0], np.ones(n[1]+1)) l2 = np.outer(np.ones(n[0]+1), vh[1]) - self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2)] + self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)] elif(self.dim == 3): - l1 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) - l2 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - l3 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2), utils.mkvc(l3)] + l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] return self._edge return locals() _edge = None @@ -410,11 +410,11 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: nF_nE = self.nFv if 'F' in locType else self.nEv - components = [utils.spzeros(loc.shape[0], n) for n in nF_nE] - components[ind] = utils.interpmat(loc, *self.getTensor(locType)) + components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) Q = sp.hstack(components) elif locType in ['CC', 'N']: - Q = utils.interpmat(loc, *self.getTensor(locType)) + Q = Utils.interpmat(loc, *self.getTensor(locType)) else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) return Q diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index a745c1cd..cbecde0a 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -2,7 +2,7 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib from mpl_toolkits.mplot3d import Axes3D -from SimPEG.utils import mkvc, animate +from SimPEG.Utils import mkvc, animate class TensorView(object): diff --git a/SimPEG/setup.py b/SimPEG/setup.py index c421cb4b..5afc57a4 100644 --- a/SimPEG/setup.py +++ b/SimPEG/setup.py @@ -1,8 +1,8 @@ import os print 'Compiling TriSolve.' -os.system('f2py -c utils/TriSolve.f -m TriSolve') +os.system('f2py -c Utils/TriSolve.f -m TriSolve') print 'TriSolve Compiled! yay.' print 'Moving TriSolve into Utils.' -os.system('mv TriSolve.so utils/TriSolve.so') +os.system('mv TriSolve.so Utils/TriSolve.so') print 'Thats it. Well Done Computer.' diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index fa4ee89a..c5490f90 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -1,9 +1,9 @@ import numpy as np import matplotlib.pyplot as plt from numpy.linalg import norm -from SimPEG.utils import mkvc, sdiag -from SimPEG import utils -from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh +from SimPEG.Utils import mkvc, sdiag +from SimPEG import Utils +from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh import numpy as np import scipy.sparse as sp import unittest @@ -112,10 +112,10 @@ class OrderTest(unittest.TestCase): else: raise Exception('Unexpected meshType') if self.meshDimension == 2: - X, Y = utils.exampleLomGird([nc, nc], kwrd) + X, Y = Utils.exampleLomGird([nc, nc], kwrd) self.M = LogicallyOrthogonalMesh([X, Y]) if self.meshDimension == 3: - X, Y, Z = utils.exampleLomGird([nc, nc, nc], kwrd) + X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd) self.M = LogicallyOrthogonalMesh([X, Y, Z]) return 1./nc @@ -212,7 +212,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole :include-source: from SimPEG.tests import checkDerivative - from SimPEG.utils import sdiag + from SimPEG.Utils import sdiag import numpy as np def simplePass(x): return np.sin(x), sdiag(np.cos(x)) diff --git a/SimPEG/tests/test_LogicallyOrthogonalMesh.py b/SimPEG/tests/test_LogicallyOrthogonalMesh.py index d760e230..241b2da8 100644 --- a/SimPEG/tests/test_LogicallyOrthogonalMesh.py +++ b/SimPEG/tests/test_LogicallyOrthogonalMesh.py @@ -1,7 +1,7 @@ import numpy as np import unittest -from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh -from SimPEG.utils import ndgrid +from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh +from SimPEG.Utils import ndgrid class BasicLOMTests(unittest.TestCase): diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py index f8c2964d..f8a8983f 100644 --- a/SimPEG/tests/test_Solver.py +++ b/SimPEG/tests/test_Solver.py @@ -1,7 +1,7 @@ import unittest from SimPEG import Solver -from SimPEG.mesh import TensorMesh -from SimPEG.utils import sdiag +from SimPEG.Mesh import TensorMesh +from SimPEG.Utils import sdiag import numpy as np import scipy.sparse as sparse diff --git a/SimPEG/tests/test_basemesh.py b/SimPEG/tests/test_basemesh.py index d81a8553..835bd29e 100644 --- a/SimPEG/tests/test_basemesh.py +++ b/SimPEG/tests/test_basemesh.py @@ -1,6 +1,6 @@ import unittest import sys -from SimPEG.mesh import BaseMesh +from SimPEG.Mesh import BaseMesh import numpy as np diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index 90312fdf..1a44b350 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -1,85 +1,85 @@ -import numpy as np -import unittest -from SimPEG.mesh import TensorMesh -from SimPEG.utils import ModelBuilder, sdiag -from SimPEG.forward import Problem -from SimPEG.examples.DC import * -from TestUtils import checkDerivative -from scipy.sparse.linalg import dsolve -from SimPEG import inverse +# import numpy as np +# import unittest +# from SimPEG.mesh import TensorMesh +# from SimPEG.Utils import ModelBuilder, sdiag +# from SimPEG.forward import Problem +# from SimPEG.examples.DC import * +# from TestUtils import checkDerivative +# from scipy.sparse.linalg import dsolve +# from SimPEG import inverse -class DCProblemTests(unittest.TestCase): +# class DCProblemTests(unittest.TestCase): - def setUp(self): - # Create the mesh - h1 = np.ones(20) - h2 = np.ones(20) - mesh = TensorMesh([h1,h2]) +# def setUp(self): +# # Create the mesh +# h1 = np.ones(20) +# h2 = np.ones(20) +# mesh = TensorMesh([h1,h2]) - # Create some parameters for the model - sig1 = 1 - sig2 = 0.01 +# # Create some parameters for the model +# sig1 = 1 +# sig2 = 0.01 - # Create a synthetic model from a block in a half-space - p0 = [2, 2] - p1 = [5, 5] - condVals = [sig1, sig2] - mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) +# # Create a synthetic model from a block in a half-space +# p0 = [2, 2] +# p1 = [5, 5] +# condVals = [sig1, sig2] +# mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) - # Set up the projection - nelec = 10 - spacelec = 2 - surfloc = 0.5 - elecini = 0.5 - elecend = 0.5+spacelec*(nelec-1) - elecLocR = np.linspace(elecini, elecend, nelec) - rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) - P = Q.T +# # Set up the projection +# nelec = 10 +# spacelec = 2 +# surfloc = 0.5 +# elecini = 0.5 +# elecend = 0.5+spacelec*(nelec-1) +# elecLocR = np.linspace(elecini, elecend, nelec) +# rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 +# q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) +# P = Q.T - # Create some data +# # Create some data - problem = DCProblem(mesh) - problem.P = P - problem.RHS = q - data = problem.createSyntheticData(mSynth, std=0.05) +# problem = DCProblem(mesh) +# problem.P = P +# problem.RHS = q +# data = problem.createSyntheticData(mSynth, std=0.05) - # Now set up the problem to do some minimization - opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - reg = inverse.Regularization(mesh) - inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4) +# # Now set up the problem to do some minimization +# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) +# reg = inverse.Regularization(mesh) +# inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4) - self.inv = inv - self.reg = reg - self.p = problem - self.mesh = mesh - self.m0 = mSynth - self.data = data +# self.inv = inv +# self.reg = reg +# self.p = problem +# self.mesh = mesh +# self.m0 = mSynth +# self.data = data - def test_misfit(self): - derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)] - passed = checkDerivative(derChk, self.m0, plotIt=False) - self.assertTrue(passed) +# def test_misfit(self): +# derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)] +# passed = checkDerivative(derChk, self.m0, plotIt=False) +# self.assertTrue(passed) - def test_adjoint(self): - # Adjoint Test - u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1]) - v = np.random.rand(self.mesh.nC) - w = np.random.rand(self.data.dobs.shape[0]) - wtJv = w.dot(self.p.J(self.m0, v, u=u)) - vtJtw = v.dot(self.p.Jt(self.m0, w, u=u)) - passed = (wtJv - vtJtw) < 1e-10 - self.assertTrue(passed) +# def test_adjoint(self): +# # Adjoint Test +# u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1]) +# v = np.random.rand(self.mesh.nC) +# w = np.random.rand(self.data.dobs.shape[0]) +# wtJv = w.dot(self.p.J(self.m0, v, u=u)) +# vtJtw = v.dot(self.p.Jt(self.m0, w, u=u)) +# passed = (wtJv - vtJtw) < 1e-10 +# self.assertTrue(passed) - def test_dataObj(self): - derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)] - checkDerivative(derChk, self.m0, plotIt=False) +# def test_dataObj(self): +# derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)] +# checkDerivative(derChk, self.m0, plotIt=False) - def test_modelObj(self): - derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] - checkDerivative(derChk, self.m0, plotIt=False) +# def test_modelObj(self): +# derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] +# checkDerivative(derChk, self.m0, plotIt=False) -if __name__ == '__main__': - unittest.main() +# if __name__ == '__main__': +# unittest.main() diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py index 3d751dde..d01dcbc3 100644 --- a/SimPEG/tests/test_interpolation.py +++ b/SimPEG/tests/test_interpolation.py @@ -1,7 +1,7 @@ import numpy as np import unittest from TestUtils import OrderTest -from SimPEG.utils import mkvc +from SimPEG.Utils import mkvc MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh'] TOLERANCES = [0.9, 0.55] diff --git a/SimPEG/tests/test_model.py b/SimPEG/tests/test_model.py new file mode 100644 index 00000000..946224f1 --- /dev/null +++ b/SimPEG/tests/test_model.py @@ -0,0 +1,27 @@ +import numpy as np +import unittest +from SimPEG import * +from TestUtils import checkDerivative +from scipy.sparse.linalg import dsolve + + +class ModelTests(unittest.TestCase): + + def setUp(self): + + a = np.array([1, 1, 1]) + b = np.array([1, 2]) + c = np.array([1, 4]) + self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5])) + + def test_modelTransforms(self): + print 'SimPEG.Model.BaseModel: Testing Model Transform' + for M in dir(Model): + if 'Model' not in M: continue + model = getattr(Model, M)() + m = model.example(self.mesh2) + passed = checkDerivative(lambda m : [model.transform(m), model.transformDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/tests/test_optimizers.py index 0253852c..12c8c9b5 100644 --- a/SimPEG/tests/test_optimizers.py +++ b/SimPEG/tests/test_optimizers.py @@ -1,11 +1,11 @@ import unittest from SimPEG import Solver -from SimPEG.mesh import TensorMesh -from SimPEG.utils import sdiag +from SimPEG.Mesh import TensorMesh +from SimPEG.Utils import sdiag import numpy as np import scipy.sparse as sp -from SimPEG import inverse -from SimPEG.tests import getQuadratic, Rosenbrock +from SimPEG import Inverse +from SimPEG.Tests import getQuadratic, Rosenbrock TOL = 1e-2 @@ -16,7 +16,7 @@ class TestOptimizers(unittest.TestCase): self.b = np.array([-5,-5]) def test_GN_Rosenbrock(self): - GN = inverse.GaussNewton() + GN = Inverse.GaussNewton() xopt = GN.minimize(Rosenbrock,np.array([0,0])) x_true = np.array([1.,1.]) print 'xopt: ', xopt @@ -24,7 +24,7 @@ class TestOptimizers(unittest.TestCase): self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) def test_GN_quadratic(self): - GN = inverse.GaussNewton() + GN = Inverse.GaussNewton() xopt = GN.minimize(getQuadratic(self.A,self.b),np.array([0,0])) x_true = np.array([5.,5.]) print 'xopt: ', xopt @@ -32,7 +32,7 @@ class TestOptimizers(unittest.TestCase): self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) def test_ProjGradient_quadraticBounded(self): - PG = inverse.ProjectedGradient(debug=True) + PG = Inverse.ProjectedGradient(debug=True) PG.lower, PG.upper = -2, 2 xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0])) x_true = np.array([2.,2.]) @@ -42,7 +42,7 @@ class TestOptimizers(unittest.TestCase): def test_ProjGradient_quadratic1Bound(self): myB = np.array([-5,1]) - PG = inverse.ProjectedGradient() + PG = Inverse.ProjectedGradient() PG.lower, PG.upper = -2, 2 xopt = PG.minimize(getQuadratic(self.A,myB),np.array([0,0])) x_true = np.array([2.,-1.]) @@ -53,7 +53,7 @@ class TestOptimizers(unittest.TestCase): def test_NewtonRoot(self): fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), sdiag( np.cos(x) ) ) x = np.array([np.pi-0.3, np.pi+0.1, 0]) - xopt = inverse.NewtonRoot(comments=False).root(fun,x) + xopt = Inverse.NewtonRoot(comments=False).root(fun,x) x_true = np.array([np.pi,np.pi,0]) print 'Newton Root Finding' print 'xopt: ', xopt diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_problem.py similarity index 51% rename from SimPEG/tests/test_forward_problem.py rename to SimPEG/tests/test_problem.py index 1d4de045..fa94e6b8 100644 --- a/SimPEG/tests/test_forward_problem.py +++ b/SimPEG/tests/test_problem.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG import mesh, forward, inverse +from SimPEG import * from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve @@ -12,15 +12,9 @@ class ProblemTests(unittest.TestCase): a = np.array([1, 1, 1]) b = np.array([1, 2]) c = np.array([1, 4]) - self.mesh2 = mesh.TensorMesh([a, b], np.array([3, 5])) - self.p2 = forward.Problem(self.mesh2) - self.reg = inverse.Regularization(self.mesh2) - - def test_modelTransform(self): - print 'SimPEG.forward.Problem: Testing Model Transform' - m = np.random.rand(self.mesh2.nC) - passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False) - self.assertTrue(passed) + self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5])) + self.p2 = Problem.BaseProblem(self.mesh2, None) + self.reg = Inverse.Regularization(self.mesh2) def test_regularization(self): derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] @@ -28,7 +22,5 @@ class ProblemTests(unittest.TestCase): checkDerivative(derChk, mSynth, plotIt=False) - - if __name__ == '__main__': unittest.main() diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index 9544dab6..3e01181b 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG.mesh import TensorMesh +from SimPEG.Mesh import TensorMesh from TestUtils import OrderTest from scipy.sparse.linalg import dsolve diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py index 058dba56..fea231f2 100644 --- a/SimPEG/tests/test_utils.py +++ b/SimPEG/tests/test_utils.py @@ -1,7 +1,7 @@ import numpy as np import unittest -from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal -from SimPEG.tests import checkDerivative +from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.Tests import checkDerivative class TestCheckDerivative(unittest.TestCase): diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py index 969551ed..be172b6b 100644 --- a/SimPEG/utils/ModelBuilder.py +++ b/SimPEG/utils/ModelBuilder.py @@ -151,7 +151,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]): .. plot:: import matplotlib.pyplot as plt - import SimPEG.utils.ModelBuilder as MB + import SimPEG.Utils.ModelBuilder as MB plt.colorbar(plt.imshow(MB.randomModel((50,50),bounds=[-4,0]))) plt.title('A very cool, yet completely random model.') plt.show() diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 25fd9b9e..55de4b4b 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -5,7 +5,7 @@ import re try: import h5py except Exception, e: - print 'Warning: SimPEG.utils.Save needs h5py to be installed.' + print 'Warning: SimPEG.Utils.Save needs h5py to be installed.' SAVEABLES = {} @@ -347,6 +347,6 @@ def loadSavable(node, pointers=None): print 'KWARGS: ', KWARGS return (cls, ARGS, KWARGS, node) else: - print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls + print 'Warning: %s Class not found in SimPEG.Utils.Save.SAVABLES' % cls return (cls, ARGS, KWARGS, node) diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 2def4332..fc1c1d0f 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,7 +3,7 @@ from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap from meshutils import exampleLomGird, meshTensors from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube from interputils import interpmat -from ipythonUtils import easyAnimate as animate +from ipythonutils import easyAnimate as animate from Solver import Solver import Save import Geophysics diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py index c8bcb4ec..43cbff6b 100644 --- a/SimPEG/utils/interputils.py +++ b/SimPEG/utils/interputils.py @@ -45,7 +45,7 @@ def interpmat(locs, x, y=None, z=None): x = np.linspace(0,1,7) dense = np.linspace(0,1,200) fun = lambda x: np.cos(2*np.pi*x) - Q = SimPEG.utils.interpmat(locs, x) + Q = SimPEG.Utils.interpmat(locs, x) plt.plot(x, fun(x), 'bs-') plt.plot(dense, fun(dense), 'y:') plt.plot(locs, Q*fun(x), 'mo') @@ -173,7 +173,7 @@ if __name__ == '__main__': x = np.linspace(0,1,7) dense = np.linspace(0,1,200) fun = lambda x: np.cos(2*np.pi*x) - Q = SimPEG.utils.interpmat(locs, x) + Q = SimPEG.Utils.interpmat(locs, x) plt.plot(x, fun(x), 'bs-') plt.plot(dense, fun(dense), 'y:') plt.plot(locs, Q*fun(x), 'mo') diff --git a/SimPEG/utils/meshutils.py b/SimPEG/utils/meshutils.py index 429db04b..0e016fb5 100644 --- a/SimPEG/utils/meshutils.py +++ b/SimPEG/utils/meshutils.py @@ -34,8 +34,8 @@ def meshTensors(*args): .. plot:: - from SimPEG import mesh, utils - M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + from SimPEG import mesh, Utils + M = mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) M.plotGrid() """ diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index c3d9f4a7..68326662 100644 --- a/SimPEG/visualize/vtk/vtkTools.py +++ b/SimPEG/visualize/vtk/vtkTools.py @@ -3,7 +3,7 @@ try: import vtk, vtk.util.numpy_support as npsup, pdb except Exception, e: print 'VTK import error. Please ensure you have VTK installed to use this visualization package.' -from SimPEG.utils import mkvc +from SimPEG.Utils import mkvc class vtkTools(object):