minor fix

Modify testings ... for TDEM
remove temporary print statemement.
2026-09-13 13:03:14 +08:00 · 2016-06-28 19:16:41 -07:00 · 2016-06-28 19:15:20 -07:00 · 2016-06-28 19:01:32 -07:00 · 2016-06-28 19:00:39 -07:00 · 2016-06-28 11:42:30 -07:00
31 changed files with 2436 additions and 1084 deletions
@@ -144,6 +144,7 @@ class BetaSchedule(InversionDirective):
            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
            self.invProb.beta /= self.coolingFactor

+
 class TargetMisfit(InversionDirective):

    chifact = 1.
@@ -242,12 +243,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
        # Save the file as a npz
        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)

-
-# class UpdateReferenceModel(Parameter):
-
-#     mref0 = None
-
-#     def nextIter(self):
 #         mref = getattr(self, 'm_prev', None)
 #         if mref is None:
 #             if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -258,56 +253,138 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 class Update_IRLS(InversionDirective):

    eps_min = None
+    eps_p = None
+    eps_q = None
+    norms = [2.,2.,2.,2.]
    factor = None
    gamma = None
    phi_m_last = None
    phi_d_last = None
+    f_old = None
+    f_min_change = 1e-2
+    beta_tol = 5e-2

+    # Solving parameter for IRLS (mode:2)
+    IRLSiter   = 0
+    minGNiter = 5
+    maxIRLSiter = 10
+    iterStart = 0
+
+    # Beta schedule
+    coolingFactor = 2.
+    coolingRate = 1
+
+    mode = 1
+
+    @property
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD*0.5
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val

    def initialize(self):

-        # Scale the regularization for changes in norm
-        if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
-            self.reg.gamma = 1.
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.gamma = self.phi_m_last / phim_new
-
-        self.reg.curModel = self.invProb.curModel
-        self.reg.gamma = self.gamma
-
-        if getattr(self, 'phi_d_last', None) is None:
-            self.phi_d_last = self.invProb.phi_d
+        if self.mode == 1:
+            self.reg.norms = [2., 2., 2., 2.]

    def endIter(self):
-        # Cool the threshold parameter if required
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor

-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
+        # After reaching target misfit with l2-norm, switch to IRLS (mode:2)
+        if self.invProb.phi_d < self.target and self.mode == 1:
+            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
+
+            self.mode = 2
+            print self.eps_p, self.eps_q, self.norms
+            self.reg.eps_p = self.eps_p
+            self.reg.eps_q = self.eps_q
+            self.reg.norms = self.norms
+            self.coolingFactor = 1.
+            self.coolingRate = 1
+            self.iterStart = self.opt.iter
+            self.phi_d_last = self.invProb.phi_d
+            self.phi_m_last = self.invProb.phi_m_last
+
+            self.reg.l2model = self.invProb.curModel
+            self.reg.curModel = self.invProb.curModel
+
+            if getattr(self, 'f_old', None) is None:
+                self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
+
+        # Beta Schedule
+        if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+            self.invProb.beta /= self.coolingFactor
+
+
+        # Only update after GN iterations
+        if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
+
+            self.IRLSiter += 1
+
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.f_change = np.abs(self.f_old - phim_new) / self.f_old
+
+            print "Regularization decrease: %6.3e" % (self.f_change)
+
+            # Check for maximum number of IRLS cycles
+            if self.IRLSiter == self.maxIRLSiter:
+                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
+                self.opt.stopNextIteration = True
+                return
+
+            # Check if the function has changed enough
+            if self.f_change < self.f_min_change and self.IRLSiter > 1:
+                print "Minimum decrease in regularization. End of IRLS"
+                self.opt.stopNextIteration = True
+                return
            else:
-                self.reg.eps = eps
+                self.f_old = phim_new

-        # Get phi_m at the end of current iteration
-        self.phi_m_last = self.invProb.phi_m_last
+            # Cool the threshold parameter if required
+            if getattr(self, 'factor', None) is not None:
+                eps = self.reg.eps / self.factor

-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
+                if getattr(self, 'eps_min', None) is not None:
+                    self.reg.eps = np.max([self.eps_min,eps])
+                else:
+                    self.reg.eps = eps

-        # Temporarely set gamma to 1. to get raw phi_m
-        self.reg.gamma = 1.
+            # Get phi_m at the end of current iteration
+            self.phi_m_last = self.invProb.phi_m_last

-        # Compute new model objective function value
-        phim_new = self.reg.eval(self.invProb.curModel)
+            # Reset the regularization matrices so that it is
+            # recalculated for current model
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None

-        # Update gamma to scale the regularization between IRLS iterations
-        self.reg.gamma = self.phi_m_last / phim_new
+             # Update the model used for the IRLS weights
+            self.reg.curModel = self.invProb.curModel

-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
+            # Temporarely set gamma to 1. to get raw phi_m
+            self.reg.gamma = 1.
+
+            # Compute new model objective function value
+            phim_new = self.reg.eval(self.invProb.curModel)
+
+            # Update gamma to scale the regularization between IRLS iterations
+            self.reg.gamma = self.phi_m_last / phim_new
+
+            # Reset the regularization matrices again for new gamma
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
+
+            # Check if misfit is within the tolerance, otherwise scale beta
+            val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+            if np.abs(1.-val) > self.beta_tol:
+                self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d

 class Update_lin_PreCond(InversionDirective):
    """
@@ -360,19 +437,3 @@ class Update_Wj(InversionDirective):
            JtJdiag = JtJdiag / max(JtJdiag)

            self.reg.wght = JtJdiag
-
-class Scale_Beta(InversionDirective):
-    """
-        Instead of a linear cooling schedule, beta is allowed to change based
-        on the ratio between the target misfit and the current data misfit. The
-        update is done only if the misfit is outside some threshold bounds.
-    """
-    tol = 0.05
-
-    def endIter(self):
-
-        # Check if misfit is within the tolerance, otherwise adjust beta
-        val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-        if np.abs(1.-val) > self.tol:
-            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
@@ -60,6 +60,20 @@ class Fields(SimPEG.Problem.Fields):

        return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)

+    def _bSecondary(self, solution, srcList):
+        """
+        Total magnetic flux density is sum of primary and secondary
+
+        :param numpy.ndarray solution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic flux density
+        """
+        if getattr(self, '_bSecondary', None) is None:
+            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
+
+        return  self._bSecondary(solution, srcList)
+
    def _h(self, solution, srcList):
        """
        Total magnetic field is sum of primary and secondary
@@ -124,6 +138,21 @@ class Fields(SimPEG.Problem.Fields):
            return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
        return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)

+    def _bSecondaryDeriv(self, src, du_dm_v, v, adjoint = False):
+        """
+        Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
+
+        :param Src src: sorce
+        :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
+        :param numpy.ndarray v: vector to take sensitivity product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative times a vector (or tuple for adjoint)
+        """
+        # TODO: modify when primary field is dependent on m
+
+        return self._bDeriv(src, du_dm_v, v, adjoint = adjoint)
+
    def _hDeriv(self, src, du_dm_v, v, adjoint = False):
        """
        Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
@@ -471,6 +500,8 @@ class Fields3D_b(Fields):
            return 'E'
        elif fieldType == 'b':
            return 'F'
+        elif fieldType == 'bSecondary':
+            return 'F'
        elif (fieldType == 'h') or (fieldType == 'j'):
            return'CCV'
        else:
@@ -97,6 +97,19 @@ class Point_b(BaseRx):
        self.projField = 'b'
        super(Point_b, self).__init__(locs, orientation, component)

+class Point_bSecondary(BaseRx):
+    """
+    Magnetic flux FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'bSecondary'
+        super(Point_bSecondary, self).__init__(locs, orientation, component)
+

 class Point_h(BaseRx):
    """
@@ -555,7 +555,7 @@ class CircularLoop(BaseSrc):
            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)

        else:
-            srcfct = MagneticDipoleVectorPotential
+            srcfct = MagneticLoopVectorPotential
            ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
            ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
            az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
@@ -0,0 +1,142 @@
+import numpy as np
+import scipy.sparse as sp
+import SimPEG
+from SimPEG import Utils
+from SimPEG.EM.Utils import omega
+from SimPEG.Utils import Zero, Identity
+
+class Fields(SimPEG.Problem.TimeFields):
+    """
+
+    Fancy Field Storage for a TDEM survey. Only one field type is stored for
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[srcList,'e']
+        b = f[srcList,'b']
+
+    If accessing all sources for a given field, use the :code:`:`
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[:,'e']
+        b = f[:,'b']
+
+    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
+    """
+
+    knownFields = {}
+    dtype = float
+
+    def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+        if adjoint is True:
+            return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
+        return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
+
+    def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+        if adjoint is True:
+            return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+        return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+class Fields_Derivs(Fields):
+    knownFields = {
+                    'bDeriv': 'F',
+                    'eDeriv': 'E',
+                    'hDeriv': 'E',
+                    'jDeriv': 'F'
+                  }
+
+
+class Fields_b(Fields):
+    """Fancy Field Storage for a TDEM survey."""
+    knownFields = {'bSolution': 'F'}
+    aliasFields = {
+                    'b': ['bSolution', 'F', '_b'],
+                    'e': ['bSolution', 'E', '_e'],
+                  }
+
+    def startup(self):
+        self.MeSigmaI      = self.survey.prob.MeSigmaI
+        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
+        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
+        self.MfMui         = self.survey.prob.MfMui
+
+    def _b(self, bSolution, srcList, tInd):
+        return bSolution
+
+    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
+        return Identity()*dun_dm_v
+
+    def _bDeriv_m(self, tInd, src, v, adjoint=False):
+        return Zero()
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+    def _e(self, bSolution, srcList, tInd):
+        e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
+        for i, src in enumerate(srcList):
+            _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
+            e[:,i] = e[:,i] - self.MeSigmaI * S_e
+        return e
+
+    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
+        if adjoint is True:
+            return self.MfMui.T * ( self.edgeCurl * ( self.MeSigmaI.T * dun_dm_v ) )
+        return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) )
+
+    def _eDeriv_m(self, tInd, src, v, adjoint = False):
+        _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
+        bSolution = self[[src],'bSolution',tInd]
+
+        _, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, adjoint=adjoint)
+
+        if adjoint is True:
+            return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv(self.MeSigmaI.T * v)
+
+        return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v)
+
+
+
+class Fields_e(Fields):
+    """Fancy Field Storage for a TDEM survey."""
+    knownFields = {'eSolution': 'E'}
+    aliasFields = {
+                    'e': ['eSolution', 'E', '_e'],
+                    'b': ['eSolution', 'F', '_b'],
+                  }
+
+    def startup(self):
+        self.MeSigmaI      = self.survey.prob.MeSigmaI
+        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
+        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
+        self.MfMui         = self.survey.prob.MfMui
+
+
+    def _e(self, eSolution, srcList, tInd):
+        return eSolution
+
+    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
+        return dun_dm_v
+
+    def _eDeriv_m(self, tInd, src, v, adjoint = False):
+        return Zero()
+
+    def _b(self, eSolution, srcList, tInd):
+        raise NotImplementedError
+
+    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
+        raise NotImplementedError
+
+    def _bDeriv_m(self, tInd, src, v, adjoint=False):
+        raise NotImplementedError
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
@@ -0,0 +1,246 @@
+import SimPEG
+from SimPEG import np, Utils
+from SimPEG.Utils import Zero, Identity
+from scipy.constants import mu_0
+from SimPEG.EM.Utils import *
+
+####################################################
+# Sources
+####################################################
+
+class BaseWaveform(object):
+
+    def __init__(self, offTime=0., hasInitialFields=False):
+        self.offTime = offTime
+        self.hasInitialFields = hasInitialFields
+
+    def _assertMatchesPair(self, pair):
+        assert (isinstance(self, pair)
+            ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__)
+
+    def eval(self, time):
+        raise NotImplementedError
+
+    def evalDeriv(self, time):
+        raise NotImplementedError # needed for E-formulation
+
+
+class StepOffWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        return 0.
+
+
+class RawWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        raise NotImplementedError('RawWaveform has not been implemented, you should write it!')
+
+
+class TriangularWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        raise NotImplementedError('TriangularWaveform has not been implemented, you should write it!')
+
+
+
+class BaseSrc(SimPEG.Survey.BaseSrc):
+
+    # rxPair = Rx
+    integrate = True
+    waveformPair = BaseWaveform
+
+    @property
+    def waveform(self):
+        "A waveform instance is not None"
+        return getattr(self, '_waveform', None)
+    @waveform.setter
+    def waveform(self, val):
+        if self.waveform is None:
+            val._assertMatchesPair(self.waveformPair)
+            self._mapping = val
+        else:
+            self._mapping = self.PropMap(val)
+
+
+    def __init__(self, rxList, waveform = StepOffWaveform(), **kwargs):
+        self.waveform = waveform
+        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
+
+
+    def bInitial(self, prob):
+        return Zero()
+
+    def bInitialDeriv(self, prob, v=None, adjoint=False):
+        return Zero()
+
+    def eInitial(self, prob):
+        return Zero()
+
+    def eInitialDeriv(self, prob, v=None, adjoint=False):
+        return Zero()
+
+    def eval(self, prob, time):
+        S_m = self.S_m(prob, time)
+        S_e = self.S_e(prob, time)
+        return S_m, S_e
+
+    def evalDeriv(self, prob, time, v=None, adjoint=False):
+        if v is not None:
+            return self.S_mDeriv(prob, time, v, adjoint), self.S_eDeriv(prob, time, v, adjoint)
+        else:
+            return lambda v: self.S_mDeriv(prob, time, v, adjoint), lambda v: self.S_eDeriv(prob, time, v, adjoint)
+
+    def S_m(self, prob, time):
+        return Zero()
+
+    def S_e(self, prob, time):
+        return Zero()
+
+    def S_mDeriv(self, prob, time, v=None, adjoint=False):
+        return Zero()
+
+    def S_eDeriv(self, prob, time, v=None, adjoint=False):
+        return Zero()
+
+
+class MagDipole(BaseSrc):
+
+    waveform = None
+    loc = None
+    orientation = 'Z'
+    moment = 1.
+    mu = mu_0
+
+    def __init__(self, rxList, **kwargs):
+        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
+        self.integrate = False
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def _bfromVectorPotential(self, prob):
+        if prob._eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif prob._eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+
+        else:
+            srcfct = MagneticDipoleVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+
+
+    def bInitial(self, prob):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        return self._bfromVectorPotential(prob)
+
+    def eInitial(self, prob):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        b = self.bInitial(prob)
+        MeSigmaI = prob.MeSigmaI
+        MfMui = prob.MfMui
+        C = prob.mesh.edgeCurl
+
+        return MeSigmaI * (C.T * (MfMui * b))
+
+    def eInitialDeriv(self, prob, v=None, adjoint=False):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        b = self.bInitial(prob)
+        MeSigmaIDeriv = prob.MeSigmaIDeriv
+        MfMui = prob.MfMui
+        C = prob.mesh.edgeCurl
+        S_e = self.S_e(prob, prob.t0)
+
+        # S_e doesn't depend on the model
+
+        if adjoint:
+            return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v
+
+        return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v
+
+
+    def S_m(self, prob, time):
+        if self.waveform.hasInitialFields is False:
+            raise NotImplementedError
+        return Zero()
+
+    def S_e(self, prob, time):
+        if self.waveform.hasInitialFields is False:
+            raise NotImplementedError
+        return Zero()
+
+class CircularLoop(MagDipole):
+
+    waveform = None
+    loc = None
+    orientation = 'Z'
+    radius = None
+    mu = mu_0
+
+    def __init__(self, rxList, **kwargs):
+        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
+        self.integrate = False
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def _bfromVectorPotential(self, prob):
+        if prob._eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif prob._eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', radius=self.radius, mu=self.mu)
+
+        else:
+            srcfct = MagneticLoopVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, radius=self.radius)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, radius=self.radius)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, radius=self.radius)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+
@@ -1,10 +1,16 @@
-from SimPEG import Utils, Survey, np
-from SimPEG.Survey import BaseSurvey
+import SimPEG
+from SimPEG import np, Utils
+from SimPEG.Utils import Zero, Identity
+from scipy.constants import mu_0
 from SimPEG.EM.Utils import *
-from BaseTDEM import FieldsTDEM
+import SrcTDEM as Src


-class RxTDEM(Survey.BaseTimeRx):
+####################################################
+# Receivers
+####################################################
+
+class Rx(SimPEG.Survey.BaseTimeRx):

    knownRxTypes = {
                    'ex':['e', 'Ex', 'N'],
@@ -21,7 +27,7 @@ class RxTDEM(Survey.BaseTimeRx):
                   }

    def __init__(self, locs, times, rxType):
-        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType)

    @property
    def projField(self):
@@ -56,144 +62,60 @@ class RxTDEM(Survey.BaseTimeRx):
        u_part = Utils.mkvc(u[src, self.projField, :])
        return P*u_part

-    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
+    def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False):
        P = self.getP(mesh, timeMesh)
-
        if not adjoint:
-            return P * Utils.mkvc(v[src, self.projField, :])
+            return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :])
        elif adjoint:
-            return P.T * v[src, self]
+            # dP_dF_T = P.T * v #[src, self]
+            # newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN )
+            return P.T * v #np.reshape(dP_dF_T, newshape, order='F')


-class SrcTDEM(Survey.BaseSrc):
-    rxPair = RxTDEM
-    radius = None
+####################################################
+# Survey
+####################################################

-    def getInitialFields(self, mesh):
-        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
-        return F0
-
-    def getJs(self, mesh, time):
-        return None
-
-
-class SrcTDEM_VMD_MVP(SrcTDEM):
-
-    def __init__(self,rxList,loc,waveformType="STEPOFF"):
-        self.loc = loc
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Vertical magnetic dipole, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-            else:
-                raise Exception('Unknown mesh for VMD')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-        else:
-            raise Exception('Unknown mesh for VMD')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
-        self.loc = loc
-        self.radius = radius
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Circular Loop, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-            else:
-                raise Exception('Unknown mesh for CircularLoop')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-        else:
-            raise Exception('Unknown mesh for CircularLoop')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SurveyTDEM(Survey.BaseSurvey):
+class Survey(SimPEG.Survey.BaseSurvey):
    """
