misc fixes

- 2to3 on setup.py to allow python 3 import
- removed uneeded circular import in Maps.py
- added missing old_div o mathutils

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.10
+current_version = 0.1.12
 files = setup.py SimPEG/__init__.py docs/conf.py
 

.gitignore

@@ -39,3 +39,5 @@ nosetests.xml
 *.sublime-workspace
 docs/_build/
 Makefile
+docs/warnings.txt
+.DS_Store

.travis.yml

@@ -1,6 +1,7 @@
 language: python
 python:
   - 2.7
+  - 3.4
 
 sudo: false
 
@@ -24,18 +25,25 @@ env:
     - TEST_DIR=tests/examples
     - TEST_DIR=tests/em/fdem/inverse/adjoint
     - TEST_DIR=tests/em/fdem/forward
+    - TEST_DIR=tests/docs;
+      GAE_PYTHONPATH=${HOME}/.cache/google_appengine;
+      PATH=$PATH:${HOME}/google-cloud-sdk/bin;
+      PYTHONPATH=${PYTHONPATH}:${GAE_PYTHONPATH};
+      CLOUDSDK_CORE_DISABLE_PROMPTS=1
 
 # Setup anaconda
 before_install:
-  - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; fi
+# Install packages
+  - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh
+    -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh
+    -O miniconda.sh; fi
   - chmod +x miniconda.sh
   - ./miniconda.sh -b
-  - export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
+  - export PATH=/home/travis/anaconda/bin:/home/travis/anaconda3/bin:/home/travis/miniconda/bin:/home/travis/miniconda3/bin:$PATH
   - conda update --yes conda
 
-# Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk sphinx
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git
@@ -46,11 +54,26 @@ install:
 
 # Run test
 script:
+  # test docs
   - nosetests $TEST_DIR --with-cov --cov SimPEG --cov-config .coveragerc -v -s
 
 # Calculate coverage
 after_success:
-  - coveralls --config_file .coveragerc
+  - bash <(curl -s https://codecov.io/bash)
+  - if [ "$TRAVIS_BRANCH" = "master" -a "$TRAVIS_PULL_REQUEST" = "false" ]; then
+      if [ ${TEST_DIR} == "tests/docs" ]; then
+        python scripts/fetch_gae_sdk.py $(dirname "${GAE_PYTHONPATH}");
+        openssl aes-256-cbc -K $encrypted_93066031461c_key -iv $encrypted_93066031461c_iv
+          -in docs/credentials.tar.gz.enc -out credentials.tar.gz -d ;
+        if [ ! -d ${HOME}/google-cloud-sdk ]; then curl https://sdk.cloud.google.com | bash; fi ;
+        tar -xzf credentials.tar.gz ;
+        gcloud auth activate-service-account --key-file client-secret.json ;
+        gcloud config set project simpegdocs;
+        gcloud -q components update gae-python;
+        gcloud -q preview app deploy ./docs/app.yaml --version ${TRAVIS_COMMIT} --promote;
+      fi;
+    fi
+
 
 notifications:
   email:

README.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
     :alt: SimPEG Logo
 
 ======
@@ -28,7 +28,11 @@ SimPEG
 .. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
     :alt: gitter chat room at https://gitter.im/simpeg/simpeg
     :target: https://gitter.im/simpeg/simpeg
-    
+
+.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
+  :target: https://codecov.io/gh/simpeg/simpeg
+
+
 Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
 
 The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

SimPEG/DCIP/BaseDC.py

@@ -1,3 +1,11 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from builtins import super
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
 
 class FieldsDC_CC(Problem.Fields):
@@ -61,7 +69,7 @@ class SrcDipole(Survey.BaseSrc):
         pts = [self.loc[0], self.loc[1]]
         inds = Utils.closestPoints(prob.mesh, pts)
         q = np.zeros(prob.mesh.nC)
-        q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
+        q[inds] = - self.current * (np.r_[1., -1.] / prob.mesh.vol[inds])
         # self._rhsDict[mesh] = q
         # return self._rhsDict[mesh]
         return q
@@ -162,8 +170,8 @@ class ProblemDC_CC(Problem.BaseProblem):
         """
             Makes the matrix A(m) for the DC resistivity problem.
 
-            :param numpy.array m: model
-            :rtype: scipy.csc_matrix
+            :param numpy.ndarray m: model
+            :rtype: scipy.sparse.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/BaseIP.py

@@ -1,5 +1,12 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
-from BaseDC import SurveyDC, FieldsDC_CC
+from .BaseDC import SurveyDC, FieldsDC_CC
 
 class SurveyIP(SurveyDC):
     """
@@ -52,7 +59,7 @@ class ProblemIP(Problem.BaseProblem):
             # sigma = self.curModel.transform
             sigma = self.sigma
             Av = self.mesh.aveF2CC
-            self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
+            self._Msig = Utils.sdiag(1//(self.mesh.dim * Av.T * (1/sigma)))
         return self._Msig
 
     @property
@@ -71,7 +78,7 @@ class ProblemIP(Problem.BaseProblem):
             Makes the matrix A(m) for the DC resistivity problem.
 
             :param numpy.array m: model
-            :rtype: scipy.csc_matrix
+            :rtype: scipy.sparse.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/DCIPUtils.py

@@ -1,6 +1,16 @@
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from builtins import open
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+from builtins import map
+from builtins import range
 from SimPEG import np, Utils
-import BaseDC as DC
-import BaseDC as IP
+from . import BaseDC as DC
+from . import BaseDC as IP
 import warnings
 
 def getActiveindfromTopo(mesh, topo):
@@ -67,7 +77,7 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
         if "!" in line.split(): continue
         elif line == '\n': continue
         elif line == ' \n': continue
-        temp =  map(float, line.split())
+        temp =  list(map(float, line.split()))
         # Read a line for the current electrode
         if len(temp) == 5: # SRC: Only X and Y are provided (assume no topography)
             #TODO consider topography and assign the closest cell center in the earth
@@ -228,10 +238,10 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
 
             elif surveyType == 'dipole-dipole':
 
-                leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
+                leg = data * 2*np.pi / (1/MA - 1/MB - 1/NB + 1/NA)
 
             else:
-                print """unitType must be 'pole-dipole' | 'dipole-dipole' """
+                print("""unitType must be 'pole-dipole' | 'dipole-dipole' """)
                 break
 
 
@@ -246,11 +256,11 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
                 rho = np.hstack([rho,leg])
 
             else:
-                print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
+                print("""unitType must be 'appResistivity' | 'appConductivity' | 'volt' """)
                 break
 
-        midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
-        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
+        midx = np.hstack([midx, (Cmid + Pmid)/2])
+        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2])
 
     # Grid points
     grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
@@ -340,11 +350,11 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
     dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
     dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
 
-    nstn = np.floor( dl_len / AM_sep )
+    nstn = np.floor(dl_len / AM_sep)
 
     # Compute discrete pole location along line
-    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
-    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep
+    stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*AM_sep
+    stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*AM_sep
 
     # Create line of P1 locations
     M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
@@ -376,15 +386,15 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
             AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
 
             # Number of receivers to fit
-            nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx])
+            nstn = np.min([(AB - MN_sep) // AM_sep, nrx])
 
             # Check if there is enough space, else break the loop
             if nstn <= 0:
                 continue
 
             # Compute discrete pole location along line
-            stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep
-            stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep
+            stn_x = N[ii,0] + dl_x*MN_sep + np.array(list(range(int(nstn))))*dl_x*AM_sep
+            stn_y = N[ii,1] + dl_y*MN_sep + np.array(list(range(int(nstn))))*dl_y*AM_sep
 
             # Create receiver poles
             # Create line of P1 locations
@@ -419,15 +429,15 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
         box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
         box_w = box_l/2.
 
-        nstn = np.floor( box_l / AM_sep )
+        nstn = box_l // AM_sep
 
         # Compute discrete pole location along line
-        stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
-        stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep
+        stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*AM_sep
+        stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*AM_sep
 
         # Define number of cross lines
-        nlin = int(np.floor( box_w / AM_sep ))
-        lind = range(-nlin,nlin+1)
+        nlin = int(box_w // AM_sep)
+        lind = list(range(-nlin,nlin+1))
 
         ngrad = nstn * len(lind)
 
@@ -449,7 +459,7 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
         srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
         SrcList.append(srcClass)
     else:
-        print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
+        print("""surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """)
 
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
@@ -668,7 +678,7 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
         obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
 
     else:
-        print "rtype must be 'DC'(default) | 'IP'"
+        print("rtype must be 'DC'(default) | 'IP'")
 
     # Pre-allocate
     srcLists = []

SimPEG/DCIP/Utils.py

@@ -1,3 +1,10 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 import numpy as np
 
 def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
@@ -5,7 +12,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
     import SimPEG.DCIP as DC
 
     elocs = np.arange(0,aSpacing*nElecs,aSpacing)
-    elocs -= (nElecs*aSpacing - aSpacing)/2
+    elocs -= (nElecs*aSpacing - aSpacing) / 2
     space = 1
     WENNER = np.zeros((0,),dtype=int)
     for ii in range(nElecs):

SimPEG/DCIP/__init__.py

@@ -1,4 +1,10 @@
-from BaseDC import *
-from BaseIP import *
-from DCIPUtils import *
-import Utils
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .BaseDC import *
+from .BaseIP import *
+from .DCIPUtils import *
+from . import Utils

SimPEG/DataMisfit.py

@@ -1,7 +1,16 @@
-import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import object
+from . import Utils, Survey, Problem
+import numpy as np, scipy.sparse as sp, gc
+from future.utils import with_metaclass
 
 
-class BaseDataMisfit(object):
+class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, object)):
     """BaseDataMisfit
 
         .. note::
@@ -9,8 +18,6 @@ class BaseDataMisfit(object):
             You should inherit from this class to create your own data misfit term.
     """
 
-    __metaclass__ = Utils.SimPEGMetaClass
-
     debug   = False  #: Print debugging information
     counter = None   #: Set this to a SimPEG.Utils.Counter() if you want to count things
 
@@ -93,11 +100,11 @@ class l2_DataMisfit(BaseDataMisfit):
             survey = self.survey
 
             if getattr(survey,'std', None) is None:
-                print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
+                print('SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%')
                 survey.std = 0.05
 
             if getattr(survey, 'eps', None) is None:
-                print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
+                print('SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||')
                 survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
 
             self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))

SimPEG/Directives.py

@@ -1,4 +1,15 @@
-import Utils, numpy as np
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from builtins import open
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+from builtins import str
+from builtins import object
+from . import Utils
+import numpy as np
 
 class InversionDirective(object):
     """InversionDirective"""
@@ -15,7 +26,7 @@ class InversionDirective(object):
     @inversion.setter
     def inversion(self, i):
         if getattr(self,'_inversion',None) is not None:
-            print 'Warning: InversionDirective %s has switched to a new inversion.' % self.__name__
+            print('Warning: InversionDirective %s has switched to a new inversion.' % self.__name__)
         self._inversion = i
 
     @property
@@ -68,7 +79,7 @@ class DirectiveList(object):
     def inversion(self, i):
         if self.inversion is i: return
         if getattr(self,'_inversion',None) is not None:
-            print 'Warning: %s has switched to a new inversion.' % self.__name__
+            print('Warning: %s has switched to a new inversion.' % self.__name__)
         for d in self.dList:
             d.inversion = i
         self._inversion = i
