Skip to content

Commit

Permalink
f for fields in data misfit, directives etc. Previously, f was used i…
Browse files Browse the repository at this point in the history
…n the InvProblem to be the function value for the objective function --> this has been renamed to phi
  • Loading branch information
lheagy committed Mar 31, 2016
1 parent 0a0cace commit d8d8915
Show file tree
Hide file tree
Showing 9 changed files with 61 additions and 61 deletions.
4 changes: 2 additions & 2 deletions SimPEG/DataMisfit.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,5 +123,5 @@ def evalDeriv(self, m, f=None):
@Utils.timeIt
def eval2Deriv(self, m, v, f=None):
"eval2Deriv(m, v, f=None)"
if f is None: f = prob.fields(m)
return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, f=f)), f=f)
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)
4 changes: 2 additions & 2 deletions SimPEG/Directives.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,10 +123,10 @@ def initialize(self):
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)

Expand Down
16 changes: 8 additions & 8 deletions SimPEG/EM/FDEM/SurveyFDEM.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ def projGLoc(self, u):
"""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.
Expand All @@ -80,27 +80,27 @@ def eval(self, src, mesh, u):
# 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
Expand Down
22 changes: 11 additions & 11 deletions SimPEG/EM/TDEM/BaseTDEM.py
Original file line number Diff line number Diff line change
Expand Up @@ -108,11 +108,11 @@ def adjoint(self, m, RHS, F=None):
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)
Expand All @@ -125,15 +125,15 @@ def Jvec(self, m, v, u=None):
"""
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)
Expand All @@ -150,15 +150,15 @@ def Jtvec(self, m, v, u=None):
"""
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)

26 changes: 13 additions & 13 deletions SimPEG/InvProblem.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,23 +82,23 @@ def warmstart(self, value):
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):
Expand All @@ -109,29 +109,29 @@ def evalFunction(self, m, return_g=True, return_H=True):
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
out += (g,)

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
Expand Down
26 changes: 13 additions & 13 deletions SimPEG/MT/BaseMT.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ def __init__(self, mesh, **kwargs):
# 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.
Expand All @@ -39,8 +39,8 @@ def Jvec(self, m, v, u=None):
"""

# 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
Expand All @@ -56,9 +56,9 @@ def Jvec(self, m, v, u=None):
# 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 )
Expand All @@ -68,13 +68,13 @@ def Jvec(self, m, v, u=None):
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.
Expand All @@ -85,8 +85,8 @@ def Jtvec(self, m, v, u=None):
:return: Data sensitivities wrt m
"""

if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)

self.curModel = m

Expand All @@ -103,15 +103,15 @@ def Jtvec(self, m, v, u=None):

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:
Expand All @@ -129,4 +129,4 @@ def Jtvec(self, m, v, u=None):
raise Exception('Must be real or imag')
# Clean the factorization, clear memory.
ATinv.clean()
return Jtv
return Jtv
6 changes: 3 additions & 3 deletions SimPEG/MT/SurveyMT.py
Original file line number Diff line number Diff line change
Expand Up @@ -427,15 +427,15 @@ def getSrcByFreq(self, freq):
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.')

#################
Expand Down
4 changes: 2 additions & 2 deletions SimPEG/Problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ def Jvec_approx(self, m, v, f=None):
: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, f=None):
Expand All @@ -142,7 +142,7 @@ def Jtvec_approx(self, m, v, f=None):
:rtype: numpy.array
:return: JTv
"""
return self.Jtvec(m, v, u)
return self.Jtvec(m, v, f)

def fields(self, m):
"""
Expand Down
14 changes: 7 additions & 7 deletions SimPEG/Survey.py
Original file line number Diff line number Diff line change
Expand Up @@ -313,20 +313,20 @@ def dpred(self, m, f=None):


@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.
Expand Down Expand Up @@ -379,8 +379,8 @@ def makeSyntheticData(self, m, std=0.05, f=None, force=False):
return self.dobs

class LinearSurvey(BaseSurvey):
def eval(self, u):
return u
def eval(self, f):
return f

@property
def nD(self):
Expand Down

0 comments on commit d8d8915

Please sign in to comment.