-        docstring for SurveyTDEM
+    Time domain electromagnetic survey
    """
-    srcPair = SrcTDEM
+
+    srcPair = Src.BaseSrc
+    rxPair = Rx

    def __init__(self, srcList, **kwargs):
        # Sort these by frequency
        self.srcList = srcList
-        Survey.BaseSurvey.__init__(self, **kwargs)
+        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)

    def eval(self, u):
-        data = Survey.Data(self)
+        data = SimPEG.Survey.Data(self)
        for src in self.srcList:
            for rx in src.rxList:
                data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u)
        return data

    def evalDeriv(self, u, v=None, adjoint=False):
-        assert v is not None, 'v to multiply must be provided.'
+        raise Exception('Use Receivers to project fields deriv.')
+        # assert v is not None, 'v to multiply must be provided.'

-        if not adjoint:
-            data = Survey.Data(self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
-            return data
-        else:
-            f = FieldsTDEM(self.mesh, self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
-                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
-                    if rx.projField not in f: # first time we are projecting
-                        f[src, rx.projField, :] = Ptv
-                    else: # there are already fields, so let's add to them!
-                        f[src, rx.projField, :] += Ptv
-            return f
+        # if not adjoint:
+        #     data = SimPEG.Survey.Data(self)
+        #     for src in self.srcList:
+        #         for rx in src.rxList:
+        #             data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+        #     return data
+        # else:
+        #     f = FieldsTDEM(self.mesh, self)
+        #     for src in self.srcList:
+        #         for rx in src.rxList:
+        #             Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+        #             Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+        #             if rx.projField not in f: # first time we are projecting
+        #                 f[src, rx.projField, :] = Ptv
+        #             else: # there are already fields, so let's add to them!
+        #                 f[src, rx.projField, :] += Ptv
+        #     return f


@@ -0,0 +1,553 @@
+from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.TDEM.SurveyTDEM import Survey as SurveyTDEM
+from SimPEG.EM.TDEM.FieldsTDEM import *
+from scipy.constants import mu_0
+import time
+
+class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
+    """
+    We start with the first order form of Maxwell's equations
+    """
+    surveyPair = SurveyTDEM
+    fieldsPair = Fields
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def fields(self, m):
+        """
+        Solve the forward problem for the fields.
+
+        :param numpy.array m: inversion model (nP,)
+        :rtype numpy.array:
+        :return F: fields
+        """
+
+        tic = time.time()
+        self.curModel = m
+
+        F = self.fieldsPair(self.mesh, self.survey)
+
+        # set initial fields
+        F[:,self._fieldType+'Solution',0] = self.getInitialFields()
+
+        # timestep to solve forward
+        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
+        Ainv = None
+        for tInd, dt in enumerate(self.timeSteps):
+            if Ainv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
+                Ainv.clean()
+                Ainv = None
+
+            if Ainv is None:
+                A = self.getAdiag(tInd)
+                if self.verbose: print 'Factoring...   (dt = %e)'%dt
+                Ainv = self.Solver(A, **self.solverOpts)
+                if self.verbose: print 'Done'
+
+            rhs = self.getRHS(tInd+1) # this is on the nodes of the time mesh
+            Asubdiag = self.getAsubdiag(tInd)
+
+            if self.verbose: print ('    Solving...   (tInd = %i)')% (tInd+1)
+            sol = Ainv * (rhs - Asubdiag * F[:,self._fieldType+'Solution',tInd]) # taking a step
+
+            if self.verbose: print '    Done...'
+
+            if sol.ndim == 1:
+                sol.shape = (sol.size,1)
+            F[:,self._fieldType+'Solution',tInd+1] = sol
+        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
+        Ainv.clean()
+        return F
+
+
+    def Jvec(self, m, v, f=None):
+        """
+        Jvec computes the sensitivity times a vector
+
+        .. math::
+            \mathbf{J} \mathbf{v} = \\frac{d\mathbf{P}}{d\mathbf{F}} \left( \\frac{d\mathbf{F}}{d\mathbf{u}} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} \\right) \mathbf{v}
+
+        where
+
+        .. math::
+            \mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}}
+        """
+
+        if f is None:
+           f = self.fields(m)
+
+        ftype = self._fieldType + 'Solution' # the thing we solved for
+        self.curModel = m
+
+        # mat to store previous time-step's solution deriv times a vector for each source
+        # size: nu x nSrc
+
+        # this is a bit silly
+
+        # if self._fieldType is 'b' or self._fieldType is 'j':
+        #     ifields = np.zeros((self.mesh.nF, len(Srcs)))
+        # elif self._fieldType is 'e' or self._fieldType is 'h':
+        #     ifields = np.zeros((self.mesh.nE, len(Srcs)))
+
+        # for i, src in enumerate(self.survey.srcList):
+        dun_dm_v = np.hstack([Utils.mkvc(self.getInitialFieldsDeriv(src,v),2) for src in self.survey.srcList]) # can over-write this at each timestep
+        #
+        df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
+
+        Adiaginv = None
+
+        for tInd, dt in zip(range(self.nT), self.timeSteps):
+            if Adiaginv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
+                Adiaginv.clean()
+                Adiaginv = None
+
+            if Adiaginv is None:
+                A = self.getAdiag(tInd)
+                Adiaginv = self.Solver(A, **self.solverOpts)
+
+            Asubdiag = self.getAsubdiag(tInd)
+
+            for i, src in enumerate(self.survey.srcList):
+
+                # here, we are lagging by a timestep, so filling in as we go
+                for projField in set([rx.projField for rx in src.rxList]):
+                    # Seogi: df_duFun?
+                    df_dmFun = getattr(f, '_%sDeriv'%projField, None)
+                    # df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
+                    # This should be called rx.footprint
+                    df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
+
+                un_src = f[src,ftype,tInd+1]
+
+                dA_dm_v   = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
+                dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
+
+                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, f[src,ftype,tInd], v)
+
+                JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
+
+                # step in time and overwrite
+                if tInd != len(self.timeSteps+1):
+                    dun_dm_v[:,i] = Adiaginv * (JRHS - Asubdiag * dun_dm_v[:,i])
+
+        # Seogi: suspcious spot
+        # Jv = self.dataPair(self.survey)
+        Jv = []
+        for src in self.survey.srcList:
+            for rx in src.rxList:
+                # Looping over data class append memory as well!!
+                # Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:]))
+                Jv.append(rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:])))
+        Adiaginv.clean()
+        # del df_dm_v, dun_dm_v, Asubdiag
+        # return Utils.mkvc(Jv)
+        return np.hstack(Jv)
+
+    def Jtvec(self, m, v, f=None):
+
+        """
+        Jvec computes the adjoint of the sensitivity times a vector
+
+        .. math::
+            \mathbf{J}^\\top \mathbf{v} =  \left( \\frac{d\mathbf{u}}{d\mathbf{m}} ^ \\top \\frac{d\mathbf{F}}{d\mathbf{u}} ^ \\top  + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} ^ \\top \\right) \\frac{d\mathbf{P}}{d\mathbf{F}} ^ \\top \mathbf{v}
+
+        where
+
+        .. math::
+            \\frac{d\mathbf{u}}{d\mathbf{m}} ^\\top \mathbf{A}^\\top  + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} ^ \\top = \\frac{d \mathbf{RHS}}{d \mathbf{m}} ^ \\top
+        """
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+        ftype = self._fieldType + 'Solution' # the thing we solved for
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        df_duT_v = Fields_Derivs(self.mesh, self.survey)
+        ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(f[self.survey.srcList[0],ftype,0])), dtype=float) # same size as fields at a single timestep
+
+        JTv = np.zeros(m.shape, dtype=float)
+
+        # Loop over sources and receivers to create a fields object: PT_v, df_duT_v, df_dmT_v
+        PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
+        for src in self.survey.srcList:
+            # Looping over initializing field class is appending memory!
+            # PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
+            # initialize size
+            df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(f[src, self._fieldType, :])
+
+            for rx in src.rxList:
+                print ('_%sDeriv')%(rx.projField)
+                PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
+
+                # PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+
+                for tInd in range(self.nT+1):
+                    cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
+                    df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
+                    JTv = cur[1] + JTv
+
+        del PT_v # no longer need this
+
+        AdiagTinv = None
+
+        # Do the back-solve through time
+        for tIndP in reversed(range(self.nT + 1)):
+            tInd = tIndP - 1
+            if AdiagTinv is not None and (tInd <= self.nT and self.timeSteps[tInd] != self.timeSteps[tInd+1]): # if the previous timestep is the same --> no need to refactor the matrix
+                AdiagTinv.clean()
+                AdiagTinv = None
+
+            # refactor if we need to
+            if AdiagTinv is None and tInd > -1:
+                Adiag = self.getAdiag(tInd)
+                AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
+
+            dAsubdiag_dm_v = Zero()
+
+            if tInd < self.nT - 1:
+                Asubdiag = self.getAsubdiag(tInd+1)
+
+
+            for isrc, src in enumerate(self.survey.srcList):
+                # solve against df_duT_v
+                if tInd >= self.nT-1:
+                    # last timestep (first to be solved)
+                    ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]
+                elif tInd > -1:
+                # else:
+                    ATinv_df_duT_v[isrc,:] = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+                else:
+                    # AdiagTinv = I
+                    ATinv_df_duT_v[isrc,:] = Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])
+                    # - Utils.mkvc(Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+                    # (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+
+                if tInd < self.nT - 1:
+                    dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
+
+                if tInd > -1:
+                    un_src = f[src,ftype,tInd+1]
+                    dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
+                    dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
+
+                    JTv = JTv + Utils.mkvc(- dAT_dm_v - dAsubdiag_dm_v + dRHST_dm_v)
+                else:
+                    # dA_dm_v = self.getInitialFieldsDeriv(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True)
+                    # print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
+                    # print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
+                    # vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
+                    # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
+                    dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
+
+                    JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
+
+
+
+                    # # dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
+                    # dRHST_dm_v0 = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
+                    # dRHST_dm_v1 = self.getInitialFieldsDeriv( Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    # JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
+
+                    # print 'here'
+                    # inFields = self.getInitialFieldsDeriv(f[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    #  # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
+                    # print inFields.shape
+                    # JTv = JTv + inFields
+                # dAsubdiag_dm_v = 0
+
+
+
+        # Missing the 0 step
+
+        # adding du_dm^T * dF_du^T * P^T vfor time 0 (no dRHS_dm_v at time 0)
+        # Asubdiag = self.getAsubdiag(0)
+        # for src in self.survey.srcList:
+        #     for projField in set(rx.projField):
+        #         v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+        #         JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, f[src, ftype, tInd], v, adjoint = True))
+        #         # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
+
+        # del df_duT_v, ATinv_df_duT_v, A, Asubdiag
+        if AdiagTinv is not None:
+            AdiagTinv.clean()
+
+        return Utils.mkvc(JTv).astype(float)
+
+
+
+    def getSourceTerm(self, tInd):
+
+        Srcs = self.survey.srcList
+
+        if self._eqLocs is 'FE':
+            S_m = np.zeros((self.mesh.nF,len(Srcs)))
+            S_e = np.zeros((self.mesh.nE,len(Srcs)))
+        elif self._eqLocs is 'EF':
+            S_m = np.zeros((self.mesh.nE,len(Srcs)))
+            S_e = np.zeros((self.mesh.nF,len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            smi, sei = src.eval(self, self.times[tInd])
+            S_m[:,i] = S_m[:,i] + smi
+            S_e[:,i] = S_e[:,i] + sei
+
+        return S_m, S_e
+
+    def getInitialFields(self):
+
+        Srcs = self.survey.srcList
+
+        if self._fieldType is 'b' or self._fieldType is 'j':
+            ifields = np.zeros((self.mesh.nF, len(Srcs)))
+        elif self._fieldType is 'e' or self._fieldType is 'h':
+            ifields = np.zeros((self.mesh.nE, len(Srcs)))
+
+        for i,src in enumerate(Srcs):
+            ifields[:,i] = ifields[:,i] + getattr(src, '%sInitial'%self._fieldType, None)(self)
+
+        return ifields
+
+    def getInitialFieldsDeriv(self, src, v, adjoint=False):
+
+        if adjoint is False:
+            if self._fieldType is 'b' or self._fieldType is 'j':
+                ifieldsDeriv = np.zeros(self.mesh.nF)
+            elif self._fieldType is 'e' or self._fieldType is 'h':
+                ifieldsDeriv = np.zeros(self.mesh.nE)
+
+        elif adjoint is True:
+            ifieldsDeriv = np.zeros(self.mapping.nP)
+
+        ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
+
+            # ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
+            # ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint)
+            # ifieldsDeriv = ifieldsDeriv.sum()
+
+        return ifieldsDeriv
+
+
+##########################################################################################
+################################ E-B Formulation #########################################
+##########################################################################################
+
+# ------------------------------- Problem_b -------------------------------------------- #
+
+class Problem_b(BaseTDEMProblem):
+    """
+    Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized)
+
+    .. math::
+
+        \mathbf{C} \mathbf{e} + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{s_m} \\\\
+        \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}
+
+    where :math:`\mathbf{s_e}` is an integrated quantity, we eliminate :math:`\mathbf{e}` using
+
+    .. math::
+        \mathbf{e} = \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}
+
+    to obtain a second order semi-discretized system in :math:`\mathbf{b}`
+
+    .. math::
+        \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}  + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
+
+    and moving everything except the time derivative to the rhs gives
+
+    .. math::
+        \\frac{\partial \mathbf{b}}{\partial t} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
+
+    For the time discretization, we use backward euler. To solve for the :math:`n+1`th time step, we have
+
+    .. math::
+        \\frac{\mathbf{b}^{n+1} - \mathbf{b}^{n}}{\mathbf{dt}} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}^{n+1} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1}
+
+    re-arranging to put :math:`\mathbf{b}^{n+1}` on the left hand side gives
+
+    .. math::
+        (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}) \mathbf{b}^{n+1} = \mathbf{b}^{n} + \mathbf{dt}(\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1})
+
+    :param Mesh mesh: mesh
+    :param Mapping mapping: mapping
+    """
+
+    _fieldType = 'b'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields_b
+    surveyPair = SurveyTDEM
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def getAdiag(self, tInd):
+        """
+        System matrix at a given time index
+
+        .. math::
+            (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f})
+
+        """
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MfMui = self.MfMui
+        I = Utils.speye(self.mesh.nF)
+
+        A = 1./dt * I + ( C * ( MeSigmaI * (C.T * MfMui ) ) )
+
+        if self._makeASymmetric is True:
+            return MfMui.T * A
+        return A
+
+    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
+        C = self.mesh.edgeCurl
+        MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x)
+        MfMui = self.MfMui
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = MfMui * v
+            return  MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v )
+
+        ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
+        if self._makeASymmetric is True:
+            return MfMui.T * ADeriv
+        return ADeriv
+
+
+    def getAsubdiag(self, tInd):
+
+        dt = self.timeSteps[tInd]
+        MfMui = self.MfMui
+        Asubdiag = - 1./dt * sp.eye(self.mesh.nF)
+
+        if self._makeASymmetric is True:
+            return MfMui.T * Asubdiag
+
+        return Asubdiag
+
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        return Zero() * v
+
+
+
+    def getRHS(self, tInd):
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MfMui = self.MfMui
+
+        S_m, S_e = self.getSourceTerm(tInd)
+
+        rhs =  (C * (MeSigmaI * S_e) + S_m)
+        if self._makeASymmetric is True:
+            return MfMui.T * rhs
+        return rhs
+
+    def getRHSDeriv(self, tInd, src, v, adjoint=False):
+
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
+        MfMui = self.MfMui
+
+        _, S_e = src.eval(tInd, self)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint)
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = self.MfMui * v
+            if isinstance(S_e, Utils.Zero):
+                MeSigmaIDerivT_v = Utils.Zero()
+            else:
+                MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v
+            RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T *  ( C.T * v ) ) + S_mDeriv(v)
+            return RHSDeriv
+
+        if isinstance(S_e, Utils.Zero):
+            MeSigmaIDeriv_v = Utils.Zero()
+        else:
+            MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v
+
+        RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v)))
+
+        if self._makeASymmetric is True:
+            return self.MfMui.T * RHSDeriv
+        return RHSDeriv
+
+
+# ------------------------------- Problem_e -------------------------------------------- #
+
+class Problem_e(BaseTDEMProblem):
+
+    _fieldType = 'e'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields_e
+    surveyPair = SurveyTDEM
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def getAdiag(self, tInd):
+        """
+        System matrix at a given time index
+
+        """
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigma = self.MeSigma
+
+        return C.T * ( MfMui * C ) + 1./dt * MeSigma
+
+
+    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigmaDeriv = self.MeSigmaDeriv(u)
+
+        if adjoint:
+            return 1./dt * MeSigmaDeriv.T * v
+
+        return 1./dt * MeSigmaDeriv * v
+
+
+    def getAsubdiag(self, tInd):
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+
+        return - 1./dt * self.MeSigma
+
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        dt = self.timeSteps[tInd]
+
+        if adjoint:
+            return - 1./dt * self.MeSigmaDeriv(u).T * v
+
+        return - 1./dt * self.MeSigmaDeriv(u) * v
+
+    def getRHS(self, tInd):
+        return Zero()
+
+    def getRHSDeriv(self, tInd, src, v, adjoint=False):
+        return Zero()
+
+
+
+
+
@@ -1,3 +1,3 @@
-from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
-from BaseTDEM import BaseTDEMProblem, FieldsTDEM
-from TDEM_b import ProblemTDEM_b
+from TDEM import BaseTDEMProblem, Problem_b, Problem_e
+from FieldsTDEM import Fields, Fields_b
+from SurveyTDEM import Survey, Src, Rx
@@ -0,0 +1,199 @@
+from SimPEG import Utils, Survey, np
+from SimPEG.Survey import BaseSurvey
+from SimPEG.EM.Utils import *
+from BaseTDEM import FieldsTDEM
+import SrcTDEM as Src
+
+class RxTDEM(Survey.BaseTimeRx):
+
+    knownRxTypes = {
+                    'ex':['e', 'Ex', 'N'],
+                    'ey':['e', 'Ey', 'N'],
+                    'ez':['e', 'Ez', 'N'],
+
+                    'bx':['b', 'Fx', 'N'],
+                    'by':['b', 'Fy', 'N'],
+                    'bz':['b', 'Fz', 'N'],
+
+                    'dbxdt':['b', 'Fx', 'CC'],
+                    'dbydt':['b', 'Fy', 'CC'],
+                    'dbzdt':['b', 'Fz', 'CC'],
+                   }
+
+    def __init__(self, locs, times, rxType):
+        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    @property
+    def projGLoc(self):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        return self.knownRxTypes[self.rxType][1]
+
+    @property
+    def projTLoc(self):
+        """Time Location projection (e.g. CC N)"""
+        return self.knownRxTypes[self.rxType][2]
+
+    def getTimeP(self, timeMesh):
+        """
+            Returns the time projection matrix.
+
+            .. note::
+
+                This is not stored in memory, but is created on demand.
+        """
+        if self.rxType in ['dbxdt','dbydt','dbzdt']:
+            return timeMesh.getInterpolationMat(self.times, self.projTLoc)*timeMesh.faceDiv
+        else:
+            return timeMesh.getInterpolationMat(self.times, self.projTLoc)
+
+    def eval(self, src, mesh, timeMesh, u):
+        P = self.getP(mesh, timeMesh)
+        u_part = Utils.mkvc(u[src, self.projField, :])
+        return P*u_part
+
+    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
+        P = self.getP(mesh, timeMesh)
+
+        if not adjoint:
+            return P * Utils.mkvc(v[src, self.projField, :])
+        elif adjoint:
+            return P.T * v[src, self]
+
+
+class SrcTDEM(Survey.BaseSrc):
+    rxPair = RxTDEM
+    radius = None
+
+    def getInitialFields(self, mesh):
+        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
+        return F0
+
+    def getJs(self, mesh, time):
+        return None
+
+
+class SrcTDEM_VMD_MVP(SrcTDEM):
+
+    def __init__(self,rxList,loc,waveformType="STEPOFF"):
+        self.loc = loc
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Vertical magnetic dipole, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+            else:
+                raise Exception('Unknown mesh for VMD')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+        else:
+            raise Exception('Unknown mesh for VMD')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SrcTDEM_CircularLoop_MVP(SrcTDEM):
+    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
+        self.loc = loc
+        self.radius = radius
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Circular Loop, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            else:
+                raise Exception('Unknown mesh for CircularLoop')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+        else:
+            raise Exception('Unknown mesh for CircularLoop')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SurveyTDEM(Survey.BaseSurvey):
+    """
+        docstring for SurveyTDEM
+    """
+    srcPair = SrcTDEM
+
+    def __init__(self, srcList, **kwargs):
+        # Sort these by frequency
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def projectFields(self, u):
+        data = Survey.Data(self)
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
+        return data
+
+    def projectFieldsDeriv(self, u, v=None, adjoint=False):
+        assert v is not None, 'v to multiply must be provided.'
+
+        if not adjoint:
+            data = Survey.Data(self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+            return data
+        else:
+            f = FieldsTDEM(self.mesh, self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+                    if rx.projField not in f: # first time we are projecting
+                        f[src, rx.projField, :] = Ptv
+                    else: # there are already fields, so let's add to them!
+                        f[src, rx.projField, :] += Ptv
+            return f
+
+
@@ -0,0 +1,3 @@
+from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
+from BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from TDEM_b import ProblemTDEM_b
@@ -42,10 +42,10 @@ def run(plotIt=True):


    rxOffset=1e-3
-    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
-    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
-    survey = EM.TDEM.SurveyTDEM([src])
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
+    survey = EM.TDEM.Survey([src])
+    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)

    prb.Solver = SolverLU
    prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
@@ -53,9 +53,9 @@ def run(plotIt=True):

    # create observed data
    std = 0.05
-    
+
    survey.dobs = survey.makeSyntheticData(mtrue,std)
-    survey.std = std 
+    survey.std = std
    survey.eps = 1e-5*np.linalg.norm(survey.dobs)

    if plotIt:
@@ -1,7 +1,7 @@
 from SimPEG import *


-def run(N=200, plotIt=True):
+def run(N=100, plotIt=True):
    """
        Inversion: Linear Problem
        =========================
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
    mesh = Mesh.TensorMesh([N])

    m0 = np.ones(mesh.nC) * 1e-4
