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

Merge branch 'dev' into em/dev

Browse Source
This commit is contained in:
Lindsey Heagy
2016-04-02 08:35:09 -07:00
parent a220c75d78 99ed4ad50e
commit f59cfa9481
13 changed files with 157 additions and 163 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/BaseDC.py
+11 -13
View File
@@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem):
return F
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: Jv
@@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem):
# Set current model; clear dependent property $\mathbf{A(m)}$
self.curModel = m
sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
if u is None:
if f is None:
# Run forward simulation if $u$ not provided
u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
else:
u = u[self.survey.srcList, 'phi_sol']
f = self.fields(self.curModel)
u = f[self.survey.srcList, 'phi_sol']
D = self.mesh.faceDiv
G = self.mesh.cellGrad
@@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem):
if self.Ainv is None:
self.Ainv = self.Solver(dA_du, **self.solverOpts)
P = self.survey.getP(self.mesh)
P = self.survey.getP(self.mesh)
Jv = - P * mkvc( self.Ainv * dCdm_x_v )
return Jv
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
self.curModel = m
sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
if u is None:
# Run forward simulation if $u$ not provided
u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
else:
u = u[self.survey.srcList, 'phi_sol']
if f is None:
# Run forward simulation if $f$ not provided
f = self.fields(self.curModel)
u = f[self.survey.srcList, 'phi_sol']
shp = u.shape
P = self.survey.getP(self.mesh)
SimPEG/DCIP/BaseIP.py
+4 -4
View File
@@ -14,12 +14,12 @@ class SurveyIP(SurveyDC):
Survey.BaseSurvey.__init__(self, **kwargs)
self._Ps = {}
def dpred(self, m, u=None):
def dpred(self, m, f=None):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
d_\\text{pred} = Pf(m)
"""
return self.prob.forward(m)
@@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem):
J_x_v = - P * mkvc( self.Ainv * dCdm_x_v )
return -J_x_v
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
return self.forward(v)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
self.curModel = m
# sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
SimPEG/DataMisfit.py
+21 -25
View File
@@ -22,11 +22,11 @@ class BaseDataMisfit(object):
Utils.setKwargs(self,**kwargs)
@Utils.timeIt
def eval(self, m, u=None):
"""eval(m, u=None)
def eval(self, m, f=None):
"""eval(m, f=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:param Fields f: fields
:rtype: float
:return: data misfit
@@ -34,11 +34,11 @@ class BaseDataMisfit(object):
raise NotImplementedError('This method should be overwritten.')
@Utils.timeIt
def evalDeriv(self, m, u=None):
"""evalDeriv(m, u=None)
def evalDeriv(self, m, f=None):
"""evalDeriv(m, f=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: data misfit derivative
@@ -47,12 +47,12 @@ class BaseDataMisfit(object):
@Utils.timeIt
def eval2Deriv(self, m, v, u=None):
"""eval2Deriv(m, v, u=None)
def eval2Deriv(self, m, v, f=None):
"""eval2Deriv(m, v, f=None)
:param numpy.array m: geophysical model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: data misfit derivative
@@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit):
self._Wd = value
@Utils.timeIt
def eval(self, m, u=None):
"eval(m, u=None)"
prob = self.prob
survey = self.survey
R = self.Wd * survey.residual(m, u)
def eval(self, m, f=None):
"eval(m, f=None)"
if f is None: f = self.prob.fields(m)
R = self.Wd * self.survey.residual(m, f)
return 0.5*np.vdot(R, R)
@Utils.timeIt
def evalDeriv(self, m, u=None):
"evalDeriv(m, u=None)"
prob = self.prob
survey = self.survey
if u is None: u = prob.fields(m)
return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u)), u)
def evalDeriv(self, m, f=None):
"evalDeriv(m, f=None)"
if f is None: f = self.prob.fields(m)
return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f)
@Utils.timeIt
def eval2Deriv(self, m, v, u=None):
"eval2Deriv(m, v, u=None)"
prob = self.prob
if u is None: u = prob.fields(m)
return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u)
def eval2Deriv(self, m, v, f=None):
"eval2Deriv(m, v, f=None)"
if f is None: f = self.prob.fields(m)
return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f)
SimPEG/Directives.py
+2 -2
View File
@@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective):
if self.debug: print 'Calculating the beta0 parameter.'
m = self.invProb.curModel
u = self.invProb.getFields(m, store=True, deleteWarmstart=False)
f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
x0 = np.random.rand(*m.shape)
t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u))
t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f))
b = x0.dot(self.reg.eval2Deriv(m, v=x0))
self.beta0 = self.beta0_ratio*(t/b)
SimPEG/EM/FDEM/SurveyFDEM.py
+8 -8
View File
@@ -67,7 +67,7 @@ class Rx(SimPEG.Survey.BaseRx):
"""Grid Location projection (e.g. Ex Fy ...)"""
return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
def eval(self, src, mesh, u):
def eval(self, src, mesh, f):
"""
Project fields to recievers to get data.
@@ -80,27 +80,27 @@ class Rx(SimPEG.Survey.BaseRx):
# projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
# projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh, self.projGLoc(u))
u_part_complex = u[src, self.projField]
P = self.getP(mesh, self.projGLoc(f))
f_part_complex = f[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
u_part = getattr(u_part_complex, real_or_imag)
f_part = getattr(f_part_complex, real_or_imag)
return P*u_part
return P*f_part
def evalDeriv(self, src, mesh, u, v, adjoint=False):
def evalDeriv(self, src, mesh, f, v, adjoint=False):
"""
Derivative of projected fields with respect to the inversion model times a vector.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields u: fields object
:param Fields f: fields object
:param numpy.ndarray v: vector to multiply
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh, self.projGLoc(u))
P = self.getP(mesh, self.projGLoc(f))
if not adjoint:
Pv_complex = P * v
SimPEG/EM/TDEM/BaseTDEM.py
+11 -11
View File
@@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
Ainv.clean()
return F
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray v: vector (model object)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
:rtype: numpy.ndarray
:return: w (data object)
@@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if u is None:
u = self.fields(m)
p = self.Gvec(m, v, u)
if f is None:
f = self.fields(m)
p = self.Gvec(m, v, f)
y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(u, v=y)
Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
return - mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
@@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
p = self.survey.evalDeriv(u, v=v, adjoint=True)
p = self.survey.evalDeriv(f, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, u)
w = self.Gtvec(m, y, f)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
return - mkvc(w)
SimPEG/FLOW/Richards/RichardsProblem.py
+29 -29
View File
@@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.count
@Utils.requires('prob')
def dpred(self, m, u=None):
def dpred(self, m, f=None):
"""
Create the projected data from a model.
The field, u, (if provided) will be used for the predicted data
The field, f, (if provided) will be used for the predicted data
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m), m)
d_\\text{pred} = P(f(m), m)
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.eval(u, m))
if f is None: f = self.prob.fields(m)
return Utils.mkvc(self.eval(f, m))
@Utils.requires('prob')
def eval(self, U, m):
@@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem):
return r, J
@Utils.timeIt
def Jfull(self, m, u=None):
if u is None:
u = self.fields(m)
def Jfull(self, m, f=None):
if f is None:
f = self.fields(m)
nn = len(u)-1
nn = len(f)-1
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
for ii in range(nn):
dt = self.timeSteps[ii]
bc = self.getBoundaryConditions(ii, u[ii])
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
bc = self.getBoundaryConditions(ii, f[ii])
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[ii+1], dt, bc)
Ad = sp.block_diag(Adiags)
zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1])
zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1])
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
B = np.array(sp.vstack(Bs).todense())
Ainv = self.Solver(A, **self.solverOpts)
P = self.survey.evalDeriv(u, m)
P = self.survey.evalDeriv(f, m)
AinvB = Ainv * B
z = np.zeros((self.mesh.nC, B.shape[1]))
zAinvB = np.vstack((z, AinvB))
@@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem):
return J
@Utils.timeIt
def Jvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
def Jvec(self, m, v, f=None):
if f is None:
f = self.fields(m)
JvC = range(len(u)-1) # Cell to hold each row of the long vector.
JvC = range(len(f)-1) # Cell to hold each row of the long vector.
# This is done via forward substitution.
bc = self.getBoundaryConditions(0, u[0])
temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
bc = self.getBoundaryConditions(0, f[0])
temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc)
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[0] = Adiaginv * (B*v)
for ii in range(1,len(u)-1):
bc = self.getBoundaryConditions(ii, u[ii])
Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
for ii in range(1,len(f)-1):
bc = self.getBoundaryConditions(ii, f[ii])
Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc)
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
P = self.survey.evalDeriv(u, m)
P = self.survey.evalDeriv(f, m)
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
@Utils.timeIt
def Jtvec(self, m, v, u=None):
if u is None:
u = self.field(m)
def Jtvec(self, m, v, f=None):
if f is None:
f = self.field(m)
P = self.survey.evalDeriv(u, m)
P = self.survey.evalDeriv(f, m)
PTv = P.T*v
# This is done via backward substitution.
minus = 0
BJtv = 0
for ii in range(len(u)-1,0,-1):
bc = self.getBoundaryConditions(ii-1, u[ii-1])
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
for ii in range(len(f)-1,0,-1):
bc = self.getBoundaryConditions(ii-1, f[ii-1])
Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
#select the correct part of v
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
SimPEG/InvProblem.py
+13 -13
View File
@@ -82,23 +82,23 @@ class BaseInvProblem(object):
self._warmstart = value
def getFields(self, m, store=False, deleteWarmstart=True):
u = None
f = None
for mtest, u_ofmtest in self.warmstart:
if m is mtest:
u = u_ofmtest
f = u_ofmtest
if self.debug: print 'InvProb is Warm Starting!'
break
if u is None:
u = self.prob.fields(m)
if f is None:
f = self.prob.fields(m)
if deleteWarmstart:
self.warmstart = []
if store:
self.warmstart += [(m,u)]
self.warmstart += [(m,f)]
return u
return f
@Utils.timeIt
def evalFunction(self, m, return_g=True, return_H=True):
@@ -109,21 +109,21 @@ class BaseInvProblem(object):
gc.collect()
# Store fields if doing a line-search
u = self.getFields(m, store=(return_g==False and return_H==False))
f = self.getFields(m, store=(return_g==False and return_H==False))
phi_d = self.dmisfit.eval(m, u=u)
phi_d = self.dmisfit.eval(m, f=f)
phi_m = self.reg.eval(m)
self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.dpred = self.survey.dpred(m, f=f) # This is a cheap matrix vector calculation.
self.phi_d, self.phi_d_last = phi_d, self.phi_d
self.phi_m, self.phi_m_last = phi_m, self.phi_m
f = phi_d + self.beta * phi_m
phi = phi_d + self.beta * phi_m
out = (f,)
out = (phi,)
if return_g:
phi_dDeriv = self.dmisfit.evalDeriv(m, u=u)
phi_dDeriv = self.dmisfit.evalDeriv(m, f=f)
phi_mDeriv = self.reg.evalDeriv(m)
g = phi_dDeriv + self.beta * phi_mDeriv
@@ -131,7 +131,7 @@ class BaseInvProblem(object):
if return_H:
def H_fun(v):
phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u)
phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f)
phi_m2Deriv = self.reg.eval2Deriv(m, v=v)
return phi_d2Deriv + self.beta * phi_m2Deriv
SimPEG/MT/BaseMT.py
+13 -13
View File
@@ -27,7 +27,7 @@ class BaseMTProblem(BaseFDEMProblem):
# Might need to add more stuff here.
## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
@@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem):
"""
# Calculate the fields
if u is None:
u = self.fields(m)
if f is None:
f= self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
@@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem):
# We need fDeriv_m = df/du*du/dm + df/dm
# Construct du/dm, it requires a solve
# NOTE: need to account for the 2 polarizations in the derivatives.
u_src = u[src,:]
f_src = f[src,:]
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
if dRHS_dm is None:
du_dm = dA_duI * ( -dA_dm )
@@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem):
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
return mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
@@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem):
:return: Data sensitivities wrt m
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
@@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, :]
f_src = f[src, :]
for rx in src.rxList:
# Get the adjoint evalDeriv
# PTv needs to be nE,
PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
# Make du_dmT
if dRHS_dmT is None:
@@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem):
raise Exception('Must be real or imag')
# Clean the factorization, clear memory.
ATinv.clean()
return Jtv
return Jtv
SimPEG/MT/SurveyMT.py
+3 -3
View File
@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def eval(self, u):
def eval(self, f):
data = Data(self)
for src in self.srcList:
sys.stdout.flush()
for rx in src.rxList:
data[src, rx] = rx.eval(src, self.mesh, u)
data[src, rx] = rx.eval(src, self.mesh, f)
return data
def evalDeriv(self, u):
def evalDeriv(self, f):
raise Exception('Use Transmitters to project fields deriv.')
#################
SimPEG/Problem.py
+16 -16
View File
@@ -88,28 +88,28 @@ class BaseProblem(object):
return self.survey is not None
@Utils.timeIt
def Jvec(self, m, v, u=None):
"""Jvec(m, v, u=None)
def Jvec(self, m, v, f=None):
"""Jvec(m, v, f=None)
Effect of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: Jv
"""
raise NotImplementedError('J is not yet implemented.')
@Utils.timeIt
def Jtvec(self, m, v, u=None):
"""Jtvec(m, v, u=None)
def Jtvec(self, m, v, f=None):
"""Jtvec(m, v, f=None)
Effect of transpose of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: JTv
"""
@@ -117,32 +117,32 @@ class BaseProblem(object):
@Utils.timeIt
def Jvec_approx(self, m, v, u=None):
"""Jvec_approx(m, v, u=None)
def Jvec_approx(self, m, v, f=None):
"""Jvec_approx(m, v, f=None)
Approximate effect of J(m) on a vector v
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: approxJv
"""
return self.Jvec(m, v, u)
return self.Jvec(m, v, f)
@Utils.timeIt
def Jtvec_approx(self, m, v, u=None):
"""Jtvec_approx(m, v, u=None)
def Jtvec_approx(self, m, v, f=None):
"""Jtvec_approx(m, v, f=None)
Approximate effect of transpose of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: JTv
"""
return self.Jtvec(m, v, u)
return self.Jtvec(m, v, f)
def fields(self, m):
"""
@@ -224,9 +224,9 @@ class LinearProblem(BaseProblem):
def fields(self, m):
return self.G.dot(m)
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
return self.G.dot(v)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
return self.G.T.dot(v)
SimPEG/Survey.py
+20 -20
View File
@@ -295,38 +295,38 @@ class BaseSurvey(object):
@Utils.count
@Utils.requires('prob')
def dpred(self, m, u=None):
"""dpred(m, u=None)
def dpred(self, m, f=None):
"""dpred(m, f=None)
Create the projected data from a model.
The field, u, (if provided) will be used for the predicted data
The fields, f, (if provided) will be used for the predicted data
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m))
d_\\text{pred} = P(f(m))
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.eval(u))
if f is None: f = self.prob.fields(m)
return Utils.mkvc(self.eval(f))
@Utils.count
def eval(self, u):
"""eval(u)
def eval(self, f):
"""eval(f)
This function projects the fields onto the data space.
.. math::
d_\\text{pred} = \mathbf{P} u(m)
d_\\text{pred} = \mathbf{P} f(m)
"""
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def evalDeriv(self, u):
"""evalDeriv(u)
def evalDeriv(self, f):
"""evalDeriv(f)
This function s the derivative of projects the fields onto the data space.
@@ -337,11 +337,11 @@ class BaseSurvey(object):
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def residual(self, m, u=None):
"""residual(m, u=None)
def residual(self, m, f=None):
"""residual(m, f=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:param numpy.array f: fields
:rtype: numpy.array
:return: data residual
@@ -352,14 +352,14 @@ class BaseSurvey(object):
\mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
"""
return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
return Utils.mkvc(self.dpred(m, f=f) - self.dobs)
@property
def isSynthetic(self):
"Check if the data is synthetic."
return self.mtrue is not None
def makeSyntheticData(self, m, std=0.05, u=None, force=False):
def makeSyntheticData(self, m, std=0.05, f=None, force=False):
"""
Make synthetic data given a model, and a standard deviation.
@@ -372,16 +372,16 @@ class BaseSurvey(object):
if getattr(self, 'dobs', None) is not None and not force:
raise Exception('Survey already has dobs. You can use force=True to override this exception.')
self.mtrue = m
self.dtrue = self.dpred(m, u=u)
self.dtrue = self.dpred(m, f=f)
noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
self.dobs = self.dtrue+noise
self.std = self.dobs*0 + std
return self.dobs
class LinearSurvey(BaseSurvey):
def eval(self, u):
return u
def eval(self, f):
return f
@property
def nD(self):
return self.prob.G.shape[0]
tests/flow/test_Richards.py
+6 -6
View File
@@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase):
v = np.random.rand(self.survey.nD)
z = np.random.rand(self.M.nC)
Hs = self.prob.fields(self.Ks)
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol
print 'Richards Adjoint Test - PressureHead'
@@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase):
v = np.random.rand(self.survey.nD)
z = np.random.rand(self.M.nC)
Hs = self.prob.fields(self.Ks)
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol
print '2D: Richards Adjoint Test - PressureHead'
@@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase):
v = np.random.rand(self.survey.nD)
z = np.random.rand(self.M.nC)
Hs = self.prob.fields(self.Ks)
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol
print '3D: Richards Adjoint Test - PressureHead'