From 5c4ec4a1cd32a9502e06fecb791492764a0e2765 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Mar 2014 09:55:34 -0700 Subject: [PATCH] Updates to data objects. --- simpegEM/FDEM/SurveyFDEM.py | 24 ++++- simpegEM/FDEM/getInterpmat.py | 156 ---------------------------- simpegEM/Tests/test_FDEM.py | 18 +++- simpegEM/Tests/test_FieldsObject.py | 15 +++ 4 files changed, 49 insertions(+), 164 deletions(-) delete mode 100644 simpegEM/FDEM/getInterpmat.py diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index a0e5a33b..82404962 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -160,9 +160,11 @@ class FieldsFDEM(object): class DataFDEM(object): """docstring for DataFDEM""" - def __init__(self, survey): + def __init__(self, survey, v=None): self.survey = survey self._dataDict = {} + if v is not None: + self.fromvec(v) def _ensureCorrectKey(self, key): if key not in self.survey.txList: @@ -176,11 +178,22 @@ class DataFDEM(object): def __getitem__(self, key): self._ensureCorrectKey(key) + if key not in self._dataDict: + raise Exception('Data for transmitter has not yet been set.') return self._dataDict[key] def tovec(self): - D = self._dataDict - return np.concatenate([D[k] for k in D]) + return np.concatenate([self[tx] for tx in self.survey.txList]) + + def fromvec(self, v): + v = Utils.mkvc(v) + assert v.size == self.survey.nD, 'v must have the correct number of data.' + indBot, indTop = 0, 0 + for tx in self.survey.txList: + indTop += tx.nD + self[tx] = v[indBot:indTop] + indBot += tx.nD + @@ -210,6 +223,11 @@ class SurveyFDEM(Survey.BaseSurvey): Survey.BaseSurvey.__init__(self, **kwargs) + @property + def nD(self): + """Number of data""" + return np.array([tx.nD for tx in self.txList]).sum() + @property def freqs(self): """Frequencies""" diff --git a/simpegEM/FDEM/getInterpmat.py b/simpegEM/FDEM/getInterpmat.py deleted file mode 100644 index 808d7716..00000000 --- a/simpegEM/FDEM/getInterpmat.py +++ /dev/null @@ -1,156 +0,0 @@ -import numpy as np -import scipy.sparse as sp -from SimPEG import utils, TensorMesh -from SimPEG.utils import spzeros, mkvc - -def interpmat(x,y,z,xr,yr,zr): - - """ Local nterpolation computed for each receiver point in turn """ - - nx = max(x.shape) - ny = max(y.shape) - nz = max(z.shape) - npts = max(xr.shape) - - Q = sp.lil_matrix((npts, nx*ny*nz)) - - for i in range(npts): - # in x-direction - im = np.argmin(abs(x-xr[i])) - if xr[i] - x[im] >= 0: # Point on the left - ind_x1 = im - ind_x2 = im+1 - elif xr[i] - x[im] < 0: # Point on the right - ind_x1 = im-1 - ind_x2 = im - dx1 = xr[i] - x[ind_x1] - dx2 = x[ind_x2] - xr[i] - # in y-direction - im = np.argmin(abs(y-yr[i])) - if yr[i] - y[im] >= 0: # Point on the left - ind_y1 = im - ind_y2 = im+1 - elif yr[i] - y[im] < 0: # Point on the right - ind_y1 = im-1 - ind_y2 = im - dy1 = yr[i] - y[ind_y1] - dy2 = y[ind_y2] - yr[i] - # in z-direction - im = np.argmin(abs(z-zr[i])) - if zr[i] - z[im] >= 0: # Point on the left - ind_z1 = im - ind_z2 = im+1 - elif zr[i] - z[im] < 0: # Point on the right - ind_z1 = im-1 - ind_z2 = im - dz1 = zr[i] - z[ind_z1] - dz2 = z[ind_z2] - zr[i] - dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) - - Dx = x[ind_x2] - x[ind_x1] - Dy = y[ind_y2] - y[ind_y1] - Dz = z[ind_z2] - z[ind_z1] - - # Get the row in the matrix - - inds = utils.sub2ind((nx,ny,nz),[ - ( ind_x1, ind_y2, ind_z1), - ( ind_x1, ind_y1, ind_z1), - ( ind_x2, ind_y1, ind_z1), - ( ind_x2, ind_y2, ind_z1), - ( ind_x1, ind_y1, ind_z2), - ( ind_x1, ind_y2, ind_z2), - ( ind_x2, ind_y1, ind_z2), - ( ind_x2, ind_y2, ind_z2)]) - - vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz), - (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz), - (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz), - (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz), - (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz), - (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz), - (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz), - (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)] - - Q[i, mkvc(inds)] = vals - Q = Q.tocsr() - return Q - -def getInterpmat(mesh, rxLoc, dataType): - """ """ - xr = rxLoc[:,0] - yr = rxLoc[:,1] - zr = rxLoc[:,2] - nrx = rxLoc.shape[0] - if dataType == 'fx': - Qx = interpmat(np.unique(mesh.gridFx[:,0]), - np.unique(mesh.gridFx[:,1]), - np.unique(mesh.gridFx[:,2]), - xr,yr,zr) - Q = sp.hstack([Qx,spzeros(nrx,mesh.nF[1]),spzeros(nrx,mesh.nF[2])]) - elif dataType == 'fy': - Qy = interpmat(np.unique(mesh.gridFy[:,0]), - np.unique(mesh.gridFy[:,1]), - np.unique(mesh.gridFy[:,2]), - xr,yr,zr) - Q = sp.hstack([spzeros(nrx,mesh.nF[0]),Qy,spzeros(nrx,mesh.nF[2])]) - elif dataType == 'fz': - Qz = interpmat(np.unique(mesh.gridFz[:,0]), - np.unique(mesh.gridFz[:,1]), - np.unique(mesh.gridFz[:,2]), - xr,yr,zr) - Q = sp.hstack([spzeros(nrx,mesh.nF[0]),spzeros(nrx,mesh.nF[1]),Qz]) - elif dataType == 'ex': - Qx = interpmat(np.unique(mesh.gridEx[:,0]), - np.unique(mesh.gridEx[:,1]), - np.unique(mesh.gridEx[:,2]), - xr, yr, zr) - Q = sp.hstack([Qx,spzeros(nrx,mesh.nE[1]),spzeros(nrx,mesh.nE[2])]) - elif dataType == 'ey': - Qy = interpmat(np.unique(mesh.gridEy[:,0]), - np.unique(mesh.gridEy[:,1]), - np.unique(mesh.gridEy[:,2]), - xr, yr, zr) - Q = sp.hstack([spzeros(nrx,mesh.nE[0]),Qy,spzeros(nrx,mesh.nE[2])]) - elif dataType == 'ez': - Qz = interpmat(np.unique(mesh.gridEz[:,0]), - np.unique(mesh.gridEz[:,1]), - np.unique(mesh.gridEz[:,2]), - xr,yr,zr) - Q = sp.hstack([spzeros(nrx,mesh.nE[0]),spzeros(nrx,mesh.nE[1]),Qz]) - else: - assert(True), "Input either face (fx, fy, fz) or edge (ex, ey, ez) option" - return Q - -if __name__ == '__main__': - pad = 1 - padfactor = 1.5 - cs = 100 - xpad = cs*(np.ones(pad)*padfactor)**np.arange(pad) - ypad = cs*(np.ones(pad)*padfactor)**np.arange(pad) - zpad = cs*(np.ones(pad)*padfactor)**np.arange(pad) - - core = 10 - xcore = cs*np.ones(core) - ycore = cs*np.ones(core) - zcore = cs*np.ones(core) - - hx = np.r_[xpad[::-1],xcore, cs, xcore,xpad] - hy = np.r_[ypad[::-1],ycore, cs, ycore, ypad] - hz = np.r_[zpad[::-1],zcore,zcore, zpad] - x0 = np.array([-np.sum(hx)/2, -np.sum(hy)/2, -np.sum(hz)/2], ) - mesh = TensorMesh([hx, hy, hz],x0) - - xr1 = np.linspace(-500,500,5) - yr1 = np.linspace(-500,500,5) - zr1 = 0 - xr, yr = np.meshgrid(xr1, yr1, indexing='ij') - zr = np.ones((xr.shape[0],xr.shape[1]))*zr1 - xr = mkvc(xr) - yr = mkvc(yr) - zr = mkvc(zr) - rxLoc = np.c_[xr, yr, zr] - Q = getInterpmat(mesh, rxLoc, 'ex') - - print Q - diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index a2378d75..836ef717 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -26,10 +26,10 @@ class FDEM_bDerivTests(unittest.TestCase): rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') Tx0 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) - dat = EM.FDEM.SurveyFDEM([Tx0]) + survey = EM.FDEM.SurveyFDEM([Tx0]) prb = EM.FDEM.ProblemFDEM_e(model) - prb.pair(dat) + prb.pair(survey) sigma = np.log(np.ones(mesh.nC)*1e-3) @@ -42,15 +42,23 @@ class FDEM_bDerivTests(unittest.TestCase): self.sigma = sigma self.prb = prb - self.dat = dat + self.survey = survey - def test_JVec(self): + def test_Jvec(self): x0 = self.sigma def fun(x): - return self.dat.dpred(x), lambda x: self.prb.Jvec(x0, x) + return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x) passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False) self.assertTrue(passed) + # def test_Jtvec(self): + # v = self.survey.nD + # x0 = self.sigma + # def fun(x): + # return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x) + # passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False) + # self.assertTrue(passed) + if __name__ == '__main__': unittest.main() diff --git a/simpegEM/Tests/test_FieldsObject.py b/simpegEM/Tests/test_FieldsObject.py index 27137c75..30ef0614 100644 --- a/simpegEM/Tests/test_FieldsObject.py +++ b/simpegEM/Tests/test_FieldsObject.py @@ -16,6 +16,7 @@ class FieldsTest(unittest.TestCase): mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) survey = EM.FDEM.SurveyFDEM(txList) self.F = EM.FDEM.FieldsFDEM(mesh, survey) + self.D = EM.FDEM.DataFDEM(survey) self.Tx0 = Tx0 self.Tx1 = Tx1 @@ -68,6 +69,20 @@ class FieldsTest(unittest.TestCase): self.assertRaises(AssertionError, EM.FDEM.SurveyFDEM, txs) + def test_dataFDEM(self): + V = [] + for tx in self.D.survey.txList: + v = np.random.rand(tx.nD) + V += [v] + self.D[tx] = v + self.assertTrue(np.all(v == self.D[tx])) + V = np.concatenate(V) + self.assertTrue(np.all(V == Utils.mkvc(self.D))) + + D2 = EM.FDEM.DataFDEM(self.D.survey, V) + self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D))) + +