+    mref = np.zeros(mesh.nC)
+
    nk = 10
    jk = np.linspace(1.,nk,nk)
    p = -2.
@@ -50,57 +52,47 @@ def run(N=200, plotIt=True):
    wr = np.sum(prob.G**2.,axis=0)**0.5
    wr = ( wr/np.max(wr) )

-    reg = Regularization.Simple(mesh)
-    reg.wght = wr
-
+#    reg = Regularization.Simple(mesh)
+#    reg.mref = mref
+#    reg.cell_weights = wr
+#
    dmis = DataMisfit.l2_DataMisfit(survey)
    dmis.Wd = 1./wd
-
-    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    invProb.curModel = m0
-
-    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-    target = Directives.TargetMisfit()
-
+#
+#    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
+#    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+#    invProb.curModel = m0
+#
+#    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
+#    target = Directives.TargetMisfit()
+#
    betaest = Directives.BetaEstimate_ByEig()
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
-
-
-    mrec = inv.run(m0)
-    ml2 = mrec
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-    # Switch regularization to sparse
-    phim = invProb.phi_m_last
-    phid =  invProb.phi_d
+#    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
+#
+#
+#    mrec = inv.run(m0)
+#    ml2 = mrec
+#    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
+#
+#    # Switch regularization to sparse
+#    phim = invProb.phi_m_last
+#    phid =  invProb.phi_d

    reg = Regularization.Sparse(mesh)