@@ -120,7 +131,7 @@ class BetaEstimate_ByEig(InversionDirective):
             :return: beta0
         """
 
-        if self.debug: print 'Calculating the beta0 parameter.'
+        if self.debug: print('Calculating the beta0 parameter.')
 
         m = self.invProb.curModel
         f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
@@ -141,7 +152,7 @@ class BetaSchedule(InversionDirective):
 
     def endIter(self):
         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
+            if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
             self.invProb.beta /= self.coolingFactor
 
 
@@ -167,7 +178,7 @@ class TargetMisfit(InversionDirective):
 
 
 
-class _SaveEveryIteration(InversionDirective):
+class SaveEveryIteration(InversionDirective):
     @property
     def name(self):
         if getattr(self, '_name', None) is None:
@@ -188,21 +199,21 @@ class _SaveEveryIteration(InversionDirective):
         self._fileName = value
 
 
-class SaveModelEveryIteration(_SaveEveryIteration):
+class SaveModelEveryIteration(SaveEveryIteration):
     """SaveModelEveryIteration"""
 
     def initialize(self):
-        print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
+        print("SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName)
 
     def endIter(self):
         np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
 
 
-class SaveOutputEveryIteration(_SaveEveryIteration):
+class SaveOutputEveryIteration(SaveEveryIteration):
     """SaveModelEveryIteration"""
 
     def initialize(self):
-        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
+        print("SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName)
         f = open(self.fileName+'.txt', 'w')
         f.write("  #     beta     phi_d     phi_m       f\n")
         f.close()
@@ -212,11 +223,11 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
         f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
         f.close()
 
-class SaveOutputDictEveryIteration(_SaveEveryIteration):
+class SaveOutputDictEveryIteration(SaveEveryIteration):
     """SaveOutputDictEveryIteration"""
 
     def initialize(self):
-        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+        print("SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName)
 
     def endIter(self):
         # Save the data.
@@ -253,8 +264,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
-    eps_p = None
-    eps_q = None
+    eps = None
     norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
@@ -263,6 +273,7 @@ class Update_IRLS(InversionDirective):
     f_old = None
     f_min_change = 1e-2
     beta_tol = 5e-2
+    prctile = 95
 
     # Solving parameter for IRLS (mode:2)
     IRLSiter   = 0
@@ -294,12 +305,25 @@ class Update_IRLS(InversionDirective):
 
         # 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..."
+            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
+
+            # Either use the supplied epsilon, or fix base on distribution of
+            # model values
+            if getattr(self, 'reg.eps', None) is None:
+                self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
+            else:
+                self.reg.eps_p = self.eps[0]
+
+            if getattr(self, 'reg.eps', None) is None:
+                self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
+            else:
+                self.reg.eps_q = self.eps[1]
+
+            print("L[p qx qy qz]-norm : " + str(self.reg.norms))
+            print("eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q))
+
             self.reg.norms = self.norms
             self.coolingFactor = 1.
             self.coolingRate = 1
@@ -315,7 +339,7 @@ class Update_IRLS(InversionDirective):
 
         # 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
+            if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
             self.invProb.beta /= self.coolingFactor
 
 
@@ -327,30 +351,30 @@ class Update_IRLS(InversionDirective):
             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)
+            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
+                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"
+                print("Minimum decrease in regularization. End of IRLS")
                 self.opt.stopNextIteration = True
                 return
             else:
                 self.f_old = phim_new
 
-            # 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])
-                else:
-                    self.reg.eps = eps
+#            # 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])
+#                else:
+#                    self.reg.eps = eps
 
             # Get phi_m at the end of current iteration
             self.phi_m_last = self.invProb.phi_m_last

SimPEG/EM/Analytics/DC.py

@@ -1,3 +1,10 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 import numpy as np
 from scipy.constants import mu_0, pi
 from scipy import special

SimPEG/EM/Analytics/FDEM.py

@@ -1,4 +1,9 @@
 from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import numpy as np
 from scipy.constants import mu_0, pi
 from scipy.special import erf

SimPEG/EM/Analytics/FDEMDipolarfields.py

@@ -0,0 +1,307 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+import numpy as np
+from scipy.constants import mu_0, pi, epsilon_0
+from scipy.special import erf
+from SimPEG import Utils
+
+omega = lambda f: 2.*np.pi*f
+# TODO:
+# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.)
+# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.):
+
+    """
+        Computing Analytic Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+    mid   = -k**2 * r**2 + 3*1j*k*r + 3
+
+    if orientation.upper() == 'X':
+        Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ey = front*(dx*dy  / r**2)*mid
+        Ez = front*(dx*dz  / r**2)*mid
+        return Ex, Ey, Ez
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ez = front*(dy*dz  / r**2)*mid
+        Ex = front*(dy*dx  / r**2)*mid
+        return Ex, Ey, Ez
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ex = front*(dz*dx  / r**2)*mid
+        Ey = front*(dz*dy  / r**2)*mid
+        return Ex, Ey, Ez
+
+
+def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+    mid   = -k**2 * r**2 + 3*1j*k*r + 3
+
+    if orientation.upper() == 'X':
+        Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.))
+        Ey_galvanic = front*(dx*dy  / r**2)*mid
+        Ez_galvanic = front*(dx*dz  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.))
+        Ez_galvanic = front*(dy*dz  / r**2)*mid
+        Ex_galvanic = front*(dy*dx  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.))
+        Ex_galvanic = front*(dz*dx  / r**2)*mid
+        Ey_galvanic = front*(dz*dy  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+
+def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+
+    if orientation.upper() == 'X':
+        Ex_inductive = front*(k**2 * r**2)
+        Ey_inductive = np.zeros_like(Ex_inductive)
+        Ez_inductive = np.zeros_like(Ex_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey_inductive = front*(k**2 * r**2)
+        Ez_inductive = np.zeros_like(Ey_inductive)
+        Ex_inductive = np.zeros_like(Ey_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez_inductive = front*(k**2 * r**2)
+        Ex_inductive = np.zeros_like(Ez_inductive)
+        Ey_inductive = np.zeros_like(Ez_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+
+def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Jx = sig*Ex
+    Jy = sig*Ey
+    Jz = sig*Ez
+    return Jx, Jy, Jz
+
+
+def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Jx_galvanic = sig*Ex_galvanic
+    Jy_galvanic = sig*Ey_galvanic
+    Jz_galvanic = sig*Ez_galvanic
+    return Jx_galvanic, Jy_galvanic, Jz_galvanic
+
+
+def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Jx_inductive = sig*Ex_inductive
+    Jy_inductive = sig*Ey_inductive
+    Jz_inductive = sig*Ez_inductive
+    return Jx_inductive, Jy_inductive, Jz_inductive
+
+
+def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Magnetic fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r)
+
+    if orientation.upper() == 'X':
+        Hy = front*(-dz  / r)
+        Hz = front*(dy  / r)
+        Hx = np.zeros_like(Hy)
+        return Hx, Hy, Hz
+
+    elif orientation.upper() == 'Y':
+        Hx = front*(dz  / r)
+        Hz = front*(-dx  / r)
+        Hy = np.zeros_like(Hx)
+        return Hx, Hy, Hz
+
+    elif orientation.upper() == 'Z':
+        Hx = front*(-dy / r)
+        Hy = front*(dx  / r)
+        Hz = np.zeros_like(Hx)
+        return Hx, Hy, Hz
+
+
+def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Magnetic flux densites from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Bx = mu*Hx
+    By = mu*Hy
+    Bz = mu*Hz
+    return Bx, By, Bz
+
+
+def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Electric vector potentials from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*r)
+
+    if orientation.upper() == 'X':
+        Ax = front*np.exp(-1j*k*r)
+        Ay = np.zeros_like(Ax)
+        Az = np.zeros_like(Ax)
+        return Ax, Ay, Az
+
+    elif orientation.upper() == 'Y':
+        Ay = front*np.exp(-1j*k*r)
+        Ax = np.zeros_like(Ay)
+        Az = np.zeros_like(Ay)
+        return Ax, Ay, Az
+
+    elif orientation.upper() == 'Z':
+        Az = front*np.exp(-1j*k*r)
+        Ax = np.zeros_like(Ay)
+        Ay = np.zeros_like(Ay)
+        return Ax, Ay, Az
+
+
+
+
+

SimPEG/EM/Analytics/FDEMcasing.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Utils, np
 from scipy.constants import mu_0, epsilon_0
 from SimPEG.EM.Utils.EMUtils import k
@@ -34,7 +40,7 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
     sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
     k2 = k(freq,sigma[2],mu[2],eps)
 
-    return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2)
+    return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1. / sqrtr2z2)
 
 
 def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):

SimPEG/EM/Analytics/TDEM.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import numpy as np
 from scipy.constants import mu_0, pi
 from scipy.special import erf

SimPEG/EM/Analytics/__init__.py

@@ -1,4 +1,11 @@
-from TDEM import hzAnalyticDipoleT
-from FDEM import hzAnalyticDipoleF
-from FDEMcasing import *
-from DC import DCAnalyticHalf, DCAnalyticSphere
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .TDEM import hzAnalyticDipoleT
+from .FDEM import hzAnalyticDipoleF
+from .FDEMcasing import *
+from .DC import DCAnalyticHalf, DCAnalyticSphere
+from .FDEMDipolarfields import *

SimPEG/EM/Base.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 
@@ -20,10 +26,10 @@ class BaseEMProblem(Problem.BaseProblem):
         Problem.BaseProblem.__init__(self, mesh, **kwargs)
 
 
-    surveyPair = Survey.BaseSurvey
-    dataPair = Survey.Data
+    surveyPair = Survey.BaseSurvey #: The survey to pair with.
+    dataPair = Survey.Data #: The data to pair with.
 
-    PropMap = EMPropMap
+    PropMap = EMPropMap #: The property mapping
 
     Solver = SimpegSolver
     solverOpts = {}
@@ -217,6 +223,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
     def eval(self, f):
         """
         Project fields to receiver locations
+
         :param Fields u: fields object
         :rtype: numpy.ndarray
         :return: data

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -1,3 +1,10 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
 import numpy as np
 import scipy.sparse as sp
 import SimPEG
@@ -6,11 +13,11 @@ from SimPEG.EM.Utils import omega
 from SimPEG.Utils import Zero, Identity, sdiag
 
 
-class Fields(SimPEG.Problem.Fields):
+class FieldsFDEM(SimPEG.Problem.Fields):
     """
 
     Fancy Field Storage for a FDEM 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
+    each problem, the rest are computed. The fields object acts like an array and is indexed by
 
     .. code-block:: python
 
@@ -42,7 +49,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total electric field
         """
         if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
-            raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting e from %s is not implemented' %list(self.knownFields.keys())[0])
 
         return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
 
@@ -56,24 +63,10 @@ class Fields(SimPEG.Problem.Fields):
         :return: total magnetic flux density
         """
         if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
-            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting b from %s is not implemented' %list(self.knownFields.keys())[0])
 
         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
@@ -84,7 +77,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total magnetic field
         """
         if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
-            raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting h from %s is not implemented' %list(self.knownFields.keys())[0])
 
         return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
 
@@ -98,7 +91,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total current density
         """
         if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
-            raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting j from %s is not implemented' %list(self.knownFields.keys())[0])
 
         return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
 
@@ -106,7 +99,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :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?
@@ -114,7 +107,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -124,7 +117,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         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 SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :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?
@@ -132,32 +125,17 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
-            raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting bDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             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
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :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?
@@ -165,7 +143,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
-            raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting hDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -175,7 +153,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :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?
@@ -183,18 +161,18 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
         return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
 
-class Fields3D_e(Fields):
+class Fields3D_e(FieldsFDEM):
     """
     Fields object for Problem3D_e.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'eSolution':'E'}
@@ -209,9 +187,6 @@ class Fields3D_e(Fields):
                     'h' : ['eSolution','CCV','_h'],
                   }
 
-    def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self, mesh, survey, **kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -317,7 +292,7 @@ class Fields3D_e(Fields):
         C = self._edgeCurl
         b = (C * eSolution)
         for i, src in enumerate(srcList):
-            b[:,i] *= - 1./(1j*omega(src.freq))
+            b[:,i] *= -1./(1j*omega(src.freq))
             s_m, _ = src.eval(self.prob)
             b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
         return b
@@ -377,7 +352,7 @@ class Fields3D_e(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the current density with respect to the field we solved for with a vector
         """
