wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-07 14:41:06 +08:00
Code Issues Packages Projects Releases Wiki Activity

Moved things around! Packages should now all be capitalized. may need to to tweak git to ensure this...

Browse Source
This commit is contained in:
rowanc1
2014-01-16 13:19:29 -08:00
parent c49ae77fbd
commit ca4932fd19
37 changed files with 546 additions and 489 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/forward/SimPEGData.py → SimPEG/Data.py
+79 -42
View File
@@ -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()
SimPEG/Model.py
+94
View File
@@ -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)))
SimPEG/forward/Problem.py → SimPEG/Problem.py
+14 -80
View File
@@ -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)
SimPEG/__init__.py
+9 -8
View File
@@ -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':
SimPEG/examples/DC.py
+52 -30
View File
@@ -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()
SimPEG/examples/Linear.py
+6 -4
View File
@@ -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)
SimPEG/forward/ModelTransforms.py
-49
View File
@@ -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)))
SimPEG/forward/__init__.py
-3
View File
@@ -1,3 +0,0 @@
from Problem import *
import ModelTransforms
from SimPEGData import *
SimPEG/inverse/Inversion.py
+17 -17
View File
@@ -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)
SimPEG/inverse/Optimize.py
+33 -33
View File
@@ -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
SimPEG/inverse/Regularization.py
+20 -20
View File
@@ -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):
"""
SimPEG/mesh/BaseMesh.py
+4 -4
View File
@@ -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]
SimPEG/mesh/Cyl1DMesh.py
+2 -2
View File
@@ -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"
SimPEG/mesh/DiffOperators.py
+1 -1
View File
@@ -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):
SimPEG/mesh/InnerProducts.py
+1 -1
View File
@@ -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
SimPEG/mesh/LogicallyOrthogonalMesh.py
+47 -47
View File
@@ -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')
SimPEG/mesh/LomView.py
+1 -1
View File
@@ -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):
SimPEG/mesh/TensorMesh.py
+30 -30
View File
@@ -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
SimPEG/mesh/TensorView.py
+1 -1
View File
@@ -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):
SimPEG/setup.py
+2 -2
View File
@@ -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.'
SimPEG/tests/TestUtils.py
+6 -6
View File
@@ -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))
SimPEG/tests/test_LogicallyOrthogonalMesh.py
+2 -2
View File
@@ -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):
SimPEG/tests/test_Solver.py
+2 -2
View File
@@ -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
SimPEG/tests/test_basemesh.py
+1 -1
View File
@@ -1,6 +1,6 @@
import unittest
import sys
from SimPEG.mesh import BaseMesh
from SimPEG.Mesh import BaseMesh
import numpy as np
SimPEG/tests/test_forward_DCproblem.py
+69 -69
View File
@@ -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()
SimPEG/tests/test_interpolation.py
+1 -1
View File
@@ -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]
SimPEG/tests/test_model.py
+27
View File
@@ -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()
SimPEG/tests/test_optimizers.py
+9 -9
View File
@@ -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
SimPEG/tests/test_forward_problem.py → SimPEG/tests/test_problem.py
+4 -12
View File
@@ -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()
SimPEG/tests/test_tensorMesh.py
+1 -1
View File
@@ -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
SimPEG/tests/test_utils.py
+2 -2
View File
@@ -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):
SimPEG/utils/ModelBuilder.py
+1 -1
View File
@@ -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()
SimPEG/utils/Save.py
+2 -2
View File
@@ -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)
SimPEG/utils/__init__.py
+1 -1
View File
@@ -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
SimPEG/utils/interputils.py
+2 -2
View File
@@ -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')
SimPEG/utils/meshutils.py
+2 -2
View File
@@ -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()
"""
SimPEG/visualize/vtk/vtkTools.py
+1 -1
View File
@@ -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):