+    reg.mref = mref
+    reg.cell_weights = wr

-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
    reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
+    eps_p = 5e-2
+    eps_q = 5e-2
+    norms   = [0., 0., 2., 2.]

-    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
+    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+    update_Jacobi = Directives.Update_lin_PreCond()
+    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)

-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
-
-    m0 = mrec
+    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])

    # Run inversion
    mrec = inv.run(m0)
@@ -117,7 +109,7 @@ def run(N=200, plotIt=True):
        axes[0].set_title('Columns of matrix G')

        axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
+        axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
        #axes[1].legend(('True Model', 'Recovered Model'))
        axes[1].set_ylim(-1.0,1.25)

@@ -1,22 +1,25 @@
 from SimPEG import Mesh, Utils, np, SolverLU

-## 2D DC forward modeling example with Tensor and Curvilinear Meshes
-
 def run(plotIt=True):
+
+    """
+        Mesh: Basic Forward 2D DC Resistivity
+        =====================================
+
+        2D DC forward modeling example with Tensor and Curvilinear Meshes
+    """
+
    # Step1: Generate Tensor and Curvilinear Mesh
    sz = [40,40]
-    # Tensor Mesh
    tM = Mesh.TensorMesh(sz)
-    # Curvilinear Mesh
    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+
    # Step2: Direct Current (DC) operator
    def DCfun(mesh, pts):
        D = mesh.faceDiv
-        G = D.T
        sigma = 1e-2*np.ones(mesh.nC)
-        Msigi = mesh.getFaceInnerProduct(1./sigma)
-        MsigI = Utils.sdInv(Msigi)
-        A = D*MsigI*G
+        MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
+        A = -D*MsigI*D.T
        A[-1,-1] /= mesh.vol[-1] # Remove null space
        rhs = np.zeros(mesh.nC)
        txind = Utils.meshutils.closestPoints(mesh, pts)
@@ -37,39 +40,17 @@ def run(plotIt=True):
    if not plotIt: return

    import matplotlib.pyplot as plt
-    import matplotlib
-    from matplotlib.mlab import griddata

    #Step4: Making Figure
    fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-    label = ["(a)", "(b)"]
-    opts = {}
    vmin, vmax = phitM.min(), phitM.max()
    dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
-
-    #TODO: At the moment Curvilinear Mesh do not have plotimage
-
-    Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
-    Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
-    PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
-    axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
-
+    dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
    cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
    cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")

-    tM.plotGrid(ax=axes[0], **opts)
    axes[0].set_title('TensorMesh')
-    rM.plotGrid(ax=axes[1], **opts)
    axes[1].set_title('CurvilinearMesh')
-    for i in range(2):
-        axes[i].set_xlim(0.025, 0.975)
-        axes[i].set_ylim(0.025, 0.975)
-        axes[i].text(0., 1.0, label[i], fontsize=20)
-        if i==0:
-            axes[i].set_ylabel("y")
-        else:
-            axes[i].set_ylabel(" ")
-        axes[i].set_xlabel("x")
    plt.show()


@@ -8,9 +8,9 @@ import EM_FDEM_Analytic_MagDipoleWholespace
 import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
-import Forward_BasicDirectCurrent
 import Inversion_IRLS
 import Inversion_Linear
+import Mesh_Basic_ForwardDC
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
 import Mesh_Operators_CahnHilliard
@@ -22,7 +22,7 @@ import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 import Utils_surface2ind_topo

-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]

 ##### AUTOIMPORTS #####

@@ -0,0 +1,157 @@
+from SimPEG import np, Mesh, Maps, Utils, DataMisfit, Regularization, Optimization, Inversion, InvProblem, Directives
+from SimPEG import SolverLU
+from SimPEG.EM import FDEM, TDEM, mu_0
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['font.size'] = 14
+
+def run(plotIt=True):
+    # Set up cylindrically symmeric mesh
+    cs, ncx, ncz, npad = 10., 15, 25, 13  # padded cyl mesh
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+
+    # Conductivity model
+    layerz = np.r_[-200., -100.]
+    layer = (mesh.vectorCCz>=layerz[0]) & (mesh.vectorCCz<=layerz[1])
+    active = mesh.vectorCCz<0.
+    sig_half  = 1e-2  # Half-space conductivity
+    sig_air   = 1e-8  # Air conductivity
+    sig_layer = 5e-2  # Layer conductivity
+    sigma = np.ones(mesh.nCz)*sig_air
+    sigma[active] = sig_half
+    sigma[layer] = sig_layer
+
+    # Mapping
+    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
+    mtrue = np.log(sigma[active])
+
+    # FDEM problem & survey
+    rxlocs = Utils.ndgrid([np.r_[50.], np.r_[0], np.r_[0.]])
+    bzi = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'real')
+    bzr = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'imag')
+
+    freqs = np.logspace(2, 3, 5)
+    srcLoc = np.array([0., 0., 0.])
+
+    print 'min skin depth = ', 500./np.sqrt(freqs.max() * sig_half), 'max skin depth = ', 500./np.sqrt(freqs.min() * sig_half)
+    print 'max x ', mesh.vectorCCx.max(), 'min z ', mesh.vectorCCz.min(), 'max z ', mesh.vectorCCz.max()
+
+    srcList = []
+    [srcList.append(FDEM.Src.MagDipole([bzr, bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+
+    surveyFD = FDEM.Survey(srcList)
+    prbFD = FDEM.Problem3D_b(mesh, mapping=mapping)
+    prbFD.pair(surveyFD)
+    std = 0.03
+    surveyFD.makeSyntheticData(mtrue, std)
+    surveyFD.eps = np.linalg.norm(surveyFD.dtrue)*1e-5
+
+    # FDEM inversion
+    np.random.seed(1)
+    dmisfit = DataMisfit.l2_DataMisfit(surveyFD)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Simple(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+    # Inversion Directives
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
+    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=10.)
+    invProb.beta = 1.
+    target = Directives.TargetMisfit()
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,target])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 5e-1
+    reg.alpha_x = 1.
+    prbFD.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.tolG = 1e-10
+    opt.eps = 1e-10
+    opt.remember('xc')
+    moptFD = inv.run(m0)
+
+    # TDEM problem
+    times = np.logspace(-4, np.log10(2e-3), 10)
+    print 'min diffusion distance ', 1.28*np.sqrt(times.min()/(sig_half*mu_0)), 'max diffusion distance ', 1.28*np.sqrt(times.max()/(sig_half*mu_0))
+    rx = TDEM.Rx(rxlocs, times, 'bz')
+    src = TDEM.Src.MagDipole([rx], waveform=TDEM.Src.StepOffWaveform(), loc=srcLoc) # same src location as FDEM problem
+
+    surveyTD = TDEM.Survey([src])
+    prbTD = TDEM.Problem_b(mesh, mapping=mapping)
+    prbTD.timeSteps = [(5e-5, 10),(1e-4, 10),(5e-4, 10)]
+    prbTD.pair(surveyTD)
+    prbTD.Solver = SolverLU
+
+    std = 0.03
+    surveyTD.makeSyntheticData(mtrue, std)
+    surveyTD.std = std
+    surveyTD.eps = np.linalg.norm(surveyTD.dtrue)*1e-5
+
+    # TDEM inversion
+    dmisfit = DataMisfit.l2_DataMisfit(surveyTD)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Simple(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
+    # Inversion Directives
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
+    invProb.beta = 1.
+    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1.)
+    target = Directives.TargetMisfit()
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 5e-1
+    reg.alpha_x = 1.
+    prbTD.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.remember('xc')
+    moptTD = inv.run(m0)
+
+    if plotIt:
+        fig, ax = plt.subplots(1,1, figsize = (4, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active], 'k-', lw=2)
+        plt.semilogx(np.exp(moptFD), mesh.vectorCCz[active], 'ko', ms=3)
+        plt.semilogx(np.exp(moptTD), mesh.vectorCCz[active], 'k*')
+        ax.set_ylim(-1000, 0)
+        ax.set_xlim(5e-3, 1e-1)
+
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+        plt.legend(['True', 'Pred (FD)', 'Pred (TD)'], fontsize=13, loc=4)
+        plt.show()
+
+        fig = plt.figure(figsize = (10*1.3, 5*1.3))
+        ax2 = plt.subplot(122)
+        ax2.plot(times, surveyTD.dobs, 'k-', lw=2)
+        ax2.plot(times, surveyTD.dpred(moptTD), 'ko', ms=4)
+        ax2.set_xscale('log')
+        ax2.set_yscale('log')
+        ax2.set_xlim(times.min(), times.max())
+        ax1 = plt.subplot(121)
+        ax1.plot(freqs, -surveyFD.dobs[::2], 'k-', lw=2)
+        ax1.plot(freqs, -surveyFD.dobs[1::2], 'k--', lw=2)
+        dpredFD = surveyFD.dpred(moptTD)
+        ax1.plot(freqs, -dpredFD[::2], 'ko', ms=4)
+        ax1.plot(freqs, -dpredFD[1::2], 'k+', markeredgewidth=2., ms=10)
+        ax1.set_xscale('log')
+        ax1.set_yscale('log')
+        ax2.set_xlabel('Time (s)', fontsize = 14)
+        ax1.set_xlabel('Frequency (Hz)', fontsize = 14)
+        ax1.set_ylabel('Vertical magnetic field (T)', fontsize = 14)
+        ax2.grid(True,which='minor')
+        ax1.grid(True,which='minor')
+        ax2.set_title("(b) TD observed vs. predicted", fontsize = 14)
+        ax1.set_title("(a) FD observed vs. predicted", fontsize = 14)
+        ax2.legend(("Obs", "Pred"), fontsize = 12)
+        ax1.legend(("Obs", "Pred (real)", "Pred (imag)"), fontsize = 12, loc=3)
+        ax1.set_xlim(freqs.max(), freqs.min())
+        plt.show()
+
+if __name__ == '__main__':
+    run()
@@ -2,6 +2,7 @@ from SimPEG import Utils, np
 from BaseMesh import BaseRectangularMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
+from View import CurvView

 # Some helper functions.
 length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
@@ -10,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))