-        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
+        n = int(self._aveE2CCV.shape[0] // self._nC) # number of components (instead of checking if cyl or not)
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         if adjoint:
@@ -414,8 +389,8 @@ class Fields3D_e(Fields):
         :rtype: numpy.ndarray
         :return: magnetic field
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
-        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
+        VI = sdiag(np.kron(np.ones(n), 1. // self.prob.mesh.vol))
 
         return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
 
@@ -429,7 +404,7 @@ class Fields3D_e(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
+        n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             v =  self._MfMui.T * (self._aveF2CCV.T  * (VI.T * du_dm_v))
@@ -446,7 +421,7 @@ class Fields3D_e(Fields):
         :rtype: numpy.ndarray
         :return: product of the magnetic field derivative with respect to the inversion model with a vector
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
+        n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             v =  self._MfMui.T * (self._aveF2CCV.T  * (VI.T * v))
@@ -455,12 +430,12 @@ class Fields3D_e(Fields):
 
 
 
-class Fields3D_b(Fields):
+class Fields3D_b(FieldsFDEM):
     """
     Fields object for Problem3D_b.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'bSolution':'F'}
@@ -475,9 +450,6 @@ class Fields3D_b(Fields):
                     'h' : ['bSolution','CCV','_h'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -500,8 +472,6 @@ 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:
@@ -644,7 +614,7 @@ class Fields3D_b(Fields):
         :return: primary current density
         """
 
-        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
@@ -661,7 +631,7 @@ class Fields3D_b(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the current density with respect to the field we solved for with a vector
         """
-        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
@@ -689,7 +659,7 @@ class Fields3D_b(Fields):
         :rtype: numpy.ndarray
         :return: magnetic field
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
 
@@ -704,7 +674,7 @@ class Fields3D_b(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         if adjoint:
@@ -724,12 +694,12 @@ class Fields3D_b(Fields):
         return Zero()
 
 
-class Fields3D_j(Fields):
+class Fields3D_j(FieldsFDEM):
     """
     Fields object for Problem3D_j.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'jSolution':'F'}
@@ -744,9 +714,6 @@ class Fields3D_j(Fields):
                     'b' : ['jSolution','CCV','_b'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -930,7 +897,7 @@ class Fields3D_j(Fields):
         :rtype: numpy.ndarray
         :return: electric field
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         return VI * (self._aveF2CCV * (self._MfRho *  self._j(jSolution, srcList)))
 
@@ -944,7 +911,7 @@ class Fields3D_j(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the electric field with respect to the field we solved for with a vector
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
@@ -961,7 +928,7 @@ class Fields3D_j(Fields):
         :return: product of the derivative of the electric field with respect to the model with a vector
         """
         jSolution = Utils.mkvc(self[src,'jSolution'])
-        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
@@ -976,7 +943,7 @@ class Fields3D_j(Fields):
         :rtype: numpy.ndarray
         :return: secondary magnetic flux density
         """
-        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
@@ -991,7 +958,7 @@ class Fields3D_j(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         if adjoint:
@@ -1009,7 +976,7 @@ class Fields3D_j(Fields):
         :return: product of the derivative of the magnetic flux density with respect to the model with a vector
         """
         jSolution = self[src,'jSolution']
-        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
 
@@ -1019,12 +986,12 @@ class Fields3D_j(Fields):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields3D_h(Fields):
+class Fields3D_h(FieldsFDEM):
     """
     Fields object for Problem3D_h.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'hSolution':'E'}
@@ -1039,9 +1006,6 @@ class Fields3D_h(Fields):
                     'b' : ['hSolution','CCV','_b'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -1194,7 +1158,7 @@ class Fields3D_h(Fields):
         :rtype: numpy.ndarray
         :return: electric field
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         return VI * (self._aveF2CCV * (self._MfRho *  self._j(hSolution, srcList)))
 
@@ -1208,7 +1172,7 @@ class Fields3D_h(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the electric field with respect to the field we solved for with a vector
         """
-        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
@@ -1225,7 +1189,7 @@ class Fields3D_h(Fields):
         :return: product of the electric field derivative with respect to the inversion model with a vector
         """
         hSolution = Utils.mkvc(self[src,'hSolution'])
-        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
@@ -1241,7 +1205,7 @@ class Fields3D_h(Fields):
         :return: magnetic flux density
         """
         h = self._h(hSolution, srcList)
-        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1256,7 +1220,7 @@ class Fields3D_h(Fields):
         :rtype: numpy.ndarray
         :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
         """
-        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
+        n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -1,7 +1,13 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
-from SurveyFDEM import Survey as SurveyFDEM
-from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
+from .SurveyFDEM import Survey as SurveyFDEM
+from .FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Utils import omega
 
@@ -31,10 +37,11 @@ class BaseFDEMProblem(BaseEMProblem):
         if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
+
     """
 
     surveyPair = SurveyFDEM
-    fieldsPair = Fields
+    fieldsPair = FieldsFDEM
 
     def fields(self, m):
         """
@@ -64,7 +71,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take sensitivity product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object
+        :param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -99,7 +106,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take adjoint product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object
+        :param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -153,8 +160,8 @@ class BaseFDEMProblem(BaseEMProblem):
         Evaluates the sources for a given frequency and puts them in matrix form
 
         :param float freq: Frequency
-        :rtype: (numpy.ndarray, numpy.ndarray)
-        :return: s_m, s_e (nE or nF, nSrc)
+        :rtype: tuple
+        :return: (s_m, s_e) (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._formulation is 'EB':
@@ -194,7 +201,7 @@ class Problem3D_e(BaseFDEMProblem):
 
     which we solve for :math:`\mathbf{e}`.
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'eSolution'
@@ -269,7 +276,7 @@ class Problem3D_e(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -305,7 +312,7 @@ class Problem3D_b(BaseFDEMProblem):
     .. note ::
         The inverse problem will not work with full anisotropy
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'bSolution'
@@ -400,7 +407,7 @@ class Problem3D_b(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -444,6 +451,7 @@ class Problem3D_j(BaseFDEMProblem):
 
         \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j}  + \mathbf{M^e} \mathbf{s_m} \\right)
 
+
     and solve for \\\(\\\mathbf{j}\\\) using
 
     .. math ::
@@ -453,7 +461,7 @@ class Problem3D_j(BaseFDEMProblem):
     .. note::
         This implementation does not yet work with full anisotropy!!
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'jSolution'
@@ -529,8 +537,8 @@ class Problem3D_j(BaseFDEMProblem):
             \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray (nE, nSrc)
-        :return: RHS
+        :rtype: numpy.ndarray
+        :return: RHS (nE, nSrc)
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -549,7 +557,7 @@ class Problem3D_j(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -591,7 +599,7 @@ class Problem3D_h(BaseFDEMProblem):
 
         \\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'hSolution'
@@ -608,9 +616,11 @@ class Problem3D_h(BaseFDEMProblem):
         .. math::
             \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
 
+
         :param float freq: Frequency
         :rtype: scipy.sparse.csr_matrix
         :return: A
+
         """
 
         MeMu = self.MeMu
@@ -653,6 +663,7 @@ class Problem3D_h(BaseFDEMProblem):
         :param float freq: Frequency
         :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
+
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -666,7 +677,7 @@ class Problem3D_h(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray

SimPEG/EM/FDEM/RxFDEM.py

@@ -1,3 +1,10 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from builtins import super
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 from SimPEG import sp
 
@@ -25,10 +32,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
 
     def eval(self, src, mesh, f):
         """
-        Project fields to recievers to get data.
+        Project fields to receivers to get data.
 
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param BaseMesh mesh: mesh used
         :param Fields f: fields object
         :rtype: numpy.ndarray
         :return: fields projected to recievers
@@ -44,8 +51,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param BaseMesh mesh: mesh used
         :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray
@@ -97,19 +104,6 @@ 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):
     """

SimPEG/EM/FDEM/SrcFDEM.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Survey, Problem, Utils, np, sp
 from scipy.constants import mu_0
 from SimPEG.EM.Utils import *
@@ -23,8 +29,8 @@ class BaseSrc(Survey.BaseSrc):
         - :math:`s_m` : magnetic source term
         - :math:`s_e` : electric source term
 
-        :param Problem prob: FDEM Problem
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :param BaseFDEMProblem prob: FDEM Problem
+        :rtype: tuple
         :return: tuple with magnetic source term and electric source term
         """
         s_m = self.s_m(prob)
@@ -37,10 +43,10 @@ class BaseSrc(Survey.BaseSrc):
         - :code:`s_mDeriv` : derivative of the magnetic source term
         - :code:`s_eDeriv` : derivative of the electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :rtype: tuple
         :return: tuple with magnetic source term and electric source term derivatives times a vector
         """
         if v is not None:
@@ -52,7 +58,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary magnetic flux density
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary magnetic flux density
         """
@@ -64,7 +70,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary magnetic field
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -76,7 +82,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary electric field
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary electric field
         """
@@ -88,7 +94,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary current density
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary current density
         """
@@ -100,7 +106,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -110,7 +116,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -120,7 +126,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Derivative of magnetic source term with respect to the inversion model
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -133,7 +139,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Derivative of electric source term with respect to the inversion model
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -162,7 +168,7 @@ class RawVec_e(BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -191,7 +197,7 @@ class RawVec_m(BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -220,7 +226,7 @@ class RawVec(BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -232,7 +238,7 @@ class RawVec(BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -301,7 +307,7 @@ class MagDipole(BaseSrc):
         """
         The primary magnetic flux density from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -339,7 +345,7 @@ class MagDipole(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -350,7 +356,7 @@ class MagDipole(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -364,7 +370,7 @@ class MagDipole(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -416,7 +422,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -455,7 +461,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -466,7 +472,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -479,7 +485,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -530,7 +536,7 @@ class CircularLoop(BaseSrc):
         """
         The primary magnetic flux density from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -555,7 +561,7 @@ class CircularLoop(BaseSrc):
             a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticLoopVectorPotential
+            srcfct = MagneticDipoleVectorPotential
             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)
@@ -567,7 +573,7 @@ class CircularLoop(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -578,7 +584,7 @@ class CircularLoop(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -591,7 +597,7 @@ class CircularLoop(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -1,10 +1,16 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 from SimPEG.EM.Utils import *
 from SimPEG.EM.Base import BaseEMSurvey
 from scipy.constants import mu_0
 from SimPEG.Utils import Zero, Identity
-import SrcFDEM as Src
-import RxFDEM as Rx
+from . import SrcFDEM as Src
+from . import RxFDEM as Rx
 from SimPEG import sp
 
 class Survey(BaseEMSurvey):

SimPEG/EM/FDEM/__init__.py

@@ -1,5 +1,11 @@
-from SurveyFDEM import Survey
-import SrcFDEM as Src
-import RxFDEM as Rx
-from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
-from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .SurveyFDEM import Survey
+from . import SrcFDEM as Src
+from . import RxFDEM as Rx
+from .ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
+from .FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h

SimPEG/EM/Static/DC/BoundaryUtils.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import numpy as np
 
 def getxBCyBC_CC(mesh, alpha, beta, gamma):

SimPEG/EM/Static/DC/FieldsDC.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 from SimPEG.Utils import Identity, Zero
 import numpy as np
@@ -9,7 +15,7 @@ class Fields(SimPEG.Problem.Fields):
 
     def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
-            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -18,7 +24,7 @@ class Fields(SimPEG.Problem.Fields):
 
     def _eDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -26,7 +32,7 @@ class Fields(SimPEG.Problem.Fields):
 
     def _jDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)

SimPEG/EM/Static/DC/FieldsDC_2D.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 from SimPEG.Utils import Identity, Zero
 import numpy as np
@@ -32,7 +38,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
 
     def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
-            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -41,7 +47,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
 
     def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -49,7 +55,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
 
     def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
 
         if adjoint:
             return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)

SimPEG/EM/Static/DC/ProblemDC.py

@@ -1,11 +1,17 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Problem, Utils
 from SimPEG.EM.Base import BaseEMProblem
-from SurveyDC import Survey
-from FieldsDC import Fields, Fields_CC, Fields_N
+from .SurveyDC import Survey
+from .FieldsDC import Fields, Fields_CC, Fields_N
 from SimPEG.Utils import sdiag
 import numpy as np
 from SimPEG.Utils import Zero
-from BoundaryUtils import getxBCyBC_CC
+from .BoundaryUtils import getxBCyBC_CC
 
 class BaseDCProblem(BaseEMProblem):
 

SimPEG/EM/Static/DC/ProblemDC_2D.py

@@ -1,11 +1,18 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import Problem, Utils
 from SimPEG.EM.Base import BaseEMProblem
-from SurveyDC import Survey, Survey_ky
-from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
+from .SurveyDC import Survey, Survey_ky
+from .FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
 from SimPEG.Utils import sdiag
 import numpy as np
 from SimPEG.Utils import Zero
-from BoundaryUtils import getxBCyBC_CC
+from .BoundaryUtils import getxBCyBC_CC
 
 class BaseDCProblem_2D(BaseEMProblem):
 

SimPEG/EM/Static/DC/RxDC.py

@@ -1,3 +1,10 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 import SimPEG
 import numpy as np
 from SimPEG.Utils import Zero, closestPoints

SimPEG/EM/Static/DC/SrcDC.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 # from SimPEG.EM.Base import BaseEMSurvey
 from SimPEG.Utils import Zero, closestPoints, mkvc

SimPEG/EM/Static/DC/SurveyDC.py

@@ -1,9 +1,15 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 from SimPEG.EM.Base import BaseEMSurvey
 from SimPEG import sp, Survey
 from SimPEG.Utils import Zero, Identity
-from RxDC import BaseRx
-from SrcDC import BaseSrc
+from .RxDC import BaseRx
+from .SrcDC import BaseSrc
 
 class Survey(BaseEMSurvey):
     rxPair  = BaseRx

SimPEG/EM/Static/DC/Utils.py

@@ -1,3 +1,10 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 import numpy as np
 
 def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):

SimPEG/EM/Static/DC/__init__.py

@@ -1,8 +1,14 @@
-from ProblemDC import Problem3D_CC, Problem3D_N
-from ProblemDC_2D import Problem2D_CC, Problem2D_N
-from SurveyDC import Survey, Survey_ky
-import SrcDC as Src #Pole
-import RxDC as Rx
-from FieldsDC import Fields_CC
-from BoundaryUtils import getxBCyBC_CC
-import Utils
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .ProblemDC import Problem3D_CC, Problem3D_N
+from .ProblemDC_2D import Problem2D_CC, Problem2D_N
+from .SurveyDC import Survey, Survey_ky
+from . import SrcDC as Src #Pole
+from . import RxDC as Rx
+from .FieldsDC import Fields_CC
+from .BoundaryUtils import getxBCyBC_CC
+from . import Utils

SimPEG/EM/Static/IP/ProblemIP.py

@@ -1,3 +1,9 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Problem, Utils, Maps, Mesh
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -5,7 +11,7 @@ from SimPEG.Utils import sdiag
 import numpy as np
 from SimPEG.Utils import Zero
 from SimPEG.EM.Static.DC import getxBCyBC_CC
-from SurveyIP import Survey
+from .SurveyIP import Survey
 
 class IPPropMap(Maps.PropMap):
     """

SimPEG/EM/Static/IP/SurveyIP.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 from SimPEG.EM.Base import BaseEMSurvey
 from SimPEG import sp, Survey

SimPEG/EM/Static/IP/__init__.py

@@ -1,2 +1,8 @@
-from ProblemIP import Problem3D_CC, Problem3D_N
-from SurveyIP import Survey
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .ProblemIP import Problem3D_CC, Problem3D_N
+from .SurveyIP import Survey

SimPEG/EM/Static/SIP/ProblemSIP.py

@@ -1,3 +1,11 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import Problem, Utils, Maps, Mesh
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -5,7 +13,7 @@ from SimPEG.Utils import sdiag
 import numpy as np
 from SimPEG.Utils import Zero
 from SimPEG.EM.Static.DC import getxBCyBC_CC
-from SurveySIP import Survey, Data
+from .SurveySIP import Survey, Data
 
 class ColeColePropMap(Maps.PropMap):
     """
@@ -105,7 +113,7 @@ class BaseSIPProblem(BaseEMProblem):
         JvAll = []
         #Assume only eta and tau (eta first then tau)
         # v = [2*Mx1]
-        v = v.reshape((int(v.size/2), 2), order='F')
+        v = v.reshape((v.size//2), 2), order='F')
 
         for tind in range(len(self.survey.times)):
             t = self.survey.times[tind]

SimPEG/EM/Static/SIP/Regularization.py

@@ -1,3 +1,10 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import Utils, Maps, Mesh, sp, np
 from SimPEG.Regularization import BaseRegularization, Simple
 

SimPEG/EM/Static/SIP/RxSIP.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 import numpy as np
 from SimPEG.Utils import Zero, closestPoints

SimPEG/EM/Static/SIP/SrcSIP.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG
 # from SimPEG.EM.Base import BaseEMSurvey
 from SimPEG.Utils import Zero, closestPoints, mkvc

SimPEG/EM/Static/SIP/SurveySIP.py

@@ -1,3 +1,10 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import str
 import SimPEG
 from SimPEG.EM.Base import BaseEMSurvey
 from SimPEG import np, sp, Survey, Utils

SimPEG/EM/Static/SIP/__init__.py

@@ -1,5 +1,11 @@
-from ProblemSIP import Problem3D_CC, Problem3D_N
-from SurveySIP import Survey, Data
-import SrcSIP as Src #Pole
-import RxSIP as Rx
-from Regularization import MultiRegularization
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .ProblemSIP import Problem3D_CC, Problem3D_N
+from .SurveySIP import Survey, Data
+from . import SrcSIP as Src #Pole
+from . import RxSIP as Rx
+from .Regularization import MultiRegularization

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -1,3 +1,11 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import np
 from SimPEG.EM.Static import DC, IP
 
@@ -85,10 +93,10 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
 
             elif stype == 'dpdp':
 
-                leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA )
-                LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA ))
+                leg = data * 2*np.pi / (1/MA - 1/MB + 1/NB - 1/NA)
+                LEG.append(1./(2*np.pi) * (1/MA - 1/MB + 1/NB - 1/NA))
             else:
-                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
+                print("""dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """)
                 break
 
 
@@ -103,7 +111,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
                 rho = np.hstack([rho,leg])
 
             else:
-                print """dtype must be 'appr' | 'appc' | 'volt' """
+                print("""dtype must be 'appr' | 'appc' | 'volt' """)
                 break
 
 
@@ -184,14 +192,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     # Mesure survey length and direction
     dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
 
-    dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
-    dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
+    dl_x = (endl[1,0] - endl[0,0]) / dl_len
+    dl_y = (endl[1,1] - endl[0,1]) / dl_len
 
-    nstn = np.floor( dl_len / a )
+    nstn = np.floor(dl_len / a)
 
     # Compute discrete pole location along line
-    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
-    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
+    stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*a
+    stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
 
     if mesh.dim==2:
         ztop = mesh.vectorNy[-1]
@@ -230,15 +238,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
             AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
 
             # Number of receivers to fit
-            nstn = np.min([np.floor( (AB - b) / a ) , n])
+            nstn = np.min([(AB - b) // a, n])
 
             # Check if there is enough space, else break the loop
             if nstn <= 0:
                 continue
 
             # Compute discrete pole location along line
-            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
-            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
+            stn_x = N[ii,0] + dl_x*b + np.array(list(range(int(nstn))))*dl_x*a
+            stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
 
             # Create receiver poles
 
@@ -275,17 +283,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
         max_y = endl[1,1] - dl_y * b
 
         box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
-        box_w = box_l/2.
+        box_w = box_l / 2.
 
-        nstn = np.floor( box_l / a )
+        nstn = np.floor(box_l / a)
 
         # Compute discrete pole location along line
-        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
-        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
+        stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*a
+        stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
 
         # Define number of cross lines
-        nlin = int(np.floor( box_w / a ))
-        lind = range(-nlin,nlin+1)
+        nlin = int(box_w // a)
+        lind = list(range(-nlin,nlin+1))
 
         ngrad = nstn * len(lind)
 
@@ -310,7 +318,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
             srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
         SrcList.append(srcClass)
     else:
-        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
+        print("""stype must be either 'pdp', 'dpdp' or 'gradient'. """)
 
 
     return SrcList

SimPEG/EM/Static/Utils/__init__.py

@@ -1 +1,7 @@
-from StaticUtils import *
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .StaticUtils import *

SimPEG/EM/Static/__init__.py

@@ -1,3 +1,9 @@
-import DC
-import IP
-import SIP
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from . import DC
+from . import IP
+from . import SIP

SimPEG/EM/TDEM/BaseTDEM.py

@@ -1,3 +1,10 @@
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import Solver, Problem
 from SimPEG.Problem import BaseTimeProblem
 from SimPEG.EM.Utils import *
@@ -47,7 +54,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         self.waveformType = "GENERAL"
 
     def fields(self, m):
-        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print('%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50))
         self.curModel = m
         # Create a fields storage object
         F = self._FieldsForward_pair(self.mesh, self.survey)
@@ -55,7 +62,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             # Set the initial conditions
             F[src,:,0] = src.getInitialFields(self.mesh)
         F = self.forward(m, self.getRHS, F=F)
-        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print('%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50))
         return F
 
     def forward(self, m, RHS, F=None):
@@ -70,13 +77,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
                 if Ainv is not None:
                     Ainv.clean()
                 A = self.getA(tInd)
-                if self.verbose: print 'Factoring...   (dt = %e)'%dt
+                if self.verbose: print('Factoring...   (dt = %e)'%dt)
                 Ainv = self.Solver(A, **self.solverOpts)
-                if self.verbose: print 'Done'
+                if self.verbose: print('Done')
             rhs = RHS(tInd, F)
-            if self.verbose: print '    Solving...   (tInd = %d)'%tInd
+            if self.verbose: print('    Solving...   (tInd = %d)'%tInd)
             sol = Ainv * rhs
-            if self.verbose: print '    Done...'
+            if self.verbose: print('    Done...')
             if sol.ndim == 1:
                 sol.shape = (sol.size,1)
             F[:,self.solType,tInd+1] = sol
@@ -95,13 +102,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
                 if Ainv is not None:
                     Ainv.clean()
                 A = self.getA(tInd)
-                if self.verbose: print 'Factoring (Adjoint)...   (dt = %e)'%dt
+                if self.verbose: print('Factoring (Adjoint)...   (dt = %e)'%dt)
                 Ainv = self.Solver(A, **self.solverOpts)
-                if self.verbose: print 'Done'
+                if self.verbose: print('Done')
             rhs = RHS(tInd, F)
-            if self.verbose: print '    Solving (Adjoint)...   (tInd = %d)'%tInd
+            if self.verbose: print('    Solving (Adjoint)...   (tInd = %d)'%tInd)
             sol = Ainv * rhs
-            if self.verbose: print '    Done...'
+            if self.verbose: print('    Done...')
             if sol.ndim == 1:
                 sol.shape = (sol.size,1)
             F[:,self.solType,tInd+1] = sol
@@ -112,7 +119,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
+            :param FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -123,21 +130,21 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
 
         """
-        if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print('%s\nCalculating J(v)\n%s'%('*'*50,'*'*50))
         self.curModel = m
         if f is None:
             f = self.fields(m)
         p = self.Gvec(m, v, f)
         y = self.solveAh(m, p)
         Jv = self.survey.evalDeriv(f, v=y)
-        if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print('%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50))
         return - mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
-            :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param numpy.ndarray v: vector (or a :class:`SimPEG.Survey.Data` object)
+            :param FieldsTDEM u: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (model object)
 
@@ -148,7 +155,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
 
         """
-        if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print('%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50))
         self.curModel = m
         if f is None:
             f = self.fields(m)
@@ -159,6 +166,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         p = self.survey.evalDeriv(f, v=v, adjoint=True)
         y = self.solveAht(m, p)
         w = self.Gtvec(m, y, f)
-        if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print('%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50))
         return - mkvc(w)
 

SimPEG/EM/TDEM/FieldsTDEM.py

@@ -1,142 +0,0 @@
-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)

SimPEG/EM/TDEM/SrcTDEM.py

@@ -1,246 +0,0 @@
-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
-

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -1,16 +1,16 @@
-import SimPEG
-from SimPEG import np, Utils
-from SimPEG.Utils import Zero, Identity
-from scipy.constants import mu_0
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from SimPEG import Utils, Survey, np
+from SimPEG.Survey import BaseSurvey
 from SimPEG.EM.Utils import *
-import SrcTDEM as Src
+from .BaseTDEM import FieldsTDEM
 
 
-####################################################
-# Receivers
-####################################################
-
-class Rx(SimPEG.Survey.BaseTimeRx):
+class RxTDEM(Survey.BaseTimeRx):
 
     knownRxTypes = {
                     'ex':['e', 'Ex', 'N'],
@@ -27,7 +27,7 @@ class Rx(SimPEG.Survey.BaseTimeRx):
                    }
 
     def __init__(self, locs, times, rxType):
-        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
 
     @property
     def projField(self):
@@ -62,60 +62,144 @@ class Rx(SimPEG.Survey.BaseTimeRx):
         u_part = Utils.mkvc(u[src, self.projField, :])
         return P*u_part
 
-    def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False):
+    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
         P = self.getP(mesh, timeMesh)
+
         if not adjoint:
-            return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :])
+            return P * Utils.mkvc(v[src, self.projField, :])
         elif adjoint:
-            # 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')
+            return P.T * v[src, self]
 
 
-####################################################
-# Survey
-####################################################
+class SrcTDEM(Survey.BaseSrc):
+    rxPair = RxTDEM
+    radius = None
 
-class Survey(SimPEG.Survey.BaseSurvey):
+    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):
     """
-    Time domain electromagnetic survey
+        docstring for SurveyTDEM
     """
-
-    srcPair = Src.BaseSrc
-    rxPair = Rx
+    srcPair = SrcTDEM
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
+        Survey.BaseSurvey.__init__(self, **kwargs)
 
     def eval(self, u):
-        data = SimPEG.Survey.Data(self)
+        data = 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):
-        raise Exception('Use Receivers to project fields deriv.')
-        # assert v is not None, 'v to multiply must be provided.'
+        assert v is not None, 'v to multiply must be provided.'
 
-        # 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
+        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
 
 

SimPEG/EM/TDEM/TDEM.py

@@ -1,553 +0,0 @@
-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()
-
-
-
-
-

SimPEG/EM/TDEM/TDEM_b.py

@@ -1,7 +1,14 @@
-from BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
+from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
 from SimPEG.Utils import mkvc, sdiag
 import numpy as np
-from SurveyTDEM import SurveyTDEM
+from .SurveyTDEM import SurveyTDEM
 
 
 class FieldsTDEM_e_from_b(FieldsTDEM):
@@ -87,8 +94,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a model)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM u: Fields resulting from m
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply G by a vector
@@ -125,9 +132,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a fields)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
-            :rtype: np.ndarray (like a model)
-            :return: p
+            :param FieldsTDEM u: Fields resulting from m
+            :rtype: numpy.ndarray
+            :return: p (like a model)
 
             Multiply G.T by a vector
         """
@@ -153,8 +160,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAh(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM p: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM p: Fields object
+            :rtype: FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
@@ -200,8 +207,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAht(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM p: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM p: Fields object
+            :rtype: FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
@@ -270,8 +277,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM vec: Fields object
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
@@ -315,8 +322,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhtVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM vec: Fields object
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where

SimPEG/EM/TDEM/__init__.py

@@ -1,3 +1,9 @@
-from TDEM import BaseTDEMProblem, Problem_b, Problem_e
-from FieldsTDEM import Fields, Fields_b
-from SurveyTDEM import Survey, Src, Rx
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
+from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from .TDEM_b import ProblemTDEM_b

SimPEG/EM/TDEM_old/SurveyTDEM.py

@@ -1,199 +0,0 @@
-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
-
-

SimPEG/EM/TDEM_old/__init__.py

@@ -1,3 +0,0 @@
-from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
-from BaseTDEM import BaseTDEMProblem, FieldsTDEM
-from TDEM_b import ProblemTDEM_b

SimPEG/EM/Utils/AnalyticUtils.py

@@ -1,3 +1,10 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
 from scipy.special import ellipk, ellipe
 from scipy.constants import mu_0, pi
@@ -17,7 +24,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
     #TODO: break this out!
 
     if type(component) in [list, tuple]:
-        out = range(len(component))
+        out = list(range(len(component)))
         for i, comp in enumerate(component):
             out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
         return np.concatenate(out)
@@ -118,7 +125,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
     """
 
     if type(component) in [list, tuple]:
-        out = range(len(component))
+        out = list(range(len(component)))
         for i, comp in enumerate(component):
             out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
         return np.concatenate(out)
@@ -158,11 +165,11 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
             # % 1/r singular at r = 0 and K(m) singular at m = 1
             Aphi = np.zeros(n)
             # % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
-            Aphi[ind] = 4e-7 / np.sqrt(m[ind])  * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
+            Aphi[ind] = 4e-7 / np.sqrt(m[ind])  * np.sqrt(radius/ r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
             if component == 'x':
-                A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
+                A[ind, i] = Aphi[ind] * (-y[ind] / r[ind])
             elif component == 'y':
-                A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
+                A[ind, i] = Aphi[ind] * (x[ind] / r[ind])
             else:
                 raise ValueError('Invalid component')
 

SimPEG/EM/Utils/EMUtils.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import numpy as np
 from scipy.constants import mu_0, epsilon_0
 
@@ -9,8 +15,8 @@ def omega(freq):
 def k(freq, sigma, mu=mu_0, eps=epsilon_0):
     """ Eq 1.47 - 1.49 in Ward and Hohmann """
     w = omega(freq)
-    alp  = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) + 1) ) 
-    beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) 
+    alp  = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) + 1)
+    beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) - 1)
     return alp - 1j*beta
 
 

SimPEG/EM/Utils/__init__.py

@@ -1,2 +1,8 @@
-from EMUtils import omega, k
-from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .EMUtils import omega, k
+from .AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential

SimPEG/EM/Utils/testingUtils.py

@@ -1,3 +1,10 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
 import unittest
 from SimPEG import *
 from SimPEG import EM
@@ -58,7 +65,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
                 Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
 
     if verbose:
-        print '  Fetching %s problem' % (fdemType)
+        print('  Fetching %s problem' % (fdemType))
 
     if fdemType == 'e':
         survey = EM.FDEM.Survey(Src)
@@ -83,7 +90,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     try:
         from pymatsolver import MumpsSolver
         prb.Solver = MumpsSolver
-    except ImportError, e:
+    except ImportError as e:
         prb.Solver = SolverLU
 
     return prb
@@ -94,7 +101,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
 
     prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
-    print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
+    print('Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp))
 
     logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
     mu = np.ones(mesh.nC)*MU
@@ -112,7 +119,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
     d1 = survey1.dpred(m)
 
     if verbose:
-        print '  Problem 1 solved'
+        print('  Problem 1 solved')
 
 
     prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
@@ -121,11 +128,11 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
     d2 = survey2.dpred(m)
 
     if verbose:
-        print '  Problem 2 solved'
+        print('  Problem 2 solved')
 
     r = d2-d1
     l2r = l2norm(r)
 
     tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
-    print l2norm(d1), l2norm(d2),  l2r , tol, l2r < tol
+    print(l2norm(d1), l2norm(d2),  l2r , tol, l2r < tol)
     return l2r < tol

SimPEG/EM/__init__.py

@@ -1,7 +1,13 @@
-import TDEM
-import FDEM
-import Static
-import Base
-import Analytics
-import Utils
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from . import TDEM
+from . import FDEM
+from . import Static
+from . import Base
+from . import Analytics
+from . import Utils
 from scipy.constants import mu_0, epsilon_0

SimPEG/Examples/DC_Analytic_Dipole.py

@@ -1,7 +1,13 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
-import SimPEG.DCIP as DC
+import SimPEG.EM.Static.DC as DC
 
-def run(plotIt=False):
+def run(plotIt=True):
     cs = 25.
     hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
     hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -21,15 +27,15 @@ def run(plotIt=False):
     #     ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
     #     ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
 
-    rx = DC.RxDipole(xyz_rxP, xyz_rxN)
-    src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
-    survey = DC.SurveyDC([src])
-    problem = DC.ProblemDC_CC(mesh)
+    rx = DC.Rx.Dipole(xyz_rxP, xyz_rxN)
+    src = DC.Src.Dipole([rx], np.r_[-200, 0, -12.5], np.r_[+200, 0, -12.5])
+    survey = DC.Survey([src])
+    problem = DC.Problem3D_CC(mesh)
     problem.pair(survey)
     try:
         from pymatsolver import MumpsSolver
         problem.Solver = MumpsSolver
-    except Exception, e:
+    except Exception as e:
         pass
     data = survey.dpred(sigma)
 
@@ -61,8 +67,8 @@ def run(plotIt=False):
         ax[0].set_title('Computed')
         plt.show()
 
-    return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
+    return np.linalg.norm(data-data_ana) / np.linalg.norm(data_ana)
 
 
 if __name__ == '__main__':
-    print run(plotIt=True)
+    print(run())

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -1,3 +1,11 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import Mesh, Utils, np, sp
 import SimPEG.DCIP as DC
 import time
@@ -57,7 +65,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
     model[ind] = sig[2]
 
     # Get index of the center
-    indy = int(mesh.nCy/2)
+    indy = int(mesh.nCy // 2)
 
     # Plot the model for reference
     # Define core mesh extent
@@ -124,10 +132,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
             tx =  np.squeeze(Tx[ii][:,0:1])
             tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
             inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
-            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )
+            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1] / mesh.vol[inds])
         else:
             inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
-            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
+            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1,1] / mesh.vol[inds])
 
         # Iterative Solve
         Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
