From a8f09ee8698465c9ae80d21d42cb3eb32f4235b9 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Mar 2014 17:09:39 -0700 Subject: [PATCH] forward problem changes in data projection. have broken the derivative test. --- simpegEM/FDEM/FDEM.py | 18 ++--- simpegEM/FDEM/SurveyFDEM.py | 106 ++++++++++++++++++++------ simpegEM/Tests/test_FDEM.py | 8 +- simpegEM/Tests/test_FDEM_analytics.py | 2 +- simpegEM/Tests/test_FieldsObject.py | 43 +++++++++-- 5 files changed, 130 insertions(+), 47 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 3fc92b07..9c5a6643 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -110,14 +110,15 @@ class BaseProblemFDEM(Problem.BaseProblem): u_tx = u[tx, self.solType] w = self.getADeriv(freq, u_tx, v) Ainvw = solver.solve(w) - fAinvw = self.calcFields(Ainvw, freq, tx.rxList.fieldType) - P = tx.projectFieldsDeriv(self.mesh, u) + fAinvw = self.calcFieldsList(Ainvw, freq, tx.rxList.fieldTypes) + P = lambda v: tx.projectFieldsDeriv(self.mesh, u, v) - df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldType, v) + df_dm = self.calcFieldsDerivList(u_tx, freq, tx.rxList.fieldTypes, v) + #TODO: this is now a list? if df_dm is None: - Jv[tx] = - P*fAinvw + Jv[tx] = - P(fAinvw) else: - Jv[tx] = - P*fAinvw + P*df_dm + Jv[tx] = - P(fAinvw) + P(df_dm) return Utils.mkvc(Jv) @@ -140,14 +141,13 @@ class BaseProblemFDEM(Problem.BaseProblem): for tx in self.survey.getTransmitters(freq): u_tx = u[tx, self.solType] - P = tx.projectFieldsDeriv(self.mesh, u) - PTv = P.T * v[tx] - fPTv = self.calcFields(PTv, freq, tx.rxList.fieldType, adjoint=True) + PTv = tx.projectFieldsDeriv(self.mesh, u, v[tx], adjoint=True) + fPTv = self.calcFields(PTv, freq, tx.rxList.fieldTypes, adjoint=True) w = solver.solve( fPTv ) Jtv_tx = self.getADeriv(freq, u_tx, w, adjoint=True) - df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldType, PTv, adjoint=True) + df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldTypes, PTv, adjoint=True) if df_dm is None: Jtv += - Jtv_tx diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 253f2aaf..e68d50fa 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -3,29 +3,71 @@ from SimPEG import Survey, Utils, np, sp class RxListFDEM(Survey.BaseRxList): knownRxTypes = { - 'ex':'Ex', - 'ey':'Ey', - 'ez':'Ez', + 'exr':['e', 'Ex', 'real'], + 'eyr':['e', 'Ey', 'real'], + 'ezr':['e', 'Ez', 'real'], + 'exi':['e', 'Ex', 'imag'], + 'eyi':['e', 'Ey', 'imag'], + 'ezi':['e', 'Ez', 'imag'], - 'bx':'Fx', - 'by':'Fy', - 'bz':'Fz', + 'bxr':['b', 'Fx', 'real'], + 'byr':['b', 'Fy', 'real'], + 'bzr':['b', 'Fz', 'real'], + 'bxi':['b', 'Fx', 'imag'], + 'byi':['b', 'Fy', 'imag'], + 'bzi':['b', 'Fz', 'imag'], } def __init__(self, locs, rxType): Survey.BaseRxList.__init__(self, locs, rxType) self._Ps = {} + for rx in self.rxTypes: + self._Ps[self._projGLoc(rx)] = {} + + def _projField(self, rx): + """Field Type projection (e.g. e b ...)""" + if type(rx) is int: rx = self.rxTypes[rx] + return self.knownRxTypes[rx][0] + + def _projGLoc(self, rx): + """Grid Location projection (e.g. Ex Fy ...)""" + if type(rx) is int: rx = self.rxTypes[rx] + return self.knownRxTypes[rx][1] + + def _projComp(self, rx): + """Component projection (real/imag)""" + if type(rx) is int: rx = self.rxTypes[rx] + return self.knownRxTypes[rx][2] @property - def fieldType(self): + def rxTypes(self): + """A list of the recieve types that are collected at this rxList locations.""" + return self.rxType.split(',') + + @property + def fieldTypes(self): #TODO: This assumes that it has the structure ex, by ... - return self.rxType[0] + return [self._projField(rx) for rx in self.rxTypes] def getP(self, mesh): - if mesh not in self._Ps: - self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.knownRxTypes[self.rxType]) - return self._Ps[mesh] + """ + Returns the projection matrices as a + list for all components collected by + the receivers. + + .. note:: + + Projection matrices are stored as a nested dict, + First gridLocation, then mesh. + """ + P = [] + for rx in self.rxTypes: + gloc = self._projGLoc(rx) + if mesh not in self._Ps[gloc]: + self._Ps[gloc][mesh] = mesh.getInterpolationMat(self.locs, gloc) + P += [self._Ps[gloc][mesh]] + return P class TxFDEM(Survey.BaseTx): @@ -43,25 +85,39 @@ class TxFDEM(Survey.BaseTx): @property def nD(self): """Number of data""" - return self.rxList.locs.shape[0] - + return self.rxList.locs.shape[0]*len(self.rxList.rxTypes) def projectFields(self, mesh, u): - if self.rxList.rxType in ['ex', 'ey', 'ez']: - u_part = u[self, 'e'] - elif self.rxList.rxType in ['bx', 'by', 'bz']: - u_part = u[self, 'b'] - else: - raise NotImplemented('Unknown receiver type.') + nRt = len(self.rxList.rxTypes) + Pu = range(nRt) + Ps = self.rxList.getP(mesh) - P = self.rxList.getP(mesh) - Pu = P*u_part - return Pu + for ii, rx in enumerate(self.rxList.rxTypes): + fieldType = self.rxList._projField(rx) + u_part_complex = u[self, fieldType] + # get the real or imag component + real_or_imag = self.rxList._projComp(rx) + u_part = getattr(u_part_complex, real_or_imag) + Pu[ii] = Ps[ii]*u_part - def projectFieldsDeriv(self, mesh, u): - P = self.rxList.getP(mesh) - return P + return np.concatenate(Pu) + + def projectFieldsDeriv(self, mesh, u, v, adjoint=False): + Ps = self.rxList.getP(mesh) + V = v.reshape((-1, len(Ps)), order='F') + Pvs = range(len(Ps)) + for ii, rx in enumerate(self.rxList.rxTypes): + Pv = Ps[ii] * V[:, ii] + real_or_imag = self.rxList._projComp(rx) + if real_or_imag == 'imag': + Pv = 1j*Pv + elif real_or_imag == 'real': + Pv = Pv.astype(complex) + else: + raise NotImplementedError('must be real or imag') + + return np.concatenate(Pvs) class FieldsFDEM(object): diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 9212d79e..eeeb9df8 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -18,9 +18,9 @@ def getProblem(fdemType): x = np.linspace(5,10,3) XYZ = Utils.ndgrid(x,x,np.r_[0]) if fdemType == 'e': - rxList = EM.FDEM.RxListFDEM(XYZ, 'bx') + rxList = EM.FDEM.RxListFDEM(XYZ, 'bxr,exi') elif fdemType == 'b': - rxList = EM.FDEM.RxListFDEM(XYZ, 'ex') + rxList = EM.FDEM.RxListFDEM(XYZ, 'exi') else: raise NotImplementedError() @@ -29,9 +29,9 @@ def getProblem(fdemType): x = np.linspace(5,10,3) XYZ = Utils.ndgrid(x,x,np.r_[0]) if fdemType == 'e': - rxList = EM.FDEM.RxListFDEM(XYZ, 'ey') + rxList = EM.FDEM.RxListFDEM(XYZ, 'eyi') elif fdemType == 'b': - rxList = EM.FDEM.RxListFDEM(XYZ, 'ey') + rxList = EM.FDEM.RxListFDEM(XYZ, 'eyr') else: raise NotImplementedError() diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index 3fdafa4d..cf6d35f2 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -21,7 +21,7 @@ class FDEM_analyticTests(unittest.TestCase): x = np.linspace(-10,10,5) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) - rxList = EM.FDEM.RxListFDEM(XYZ, 'ex') + rxList = EM.FDEM.RxListFDEM(XYZ, 'exi') Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList) survey = EM.FDEM.SurveyFDEM([Tx0]) diff --git a/simpegEM/Tests/test_FieldsObject.py b/simpegEM/Tests/test_FieldsObject.py index 9eebb4a5..27ee18f3 100644 --- a/simpegEM/Tests/test_FieldsObject.py +++ b/simpegEM/Tests/test_FieldsObject.py @@ -5,20 +5,27 @@ import simpegEM as EM class FieldsTest(unittest.TestCase): def setUp(self): - x = np.linspace(5,10,3) - XYZ = Utils.ndgrid(x,x,np.r_[0]) - rxList = EM.FDEM.RxListFDEM(XYZ, 'ex') - Tx0 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) - Tx1 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) - Tx2 = EM.FDEM.TxFDEM(None, 'VMD', 2, rxList) - Tx3 = EM.FDEM.TxFDEM(None, 'VMD', 1, rxList) - txList = [Tx0,Tx1,Tx2,Tx3] mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) + x = np.linspace(5,10,3) + XYZ = Utils.ndgrid(x,x,np.r_[0.]) + txLoc = np.r_[0,0,0.] + rxList0 = EM.FDEM.RxListFDEM(XYZ, 'exi,exr,eyi,eyr,ezi,ezr') + Tx0 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., rxList0) + rxList1 = EM.FDEM.RxListFDEM(XYZ, 'bxi,bxr,byi,byr,bzi,bzr') + Tx1 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., rxList1) + rxList2 = EM.FDEM.RxListFDEM(XYZ, 'bxi,eyr') + Tx2 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., rxList2) + rxList3 = EM.FDEM.RxListFDEM(XYZ, 'bxi') + Tx3 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., rxList3) + Tx4 = EM.FDEM.TxFDEM(txLoc, 'VMD', 1., rxList0) + txList = [Tx0,Tx1,Tx2,Tx3,Tx4] survey = EM.FDEM.SurveyFDEM(txList) self.F = EM.FDEM.FieldsFDEM(mesh, survey) self.D = EM.FDEM.DataFDEM(survey) self.Tx0 = Tx0 self.Tx1 = Tx1 + self.mesh = mesh + self.XYZ = XYZ def test_SetGet(self): F = self.F @@ -82,7 +89,27 @@ class FieldsTest(unittest.TestCase): D2 = EM.FDEM.DataFDEM(self.D.survey, V) self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D))) + def test_FieldProjections(self): + F = self.F + for freq in F.survey.freqs: + nFreq = F.survey.nTx[freq] + e = np.random.rand(F.mesh.nE, nFreq) + b = np.random.rand(F.mesh.nF, nFreq) + F[freq] = {'b':b,'e':e} + Txs = F.survey.getTransmitters(freq) + for ii, tx in enumerate(Txs): + dat = tx.projectFields(self.mesh, F) + self.assertTrue(dat.dtype == float) + dat = dat.reshape((self.XYZ.shape[0], len(tx.rxList.rxTypes)), order='F') + for jj, rx in enumerate(tx.rxList.rxTypes): + fieldType = tx.rxList._projField(rx) + u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType] + real_or_imag = tx.rxList._projComp(rx) + u = getattr(u, real_or_imag) + gloc = tx.rxList._projGLoc(rx) + d = self.mesh.getInterpolationMat(self.XYZ, gloc)*u + self.assertTrue(np.all(dat[:, jj] == d))