-class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
    """
    CurvilinearMesh is a mesh class that deals with curvilinear meshes.

@@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):



-    #############################################
-    #            Plotting Functions             #
-    #############################################
-
-    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
-        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
-
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.CurvilinearMesh([X, Y])
-            M.plotGrid(showIt=True)
-
-        """
-        import matplotlib.pyplot as plt
-        import matplotlib
-        from mpl_toolkits.mplot3d import Axes3D
-        mkvc = Utils.mkvc
-
-        axOpts = {'projection':'3d'} if self.dim == 3 else {}
-        if ax is None: ax = plt.subplot(111, **axOpts)
-
-        NN = self.r(self.gridN, 'N', 'N', 'M')
-        if self.dim == 2:
-
-            if lines:
-                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
-
-                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
-
-                X = np.r_[X1, X2]
-                Y = np.r_[Y1, Y2]
-
-                ax.plot(X, Y, 'b-')
-            if centers:
-                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
-
-            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')
-
-            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
-
-            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            # ax.plot(nX, nY, 'r-')
-
-            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            # ax.plot(nX, nY, 'g-')
-
-            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            # ax.plot(tX, tY, 'r-')
-
-            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            # ax.plot(nX, nY, 'g-')
-
-        elif self.dim == 3:
-            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
-            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
-            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
-
-            X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
-            Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
-            Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
-
-            X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
-            Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
-            Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
-
-            X = np.r_[X1, X2, X3]
-            Y = np.r_[Y1, Y2, Y3]
-            Z = np.r_[Z1, Z2, Z3]
-
-            ax.plot(X, Y, 'b', zs=Z)
-            ax.set_zlabel('x3')
-
-        ax.grid(True)
-        ax.set_xlabel('x1')
-        ax.set_ylabel('x2')
-
-        if showIt: plt.show()
-
-
 if __name__ == '__main__':
    nc = 5
    h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
@@ -552,7 +552,8 @@ class CurvView(object):
    def __init__(self):
        pass

-    def plotGrid(self, length=0.05, showIt=False):
+
+    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.


@@ -560,60 +561,63 @@ class CurvView(object):
            :include-source:

            from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleCurvGird([3,3],'rotate')
+            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
            M = Mesh.CurvilinearMesh([X, Y])
            M.plotGrid(showIt=True)

        """
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+
+        axOpts = {'projection':'3d'} if self.dim == 3 else {}
+        if ax is None: ax = plt.subplot(111, **axOpts)
+
        NN = self.r(self.gridN, 'N', 'N', 'M')
        if self.dim == 2:
-            fig = plt.figure(2)
-            fig.clf()
-            ax = plt.subplot(111)
-            X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-            Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()

-            X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-            Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
+            if lines:
+                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
+                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()

-            X = np.r_[X1, X2]
-            Y = np.r_[Y1, Y2]
+                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
+                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()

-            plt.plot(X, Y)
+                X = np.r_[X1, X2]
+                Y = np.r_[Y1, Y2]

-            plt.hold(True)
-            Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            Ty = self.r(self.tangents, 'E', 'Ey', 'V')
+                ax.plot(X, Y, 'b-')
+            if centers:
+                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')

-            plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
+            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
+            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
+            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
+            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')

-            nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            plt.plot(nX, nY, 'r-')
+            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')

-            nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            plt.plot(nX, nY, 'g-')
+            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
+            # ax.plot(nX, nY, 'r-')

-            tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            plt.plot(tX, tY, 'r-')
+            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
+            # ax.plot(nX, nY, 'g-')

-            nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            plt.plot(nX, nY, 'g-')
-            plt.axis('equal')
+            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
+            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
+            # ax.plot(tX, tY, 'r-')
+
+            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
+            # ax.plot(nX, nY, 'g-')

        elif self.dim == 3:
-            fig = plt.figure(3)
-            fig.clf()
-            ax = fig.add_subplot(111, projection='3d')
            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
@@ -630,16 +634,50 @@ class CurvView(object):
            Y = np.r_[Y1, Y2, Y3]
            Z = np.r_[Z1, Z2, Z3]

-            plt.plot(X, Y, 'b', zs=Z)
+            ax.plot(X, Y, 'b', zs=Z)
            ax.set_zlabel('x3')

        ax.grid(True)
-        ax.hold(False)
        ax.set_xlabel('x1')
        ax.set_ylabel('x2')

        if showIt: plt.show()

+    def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
+        if self.dim == 3: raise NotImplementedError('This is not yet done!')
+
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        import matplotlib.colors as colors
+        import matplotlib.cm as cmx
+
+        if ax is None: ax = plt.subplot(111)
+        jet = cm = plt.get_cmap('jet')
+        cNorm  = colors.Normalize(
+            vmin=I.min() if clim is None else clim[0],
+            vmax=I.max() if clim is None else clim[1])
+
+        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
+        # ax.set_xlim((self.x0[0], self.h[0].sum()))
+        # ax.set_ylim((self.x0[1], self.h[1].sum()))
+
+        Nx = self.r(self.gridN[:,0],'N','N','M')
+        Ny = self.r(self.gridN[:,1],'N','N','M')
+        cell = self.r(I,'CC','CC','M')
+
+        for ii in range(self.nCx):
+            for jj in range(self.nCy):
+                I = [ii,ii+1,ii+1,ii]
+                J = [jj,jj,jj+1,jj+1]
+                ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none'))
+
+        scalarMap._A = []  # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        if showIt: plt.show()
+        return [scalarMap]
+

 if __name__ == '__main__':
    from SimPEG import *
@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
            delx[indx] = 0.

-            return delx
+            return delx
@@ -1,4 +1,6 @@
-import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
+import Utils, Maps, Mesh
+import numpy as np
+import scipy.sparse as sp

 class RegularizationMesh(object):
    """
@@ -403,7 +405,238 @@ class BaseRegularization(object):

        return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )

-class Tikhonov(BaseRegularization):
+class Simple(BaseRegularization):
+    """
+    Simple regularization that does not include length scales in the derivatives.
+    """
+
+    mrefInSmooth = False  #: include mref in the smoothness?
+    alpha_s      = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
+    alpha_x      = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
+    alpha_y      = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
+    alpha_z      = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
+    cell_weights = 1.
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
+
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self,'_Wsmall', None) is None:
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
+        return self._Wsmall
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
+        return self._Wz
+
+#    @property
+#    def Wsmooth(self):
+#        """Full smoothness regularization matrix W"""
+#        print 'wtf why are we using Wsmooth'
+#        raise NotImplementedError
+#        if getattr(self, '_Wsmooth', None) is None:
+#            wlist = (self.Wx,)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._Wsmooth = sp.vstack(wlist)
+#        return self._Wsmooth
+#
+#    @property
+#    def W(self):
+#        """Full regularization matrix W"""
+#        print 'wtf why are we using W'
+#        if getattr(self, '_W', None) is None:
+#            wlist = (self.Wsmall, self.Wx)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._W = sp.vstack(wlist)
+#        return self._W
+
+
+    @Utils.timeIt
+    def _evalSmall(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmallDeriv(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
+
+    @Utils.timeIt
+    def _evalSmall2Deriv(self, m, v = None):
+        rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothx(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth(self, m):
+        phiSmooth = self._evalSmoothx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth += self._evalSmoothy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth += self._evalSmoothz(m)
+        return phiSmooth
+
+    @Utils.timeIt
+    def _evalSmoothxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * m )
+            return r.T * ( self.Wx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wx * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * m )
+            return r.T * ( self.Wy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wy * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * m )
+            return r.T * ( self.Wz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wz * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv(self, m):
+        deriv = self._evalSmoothxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv(self, m, v=None):
+        deriv = self._evalSmoothx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothz2Deriv(m, v)
+        return deriv
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        return self._evalSmall(m) + self._evalSmooth(m)
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+
+    @Utils.timeIt
+    def eval2Deriv(self, m, v=None):
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
+
+
+
+class Tikhonov(Simple):
    """
    L2 Tikhonov regularization with both smallness and smoothness (first order
    derivative) contributions.
@@ -493,56 +726,131 @@ class Tikhonov(BaseRegularization):
            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
        return self._Wzz

+
    @property
-    def Wsmooth(self):
+    def Wsmooth2(self):
        """Full smoothness regularization matrix W"""
        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
+            wlist = (self.Wxx)
            if self.regmesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
+                wlist += (self.Wyy)
            if self.regmesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
+                wlist += (self.Wzz)
            self._Wsmooth = sp.vstack(wlist)
        return self._Wsmooth

-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
+    def _evalSmoothxx(self, m):
        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+            r = self.Wxx * ( self.mapping * (m - self.mref) )
        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * (m) )
+            r = self.Wxx * ( self.mapping * (m) )
        return 0.5 * r.dot(r)

+    @Utils.timeIt
+    def _evalSmoothyy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothzz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth2(self, m):
+        phiSmooth2 = self._evalSmoothxx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth2 += self._evalSmoothyy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth2 += self._evalSmoothzz(m)
+        return phiSmooth2
+
+    @Utils.timeIt
+    def _evalSmoothxxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wxx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wxx * ( self.mapping * m )
+            return r.T * ( self.Wxx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothyyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * m )
+            return r.T * ( self.Wyy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothzzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * m )
+            return r.T * ( self.Wzz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothxx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wxx * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wyy * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wzz * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv2(self, m):
+        deriv = self._evalSmoothxxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv2(self, m, v=None):
+        deriv = self._evalSmoothxx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzz2Deriv(m, v)
+        return deriv
+
+
    @Utils.timeIt
    def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
-
-    @Utils.timeIt
-    def _evalSmallDeriv(self,m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @Utils.timeIt
-    def _evalSmoothDeriv(self,m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
+        return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)

    @Utils.timeIt
    def evalDeriv(self, m):
@@ -560,184 +868,134 @@ class Tikhonov(BaseRegularization):
            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}

        """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
+
+    def eval2Deriv(self, m, v=None):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)


-class Simple(Tikhonov):
+
+class Sparse(Simple):
    """
-    Simple regularization that does not include length scales in the derivatives.
+        The regularization is:
+    
+        .. math::
+    
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
+    
+        where the IRLS weight
+    
+        .. math::
+    
+            R = \eta TO FINISH LATER!!!
+    
+        So the derivative is straight forward:
+    
+        .. math::
+    
+            R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
+    
+        The IRLS weights are recomputed after each beta solves.
+        It is strongly recommended to do a few Gauss-Newton iterations
+        before updating.
    """
-
-    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-    wght        = 1.
+        
+    # set default values
+    eps_p = 1e-1        # Threshold value for the model norm
+    eps_q = 1e-1        # Threshold value for the model gradient norm
+    curModel = None     # Requires model to compute the weights
+    l2model = None
+    gamma = 1.          # Model norm scaling to smooth out convergence
+    norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
+    cell_weights = 1.        # Consider overwriting with sensitivity weights

    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)

-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights

    @property
    def Wsmall(self):
        """Regularization matrix Wsmall"""
        if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+            if getattr(self, 'curModel', None) is None:
+                self.Rs = Utils.speye(self.regmesh.nC)
+
+            else:
+                f_m = self.mapping * (self.curModel - self.reg.mref)
+                self.rs = self.R(f_m , self.eps_p, self.norms[0])
+                self.Rs = Utils.sdiag( self.rs )
+
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
+
        return self._Wsmall

    @property
    def Wx(self):
        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
+        if getattr(self,'_Wx', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
+                self.rx = self.R( f_m , self.eps_q, self.norms[1])
+                self.Rx = Utils.sdiag( self.rx )
+
+            self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
+
        return self._Wx

    @property
    def Wy(self):
        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
+        if getattr(self,'_Wy', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
+                self.ry = self.R( f_m , self.eps_q, self.norms[2])
+                self.Ry = Utils.sdiag( self.ry )
+
+            self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
+
        return self._Wy

    @property
    def Wz(self):
        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
+        if getattr(self,'_Wz', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
+                self.rz = self.R( f_m , self.eps_q, self.norms[3])
+                self.Rz = Utils.sdiag( self.rz )
+
+            self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
+
        return self._Wz

-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx,)
-            if self.regmesh.dim > 1:
-                wlist += (self.Wy,)
-            if self.regmesh.dim > 2:
-                wlist += (self.Wz,)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m)
-        return 0.5 * r.dot(r)
-
-
-class Sparse(Simple):
-
-    # set default values
-    eps_p      = 1e-1
-    eps_q      = 1e-1
-    curModel = None # use a model to compute the weights
-    gamma    = 1.
-    norms    = [0., 2., 2., 2.]
-    wght     = 1.
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
-
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
-
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self, 'curModel', None) is None:
-            self.Rs = Utils.speye(self.regmesh.nC)
-
-        else:
-            f_m = self.curModel - self.reg.mref
-            self.rs = self.R(f_m , self.eps_p, self.norms[0])
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
-
-        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
-
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
-            self.rx = self.R( f_m , self.eps_q, self.norms[1])
-            self.Rx = Utils.sdiag( self.rx )
-
-        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.eps_q, self.norms[2])
-            self.Ry = Utils.sdiag( self.ry )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.eps_q, self.norms[3])
-            self.Rz = Utils.sdiag( self.rz )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
-
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        #if getattr(self, '_Wsmooth', None) is None:
-        wlist = (self.Wx,)
-        if self.regmesh.dim > 1:
-            wlist += (self.Wy,)
-        if self.regmesh.dim > 2:
-            wlist += (self.Wz,)
-        #self._Wsmooth = sp.vstack(wlist)
-        return sp.vstack(wlist)
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
    def R(self, f_m , eps, exponent):

+        # Eta scaling is important for mix-norms...do not mess with it
        eta = (eps**(1.-exponent/2.))**0.5
        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)

@@ -1,4 +1,5 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
+import gc

 class BaseRx(object):
    """SimPEG Receiver Object"""
@@ -347,10 +347,10 @@ and



-TDEM - B formulation
-====================
+TDEM Problem
+============

-.. automodule:: SimPEG.EM.TDEM.TDEM_b
+.. automodule:: SimPEG.EM.TDEM.TDEM
    :show-inheritance:
    :members:
    :undoc-members:
@@ -359,7 +359,7 @@ TDEM - B formulation
 Field Storage
 =============

-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Fields
    :show-inheritance:
    :members:
    :undoc-members:
@@ -369,19 +369,19 @@ Field Storage
 TDEM Survey Classes
 ===================

-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.SurveyTDEM
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Survey
    :show-inheritance:
    :members:
    :undoc-members:
    :inherited-members:


-Base Classes
-============
+.. Base Classes
+.. ============

-.. automodule:: SimPEG.EM.TDEM.BaseTDEM
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:
+.. .. automodule:: SimPEG.EM.TDEM.BaseTDEM
+..     :show-inheritance:
+..     :members:
+..     :undoc-members:
+..     :inherited-members:

@@ -22,7 +22,6 @@ radi    = Radius of spheres [r1,r2]
 param   = Conductivity of background and two spheres [m0,m1,m2]
 surveyType = survey type 'pole-dipole' or 'dipole-dipole'
 unitType = Data type "appResistivity" | "appConductivity"  | "volt"
-
 Created by @fourndo


@@ -1,4 +1,4 @@
-.. _examples_Forward_BasicDirectCurrent:
+.. _examples_Mesh_Basic_ForwardDC:

 .. --------------------------------- ..
 ..                                   ..
@@ -8,14 +8,18 @@
 ..                                   ..
 .. --------------------------------- ..

-Forward BasicDirectCurrent
-==========================
+
+Mesh: Basic Forward 2D DC Resistivity
+=====================================
+
+2D DC forward modeling example with Tensor and Curvilinear Meshes
+

 .. plot::

    from SimPEG import Examples
-    Examples.Forward_BasicDirectCurrent.run()
+    Examples.Mesh_Basic_ForwardDC.run()

-.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_ForwardDC.py
    :language: python
    :linenos:
@@ -3,310 +3,218 @@ from SimPEG import *
 from SimPEG import EM

 plotIt = False
-tol = 1e-6

-class TDEM_bDerivTests(unittest.TestCase):
+testDeriv   = True
+testAdjoint = True

-    def setUp(self):
+TOL = 1e-5

-        cs = 5.
-        ncx = 20
-        ncy = 6
-        npad = 20
-        hx = [(cs,ncx), (cs,npad,1.3)]
-        hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-        mesh = Mesh.CylMesh([hx,1,hy], '00C')
+def setUp(prbtype='b', rxcomp='bz'):
+    cs = 5.
+    ncx = 20
+    ncy = 15
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+#
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap

-        active = mesh.vectorCCz<0.
-        activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-        mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+    rxOffset = 10.
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))

-        rxOffset = 40.
-        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
-        src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
+    survey = EM.TDEM.Survey([src])

-        survey = EM.TDEM.SurveyTDEM([src])
+    if prbtype == 'b':
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    elif prbtype == 'e':
+        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)

-        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-        # self.prb.timeSteps = [1e-5]
-        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # self.prb.timeSteps = [(1e-05, 100)]
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]

-        try:
-            from pymatsolver import MumpsSolver
-            self.prb.Solver = MumpsSolver
-        except ImportError, e:
-            self.prb.Solver = SolverLU
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU

-        self.sigma = np.ones(mesh.nCz)*1e-8
-        self.sigma[mesh.vectorCCz<0] = 1e-1
-        self.sigma = np.log(self.sigma[active])
+    m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)

-        self.prb.pair(survey)
-        self.mesh = mesh
+    prb.pair(survey)
+    mesh = mesh

-    def test_AhVec(self):
-        """
-            Test that fields and AhVec produce consistent results
-        """
+    return prb, m, mesh

-        prb = self.prb
-        sigma = self.sigma

-        u = prb.fields(sigma)
-        Ahu = prb._AhVec(sigma, u)
+class TDEM_DerivTests(unittest.TestCase):

-        V1 = Ahu[:,'b',1]
-        V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
-        self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
+# ====== TEST A ========== #

-        V1 = Ahu[:,'e',1]
-        return np.linalg.norm(V1) < 1.e-6
+    def AderivTest(self, prbtype):
+        prb, m0, mesh = setUp(prbtype)
+        tInd = 2
+        if prbtype == 'b':
+            nu = mesh.nF
+        elif prbtype == 'e':
+            nu = mesh.nE
+        v = np.random.rand(nu)

-        for i in range(2,prb.nT):
+        def AderivFun(m):
+            prb.curModel = m
+            A = prb.getAdiag(tInd)
+            Av = A*v
+            prb.curModel = m0
+            ADeriv_dm = lambda dm: prb.getAdiagDeriv(tInd, v, dm)

-            dt = prb.timeSteps[i]
+            return Av, ADeriv_dm

-            V1 = Ahu[:,'b',i]
-            V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
-            # print np.linalg.norm(V1), np.linalg.norm(V2)
-            self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
+        print '\n Testing ADeriv %s'%(prbtype)
+        Tests.checkDerivative(AderivFun, m0, plotIt=False, num=4, eps=1e-20)

-            V1 = Ahu[:,'e',i]
-            V2 = prb.MeSigma*u[:,'e',i]
-            # print np.linalg.norm(V1), np.linalg.norm(V2)
-            return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
-
-    def test_AhVecVSMat_OneTS(self):
-
-        prb = self.prb
-        prb.timeSteps = [1e-05]
-        sigma = self.sigma
-        prb.curModel = sigma
-
-        dt = prb.timeSteps[0]
-        a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF)
-        a12 = prb.MfMui*prb.mesh.edgeCurl
-        a21 = prb.mesh.edgeCurl.T*prb.MfMui
-        a22 = -prb.MeSigma
-        A = sp.bmat([[a11,a12],[a21,a22]])
-
-        f = prb.fields(sigma)
-        u1 = A*f.tovec()
-        u2 = prb._AhVec(sigma,f).tovec()
-
-        self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
-
-    def test_solveAhVSMat_OneTS(self):
-        prb = self.prb
-
-        prb.timeSteps = [1e-05]
-
-        sigma = self.sigma
-        prb.curModel = sigma
-
-        dt = prb.timeSteps[0]
-        a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF)
-        a12 = prb.MfMui*prb.mesh.edgeCurl
-        a21 = prb.mesh.edgeCurl.T*prb.MfMui
-        a22 = -prb.MeSigma
-        A = sp.bmat([[a11,a12],[a21,a22]])
-
-        f = prb.fields(sigma)
-        f[:,:,0] = {'b':0}
-        f[:,'b',1] = 0
-
-        self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
-
-        u1 = prb.solveAh(sigma,f).tovec().flatten()
-        u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
-
-        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
-
-    def test_solveAhVsAhVec(self):
-
-        prb = self.prb
-        mesh = self.prb.mesh
-        sigma = self.sigma
-        self.prb.curModel = sigma
-
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        f[:,'b',:] = 0.0
-        for i in range(prb.nT):
-            f[:,'e', i] = np.random.rand(mesh.nE, 1)
-
-        Ahf = prb._AhVec(sigma, f)
-        f_test = prb.solveAh(sigma, Ahf)
-
-        u1 = f.tovec()
-        u2 = f_test.tovec()
-        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
-
-    def test_DerivG(self):
-        """
-            Test the derivative of c with respect to sigma
-        """
-
-        # Random model and perturbation
-        sigma = np.random.rand(self.prb.mapping.nP)
-
-        f = self.prb.fields(sigma)
-        dm = 1000*np.random.rand(self.prb.mapping.nP)
-        h = 0.01
-
-        derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
-        print '\ntest_DerivG'
-        passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-        return passed
-
-    def test_Deriv_dUdM(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        dm = 10*np.random.rand(prb.mapping.nP)
-        f = prb.fields(sigma)
-
-        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
-        print '\n'
-        print 'test_Deriv_dUdM'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_J(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-        d_sig = 10*np.random.rand(prb.mapping.nP)
-
-
-        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-        print '\n'
-        print 'test_Deriv_J'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
-
-    def test_projectAdjoint(self):
-        prb = self.prb
-        survey = prb.survey
-        mesh = self.mesh
-
-        # Generate random fields and data
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(prb.nT):
-            f[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f[:,'e',i] = np.random.rand(mesh.nE, 1)
-        d_vec = np.random.rand(survey.nD)
-        d = Survey.Data(survey,v=d_vec)
-
-        # Check that d.T*Q*f = f.T*Q.T*d
-        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
-        V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())
-
-        self.assertTrue((V1-V2)/np.abs(V1) < tol)
-
-    def test_adjointAhVsAht(self):
-        prb = self.prb
-        mesh = self.mesh
-        sigma = self.sigma
-
-        f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            f2[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f2[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
-        V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
-
-    # def test_solveAhtVsAhtVec(self):
-    #     prb = self.prb
-    #     mesh = self.mesh
-    #     sigma = np.random.rand(prb.mapping.nP)
-
-    #     f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     f2 = prb.solveAht(sigma, f1)
-    #     f3 = prb._AhtVec(sigma, f2)
-
-    #     if True:
-    #         import matplotlib.pyplot as plt
-    #         plt.plot(f3.tovec(),'b')
-    #         plt.plot(f1.tovec(),'r')
-    #         plt.show()
-    #     V1 = np.linalg.norm(f3.tovec()-f1.tovec())
-    #     V2 = np.linalg.norm(f1.tovec())
-    #     print 'AhtVsAhtVec', V1, V2, f1.tovec()
-    #     print 'I am gunna fail this one: boo. :('
-    #     self.assertLess(V1/V2, 1e-6)
-
-    # def test_adjointsolveAhVssolveAht(self):
-    #     prb = self.prb
-    #     mesh = self.mesh
-    #     sigma = self.sigma
-
-    #     f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f2[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f2[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
-    #     V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
-    #     print V1, V2
-    #     self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
-
-    def test_adjointGvecVsGtvec(self):
-        mesh = self.mesh
-        prb = self.prb
+    def A_adjointTest(self,prbtype):
+        prb, m0, mesh = setUp(prbtype)
+        tInd = 2

+        print '\n Testing A_adjoint'
        m = np.random.rand(prb.mapping.nP)
-        sigma = np.random.rand(prb.mapping.nP)
+        if prbtype == 'b':
+            nu = prb.mesh.nF
+        elif prbtype == 'e':
+            nu = prb.mesh.nE

-        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            u[:,'b',i] = np.random.rand(mesh.nF, 1)
-            u[:,'e',i] = np.random.rand(mesh.nE, 1)
+        v = np.random.rand(nu)
+        u = np.random.rand(nu)
+        prb.curModel = m0

-        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            v[:,'b',i] = np.random.rand(mesh.nF, 1)
-            v[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        V1 = m.dot(prb.Gtvec(sigma, v, u))
-        V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
-
-    def test_adjointJvecVsJtvec(self):
-        mesh = self.mesh
-        prb = self.prb
-        sigma = self.sigma
-
-        m = np.random.rand(prb.mapping.nP)
-        d = np.random.rand(prb.survey.nD)
-
-        V1 = d.dot(prb.Jvec(sigma, m))
-        V2 = m.dot(prb.Jtvec(sigma, d))
-        passed = np.abs(V1-V2)/np.abs(V1) < tol
-        print 'AdjointTest', V1, V2, passed
+        tInd = 2 # not actually used
+        V1 = v.dot(prb.getAdiagDeriv(tInd, u, m))
+        V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True))
+        passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2.
+        print 'AdjointTest %s'%(prbtype), V1, V2, passed
        self.assertTrue(passed)

+    def test_Aderiv_b(self):
+        self.AderivTest('b')
+    def test_Aderiv_e(self):
+        self.AderivTest('e')
+
+    def test_Aadjoint_b(self):
+        self.A_adjointTest('b')
+    def test_Aadjoint_e(self):
+        self.A_adjointTest('e')
+
+# ====== TEST Fields Deriv Pieces ========== #
+
+    def test_eDeriv_m_adjoint(self):
+        prb, m0, mesh = setUp()
+        tInd = 0
+
+        v = np.random.rand(mesh.nF)
+
+        print '\n Testing eDeriv_m Adjoint'
+
+        prb, m0, mesh = setUp()
+        f = prb.fields(m0)
+
+        m = np.random.rand(prb.mapping.nP)
+        e = np.random.randn(prb.mesh.nE)
+        V1 = e.dot(f._eDeriv_m(1, prb.survey.srcList[0], m))
+        V2 = m.dot(f._eDeriv_m(1, prb.survey.srcList[0], e, adjoint=True))
+        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+        passed = np.abs(V1-V2) < tol
+
+        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        self.assertTrue(passed)
+
+    def test_eDeriv_u_adjoint(self):
+        print '\n Testing eDeriv_u Adjoint'
+
+        prb, m0, mesh = setUp()
+        f = prb.fields(m0)
+
+        b = np.random.rand(prb.mesh.nF)
+        e = np.random.randn(prb.mesh.nE)
+        V1 = e.dot(f._eDeriv_u(1, prb.survey.srcList[0], b))
+        V2 = b.dot(f._eDeriv_u(1, prb.survey.srcList[0], e, adjoint=True))
+        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+        passed = np.abs(V1-V2) < tol
+
+        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        self.assertTrue(passed)
+
+
+# ====== TEST Jvec ========== #
+
+    if testDeriv:
+
+        def JvecTest(self, prbtype, rxcomp):
+            prb, m, mesh = setUp(prbtype, rxcomp)
+
+            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
+            print '\n'
+            print 'test_Jvec_%s_%s' %(prbtype, rxcomp)
+            Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
+
+        def test_Jvec_b_bx(self):
+            self.JvecTest('b','bx')
+
+        def test_Jvec_b_bz(self):
+            self.JvecTest('b','bz')
+
+        def test_Jvec_b_dbxdt(self):
+            self.JvecTest('b','dbxdt')
+
+        def test_Jvec_b_dbzdt(self):
+            self.JvecTest('b','dbzdt')
+
+        def test_Jvec_b_ey(self):
+            self.JvecTest('b','ey')
+
+        def test_Jvec_e_ey(self):
+            self.JvecTest('e','ey')
+
+
+# ====== TEST Jtvec ========== #
+
+    if testAdjoint:
+
+        def JvecVsJtvecTest(self, prbtype='b', rxcomp='bz'):
+
+            print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp)
+
+            prb, m0, mesh = setUp(prbtype, rxcomp)
+            m = np.random.rand(prb.mapping.nP)
+            d = np.random.randn(prb.survey.nD)
+            V1 = d.dot(prb.Jvec(m0, m))
+            V2 = m.dot(prb.Jtvec(m0, d))
+            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+            passed = np.abs(V1-V2) < tol
+
+            print '    ', V1, V2, np.abs(V1-V2), tol, passed
+            self.assertTrue(passed)
+
+        def test_Jvec_adjoint_b_bx(self):
+            self.JvecVsJtvecTest('b', 'bx')
+
+        def test_Jvec_adjoint_b_bz(self):
+            self.JvecVsJtvecTest('b', 'bz')
+
+        def test_Jvec_adjoint_b_dbxdt(self):
+            self.JvecVsJtvecTest('b', 'bx')
+
+        def test_Jvec_adjoint_b_dbzdt(self):
+            self.JvecVsJtvecTest('b', 'bz')
+
+        def test_Jvec_adjoint_b_ey(self):
+            self.JvecVsJtvecTest('b', 'ey')
+
+        # This is not working because Problem_e has not done
+        # def test_Jvec_adjoint_e_ey(self):
+        #     self.JvecVsJtvecTest('e', 'ey')



@@ -3,10 +3,12 @@ from SimPEG import *
 from SimPEG import EM

 plotIt = False
+testDeriv = True
+testAdjoint = True

-class TDEM_bDerivTests(unittest.TestCase):
+TOL = 1e-5

-    def setUp(self):
+def setUp(self, rxcomp='bz'):

        cs = 5.
        ncx = 20
@@ -21,131 +23,78 @@ class TDEM_bDerivTests(unittest.TestCase):
        mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap

        rxOffset = 40.
-        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
-        src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
-        rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
-        src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
+        rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxcomp)
+        src = EM.TDEM.Src.MagDipole( [rx], loc=np.array([0., 0., 0.]))
+        rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), rxcomp)
+        src2 = EM.TDEM.Src.MagDipole( [rx2], loc=np.array([0., 0., 0.]))