@@ -143,10 +151,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
         dtemp = (P1*phi - P2*phi)*np.pi
 
         data.append( dtemp )
-        print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
+        print('\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time), end=' ')
 
-    print 'Transmitter {0} of {1}'.format(ii,len(Tx))
-    print 'Forward completed'
+    print('Transmitter {0} of {1}'.format(ii,len(Tx)))
+    print('Forward completed')
 
     # Let's just convert the 3D format into 2D (distance along line) and plot
     survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 import SimPEG.EM as EM
 from SimPEG.EM import mu_0
@@ -56,7 +62,7 @@ def run(plotIt=True):
     try:
         from pymatsolver import MumpsSolver
         prb.Solver = MumpsSolver
-    except ImportError, e:
+    except ImportError as e:
         prb.Solver = SolverLU
 
     prb.pair(survey)

SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 import SimPEG.EM as EM
 

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -1,3 +1,10 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
 from SimPEG.EM import FDEM, Analytics, mu_0
 import time
@@ -19,10 +26,13 @@ def run(plotIt=True):
         Morrison Casing Model, and the results are used in a 2016 SEG abstract by
         Yang et al.
 
-        - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
+        .. code-block:: text
+
+            Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
 
 
         The model consists of:
+
         - Air: Conductivity 1e-8 S/m, above z = 0
         - Background: conductivity 1e-2 S/m, below z = 0
         - Casing: conductivity 1e6 S/m
@@ -64,8 +74,8 @@ def run(plotIt=True):
     casing_l = 300   # length of the casing
 
     casing_r = 0.1
-    casing_a = casing_r - casing_t/2. # inner radius
-    casing_b = casing_r + casing_t/2. # outer radius
+    casing_a = casing_r - casing_t / 2. # inner radius
+    casing_b = casing_r + casing_t / 2. # outer radius
     casing_z = np.r_[-casing_l,0.]
 
 
@@ -75,25 +85,25 @@ def run(plotIt=True):
     src_loc = np.r_[0.,0.,dsz]
     inf_loc = np.r_[0.,0.,1e4]
 
-    print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
+    print('Skin Depth: ', [(500. / np.sqrt(sigmaback*_)) for _ in freqs])
 
 
     # ------------------ MESH ------------------
     # fine cells near well bore
     csx1, csx2 = 2e-3, 60.
     pfx1, pfx2 = 1.3, 1.3
-    ncx1 = np.ceil(casing_b/csx1+2)
+    ncx1 = np.ceil(casing_b/csx1)+2
 
     # pad nicely to second cell size
-    npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
+    npadx1 = np.log(csx2/csx1) // np.log(pfx1)
     hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
     dx1 = sum(hx1a)+sum(hx1b)
-    dx1 = np.floor(dx1/csx2)
-    hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
+    dx1 = dx1 // csx2
+    hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
 
     # second chunk of mesh
     dx2 = 300. # uniform mesh out to here
-    ncx2 = np.ceil((dx2 - dx1)/csx2)
+    ncx2 = np.ceil((dx2 - dx1) / csx2)
     npadx2 = 45
     hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
     hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
@@ -107,8 +117,8 @@ def run(plotIt=True):
     # Mesh
     mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
 
-    print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
-    print 'Number of cells', mesh.nC
+    print('Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()))
+    print('Number of cells', mesh.nC)
 
     if plotIt is True:
         fig, ax = plt.subplots(1, 1, figsize=(6, 4))
@@ -215,13 +225,13 @@ def run(plotIt=True):
     # ------------ Problem and Survey ---------------
     survey = FDEM.Survey(sg_p + dg_p)
     mapping = [('sigma', Maps.IdentityMap(mesh))]
-    problem = FDEM.Problem3D_h(mesh, mapping=mapping)
+    problem = FDEM.Problem3D_h(mesh, mapping=mapping, Solver=solver)
     problem.pair(survey)
 
     # ------------- Solve ---------------------------
     t0 = time.time()
     fieldsCasing = problem.fields(sigCasing)
-    print 'Time to solve 2 sources', time.time() - t0
+    print('Time to solve 2 sources', time.time() - t0)
 
     # Plot current
 
@@ -248,9 +258,9 @@ def run(plotIt=True):
     in1_in = in1[np.r_[inds]]
     z_in = mesh.gridFz[inds_fz,2]
 
-    in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
-    in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
-    z_in = z_in.reshape([z_in.shape[0]/3,3])
+    in0_in = in0_in.reshape([in0_in.shape[0]//3,3])
+    in1_in = in1_in.reshape([in1_in.shape[0]//3,3])
+    z_in = z_in.reshape([z_in.shape[0]//3,3])
 
     I0 = in0_in.sum(1).real
     I1 = in1_in.sum(1).real

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 import SimPEG.EM as EM
 from SimPEG.EM import mu_0
@@ -42,10 +48,10 @@ def run(plotIt=True):
 
 
     rxOffset=1e-3
-    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)
+    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)
 
     prb.Solver = SolverLU
     prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
@@ -53,9 +59,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:

SimPEG/Examples/FLOW_Richards_1D_Celia1990.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 from SimPEG.FLOW import Richards
 

SimPEG/Examples/Inversion_IRLS.py

@@ -1,3 +1,11 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import str
+from builtins import range
 from SimPEG import *
 
 
@@ -42,62 +50,40 @@ def run(N=100, plotIt=True):
     survey = Survey.LinearSurvey()
     survey.pair(prob)
     survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
-    #survey.makeSyntheticData(mtrue, std=std_noise)
 
     wd = np.ones(nk) * std_noise
 
-    #print survey.std[0]
-    #M = prob.mesh
     # Distance weighting
     wr = np.sum(prob.G**2.,axis=0)**0.5
-    wr = ( wr/np.max(wr) )
+    wr = ( wr/np.max(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=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
 
     reg = Regularization.Sparse(mesh)
     reg.mref = mref
     reg.cell_weights = wr
 
     reg.mref = np.zeros(mesh.nC)
-    eps_p = 5e-2
-    eps_q = 5e-2
-    norms   = [0., 0., 2., 2.]
+
 
     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)
+
+    # Set the IRLS directive, penalize the lowest 25 percentile of model values
+    # Start with an l2-l2, then switch to lp-norms
+    norms   = [0., 0., 2., 2.]
+    IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
 
     inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 
     # Run inversion
     mrec = inv.run(m0)
 
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
+    print("Final misfit:" + str(invProb.dmisfit.eval(mrec)))
 
 
     if plotIt:

SimPEG/Examples/Inversion_Linear.py

@@ -1,3 +1,10 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
 
 

SimPEG/Examples/MT_1D_ForwardAndInversion.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 import SimPEG as simpeg
 import numpy as np
 import SimPEG.MT as MT
@@ -7,7 +13,7 @@ import matplotlib.pyplot as plt
 def run(plotIt=True):
     """
         MT: 1D: Inversion
-        =======================
+        =================
 
         Forward model 1D MT data.
         Setup and run a MT 1D inversion.
@@ -50,7 +56,7 @@ def run(plotIt=True):
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
-    actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
+    actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
     mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
 
     ## Setup the layout of the survey, set the sources and the connected receivers
@@ -76,7 +82,7 @@ def run(plotIt=True):
     survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
         fig.suptitle('Target - smooth true')
 
 

SimPEG/Examples/MT_3D_Foward.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 # Test script to use SimPEG.MT platform to forward model synthetic data.
 
 # Import
@@ -12,7 +18,7 @@ except:
 def run(plotIt=True, nFreq=1):
     """
         MT: 3D: Forward
-        =======================
+        ===============
 
         Forward model 3D MT data.
 
@@ -46,16 +52,15 @@ def run(plotIt=True, nFreq=1):
     survey = MT.Survey(srcList)
 
     ## Setup the problem object
-    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
+    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG, Solver=Solver)
     problem.pair(survey)
-    problem.Solver = Solver
 
     # Calculate the data
     fields = problem.fields(sig)
     dataVec = survey.eval(fields)
 
     # Make the data
-    mtData = MT.Data(survey,dataVec)
+    mtData = MT.Data(survey, dataVec)
     # Add plots
     if plotIt:
         pass

SimPEG/Examples/Maps_ComboMaps.py

@@ -0,0 +1,69 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from builtins import dict
+from future import standard_library
+standard_library.install_aliases()
+from SimPEG import Mesh, Maps, np
+
+def run(plotIt=True):
+    """
+
+        Maps: ComboMaps
+        ===============
+
+        We will use an example where we want a 1D layered earth as
+        our model, but we want to map this to a 2D discretization to do our forward
+        modeling. We will also assume that we are working in log conductivity still,
+        so after the transformation we want to map to conductivity space.
+        To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
+        which does the first part of what we just described. The second part will be
+        done by the :class:`SimPEG.Maps.ExpMap` described above.
+
+        .. code-block:: python
+            :linenos:
+
+            M = Mesh.TensorMesh([7,5])
+            v1dMap = Maps.SurjectVertical1D(M)
+            expMap = Maps.ExpMap(M)
+            myMap = expMap * v1dMap
+            m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+            sig = myMap * m
+
+        If you noticed, it was pretty easy to combine maps. What is even cooler is
+        that the derivatives also are made for you (if everything goes right).
+        Just to be sure that the derivative is correct, you should always run the test
+        on the mapping that you create.
+
+    """
+
+
+    M = Mesh.TensorMesh([7,5])
+    v1dMap = Maps.SurjectVertical1D(M)
+    expMap = Maps.ExpMap(M)
+    myMap = expMap * v1dMap
+    m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+    sig = myMap * m
+
+    if not plotIt: return
+
+    import matplotlib.pyplot as plt
+    figs, axs = plt.subplots(1,2)
+    axs[0].plot(m, M.vectorCCy, 'b-o')
+    axs[0].set_title('Model')
+    axs[0].set_ylabel('Depth, y')
+    axs[0].set_xlabel('Value, $m_i$')
+    axs[0].set_xlim(0,3)
+    axs[0].set_ylim(0,1)
+    clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
+    axs[1].set_title('Physical Property')
+    axs[1].set_ylabel('Depth, y')
+    clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
+    plt.tight_layout()
+    plt.show()
+
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/Maps_Mesh2Mesh.py

@@ -0,0 +1,47 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from SimPEG import Mesh, Maps, Utils
+
+def run(plotIt=True):
+    """
+
+        Maps: Mesh2Mesh
+        ===============
+
+        This mapping allows you to go from one mesh to another.
+
+    """
+
+    M = Mesh.TensorMesh([100,100])
+    h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
+    h1 = h1/h1.sum()
+    M2 = Mesh.TensorMesh([h1,h1])
+    V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
+    v = Utils.mkvc(V)
+    modh = Maps.Mesh2Mesh([M,M2])
+    modH = Maps.Mesh2Mesh([M2,M])
+    H = modH * v
+    h = modh * H
+
+    if not plotIt: return
+
+    import matplotlib.pyplot as plt
+    ax = plt.subplot(131)
+    M.plotImage(v, ax=ax)
+    ax.set_title('Fine Mesh (Original)')
+    ax = plt.subplot(132)
+    M2.plotImage(H,clim=[0,1],ax=ax)
+    ax.set_title('Course Mesh')
+    ax = plt.subplot(133)
+    M.plotImage(h,clim=[0,1],ax=ax)
+    ax.set_title('Fine Mesh (Interpolated)')
+    plt.show()
+
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/Mesh_Basic_ForwardDC.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Mesh, Utils, np, SolverLU
 
 def run(plotIt=True):

SimPEG/Examples/Mesh_Basic_PlotImage.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 
 def run(plotIt=True):

SimPEG/Examples/Mesh_Basic_Types.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 
 def run(plotIt=True):

SimPEG/Examples/Mesh_Operators_CahnHilliard.py

@@ -1,3 +1,11 @@
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+from __future__ import absolute_import
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+from builtins import zip
 from SimPEG import *
 
 def run(plotIt=True, n=60):
@@ -87,7 +95,7 @@ def run(plotIt=True, n=60):
         if elapsed > capture[jj]:
             PHIS += [(elapsed, phi.copy())]
             jj += 1
-        if ii % 10 == 0: print ii, elapsed
+        if ii % 10 == 0: print(ii, elapsed)
         ii += 1
 
     if plotIt:

SimPEG/Examples/Mesh_QuadTree_Creation.py

@@ -1,3 +1,10 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
 
 def run(plotIt=True):

SimPEG/Examples/Mesh_QuadTree_FaceDiv.py

@@ -1,3 +1,11 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import zip
+from builtins import range
 from SimPEG import *
 
 def run(plotIt=True, n=60):
@@ -28,15 +36,15 @@ def run(plotIt=True, n=60):
         axes[0].set_xlim([-1,17])
         axes[0].set_ylim([-1,17])
 
-        for ii, loc in zip(range(M.nC),M.gridCC):
+        for ii, loc in zip(list(range(M.nC)),M.gridCC):
             axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
 
         axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
-        for ii, loc in zip(range(M.nFx),M.gridFx):
+        for ii, loc in zip(list(range(M.nFx)),M.gridFx):
             axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
 
         axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
-        for ii, loc in zip(range(M.nFy),M.gridFy):
+        for ii, loc in zip(list(range(M.nFy)),M.gridFy):
             axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
 
         axes[1].spy(M.faceDiv)

SimPEG/Examples/Mesh_QuadTree_HangingNodes.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 
 def run(plotIt=True):

SimPEG/Examples/Mesh_Tensor_Creation.py

@@ -1,3 +1,9 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 
 def run(plotIt=True):

SimPEG/Examples/Utils_surface2ind_topo.py

@@ -1,9 +1,19 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import *
 from SimPEG.Utils import surface2ind_topo
 
 
-def run(plotIt=False, nx = 5, ny = 5):
+def run(plotIt=True, nx=5, ny=5):
     """
+
+        Utils: surface2ind_topo
+        =======================
+
         Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
         a topographic surface.
 
@@ -13,27 +23,25 @@ def run(plotIt=False, nx = 5, ny = 5):
     xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max())
     topo = 0.4*np.sin(xtopo*5) # define a topographic surface
 
-    Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array
+    Topo = np.hstack([Utils.mkvc(xtopo,2), Utils.mkvc(topo,2)]) #make it an array
 
-    indcc = surface2ind_topo(mesh, Topo,'CC')
+    indcc = surface2ind_topo(mesh, Topo, 'CC')
 
     if plotIt:
         from matplotlib.pylab import plt
         from scipy.interpolate import interp1d
-        fig, ax = plt.subplots(1,1,figsize=(6,6))
+        fig, ax = plt.subplots(1,1, figsize=(6,6))
         mesh.plotGrid(ax=ax, nodes=True, centers=True)
         ax.plot(xtopo,topo,'k',linewidth=1)
-        # ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3)
         ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3)
 
-
         aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T
         a = aveN2CC * indcc
         a[a > 0] = 1.
         a[a < 0.25] = np.nan
         a = a.reshape(mesh.vnN, order='F')
         masked_array = np.ma.array(a, mask=np.isnan(a))
-        ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2)
+        ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap=plt.cm.gray, alpha=0.2)
         plt.show()
 
 

SimPEG/Examples/__init__.py

@@ -1,28 +1,37 @@
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import division
+from builtins import open
+from future import standard_library
+standard_library.install_aliases()
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import DC_Analytic_Dipole
-import DC_Forward_PseudoSection
-import EM_FDEM_1D_Inversion
-import EM_FDEM_Analytic_MagDipoleWholespace
-import EM_Schenkel_Morrison_Casing
-import EM_TDEM_1D_Inversion
-import FLOW_Richards_1D_Celia1990
-import Inversion_IRLS
-import Inversion_Linear
-import Mesh_Basic_ForwardDC
-import Mesh_Basic_PlotImage
-import Mesh_Basic_Types
-import Mesh_Operators_CahnHilliard
-import Mesh_QuadTree_Creation
-import Mesh_QuadTree_FaceDiv
-import Mesh_QuadTree_HangingNodes
-import Mesh_Tensor_Creation
-import MT_1D_ForwardAndInversion
-import MT_3D_Foward
-import Utils_surface2ind_topo
+from . import DC_Analytic_Dipole
+from . import DC_Forward_PseudoSection
+from . import EM_FDEM_1D_Inversion
+from . import EM_FDEM_Analytic_MagDipoleWholespace
+from . import EM_Schenkel_Morrison_Casing
+from . import EM_TDEM_1D_Inversion
+from . import FLOW_Richards_1D_Celia1990
+from . import Inversion_IRLS
+from . import Inversion_Linear
+from . import Maps_ComboMaps
+from . import Maps_Mesh2Mesh
+from . import Mesh_Basic_ForwardDC
+from . import Mesh_Basic_PlotImage
+from . import Mesh_Basic_Types
+from . import Mesh_Operators_CahnHilliard
+from . import Mesh_QuadTree_Creation
+from . import Mesh_QuadTree_FaceDiv
+from . import Mesh_QuadTree_HangingNodes
+from . import Mesh_Tensor_Creation
+from . import MT_1D_ForwardAndInversion
+from . import MT_3D_Foward
+from . 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", "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"]
+__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", "Maps_ComboMaps", "Maps_Mesh2Mesh", "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 #####
 
@@ -38,7 +47,7 @@ if __name__ == '__main__':
 
     # Create the examples dir in the docs folder.
     fName = os.path.realpath(__file__)
-    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
+    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'content', 'examples'])
     shutil.rmtree(docExamplesDir)
     os.makedirs(docExamplesDir)
 
@@ -95,14 +104,14 @@ if __name__ == '__main__':
     from SimPEG import Examples
     Examples.%s.run()
 
-.. literalinclude:: ../../SimPEG/Examples/%s.py
+.. literalinclude:: ../../../SimPEG/Examples/%s.py
     :language: python
     :linenos:
 """%(name,doc,name,name)
 
-        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
+        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
 
-        print 'Creating: %s.rst'%name
+        print('Creating: %s.rst'%name)
         f = open(rst, 'w')
         f.write(out)
         f.close()

SimPEG/Examples/simpegEMpaper_example_1.py

@@ -1,157 +0,0 @@
-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()

SimPEG/FLOW/Richards/Empirical.py

@@ -1,14 +1,20 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import object
 from SimPEG import Mesh, Maps, Utils, np
+from future.utils import with_metaclass
 
 
-class NonLinearMap(object):
+class NonLinearMap(with_metaclass(Utils.SimPEGMetaClass, object)):
     """
     SimPEG NonLinearMap
 
     """
 
-    __metaclass__ = Utils.SimPEGMetaClass
-
     counter = None   #: A SimPEG.Utils.Counter object
     mesh = None      #: A SimPEG Mesh
 
@@ -31,7 +37,7 @@ class NonLinearMap(object):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -44,7 +50,7 @@ class NonLinearMap(object):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -186,7 +192,7 @@ class _haverkamp_theta(NonLinearMap):
 
     def transformDerivU(self, u, m):
         self.setModel(m)
-        g = (self.alpha*((self.theta_s - self.theta_r)/
+        g = (self.alpha*((self.theta_s - self.theta_r) /
              (self.alpha + abs(u)**self.beta)**2)
              *(-self.beta*abs(u)**(self.beta-1)*np.sign(u)))
         g[u >= 0] = 0
@@ -273,7 +279,7 @@ class _vangenuchten_theta(NonLinearMap):
     def transform(self, u, m):
         self.setModel(m)
         m = 1 - 1.0/self.n
-        f = ((  self.theta_s  -  self.theta_r  )/
+        f = ((  self.theta_s  -  self.theta_r  ) /
              ((1+abs(self.alpha*u)**self.n)**m)   +  self.theta_r)
         if Utils.isScalar(self.theta_s):
             f[u >= 0] = self.theta_s
@@ -343,7 +349,7 @@ class _vangenuchten_k(NonLinearMap):
         Ks = self.Ks
         m = 1.0 - 1.0/n
 
-        g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))
+        g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1))+ 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)
         g[u >= 0] = 0
         g = Utils.sdiag(g)
         return g

SimPEG/FLOW/Richards/RichardsProblem.py

@@ -1,5 +1,12 @@
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
 from SimPEG import *
-from Empirical import RichardsMap
+from .Empirical import RichardsMap
 import time
 
 
@@ -61,7 +68,7 @@ class RichardsSurvey(Survey.BaseSurvey):
 
     @Utils.requires('prob')
     def eval(self, U, m):
-        Ds = range(len(self.rxList))
+        Ds = list(range(len(self.rxList)))
         for ii, rx in enumerate(self.rxList):
             Ds[ii] = rx.eval(U, m,
                                 self.prob.mapping,
@@ -73,7 +80,7 @@ class RichardsSurvey(Survey.BaseSurvey):
     @Utils.requires('prob')
     def evalDeriv(self, U, m):
         """The Derivative with respect to the fields."""
-        Ds = range(len(self.rxList))
+        Ds = list(range(len(self.rxList)))
         for ii, rx in enumerate(self.rxList):
             Ds[ii] = rx.evalDeriv(U, m,
                                 self.prob.mapping,
@@ -135,12 +142,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
     @Utils.timeIt
     def fields(self, m):
         tic = time.time()
-        u = range(self.nT+1)
+        u = list(range(self.nT+1))
         u[0] = self.initialConditions
         for ii, dt in enumerate(self.timeSteps):
             bc = self.getBoundaryConditions(ii, u[ii])
             u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
-            if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
+            if self.debug: print("Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic))
         return u
 
     @Utils.timeIt
@@ -238,7 +245,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
             f = self.fields(m)
 
         nn = len(f)-1
-        Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
+        Asubs, Adiags, Bs = list(range(nn)), list(range(nn)), list(range(nn))
         for ii in range(nn):
             dt = self.timeSteps[ii]
             bc = self.getBoundaryConditions(ii, f[ii])
@@ -263,7 +270,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
         if f is None:
             f = self.fields(m)
 
-        JvC = range(len(f)-1) # Cell to hold each row of the long vector.
+        JvC = list(range(len(f)-1)) # Cell to hold each row of the long vector.
 
         # This is done via forward substitution.
         bc = self.getBoundaryConditions(0, f[0])
@@ -295,7 +302,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
             bc = self.getBoundaryConditions(ii-1, f[ii-1])
             Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
             #select the correct part of v
-            vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
+            vpart = list(range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]))
             AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
             JTvC = AdiaginvT * (PTv[vpart] - minus)
             minus = Asub.T*JTvC  # this is now the super diagonal.

SimPEG/FLOW/Richards/__init__.py

@@ -1,2 +1,8 @@
-import Empirical
-from RichardsProblem import *
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from . import Empirical
+from .RichardsProblem import *

SimPEG/FLOW/__init__.py

@@ -1 +1,7 @@
-import Richards
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from . import Richards

SimPEG/Fields.py

@@ -1,4 +1,13 @@
-import Utils, numpy as np, scipy.sparse as sp
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
+from builtins import range
+from builtins import object
+from . import Utils
+import numpy as np, scipy.sparse as sp
 
 class Fields(object):
     """Fancy Field Storage
@@ -244,7 +253,7 @@ class TimeFields(Fields):
                 out = func(pointerFields, srcII, timeII)
             else: #loop over the time steps
                 nT = pointerShape[2]
-                out = range(nT)
+                out = list(range(nT))
                 for i, TIND_i in enumerate(timeII):
                     fieldI = pointerFields[:,:,i]
                     if fieldI.shape[0] == fieldI.size:

SimPEG/InvProblem.py

@@ -1,14 +1,21 @@
-import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
-from Utils.SolverUtils import *
-import DataMisfit
-import Regularization
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from builtins import object
+from . import Utils, Survey, Problem
+import numpy as np, scipy.sparse as sp, gc
+from .Utils.SolverUtils import *
+from . import DataMisfit
+from . import Regularization
+from future.utils import with_metaclass
 
 
-class BaseInvProblem(object):
+class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
     """BaseInvProblem(dmisfit, reg, opt)"""
 
-    __metaclass__ = Utils.SimPEGMetaClass
-
     beta    = 1.0    #: Trade-off parameter
 
     debug   = False  #: Print debugging information
@@ -54,10 +61,10 @@ class BaseInvProblem(object):
 
             Called when inversion is first starting.
         """
-        if self.debug: print 'Calling InvProblem.startup'
+        if self.debug: print('Calling InvProblem.startup')
 
         if self.reg.mref is None:
-            print 'SimPEG.InvProblem will set Regularization.mref to m0.'
+            print('SimPEG.InvProblem will set Regularization.mref to m0.')
             self.reg.mref = m0
 
         self.phi_d = np.nan
@@ -65,8 +72,8 @@ class BaseInvProblem(object):
 
         self.curModel = m0
 
-        print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
-                    ***Done using same Solver and solverOpts as the problem***"""
+        print("""SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
+                    ***Done using same Solver and solverOpts as the problem***""")
         self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts)
 
     @property
@@ -87,7 +94,7 @@ class BaseInvProblem(object):
         for mtest, u_ofmtest in self.warmstart:
             if m is mtest:
                 f = u_ofmtest
-                if self.debug: print 'InvProb is Warm Starting!'
+                if self.debug: print('InvProb is Warm Starting!')
                 break
 
         if f is None:

SimPEG/Inversion.py

@@ -1,18 +1,24 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from builtins import object
 import SimPEG
 from SimPEG import Utils, sp, np
-from Optimization import Remember, IterationPrinters, StoppingCriteria
-import Directives
+from .Optimization import Remember, IterationPrinters, StoppingCriteria
+from . import Directives
+from future.utils import with_metaclass
 
 
-class BaseInversion(object):
+class BaseInversion(with_metaclass(Utils.SimPEGMetaClass, object)):
     """
 
         Inversion Class.
 
     """
 
-    __metaclass__ = Utils.SimPEGMetaClass
-
     name = 'BaseInversion'
 
     debug   = False    #: Print debugging information

SimPEG/MT/BaseMT.py

@@ -1,7 +1,13 @@
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
 from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
-from SurveyMT import Survey, Data
-from FieldsMT import BaseMTFields
+from .SurveyMT import Survey, Data
+from .FieldsMT import BaseMTFields
 
 
 class BaseMTProblem(BaseFDEMProblem):

SimPEG/MT/FieldsMT.py

@@ -1,3 +1,9 @@
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Survey, Utils, Problem, np, sp, mkvc
 from scipy.constants import mu_0
 import sys
@@ -63,7 +69,7 @@ class Fields1D_e(BaseMTFields):
         C = self.mesh.nodalGrad
         b = (C * eSolution)
         for i, src in enumerate(srcList):
-            b[:,i] *= - 1./(1j*omega(src.freq))
+            b[:,i] *= -1./(1j*omega(src.freq))
             # There is no magnetic source in the MT problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
@@ -188,7 +194,7 @@ class Fields3D_e(BaseMTFields):
             # adjoint: returns a 2*nE long vector with zero's for py
             return np.vstack((v,np.zeros_like(v)))
         # Not adjoint: return only the px part of the vector
-        return v[:len(v)/2]
+        return v[:len(v)//2]
 
     def _e_pyDeriv_u(self, src, v, adjoint = False):
         '''
@@ -198,7 +204,7 @@ class Fields3D_e(BaseMTFields):
             # adjoint: returns a 2*nE long vector with zero's for px
             return np.vstack((np.zeros_like(v),v))
         # Not adjoint: return only the px part of the vector
-        return v[len(v)/2::]
+        return v[len(v)//2::]
 
     def _e_pxDeriv_m(self, src, v, adjoint = False):
         # assuming primary does not depend on the model
@@ -227,7 +233,7 @@ class Fields3D_e(BaseMTFields):
         C = self.mesh.edgeCurl
         b = (C * e_pxSolution)
         for i, src in enumerate(srcList):
-            b[:,i] *= - 1./(1j*omega(src.freq))
+            b[:,i] *= -1./(1j*omega(src.freq))
             # There is no magnetic source in the MT problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
@@ -238,7 +244,7 @@ class Fields3D_e(BaseMTFields):
         C = self.mesh.edgeCurl
         b = (C * e_pySolution)
         for i, src in enumerate(srcList):
-            b[:,i] *= - 1./(1j*omega(src.freq))
+            b[:,i] *= -1./(1j*omega(src.freq))
             # There is no magnetic source in the MT problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
@@ -348,4 +354,4 @@ class Fields3D_e(BaseMTFields):
         This function stacks the fields derivatives appropriately
         """
         # The fields have no dependance to the model.
-        return None
+        return None

SimPEG/MT/Problem1D/Probs.py

@@ -1,3 +1,9 @@
+from __future__ import print_function
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG.EM.Utils import omega
 from SimPEG import mkvc
 from scipy.constants import mu_0
@@ -46,7 +52,7 @@ class eForm_psField(BaseMTProblem):
             Edge inner product matrix
         """
         if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+            self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
         return self._MeMui
 
     @property
@@ -142,7 +148,7 @@ class eForm_psField(BaseMTProblem):
         for freq in self.survey.freqs:
             if self.verbose:
                 startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
+                print('Starting work for {:.3e}'.format(freq))
                 sys.stdout.flush()
             A = self.getA(freq)
             rhs  = self.getRHS(freq)
@@ -158,7 +164,7 @@ class eForm_psField(BaseMTProblem):
             # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
             # F[Src, 'b_1d'] = b[:,1]
             if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                print('Ran for {:f} seconds'.format(time.time()-startTime))
                 sys.stdout.flush()
         return F
 
@@ -191,7 +197,7 @@ class eForm_TotalField(BaseMTProblem):
             Edge inner product matrix
         """
         if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+            self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
         return self._MeMui
 
     @property
@@ -249,7 +255,7 @@ class eForm_TotalField(BaseMTProblem):
         Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
         Etot = (Ed + Eu)
         sourceAmp = 1.0
-        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
+        Etot = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
         ## Note: The analytic solution is derived with e^iwt
         eBC = np.r_[Etot[0],Etot[-1]]
         # The right hand side
@@ -274,7 +280,7 @@ class eForm_TotalField(BaseMTProblem):
         for freq in self.survey.freqs:
             if self.verbose:
                 startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
+                print('Starting work for {:.3e}'.format(freq))
                 sys.stdout.flush()
             A = self.getA(freq)
             rhs, e_o = self.getRHS(freq)
@@ -286,6 +292,6 @@ class eForm_TotalField(BaseMTProblem):
             # NOTE: only store e fields
             F[Src, 'e_1dSolution'] = e[:,0]
             if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                print('Ran for {:f} seconds'.format(time.time()-startTime))
                 sys.stdout.flush()
         return F

SimPEG/MT/Problem1D/__init__.py

@@ -1 +1,7 @@
-from Probs import eForm_TotalField, eForm_psField
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .Probs import eForm_TotalField, eForm_psField

SimPEG/MT/Problem2D/__init__.py

@@ -1 +1,7 @@
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 pass

SimPEG/MT/Problem3D/Probs.py

@@ -1,3 +1,9 @@
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+from __future__ import absolute_import
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
 from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
@@ -115,7 +121,7 @@ class eForm_ps(BaseMTProblem):
         for freq in self.survey.freqs:
             if self.verbose:
                 startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
+                print('Starting work for {:.3e}'.format(freq))
                 sys.stdout.flush()
             A = self.getA(freq)
             rhs  = self.getRHS(freq)
@@ -131,7 +137,7 @@ class eForm_ps(BaseMTProblem):
             # Note curl e = -iwb so b = -curl/iw
 
             if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                print('Ran for {:f} seconds'.format(time.time()-startTime))
                 sys.stdout.flush()
             Ainv.clean()
         return F

SimPEG/MT/Problem3D/__init__.py

@@ -1 +1,7 @@
-from Probs import eForm_ps
+from __future__ import absolute_import
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import division
+from future import standard_library
+standard_library.install_aliases()
+from .Probs import eForm_ps

SimPEG/MT/SrcMT.py

@@ -1,10 +1,16 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Utils, Problem, Maps, np, sp, mkvc
 from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
 from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
-from Utils.sourceUtils import homo1DModelSource
-from Utils import rec2ndarr
+from .Utils.sourceUtils import homo1DModelSource
+from .Utils import rec2ndarr
 import sys
 
 #################
@@ -78,7 +84,7 @@ class polxy_1Dprimary(BaseMTSrc):
             C = problem.mesh.nodalGrad
         elif problem.mesh.dim == 3:
             C = problem.mesh.edgeCurl
-        bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
+        bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
         return bBG_bp
 
     def S_e(self,problem):
@@ -86,7 +92,7 @@ class polxy_1Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -155,7 +161,7 @@ class polxy_3Dprimary(BaseMTSrc):
             C = problem.mesh.nodalGrad
         elif problem.mesh.dim == 3:
             C = problem.mesh.edgeCurl
-        bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
+        bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
         return bBG_bp
 
     def S_e(self,problem):
@@ -163,7 +169,7 @@ class polxy_3Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)

SimPEG/MT/SurveyMT.py

@@ -1,10 +1,16 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from future import standard_library
+standard_library.install_aliases()
 from SimPEG import Survey as SimPEGsurvey, Utils, Problem, Maps, np, sp, mkvc
 from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
 from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
-from Utils import rec2ndarr
-import SrcMT
+from .Utils import rec2ndarr
+from . import SrcMT
 import sys
 
 #################
@@ -76,7 +82,7 @@ class Rx(SimPEGsurvey.BaseRx):
             bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
             # Note: Has a minus sign in front, to comply with quadrant calculations.
             # Can be derived from zyx case for the 3D case.
-            f_part_complex = -ex/bx
+            f_part_complex = old_div(-ex,bx)
         # elif self.projType is 'Z2D':
         elif self.projType is 'Z3D':
             ## NOTE: Assumes that e is on edges and b on the faces. Need to generalize that or use a prop of fields to determine that.
@@ -103,13 +109,13 @@ class Rx(SimPEGsurvey.BaseRx):
             hy_py = Pby*f[src,'b_py']/mu_0
             # Make the complex data
             if 'zxx' in self.rxType:
-                f_part_complex = ( ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
+                f_part_complex = old_div(( ex_px*hy_py - ex_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
             elif 'zxy' in self.rxType:
-                f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
+                f_part_complex = old_div((-ex_px*hx_py + ex_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
             elif 'zyx' in self.rxType:
-                f_part_complex = ( ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
+                f_part_complex = old_div(( ey_px*hy_py - ey_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
             elif 'zyy' in self.rxType:
-                f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
+                f_part_complex = old_div((-ey_px*hx_py + ey_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
         elif self.projType is 'T3D':
             if self.locs.ndim == 3:
                 horLoc = self.locs[:,:,0]
@@ -127,9 +133,9 @@ class Rx(SimPEGsurvey.BaseRx):
             by_py = Pby*f[src,'b_py']
             bz_py = Pbz*f[src,'b_py']
             if 'tzx' in self.rxType:
-                f_part_complex = (- by_px*bz_py + by_py*bz_px)/(bx_px*by_py - bx_py*by_px)
+                f_part_complex = old_div((- by_px*bz_py + by_py*bz_px),(bx_px*by_py - bx_py*by_px))
             if 'tzy' in self.rxType:
-                f_part_complex = (  bx_px*bz_py - bx_py*bz_px)/(bx_px*by_py - bx_py*by_px)
+                f_part_complex = old_div((  bx_px*bz_py - bx_py*bz_px),(bx_px*by_py - bx_py*by_px))
 
         else:
             NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType))
@@ -157,8 +163,8 @@ class Rx(SimPEGsurvey.BaseRx):
                 Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
                 # ex = Pex*mkvc(f[src,'e_1d'],2)
                 # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
-                dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
-                dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
+                dP_de = -mkvc(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pex*v),2)
+                dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
                 PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
             elif self.projType is 'Z2D':
                 raise NotImplementedError('Has not been implement for 2D impedance tensor')
@@ -198,7 +204,7 @@ class Rx(SimPEGsurvey.BaseRx):
                 # Update the input vector
                 sDiag = lambda t: Utils.sdiag(mkvc(t,2))
                 # Define the components of the derivative
-                Hd = sDiag(1./(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px))
+                Hd = sDiag(old_div(1.,(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px)))
                 Hd_uV = sDiag(hy_py)*hx_px_u(v) + sDiag(hx_px)*hy_py_u(v) - sDiag(hx_py)*hy_px_u(v) - sDiag(hy_px)*hx_py_u(v)
                 # Calculate components
                 if 'zxx' in self.rxType:
@@ -247,7 +253,7 @@ class Rx(SimPEGsurvey.BaseRx):
                 # Update the input vector
                 sDiag = lambda t: Utils.sdiag(mkvc(t,2))
                 # Define the components of the derivative
-                Hd = sDiag(1./(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px))
+                Hd = sDiag(old_div(1.,(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px)))
                 Hd_uV = sDiag(by_py)*bx_px_u(v) + sDiag(bx_px)*by_py_u(v) - sDiag(bx_py)*by_px_u(v) - sDiag(by_px)*bx_py_u(v)
                 if 'tzx' in self.rxType:
                     Tij = sDiag(Hd*( - sDiag(by_px)*bz_py + sDiag(by_py)*bz_px ))
@@ -267,8 +273,8 @@ class Rx(SimPEGsurvey.BaseRx):
                 Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
                 # ex = Pex*mkvc(f[src,'e_1d'],2)
                 # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
-                dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
-                db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
+                dP_deTv = -mkvc(Pex.T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*v,2)
+                db_duv = Pbx.T/mu_0*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
                 dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
                 PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
             elif self.projType is 'Z2D':
@@ -300,17 +306,17 @@ class Rx(SimPEGsurvey.BaseRx):
                 aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True)
                 aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True)
                 aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True)
-                ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
-                ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
-                ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
-                ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
+                ahx_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
+                ahy_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
+                ahx_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
+                ahy_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
 
                 # Update the input vector
                 # Define shortcuts
                 sDiag = lambda t: Utils.sdiag(mkvc(t,2))
                 sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
                 # Define the components of the derivative
-                aHd = sDiag(1./(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px))
+                aHd = sDiag(old_div(1.,(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px)))
                 aHd_uV = lambda x: ahx_px_u(sDiag(ahy_py)*x) + ahx_px_u(sDiag(ahy_py)*x) - ahy_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(ahy_px)*x)
                 # Need to fix this to reflect the adjoint
                 if 'zxx' in self.rxType:
@@ -362,7 +368,7 @@ class Rx(SimPEGsurvey.BaseRx):
                 sDiag = lambda t: Utils.sdiag(mkvc(t,2))
                 sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
                 # Define the components of the derivative
-                aHd = sDiag(1./(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px))
+                aHd = sDiag(old_div(1.,(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px)))
                 aHd_uV = lambda x: abx_px_u(sDiag(aby_py)*x) + abx_px_u(sDiag(aby_py)*x) - aby_px_u(sDiag(abx_py)*x) - abx_py_u(sDiag(aby_px)*x)
                 # Need to fix this to reflect the adjoint
                 if 'tzx' in self.rxType: