From cc3f2d867cc1daaf5bc3df9ff74089e6c6c3f878 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 13 Feb 2014 17:42:51 -0800
Subject: [PATCH] Data projections

---
 docs/api_TDEM_derivation.rst | 16 ++++++++++++++--
 simpegEM/TDEM/BaseTDEM.py    |  5 +++++
 simpegEM/TDEM/DataTDEM.py    | 24 ++++++++++++++++++++++--
 simpegEM/TDEM/FieldsTDEM.py  |  4 ++++
 4 files changed, 45 insertions(+), 4 deletions(-)

diff --git a/docs/api_TDEM_derivation.rst b/docs/api_TDEM_derivation.rst
index dc247988..e38be920 100644
--- a/docs/api_TDEM_derivation.rst
+++ b/docs/api_TDEM_derivation.rst
@@ -1,4 +1,4 @@
-.. _api_TDEM:
+.. _api_TDEM_derivation:
 
 
 .. math::
@@ -229,7 +229,7 @@ with
 
 
 Implementing **J** times a vector
-****************************************
+*********************************
 
 Multiplying **J** onto a vector can be broken into three steps
 
@@ -290,3 +290,15 @@ and
         + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t+1)} + \MfMui \vec{p}_b^{(t+1)} \\
         \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
     \end{align}
+
+
+
+Implementing \\(\\mathbf{J}^\\top\\) times a vector
+*********************************
+
+Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps
+
+
+* Compute \\(\\vec{u} = \\mathbf{Q}^\\top \\vec{v}\\)
+* Solve \\(\\hat{\\mathbf{A}}^\\top \\vec{y} = \\vec{u}\\)
+* Compute \\(\\vec{w} = -\\mathbf{G}^\\top y\\)
diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py
index 4e9f905a..943c2611 100644
--- a/simpegEM/TDEM/BaseTDEM.py
+++ b/simpegEM/TDEM/BaseTDEM.py
@@ -77,6 +77,11 @@ class MixinTimeStuff(object):
         self._dt = dt
         self._nsteps = nsteps
 
+    @property
+    def nTimes(self):
+        return self.times.size
+
+
 class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
     """docstring for ProblemTDEM1D"""
     def __init__(self, mesh, model, **kwargs):
diff --git a/simpegEM/TDEM/DataTDEM.py b/simpegEM/TDEM/DataTDEM.py
index b69f0c4a..73946ac1 100644
--- a/simpegEM/TDEM/DataTDEM.py
+++ b/simpegEM/TDEM/DataTDEM.py
@@ -1,5 +1,6 @@
-from SimPEG import Utils
+from SimPEG import Utils, np
 from SimPEG.Data import BaseData
+from FieldsTDEM import FieldsTDEM
 
 class DataTDEM1D(BaseData):
     """
@@ -11,13 +12,32 @@ class DataTDEM1D(BaseData):
     rxLoc = None #: rxLoc
     rxType = None #: rxType
     timeCh = None #: timeCh
+    nTx    = 1 #: Number of transmitters
+
+    @property
+    def nTimeCh(self):
+        """Number of time channels"""
+        return self.timeCh.size
 
     def __init__(self, **kwargs):
         BaseData.__init__(self, **kwargs)
         Utils.setKwargs(self, **kwargs)
 
     def projectFields(self, u):
-        return self.Qrx.dot(u.b[:,:,0].T)
+        #TODO: this is hardcoded to 1Tx
+        return self.Qrx.dot(u.b[:,:,0].T).T
+
+    def projectFieldsAdjoint(self, d):
+        # TODO: make the following self.nTimeCh
+        d = d.reshape((self.prob.nTimes, self.nTx), order='F')
+        #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??)
+        ii = 0
+        F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nTimes, 'b')
+        for ii in range(self.prob.nTimes):
+            b = self.Qrx.T*d[ii,:]
+            F.set_b(b, ii)
+            F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii)
+        return F
 
     ####################################################
     # Interpolation Matrices
diff --git a/simpegEM/TDEM/FieldsTDEM.py b/simpegEM/TDEM/FieldsTDEM.py
index c6d59d58..e1bb482b 100644
--- a/simpegEM/TDEM/FieldsTDEM.py
+++ b/simpegEM/TDEM/FieldsTDEM.py
@@ -60,10 +60,14 @@ class FieldsTDEM(object):
         if self.b is None:
             self.b = np.zeros((self.nTimes, np.sum(self.mesh.nF), self.nTx))
             self.b[:] = np.nan
+        if len(b.shape) == 1:
+            b = b[:, np.newaxis]
         self.b[ind,:,:] = b
 
     def set_e(self, e, ind):
         if self.e is None:
             self.e = np.zeros((self.nTimes, np.sum(self.mesh.nE), self.nTx))
             self.e[:] = np.nan
+        if len(e.shape) == 1:
+            e = e[:, np.newaxis]
         self.e[ind,:,:] = e