-        survey = EM.TDEM.SurveyTDEM([src,src2])
+        survey = EM.TDEM.Survey([src,src2])

-        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-        # self.prb.timeSteps = [1e-5]
-        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # self.prb.timeSteps = [(1e-05, 100)]
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+        # prb.timeSteps = [1e-5]
+        prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+        # prb.timeSteps = [(1e-05, 100)]

        try:
            from pymatsolver import MumpsSolver
-            self.prb.Solver = MumpsSolver
+            prb.Solver = MumpsSolver
        except ImportError, e:
-            self.prb.Solver  = SolverLU
+            prb.Solver  = SolverLU

-        self.sigma = np.ones(mesh.nCz)*1e-8
-        self.sigma[mesh.vectorCCz<0] = 1e-1
-        self.sigma = np.log(self.sigma[active])
+        m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)

-        self.prb.pair(survey)
-        self.mesh = mesh
+        prb.pair(survey)

-    def test_DerivG(self):
-        """
-            Test the derivative of c with respect to sigma
-        """
+        return mesh, prb, m

-        # Random model and perturbation
-        sigma = np.random.rand(self.prb.mapping.nP)
-
-        f = self.prb.fields(sigma)
-        dm = 1000*np.random.rand(self.prb.mapping.nP)
-        h = 0.01
-
-        derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
-        print '\ntest_DerivG'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_dUdM(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        dm = 10*np.random.rand(prb.mapping.nP)
-        f = prb.fields(sigma)
-
-        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
-        print '\n'
-        print 'test_Deriv_dUdM'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_J(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-        d_sig = 10*np.random.rand(prb.mapping.nP)
+class TDEM_bDerivTests(unittest.TestCase):


-        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-        print '\n'
-        print 'test_Deriv_J'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
+    if testDeriv:
+        def Deriv_J(self, rxcomp='bz'):

-    def test_projectAdjoint(self):
-        prb = self.prb
-        survey = prb.survey
-        nSrc = survey.nSrc
-        mesh = self.mesh
+            mesh, prb, m0 = setUp(rxcomp)

-        # Generate random fields and data
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(prb.nT):
-            f[:,'b',i] = np.random.rand(mesh.nF, nSrc)
-            f[:,'e',i] = np.random.rand(mesh.nE, nSrc)
-        d_vec = np.random.rand(survey.nD)
-        d = Survey.Data(survey,v=d_vec)
+            prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]

-        # Check that d.T*Q*f = f.T*Q.T*d
-        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
-        V2 = np.sum((f.tovec())*(survey.evalDeriv(None, v=d, adjoint=True).tovec()))
+            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m0, mx)]
+            print '\n'
+            print 'test_Deriv_J %s'%rxcomp
+            Tests.checkDerivative(derChk, m0, plotIt=False, num=3, eps=1e-20)

-        self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
+        def test_Jvec_bx(self):
+            self.Deriv_J('bx')

-    def test_adjointGvecVsGtvec(self):
-        mesh = self.mesh
-        prb = self.prb
+        def test_Jvec_bz(self):
+            self.Deriv_J('bz')

-        m = np.random.rand(prb.mapping.nP)
-        sigma = np.random.rand(prb.mapping.nP)
+        def test_Jvec_ey(self):
+            self.Deriv_J('ey')

-        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            u[:,'b',i] = np.random.rand(mesh.nF, 2)
-            u[:,'e',i] = np.random.rand(mesh.nE, 2)
+    if testAdjoint:
+        def adjointJvecVsJtvec(self, rxcomp='bz'):
+            print ' \n Testing Adjoint %s' %rxcomp
+            mesh, prb, m0 = setUp(rxcomp)

-        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            v[:,'b',i] = np.random.rand(mesh.nF, 2)
-            v[:,'e',i] = np.random.rand(mesh.nE, 2)
+            m = np.random.rand(prb.mapping.nP)
+            d = np.random.rand(prb.survey.nD)

-        V1 = m.dot(prb.Gtvec(sigma, v, u))
-        V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6)
+            V1 = d.dot(prb.Jvec(m0, m))
+            V2 = m.dot(prb.Jtvec(m0, d))

-    def test_adjointJvecVsJtvec(self):
-        mesh = self.mesh
-        prb = self.prb
-        sigma = self.sigma
+            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+            passed = np.abs(V1-V2) < tol
+            print '    ', V1, V2, np.abs(V1-V2), tol, passed
+            self.assertTrue(passed)

-        m = np.random.rand(prb.mapping.nP)
-        d = np.random.rand(prb.survey.nD)
+        def test_JvecVsJtvec_bx(self):
+            self.adjointJvecVsJtvec('bx')

-        V1 = d.dot(prb.Jvec(sigma, m))
-        V2 = m.dot(prb.Jtvec(sigma, d))
-        print 'AdjointTest', V1, V2
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
+        def test_JvecVsJtvec_bz(self):
+            self.adjointJvecVsJtvec('bz')
+
+        def test_JvecVsJtvec_ey(self):
+            self.adjointJvecVsJtvec('ey')



@@ -1,94 +0,0 @@
-import unittest
-from SimPEG import *
-from SimPEG import EM
-
-plotIt = False
-
-def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
-    cs = 5.
-    ncx = 20
-    ncy = 6
-    npad = 20
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hy], '00C')
-
-    active = mesh.vectorCCz<0.
-    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
-
-    rxOffset = 40.
-
-    srcs = []
-    for ii in range(nSrc):
-        rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
-        srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
-
-    survey = EM.TDEM.SurveyTDEM(srcs)
-
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-    # prb.timeSteps = [1e-5]
-    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-    # prb.timeSteps = [(1e-05, 100)]
-
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-        prb.Solver  = SolverLU
-
-    sigma = np.ones(mesh.nCz)*1e-8
-    sigma[mesh.vectorCCz<0] = 1e-1
-    sigma = np.log(sigma[active])
-
-    prb.pair(survey)
-    return prb, mesh, sigma
-
-def dotestJvec(prb, mesh, sigma):
-    prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-    # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-    d_sig = 10*np.random.rand(prb.mapping.nP)
-    derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-    return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
-
-def dotestAdjoint(prb, mesh, sigma):
-    m = np.random.rand(prb.mapping.nP)
-    d = np.random.rand(prb.survey.nD)
-
-    V1 = d.dot(prb.Jvec(sigma, m))
-    V2 = m.dot(prb.Jtvec(sigma, d))
-    print 'AdjointTest', V1, V2
-    return np.abs(V1-V2)/np.abs(V1), 1e-6
-
-class TDEM_bDerivTests(unittest.TestCase):
-
-    def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
-    def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
-
-    def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
-    def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
-
-    def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
-    def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
-
-    def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
-    def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
-
-    def test_Jvec_dbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt')))
-    def test_Adjoint_dbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt')))
-
-    def test_Jvec_dbzdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbzdt')))
-    def test_Adjoint_dbzdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbzdt')))
-
-    def test_Jvec_dbxdtbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt,bz')))
-    def test_Adjoint_dbxdtbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt,bz')))
-
-    def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey')))
-    def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey')))
-
-    def test_Jvec_eybzdbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey,bz,dbxdt')))
-    def test_Adjoint_eybzdbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey,bz,dbxdt')))
-
-
-if __name__ == '__main__':
-    unittest.main()
@@ -0,0 +1,76 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+
+TOL = 1e-5
+FLR = 1e-20
+
+np.random.seed(seed=25) # set a seed so that the same conductivity model is used for all runs
+
+def setUp(prbtype = 'b', rxcomp='bz'):
+    cs = 5.
+    ncx = 20
+    ncy = 15
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+#
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+
+    rxOffset = 10.
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
+
+    survey = EM.TDEM.Survey([src])
+
+    if prbtype == 'b':
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    elif prbtype == 'e':
+        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
+
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
+
+    m = np.log(1e-1)*np.ones(prb.mapping.nP) #+ 1e-2*np.random.randn(prb.mapping.nP)
+
+    prb.pair(survey)
+    mesh = mesh
+
+    return prb, m, mesh
+
+def CrossCheck(prbtype1='b', prbtype2='e', rxcomp='bz'):
+
+    prb1,m1,mesh1 = setUp(prbtype1, rxcomp)
+    prb2,m2,mesh2 = setUp(prbtype2, rxcomp)
+
+    assert (m1 == m2).all(), 'Models for two formulations are different'
+
+    d1 = prb1.survey.dpred(m1)
+    d2 = prb2.survey.dpred(m2)
+
+
+    check = np.linalg.norm(d1 - d2)
+    tol   = 0.5 * (np.linalg.norm(d1) + np.linalg.norm(d2)) * TOL
+    passed = check < tol
+
+    print 'Checking %s, %s for %s data'%(prbtype1, prbtype2, rxcomp)
+    print '  ', np.linalg.norm(d1), np.linalg.norm(d2), np.linalg.norm(check), tol, passed
+
+    assert passed
+
+class TDEM_cross_check_EB(unittest.TestCase):
+    def test_EB_ey(self):
+        CrossCheck('b','e','ey')
+
+if __name__ == '__main__':
+    unittest.main()
+
@@ -29,12 +29,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap

-    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
-    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
+    src = EM.TDEM.Src.MagDipole([rx], waveform= EM.TDEM.Src.StepOffWaveform(), loc=np.array([0., 0., 0.]))
    # src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))

-    survey = EM.TDEM.SurveyTDEM([src])
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    survey = EM.TDEM.Survey([src])
+    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
    prb.Solver = MumpsSolver

    prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
@@ -50,6 +50,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,

    ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1])
    log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
+
+    print ' |bz_ana| = ',np.linalg.norm(bz_ana), ' |bz_num| = ', np.linalg.norm(bz_calc), ' |bz_ana - bz_num| =', np.linalg.norm(bz_ana-bz_calc)
    print 'Difference: ', log10diff

    if showIt == True:
@@ -61,6 +63,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
    return log10diff


+class TDEM_SimpleSrcTests(unittest.TestCase):
+    def test_source(self):
+        waveform = EM.TDEM.Src.StepOffWaveform()
+        assert waveform.eval(0.) == 0.
+
+
 class TDEM_bTests(unittest.TestCase):

    def test_analytic_p2_CYL_50m(self):
Author	SHA1	Message	Date
seogi_macbook	d010886c7e	minor fix	2016-06-28 19:16:41 -07:00
seogi_macbook	222035d0b8	Modify testings ... for TDEM	2016-06-28 19:15:20 -07:00
seogi_macbook	f4f10ff027	remove temporary print statemement.	2016-06-28 19:01:32 -07:00
seogi_macbook	53adb3abd7	Tracking memory leakages in TDEM code - make PTv out from the loop - do not use Data class in the loop for Jv	2016-06-28 19:00:39 -07:00
seogi_macbook	61a1d3207b	delete fields_deriv	2016-06-28 11:42:30 -07:00
Lindsey Heagy	875885c684	typo fix	2016-06-28 11:28:44 -07:00
Lindsey Heagy	e179b5e71d	make sure factors of JtVec get cleaned	2016-06-28 11:25:17 -07:00
seogi_macbook	43d0b0369f	Add printing for TDEM fwd problem.	2016-06-28 01:02:03 -07:00
seogi_macbook	b7546d332b	- Add CircularLoop for TDEM source - Fix bug for verbose option in Problem_b (TDEM)	2016-06-27 22:59:57 -07:00
Lindsey Heagy	4e02ff3561	Merge branch 'em/dev' into em/ref/tdem	2016-06-20 11:45:44 -07:00
Lindsey Heagy	39ddec8702	Merge pull request #339 from simpeg/em/secondary_Rx Implementation for Inverting Secondary B field	2016-06-20 12:45:00 -06:00
seogi_macbook	e30a657abb	update 1D example for simpegEMpaper	2016-06-17 18:22:57 -07:00
seogi_macbook	71213b9f91	Merge branch 'em/secondary_Rx' of https://github.com/simpeg/simpeg into em/ref/tdem	2016-06-17 17:39:01 -07:00
seogi_macbook	2fb0f3fbbb	Implementation for Inverting Secondary B field	2016-06-17 07:50:13 -07:00
seogi_macbook	fcac3e9bc7	simpegEM paper example	2016-06-12 10:33:42 -07:00
Lindsey Heagy	9c7c89b8ba	Merge branch 'em/dev' into em/ref/tdem # Conflicts: # SimPEG/EM/TDEM/SurveyTDEM.py	2016-06-02 07:41:03 -07:00
seogi_macbook	6f8f065643	Update current changes ... about u -> f	2016-06-01 23:45:52 -07:00
seogi_macbook	8e6536b2ed	Merge branch 'em/ref/tdem' of https://github.com/simpeg/simpeg into em/ref/tdem Conflicts: SimPEG/EM/FDEM/FieldsFDEM.py SimPEG/EM/FDEM/SurveyFDEM.py	2016-06-01 22:11:13 -07:00
Lindsey	382b31bd12	Merge pull request #310 from simpeg/ref/regularization modularizing regularization	2016-05-31 14:09:54 -07:00
Lindsey	6cc509020a	Merge pull request #323 from simpeg/feat/plotImage-curvilinear plotImage for curvilinear mesh	2016-05-30 08:27:56 -07:00
Lindsey Heagy	f2e13182bf	Merge branch 'dev' into feat/plotImage-curvilinear # Conflicts: # SimPEG/Examples/__init__.py	2016-05-29 16:43:40 -07:00
D Fournier	09cd9c7fa3	Merge branch 'dev' into ref/regularization	2016-05-29 16:26:58 -07:00
Lindsey Heagy	62eb4541cb	update example name --> based on mesh	2016-05-29 15:48:15 -07:00
D Fournier	f8b86abd5a	Merge branch 'ref/dev' into ref/regularization	2016-05-29 14:10:46 -07:00
D Fournier	e476bf0059	Cleanup Sparse Reg and Directives	2016-05-29 14:09:29 -07:00
Rowan Cockett	825511e9d3	Add a curvilinear plotImage function, update example.	2016-05-29 13:17:22 -07:00
D Fournier	3b4bec9c0b	Refactor IRLS iterations, full solves from l2->lp Adapt Example	2016-05-28 11:27:09 -07:00
D Fournier	022e1f7660	Update IRLS directive to allow multiple GN iterations. Remove modifications to the ProjGN solver. Update IRLS example.	2016-05-27 13:11:31 -07:00
D Fournier	fd3bde787f	Propose change to the Projected_GNCG solver. Add inner GN iterations. Nice improvement to the convergence of IRLS	2016-05-12 14:58:16 -07:00
D Fournier	3cc46131a3	Temporary change ... comment out W and Wsmooth	2016-05-12 08:31:01 -07:00
D Fournier	cd2360b815	Stash the regularization between each beta	2016-05-11 23:04:14 -07:00
Lindsey Heagy	6de786c972	Merge branch 'em/dev' into em/ref/tdem # Conflicts: # SimPEG/EM/FDEM/FDEM.py # SimPEG/EM/FDEM/SrcFDEM.py # SimPEG/EM/FDEM/SurveyFDEM.py # SimPEG/Examples/EM_FDEM_1D_Inversion.py	2016-05-11 09:21:59 -07:00
D Fournier	e10d6878fb	Remove Wsmooth from def W and replace by parts	2016-05-11 07:58:06 -07:00
Lindsey Heagy	3dd9ecc9cd	fix tikhonov 2Deriv	2016-05-10 22:10:19 -07:00
Lindsey Heagy	90a3030796	fixed 2 deriv	2016-05-10 21:39:15 -07:00
D Fournier	7964ebce50	Update directive to None the Wsmooth after iteration.	2016-05-10 17:20:46 -07:00
Lindsey Heagy	955bd54019	notation cleanup in Regularization	2016-05-10 16:46:11 -07:00
Lindsey Heagy	2a802c1aa3	weights --> cell_weights, removed vol term from simple regularization	2016-05-10 16:39:47 -07:00
Lindsey Heagy	3f0c89f10b	remove extra Ws	2016-05-10 14:53:43 -07:00
Lindsey Heagy	eaa37f42e4	remove duplicate evalSmall	2016-05-10 14:51:56 -07:00
D Fournier	fb5434695f	Alpha_s default to 1.0	2016-05-10 14:31:45 -07:00
D Fournier	e037597ecd	Merge branch 'feat/sparse-regularization' into ref/regularization	2016-05-10 13:37:16 -07:00
D Fournier	b4ab60c260	Add model mapping to sparse regularization	2016-05-05 11:55:56 -07:00
Lindsey Heagy	fbb8cf2731	modularizing regularization	2016-05-04 23:17:01 -07:00
seogi_macbook	20d75a9704	simple fix for inverting secondary magnetic fields	2016-04-05 09:38:49 -07:00
Lindsey Heagy	1be4082ea3	parse out SrcTDEM	2016-03-20 22:41:40 -07:00
Lindsey Heagy	2d8bbdce45	working to debug JTv for e-formulation... still some work to do	2016-03-20 11:42:31 -07:00
Lindsey Heagy	0be942730a	start of Problem_e (e Jvec working)	2016-03-18 15:47:31 -07:00
Lindsey Heagy	74f4705048	cleanup of TDEM example	2016-03-14 13:25:09 -07:00
Lindsey Heagy	a9efb2fc8a	add dbdt to testing	2016-03-14 13:05:09 -07:00
Lindsey Heagy	c91815d14f	adjoint hooked up for b formulation	2016-03-14 12:28:05 -07:00
Lindsey Heagy	fe91312917	bx, bz running and passing, ey failing adjoint --> I think the initial fields are not being taken care of correctly in the deriv	2016-03-13 14:02:32 -07:00
Lindsey Heagy	605e19eb22	combos will be tested in TDEM_b_DerivAdjoint	2016-03-13 12:35:44 -07:00
Lindsey Heagy	f549756208	create only 1 fields object in Jtvec	2016-03-13 12:28:38 -07:00
Lindsey Heagy	576459d17c	only save previous tilmestep for back solve (don't need all times)	2016-03-13 12:17:38 -07:00
Lindsey Heagy	ea4721a941	first pass at multisrc Jtvec (will be hugely memory inefficient at the moment)	2016-03-13 12:07:26 -07:00
Lindsey Heagy	c708ceb53d	cleanup and minimal docs for Jvec, JTvec	2016-03-13 11:09:11 -07:00
Lindsey Heagy	a1ecef0709	first shot through of passing Jtvec for TDEM problem (code will need to be cleaned up, but it passes!)	2016-03-12 15:13:17 -08:00
Lindsey Heagy	fb66acea11	things in the adjoint are the right sizes, but not passing... +1,-1 somewhere??	2016-03-10 13:08:39 -08:00
Lindsey Heagy	1b401feb54	cleaned up Jvec	2016-03-08 19:56:14 -08:00
Lindsey Heagy	ceff861413	forward and Jvec using Adiag, Asubdiag, RHS	2016-03-08 19:38:44 -08:00
Lindsey Heagy	705cdd0c52	jtvec runs, fails	2016-03-08 16:38:45 -08:00
Lindsey Heagy	1d2eac62a3	e hooked up with Jvec	2016-03-06 21:42:30 -08:00
Lindsey Heagy	5cf0acd153	TDEM bderiv from b formulation working	2016-03-06 16:06:14 -08:00
Lindsey Heagy	664adb04ac	light notation cleanup in Jvec, testing ADeriv --> passes, Jvec is still first order	2016-03-06 15:10:13 -08:00
Lindsey Heagy	4f31e4e002	tdem deriv runs but is first order at the moment	2016-03-06 11:13:48 -08:00
Lindsey Heagy	0bfc816ecc	starting sensitivities	2016-03-04 10:43:19 -08:00
Lindsey Heagy	7ece7c3edb	sketching out code	2016-03-04 08:45:36 -08:00
Lindsey Heagy	8bd027c2b2	sketch of derivs	2016-02-24 15:24:56 -08:00
Lindsey Heagy	617241ad4e	TDEM forward refactor (no derive yet)	2016-02-22 18:15:29 -08:00
Lindsey Heagy	d5967d20b9	forward is running, but not passing	2016-02-22 17:42:22 -08:00
Lindsey Heagy	ecbd5c21f5	sketch of sources and waveforms	2016-02-22 15:13:25 -08:00
Lindsey Heagy	341e902469	merged in em/dev	2016-02-22 11:26:30 -08:00
Lindsey Heagy	cd51ab8be7	sketching out TDEM problem	2016-02-22 11:04:18 -08:00