Fixed directive to save iteration dictionary.

Fixing MT namespace to NSEM.

.bumpversion.cfg

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

.gitignore

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

.travis.yml

@@ -1,7 +1,6 @@
 language: python
 python:
   - 2.7
-  - 3.4
 
 sudo: false
 
@@ -25,25 +24,18 @@ 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:
-# 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
+  - 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/anaconda3/bin:/home/travis/miniconda/bin:/home/travis/miniconda3/bin:$PATH
+  - export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
   - conda update --yes conda
 
+# Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk sphinx
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython ipywidgets nose vtk
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git
@@ -54,26 +46,11 @@ install:
 
 # Run test
 script:
-  # test docs
   - nosetests $TEST_DIR --with-cov --cov SimPEG --cov-config .coveragerc -v -s
 
 # Calculate coverage
 after_success:
-  - 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
-
+  - coveralls --config_file .coveragerc
 
 notifications:
   email:

README.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
     :alt: SimPEG Logo
 
 ======
@@ -28,11 +28,7 @@ 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,11 +1,3 @@
-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):
@@ -69,7 +61,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
@@ -170,8 +162,8 @@ class ProblemDC_CC(Problem.BaseProblem):
         """
             Makes the matrix A(m) for the DC resistivity problem.
 
-            :param numpy.ndarray m: model
-            :rtype: scipy.sparse.csc_matrix
+            :param numpy.array m: model
+            :rtype: scipy.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/BaseIP.py

@@ -1,12 +1,5 @@
-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):
     """
@@ -59,7 +52,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
@@ -78,7 +71,7 @@ class ProblemIP(Problem.BaseProblem):
             Makes the matrix A(m) for the DC resistivity problem.
 
             :param numpy.array m: model
-            :rtype: scipy.sparse.csc_matrix
+            :rtype: scipy.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/DCIPUtils.py

@@ -1,16 +1,6 @@
-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
-from . import BaseDC as DC
-from . import BaseDC as IP
+import BaseDC as DC
+import BaseDC as IP
 import warnings
 
 def getActiveindfromTopo(mesh, topo):
@@ -77,7 +67,7 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
         if "!" in line.split(): continue
         elif line == '\n': continue
         elif line == ' \n': continue
-        temp =  list(map(float, line.split()))
+        temp =  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
@@ -238,10 +228,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
 
 
@@ -256,11 +246,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)]
@@ -350,11 +340,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(list(range(int(nstn))))*dl_x*AM_sep
-    stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*AM_sep
+    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
 
     # Create line of P1 locations
     M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
@@ -386,15 +376,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([(AB - MN_sep) // AM_sep, nrx])
+            nstn = np.min([np.floor( (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(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
+            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
 
             # Create receiver poles
             # Create line of P1 locations
@@ -429,15 +419,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 = box_l // AM_sep
+        nstn = np.floor( box_l / AM_sep )
 
         # Compute discrete pole location along line
-        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
+        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
 
         # Define number of cross lines
-        nlin = int(box_w // AM_sep)
-        lind = list(range(-nlin,nlin+1))
+        nlin = int(np.floor( box_w / AM_sep ))
+        lind = range(-nlin,nlin+1)
 
         ngrad = nstn * len(lind)
 
@@ -459,7 +449,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
@@ -678,7 +668,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,10 +1,3 @@
-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):
@@ -12,7 +5,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,10 +1,4 @@
-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
+from BaseDC import *
+from BaseIP import *
+from DCIPUtils import *
+import Utils

SimPEG/DataMisfit.py

@@ -1,16 +1,7 @@
-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
+import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
 
 
-class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, object)):
+class BaseDataMisfit(object):
     """BaseDataMisfit
 
         .. note::
@@ -18,6 +9,8 @@ class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, 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
 
@@ -100,11 +93,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,15 +1,4 @@
-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
+import Utils, numpy as np
 
 class InversionDirective(object):
     """InversionDirective"""
@@ -26,7 +15,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
@@ -79,7 +68,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
@@ -131,7 +120,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)
@@ -152,7 +141,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
 
 
@@ -178,7 +167,7 @@ class TargetMisfit(InversionDirective):
 
 
 
-class SaveEveryIteration(InversionDirective):
+class _SaveEveryIteration(InversionDirective):
     @property
     def name(self):
         if getattr(self, '_name', None) is None:
@@ -199,21 +188,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()
@@ -223,48 +212,40 @@ 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):
-    """SaveOutputDictEveryIteration"""
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """
+    Saves inversion parameters at every iteraion.
+
+
+    """
 
     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.
-        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
-        phi_ms = 0.5*ms.dot(ms)
-        if self.reg.mrefInSmooth == True:
-            mref = self.reg.mref
-        else:
-            mref = 0
-        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
-        phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim >= 2:
-            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
-            phi_my = 0.5 * my.dot(my)
-        else:
-            phi_my = 'NaN'
-        if self.prob.mesh.dim==3:
-            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
-            phi_mz = 0.5 * mz.dot(mz)
-        else:
-            phi_mz = 'NaN'
 
+        # Initialize the output dict
+        outDict = {}
+        # Save the data.
+        outDict['iter'] = self.opt.iter
+        outDict['beta'] = self.invProb.beta
+        outDict['phi_d'] = self.invProb.phi_d
+        outDict['phi_ms'] = self.reg._evalSmall(self.invProb.curModel)
+        outDict['phi_mx'] = self.reg._evalSmoothx(self.invProb.curModel)
+        outDict['phi_my'] = self.reg._evalSmoothy(self.invProb.curModel) if self.prob.mesh.dim >= 2 else 'NaN'
+        outDict['phi_mz'] = self.reg._evalSmoothz(self.invProb.curModel) if self.prob.mesh.dim==3 else 'NaN'
+        outDict['f'] = self.opt.f
+        outDict['m'] = self.invProb.curModel
+        outDict['dpred'] = self.invProb.dpred
 
         # Save the file as a npz
-        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
-
-#         mref = getattr(self, 'm_prev', None)
-#         if mref is None:
-#             if self.debug: print 'UpdateReferenceModel is using mref0'
-#             mref = self.mref0
-#         self.m_prev = self.invProb.m_current
-#         return mref
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), outDict)
 
 class Update_IRLS(InversionDirective):
 
     eps_min = None
-    eps = None
+    eps_p = None
+    eps_q = None
     norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
@@ -273,7 +254,6 @@ 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
@@ -305,25 +285,12 @@ 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
-
-            # 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))
-
+            print self.eps_p, self.eps_q, self.norms
+            self.reg.eps_p = self.eps_p
+            self.reg.eps_q = self.eps_q
             self.reg.norms = self.norms
             self.coolingFactor = 1.
             self.coolingRate = 1
@@ -339,7 +306,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
 
 
@@ -351,30 +318,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,10 +1,3 @@
-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,9 +1,4 @@
 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

@@ -1,307 +0,0 @@
-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,9 +1,3 @@
-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
@@ -40,7 +34,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,9 +1,3 @@
-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,11 +1,4 @@
-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 *
+from TDEM import hzAnalyticDipoleT
+from FDEM import hzAnalyticDipoleF
+from FDEMcasing import *
+from DC import DCAnalyticHalf, DCAnalyticSphere

SimPEG/EM/Base.py

@@ -1,9 +1,3 @@
-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
 
@@ -26,10 +20,10 @@ class BaseEMProblem(Problem.BaseProblem):
         Problem.BaseProblem.__init__(self, mesh, **kwargs)
 
 
-    surveyPair = Survey.BaseSurvey #: The survey to pair with.
-    dataPair = Survey.Data #: The data to pair with.
+    surveyPair = Survey.BaseSurvey
+    dataPair = Survey.Data
 
-    PropMap = EMPropMap #: The property mapping
+    PropMap = EMPropMap
 
     Solver = SimpegSolver
     solverOpts = {}
@@ -223,7 +217,6 @@ 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,10 +1,3 @@
-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
@@ -13,11 +6,11 @@ from SimPEG.EM.Utils import omega
 from SimPEG.Utils import Zero, Identity, sdiag
 
 
-class FieldsFDEM(SimPEG.Problem.Fields):
+class Fields(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 object acts like an array and is indexed by
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
 
     .. code-block:: python
 
@@ -49,7 +42,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
 
@@ -63,7 +56,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
@@ -77,7 +70,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
 
@@ -91,7 +84,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
 
@@ -99,7 +92,7 @@ class FieldsFDEM(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 SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
+        :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?
@@ -107,7 +100,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -117,7 +110,7 @@ class FieldsFDEM(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 SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
+        :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?
@@ -125,7 +118,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -135,7 +128,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
         """
         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 SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
+        :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?
@@ -143,7 +136,7 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -153,7 +146,7 @@ class FieldsFDEM(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 SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
+        :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?
@@ -161,18 +154,18 @@ class FieldsFDEM(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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %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(FieldsFDEM):
+class Fields3D_e(Fields):
     """
     Fields object for Problem3D_e.
 
-    :param BaseMesh mesh: mesh
-    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
+    :param Mesh mesh: mesh
+    :param Survey survey: survey
     """
 
     knownFields = {'eSolution':'E'}
@@ -187,6 +180,9 @@ class Fields3D_e(FieldsFDEM):
                     '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
@@ -292,7 +288,7 @@ class Fields3D_e(FieldsFDEM):
         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
@@ -352,7 +348,7 @@ class Fields3D_e(FieldsFDEM):
         :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:
@@ -389,8 +385,8 @@ class Fields3D_e(FieldsFDEM):
         :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)))
 
@@ -404,7 +400,7 @@ class Fields3D_e(FieldsFDEM):
         :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))
@@ -421,7 +417,7 @@ class Fields3D_e(FieldsFDEM):
         :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))
@@ -430,12 +426,12 @@ class Fields3D_e(FieldsFDEM):
 
 
 
-class Fields3D_b(FieldsFDEM):
+class Fields3D_b(Fields):
     """
     Fields object for Problem3D_b.
 
-    :param BaseMesh mesh: mesh
-    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
+    :param Mesh mesh: mesh
+    :param Survey survey: survey
     """
 
     knownFields = {'bSolution':'F'}
@@ -450,6 +446,9 @@ class Fields3D_b(FieldsFDEM):
                     '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
@@ -614,7 +613,7 @@ class Fields3D_b(FieldsFDEM):
         :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 ) ) )
@@ -631,7 +630,7 @@ class Fields3D_b(FieldsFDEM):
         :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) ) )
@@ -659,7 +658,7 @@ class Fields3D_b(FieldsFDEM):
         :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)))
 
@@ -674,7 +673,7 @@ class Fields3D_b(FieldsFDEM):
         :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:
@@ -694,12 +693,12 @@ class Fields3D_b(FieldsFDEM):
         return Zero()
 
 
-class Fields3D_j(FieldsFDEM):
+class Fields3D_j(Fields):
     """
     Fields object for Problem3D_j.
 
-    :param BaseMesh mesh: mesh
-    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
+    :param Mesh mesh: mesh
+    :param Survey survey: survey
     """
 
     knownFields = {'jSolution':'F'}
@@ -714,6 +713,9 @@ class Fields3D_j(FieldsFDEM):
                     '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
@@ -897,7 +899,7 @@ class Fields3D_j(FieldsFDEM):
         :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)))
 
@@ -911,7 +913,7 @@ class Fields3D_j(FieldsFDEM):
         :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 ) )
@@ -928,7 +930,7 @@ class Fields3D_j(FieldsFDEM):
         :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 ) )
@@ -943,7 +945,7 @@ class Fields3D_j(FieldsFDEM):
         :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)) )
@@ -958,7 +960,7 @@ class Fields3D_j(FieldsFDEM):
         :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:
@@ -976,7 +978,7 @@ class Fields3D_j(FieldsFDEM):
         :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)
 
@@ -986,12 +988,12 @@ class Fields3D_j(FieldsFDEM):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields3D_h(FieldsFDEM):
+class Fields3D_h(Fields):
     """
     Fields object for Problem3D_h.
 
-    :param BaseMesh mesh: mesh
-    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
+    :param Mesh mesh: mesh
+    :param Survey survey: survey
     """
 
     knownFields = {'hSolution':'E'}
@@ -1006,6 +1008,9 @@ class Fields3D_h(FieldsFDEM):
                     '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
@@ -1158,7 +1163,7 @@ class Fields3D_h(FieldsFDEM):
         :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)))
 
@@ -1172,7 +1177,7 @@ class Fields3D_h(FieldsFDEM):
         :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 ) ) )
@@ -1189,7 +1194,7 @@ class Fields3D_h(FieldsFDEM):
         :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) ) )
@@ -1205,7 +1210,7 @@ class Fields3D_h(FieldsFDEM):
         :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))
@@ -1220,7 +1225,7 @@ class Fields3D_h(FieldsFDEM):
         :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,13 +1,7 @@
-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 FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
+from SurveyFDEM import Survey as SurveyFDEM
+from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Utils import omega
 
@@ -37,11 +31,10 @@ 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 = FieldsFDEM
+    fieldsPair = Fields
 
     def fields(self, m):
         """
@@ -71,7 +64,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.FieldsFDEM.FieldsFDEM u: fields object
+        :param SimPEG.EM.FDEM.Fields u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -106,7 +99,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.FieldsFDEM.FieldsFDEM u: fields object
+        :param SimPEG.EM.FDEM.Fields u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -160,8 +153,8 @@ class BaseFDEMProblem(BaseEMProblem):
         Evaluates the sources for a given frequency and puts them in matrix form
 
         :param float freq: Frequency
-        :rtype: tuple
-        :return: (s_m, s_e) (nE or nF, nSrc)
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: s_m, s_e (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._formulation is 'EB':
@@ -201,7 +194,7 @@ class Problem3D_e(BaseFDEMProblem):
 
     which we solve for :math:`\mathbf{e}`.
 
-    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _solutionType = 'eSolution'
@@ -276,7 +269,7 @@ class Problem3D_e(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param SimPEG.EM.FDEM.Src src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -312,7 +305,7 @@ class Problem3D_b(BaseFDEMProblem):
     .. note ::
         The inverse problem will not work with full anisotropy
 
-    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _solutionType = 'bSolution'
@@ -407,7 +400,7 @@ class Problem3D_b(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param SimPEG.EM.FDEM.Src src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -451,7 +444,6 @@ 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 ::
@@ -461,7 +453,7 @@ class Problem3D_j(BaseFDEMProblem):
     .. note::
         This implementation does not yet work with full anisotropy!!
 
-    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _solutionType = 'jSolution'
@@ -537,8 +529,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
-        :return: RHS (nE, nSrc)
+        :rtype: numpy.ndarray (nE, nSrc)
+        :return: RHS
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -557,7 +549,7 @@ class Problem3D_j(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param SimPEG.EM.FDEM.Src src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -599,7 +591,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.BaseMesh.BaseMesh mesh: mesh
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _solutionType = 'hSolution'
@@ -616,11 +608,9 @@ 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
@@ -663,7 +653,6 @@ class Problem3D_h(BaseFDEMProblem):
         :param float freq: Frequency
         :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
-
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -677,7 +666,7 @@ class Problem3D_h(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param SimPEG.EM.FDEM.Src 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,10 +1,3 @@
-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
 
@@ -32,10 +25,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
 
     def eval(self, src, mesh, f):
         """
-        Project fields to receivers to get data.
+        Project fields to recievers to get data.
 
-        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
-        :param BaseMesh mesh: mesh used
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
         :param Fields f: fields object
         :rtype: numpy.ndarray
         :return: fields projected to recievers
@@ -51,8 +44,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
-        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
-        :param BaseMesh mesh: mesh used
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
         :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray

SimPEG/EM/FDEM/SrcFDEM.py

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

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -1,16 +1,10 @@
-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
-from . import SrcFDEM as Src
-from . import RxFDEM as Rx
+import SrcFDEM as Src
+import RxFDEM as Rx
 from SimPEG import sp
 
 class Survey(BaseEMSurvey):

SimPEG/EM/FDEM/__init__.py

@@ -1,11 +1,5 @@
-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
+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

SimPEG/EM/Static/DC/BoundaryUtils.py

@@ -1,9 +1,3 @@
-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,9 +1,3 @@
-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
@@ -15,7 +9,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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -24,7 +18,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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -32,7 +26,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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %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,9 +1,3 @@
-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
@@ -38,7 +32,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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -47,7 +41,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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -55,7 +49,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' %list(self.knownFields.keys())[0])
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %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,17 +1,11 @@
-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,18 +1,11 @@
-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,10 +1,3 @@
-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,9 +1,3 @@
-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,15 +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()
 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,10 +1,3 @@
-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,14 +1,8 @@
-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
+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

SimPEG/EM/Static/IP/ProblemIP.py

@@ -1,9 +1,3 @@
-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
@@ -11,7 +5,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,9 +1,3 @@
-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,8 +1,2 @@
-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
+from ProblemIP import Problem3D_CC, Problem3D_N
+from SurveyIP import Survey

SimPEG/EM/Static/SIP/ProblemSIP.py

@@ -1,11 +1,3 @@
-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
@@ -13,7 +5,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):
     """
@@ -113,7 +105,7 @@ class BaseSIPProblem(BaseEMProblem):
         JvAll = []
         #Assume only eta and tau (eta first then tau)
         # v = [2*Mx1]
-        v = v.reshape((v.size//2), 2), order='F')
+        v = v.reshape((int(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,10 +1,3 @@
-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,9 +1,3 @@
-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,9 +1,3 @@
-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,10 +1,3 @@
-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,11 +1,5 @@
-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
+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

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -1,11 +1,3 @@
-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
 
@@ -93,10 +85,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
 
 
@@ -111,7 +103,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
 
 
@@ -192,14 +184,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(list(range(int(nstn))))*dl_x*a
-    stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
+    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
 
     if mesh.dim==2:
         ztop = mesh.vectorNy[-1]
@@ -238,15 +230,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([(AB - b) // a, n])
+            nstn = np.min([np.floor( (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(list(range(int(nstn))))*dl_x*a
-            stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
+            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
 
             # Create receiver poles
 
@@ -283,17 +275,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(list(range(int(nstn))))*dl_x*a
-        stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
+        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
+        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
 
         # Define number of cross lines
-        nlin = int(box_w // a)
-        lind = list(range(-nlin,nlin+1))
+        nlin = int(np.floor( box_w / a ))
+        lind = range(-nlin,nlin+1)
 
         ngrad = nstn * len(lind)
 
@@ -318,7 +310,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,7 +1 @@
-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 *
+from StaticUtils import *

SimPEG/EM/Static/__init__.py

@@ -1,9 +1,3 @@
-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
+import DC
+import IP
+import SIP

SimPEG/EM/TDEM/BaseTDEM.py

@@ -1,10 +1,3 @@
-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 *
@@ -54,7 +47,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)
@@ -62,7 +55,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):
@@ -77,13 +70,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
@@ -102,13 +95,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
@@ -119,7 +112,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param FieldsTDEM f: Fields resulting from m
+            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -130,21 +123,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 v: vector (or a :class:`SimPEG.Survey.Data` object)
-            :param FieldsTDEM u: Fields resulting from m
+            :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
+            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (model object)
 
@@ -155,7 +148,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)
@@ -166,6 +159,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/SurveyTDEM.py

@@ -1,13 +1,7 @@
-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 *
-from .BaseTDEM import FieldsTDEM
+from BaseTDEM import FieldsTDEM
 
 
 class RxTDEM(Survey.BaseTimeRx):
@@ -93,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
     def getInitialFields(self, mesh):
         """Vertical magnetic dipole, magnetic vector potential"""
         if self.waveformType == "STEPOFF":
-            print(">> Step waveform: Non-zero initial condition")
+            print ">> Step waveform: Non-zero initial condition"
             if mesh._meshType is 'CYL':
                 if mesh.isSymmetric:
                     MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -105,7 +99,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
                 raise Exception('Unknown mesh for VMD')
             return {"b": mesh.edgeCurl*MVP}
         elif self.waveformType == "GENERAL":
-            print(">> General waveform: Zero initial condition")
+            print ">> General waveform: Zero initial condition"
             return {"b": np.zeros(mesh.nF)}
         else:
             raise NotImplementedError("Only use STEPOFF or GENERAL")
@@ -133,7 +127,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
     def getInitialFields(self, mesh):
         """Circular Loop, magnetic vector potential"""
         if self.waveformType == "STEPOFF":
-            print(">> Step waveform: Non-zero initial condition")
+            print ">> Step waveform: Non-zero initial condition"
             if mesh._meshType is 'CYL':
                 if mesh.isSymmetric:
                     MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -145,7 +139,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
                 raise Exception('Unknown mesh for CircularLoop')
             return {"b": mesh.edgeCurl*MVP}
         elif self.waveformType == "GENERAL":
-            print(">> General waveform: Zero initial condition")
+            print ">> General waveform: Zero initial condition"
             return {"b": np.zeros(mesh.nF)}
         else:
             raise NotImplementedError("Only use STEPOFF or GENERAL")

SimPEG/EM/TDEM/TDEM_b.py

@@ -1,14 +1,7 @@
-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 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):
@@ -94,8 +87,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a model)
-            :param FieldsTDEM u: Fields resulting from m
-            :rtype: FieldsTDEM
+            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :rtype: simpegEM.TDEM.FieldsTDEM
             :return: f
 
             Multiply G by a vector
@@ -132,9 +125,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a fields)
-            :param FieldsTDEM u: Fields resulting from m
-            :rtype: numpy.ndarray
-            :return: p (like a model)
+            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :rtype: np.ndarray (like a model)
+            :return: p
 
             Multiply G.T by a vector
         """
@@ -160,8 +153,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAh(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param FieldsTDEM p: Fields object
-            :rtype: FieldsTDEM
+            :param simpegEM.TDEM.FieldsTDEM p: Fields object
+            :rtype: simpegEM.TDEM.FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
@@ -207,8 +200,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAht(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param FieldsTDEM p: Fields object
-            :rtype: FieldsTDEM
+            :param simpegEM.TDEM.FieldsTDEM p: Fields object
+            :rtype: simpegEM.TDEM.FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
@@ -277,8 +270,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param FieldsTDEM vec: Fields object
-            :rtype: FieldsTDEM
+            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
+            :rtype: simpegEM.TDEM.FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
@@ -322,8 +315,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhtVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param FieldsTDEM vec: Fields object
-            :rtype: FieldsTDEM
+            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
+            :rtype: simpegEM.TDEM.FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where

SimPEG/EM/TDEM/__init__.py

@@ -1,9 +1,3 @@
-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
+from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
+from BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from TDEM_b import ProblemTDEM_b

SimPEG/EM/Utils/AnalyticUtils.py

@@ -1,10 +1,3 @@
-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
@@ -24,7 +17,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
     #TODO: break this out!
 
     if type(component) in [list, tuple]:
-        out = list(range(len(component)))
+        out = range(len(component))
         for i, comp in enumerate(component):
             out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
         return np.concatenate(out)
@@ -125,7 +118,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
     """
 
     if type(component) in [list, tuple]:
-        out = list(range(len(component)))
+        out = range(len(component))
         for i, comp in enumerate(component):
             out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
         return np.concatenate(out)
@@ -165,11 +158,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,9 +1,3 @@
-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
 
@@ -15,8 +9,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,8 +1,2 @@
-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
+from EMUtils import omega, k
+from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential

SimPEG/EM/Utils/testingUtils.py

@@ -1,10 +1,3 @@
-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
@@ -65,7 +58,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)
@@ -90,7 +83,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     try:
         from pymatsolver import MumpsSolver
         prb.Solver = MumpsSolver
-    except ImportError as e:
+    except ImportError, e:
         prb.Solver = SolverLU
 
     return prb
@@ -101,7 +94,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
@@ -119,7 +112,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)
@@ -128,11 +121,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,13 +1,7 @@
-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
+import TDEM
+import FDEM
+import Static
+import Base
+import Analytics
+import Utils
 from scipy.constants import mu_0, epsilon_0

SimPEG/Examples/DC_Analytic_Dipole.py

@@ -1,13 +1,7 @@
-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.EM.Static.DC as DC
+import SimPEG.DCIP as DC
 
-def run(plotIt=True):
+def run(plotIt=False):
     cs = 25.
     hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
     hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -27,15 +21,15 @@ def run(plotIt=True):
     #     ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
     #     ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
 
-    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)
+    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)
     problem.pair(survey)
     try:
         from pymatsolver import MumpsSolver
         problem.Solver = MumpsSolver
-    except Exception as e:
+    except Exception, e:
         pass
     data = survey.dpred(sigma)
 
@@ -67,8 +61,8 @@ def run(plotIt=True):
         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())
+    print run(plotIt=True)

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -1,11 +1,3 @@
-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
@@ -65,7 +57,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
@@ -132,10 +124,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)
@@ -151,10 +143,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), end=' ')
+        print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
 
-    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,9 +1,3 @@
-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
@@ -62,7 +56,7 @@ def run(plotIt=True):
     try:
         from pymatsolver import MumpsSolver
         prb.Solver = MumpsSolver
-    except ImportError as e:
+    except ImportError, e:
         prb.Solver = SolverLU
 
     prb.pair(survey)

SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py

@@ -1,9 +1,3 @@
-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,10 +1,3 @@
-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
@@ -26,13 +19,10 @@ def run(plotIt=True):
         Morrison Casing Model, and the results are used in a 2016 SEG abstract by
         Yang et al.
 
-        .. 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.
+        - 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
@@ -74,8 +64,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.]
 
 
@@ -85,25 +75,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.log(csx2/csx1) // np.log(pfx1)
+    npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
     hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
     dx1 = sum(hx1a)+sum(hx1b)
-    dx1 = dx1 // csx2
-    hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
+    dx1 = np.floor(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])
@@ -117,8 +107,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))
@@ -225,13 +215,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, Solver=solver)
+    problem = FDEM.Problem3D_h(mesh, mapping=mapping)
     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
 
@@ -258,9 +248,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,9 +1,3 @@
-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

SimPEG/Examples/FLOW_Richards_1D_Celia1990.py

@@ -1,9 +1,3 @@
-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,11 +1,3 @@
-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 *
 
 
@@ -50,40 +42,62 @@ 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()
-
-    # 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)
+    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)
 
     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,10 +1,3 @@
-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,19 +1,15 @@
-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
+from SimPEG import NSEM
 from scipy.constants import mu_0
 import matplotlib.pyplot as plt
 
+np.random.seed(1983)
+
 def run(plotIt=True):
     """
         MT: 1D: Inversion
-        =================
+        =======================
 
         Forward model 1D MT data.
         Setup and run a MT 1D inversion.
@@ -23,13 +19,13 @@ def run(plotIt=True):
     ## Setup the forward modeling
     # Setting up 1D mesh and conductivity models to forward model data.
     # Frequency
-    nFreq = 31
-    freqs = np.logspace(3,-3,nFreq)
+    nFreq = 26
+    freqs = np.logspace(2,-3,nFreq)
     # Set mesh parameters
-    ct = 20
-    air = simpeg.Utils.meshTensor([(ct,16,1.4)])
+    ct = 10
+    air = simpeg.Utils.meshTensor([(ct,25,1.4)])
     core = np.concatenate( (  np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
-    bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
+    bot = simpeg.Utils.meshTensor([(core[0],25,-1.4)])
     x0 = -np.array([np.sum(np.concatenate((core,bot)))])
     # Make the model
     m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
@@ -39,7 +35,7 @@ def run(plotIt=True):
     layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
     layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
     # Set the conductivity values
-    sig_half = 2e-3
+    sig_half = 1e-2
     sig_air = 1e-8
     sig_layer1 = .2
     sig_layer2 = .2
@@ -63,31 +59,31 @@ def run(plotIt=True):
     # Receivers
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([-0.5]),2).T,rxType))
     # Source list
     srcList =[]
     for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
     # Make the survey
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
     survey.mtrue = m_true
 
     ## Set the problem
-    problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
+    problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
     problem.pair(survey)
 
     ## Forward model data
     # Project the data
     survey.dtrue = survey.dpred(m_true)
-    survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
+    survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
+        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
         fig.suptitle('Target - smooth true')
 
 
     # Assign uncertainties
-    std = 0.05 # 5% std
+    std = 0.025 # 5% std
     survey.std = np.abs(survey.dobs*std)
     # Assign the data weight
     Wd = 1./survey.std
@@ -96,30 +92,33 @@ def run(plotIt=True):
     # Define a counter
     C =  simpeg.Utils.Counter()
     # Set the optimization
-    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
+    opt = simpeg.Optimization.ProjectedGNCG(maxIter = 25)
     opt.counter = C
-    opt.LSshorten = 0.5
+    opt.lower = np.log(1e-4)
+    opt.upper = np.log(5)
+    opt.LSshorten = 0.1
     opt.remember('xc')
     # Data misfit
     dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = Wd
     # Regularization - with a regularization mesh
-    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
+    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
     reg.mrefInSmooth = True
-    reg.alpha_s = 1e-7
+    reg.alpha_s = 1e-1
     reg.alpha_x = 1.
+
     # Inversion problem
     invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
     invProb.counter = C
     # Beta cooling
     beta = simpeg.Directives.BetaSchedule()
-    beta.coolingRate = 4
-    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
+    beta.coolingRate = 4.
+    beta.coolingFactor = 4.
+    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
+    betaest.beta0 = 1.
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
-    saveModel = simpeg.Directives.SaveModelEveryIteration()
-    saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
     # Create an inversion object
     inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
 
@@ -127,8 +126,9 @@ def run(plotIt=True):
     mopt = inv.run(m_0)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
+        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
         fig.suptitle('Target - smooth true')
+        fig.axes[0].set_ylim([-10000,500])
         plt.show()
 
 if __name__ == '__main__':

SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

@@ -0,0 +1,428 @@
+from scipy.constants import epsilon_0, mu_0
+import matplotlib.pyplot as plt
+import numpy as np
+from ipywidgets import *
+from SimPEG.EM.Utils import k, omega
+
+"""
+MT1D: n layered earth problem
+*****************************
+
+Author: Thibaut Astic
+Contact: thast@eos.ubc.ca
+Date: January 2016
+
+This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
+
+We start by looking at Maxwell's equations in the electric
+field \\\(\\\mathbf{E}\\) and the magnetic flux
+\\\(\\\mathbf{H}\\) to write the wave equations
+\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
+\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
+
+Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
+\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
+\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
+
+With U and D the Up and Down components of the E-field.
+
+The iteration from one layer to another is ensure by:
+
+\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
+    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
+
+And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
+and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
+
+The layer 0 is assumed to be the air layer.
+
+"""
+
+#Define a frquency range for a survey
+frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
+
+#Functions to create random physical Perties for a n-layered earth
+thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
+                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
+                                                            ,decimals =1))
+
+sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
+                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
+                                                      ,decimals=3))
+
+mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
+                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
+                                                    ,decimals=1))
+
+eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
+                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
+                                                                ,decimals=1))
+
+#Evaluate Impedance Z of a layer
+ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
+
+#Complex Cole-Cole Conductivity - EM utils
+PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
+
+#Converted thickness array into top of layer array
+top = lambda thick: np.cumsum(thick)
+
+#Propagation Matrix and theirs inverses
+
+#matrix T for transition of Up and Down components accross a layer
+T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
+
+Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
+
+#transition of Up and Down components accross a layer
+UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
+
+
+#matrix P relating Up and Down components with E and H fields
+P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
+
+Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
+
+
+#Time Variation of E and H
+E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
+H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
+
+#Plot the configuration of the problem
+def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
+
+    topn = top(thick)
+    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
+
+    ax.set_ylim([z.min(),z.max()])
+    ax.set_xlim([-widthg,widthg])
+
+    ax.set_ylabel("Depth (m)", fontsize=16.)
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position("right")
+
+    #define filling for the different layers
+    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
+
+    #Write the physical properties of air
+    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
+            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
+    for i in range(1,len(topn)-1,1):
+        if topn[i] == topn[i+1]:
+            pass
+        else:
+            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
+                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
+                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
+                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
+            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
+
+    #Write the physical properties of the n-th layer and fill it with pattern
+    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
+    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
+
+    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
+            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
+            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
+            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    #plot Trees!
+    ax.annotate("",
+            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+
+    ax.invert_yaxis()
+
+    return ax
+
+#Propagate Up and Down component for a certain frequency & evaluate E and H field
+
+def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
+
+    sigcm = np.zeros_like(sig,dtype='complex_')
+
+    for j in range(1,len(sig)):
+        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
+
+    K = k(f, sigcm, mu, eps)
+    Z = ImpZ(f,mu,K)
+
+    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+
+    UD[1,-1] = 1.
+
+    for i in range(-2,-(n+2),-1):
+
+        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
+        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
+
+    for j in range(0,n+1):
+        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
+
+    return UD, EH, Z ,K
+
+
+#Evaluate the apparent resistivity and phase for a frequency range
+def appres(F,H,sig,chg,taux,c,mu,eps,n):
+
+    Res = np.zeros_like(F)
+    Phase = np.zeros_like(F)
+    App_ImpZ= np.zeros_like(F,dtype='complex_')
+
+    for i in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
+
+        App_ImpZ[i] = EH[0,1]/EH[1,1]
+
+        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
+        Phase[i] = np.angle(App_ImpZ[i], deg = True)
+
+    return Res,Phase
+
+#Evaluate Up, Down components, E and H field, for a frequency range,
+#a discretized depth range and a time range (use to calculate envelope)
+def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
+
+    topc = top(H)
+
+    layer = np.zeros(len(zsample),dtype=np.int)-1
+
+    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
+
+    for i in range(0,n+1,1):
+        layer = layer+(zsample>=topc[i])*1
+
+    for j in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
+
+        for p in range(0,len(zsample)):
+
+            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
+
+            for q in range(0,len(tsample)):
+
+                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
+                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
+                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+
+    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
+
+
+#Function to Plot Apparent Resistivity and Phase
+def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
+
+    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
+
+    fig,ax = plt.subplots(1,2,figsize=(16,10))
+
+    ax[0].scatter(Res,F,color='black')
+    ax[0].set_xscale('Log')
+    ax[0].set_yscale('Log')
+    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
+    ax[0].set_ylim([F.min(),F.max()])
+    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
+    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
+    ax[0].grid(which='major')
+
+    ax0 = ax[0].twiny()
+
+    ax0.set_xlim([0.,90.])
+    ax0.set_ylim([F.min(),F.max()])
+    ax0.scatter(Phase,F,color='purple')
+    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
+
+    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
+
+    ax[0].tick_params(labelsize=16)
+    ax[1].tick_params(labelsize=16)
+    ax0.tick_params(labelsize=16)
+
+    if PlotEnvelope:
+
+        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
+        fenvelope1n=np.ones(100)*fenvelope
+        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
+
+        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
+        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
+
+        ax1=ax[1].twiny()
+
+        ax[1].tick_params(labelsize=16)
+        ax1.tick_params(labelsize=16)
+
+        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
+
+        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
+
+        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
+                      color='blue', alpha=0.1)
+
+        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
+                      color='red', alpha=0.1)
+
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+    else:
+        print 'No envelop (if True, might be slow)'
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
+        ax[1].get_xaxis().set_ticks([])
+
+    plt.show()
+
+#Interactive MT for Notebook
+def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
+
+    frangn=frange(-5,5,100.)
+    sig3= np.array([0.,0.001,0.1, 0.001])
+    thick3 = np.array([120000.,50.,50.])
+    eps3=np.array([1.,1.,1.,1])
+    mu3=np.array([1.,1.,1.,1])
+    chg3=np.array([0.,0.1,0.,0.2])
+    chg3_0=np.array([0.,0.1,0.,0.])
+    taux3=np.array([0.,0.1,0.,0.1])
+    c3=np.array([1.,1.,1.,1.])
+
+    sig3[1]=sigl1
+    sig3[1]=10.**sig3[1]
+    sig3[2]=sigl2
+    sig3[2]=10.**sig3[2]
+    sig3[3]=sigl3
+    sig3[3]=10.**sig3[3]
+    mu3[1]=mul1
+    mu3[2]=mul2
+    mu3[3]=mul3
+    eps3[1]=epsl1
+    eps3[2]=epsl2
+    eps3[3]=epsl3
+    thick3[1]=h1
+    thick3[2]=h2
+
+    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
+
+
+def run(n=3,plotIt=True):
+    # something to make a plot
+
+    F = frange(-5.,5.,20)
+    H = thick(50.,100.,n)
+    sign = sig(-5.,0.,n)
+    mun = mu(1.,2.,n)
+    epsn = eps(1.,9.,n)
+    chg = np.zeros_like(sign)
+    taux = np.zeros_like(sign)
+    c = np.zeros_like(sign)
+
+    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
+
+    if plotIt:
+
+        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
+
+    return Res, Phase
+
+if __name__ == '__main__':
+    run()
+
+
+
+
+

SimPEG/Examples/MT_3D_Foward.py

@@ -1,14 +1,8 @@
-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
 import SimPEG as simpeg
-from SimPEG import MT
+from SimPEG import NSEM
 import numpy as np
 try:
     from pymatsolver import MumpsSolver as Solver
@@ -18,7 +12,7 @@ except:
 def run(plotIt=True, nFreq=1):
     """
         MT: 3D: Forward
-        ===============
+        =======================
 
         Forward model 3D MT data.
 
@@ -43,24 +37,25 @@ def run(plotIt=True, nFreq=1):
     for loc in rx_loc:
         # NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
         for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rxList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
     # Source list
     srcList =[]
     for freq in np.logspace(3,-3,nFreq):
-        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
     # Survey MT
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
 
     ## Setup the problem object
-    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG, Solver=Solver)
+    problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
     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 = NSEM.Data(survey,dataVec)
     # Add plots
     if plotIt:
         pass

SimPEG/Examples/Maps_ComboMaps.py

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

@@ -1,47 +0,0 @@
-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,9 +1,3 @@
-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,9 +1,3 @@
-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,9 +1,3 @@
-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,11 +1,3 @@
-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):
@@ -95,7 +87,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,10 +1,3 @@
-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,11 +1,3 @@
-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):
@@ -36,15 +28,15 @@ def run(plotIt=True, n=60):
         axes[0].set_xlim([-1,17])
         axes[0].set_ylim([-1,17])
 
-        for ii, loc in zip(list(range(M.nC)),M.gridCC):
+        for ii, loc in zip(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(list(range(M.nFx)),M.gridFx):
+        for ii, loc in zip(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(list(range(M.nFy)),M.gridFy):
+        for ii, loc in zip(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,9 +1,3 @@
-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,9 +1,3 @@
-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,19 +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 *
 from SimPEG.Utils import surface2ind_topo
 
 
-def run(plotIt=True, nx=5, ny=5):
+def run(plotIt=False, nx = 5, ny = 5):
     """
-
-        Utils: surface2ind_topo
-        =======================
-
         Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
         a topographic surface.
 
@@ -23,25 +13,27 @@ def run(plotIt=True, 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,37 +1,29 @@
-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 #####
-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
+import EM_FDEM_1D_Inversion
+import Mesh_QuadTree_Creation
+import EM_TDEM_1D_Inversion
+import Mesh_QuadTree_FaceDiv
+import Mesh_Tensor_Creation
+import FLOW_Richards_1D_Celia1990
+import DC_Forward_PseudoSection
+import Mesh_Operators_CahnHilliard
+import Mesh_Basic_Types
+import Inversion_IRLS
+import Inversion_Linear
+import EM_Schenkel_Morrison_Casing
+import MT_3D_Foward
+import Mesh_Basic_ForwardDC
+import MT_1D_ForwardAndInversion
+import Utils_surface2ind_topo
+import MT_1D_analytic_nlayer_Earth
+import EM_FDEM_Analytic_MagDipoleWholespace
+import Mesh_Basic_PlotImage
+import DC_Analytic_Dipole
+import Mesh_QuadTree_HangingNodes
 
-__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"]
+__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "Mesh_Basic_ForwardDC", "MT_1D_ForwardAndInversion", "Utils_surface2ind_topo", "MT_1D_analytic_nlayer_Earth", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"]
 
 ##### AUTOIMPORTS #####
 
@@ -47,7 +39,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', 'content', 'examples'])
+    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
     shutil.rmtree(docExamplesDir)
     os.makedirs(docExamplesDir)
 
@@ -104,14 +96,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', 'content', 'examples', name + '.rst']))
+        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
 
-        print('Creating: %s.rst'%name)
+        print 'Creating: %s.rst'%name
         f = open(rst, 'w')
         f.write(out)
         f.close()

SimPEG/FLOW/Richards/Empirical.py

@@ -1,20 +1,14 @@
-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(with_metaclass(Utils.SimPEGMetaClass, object)):
+class NonLinearMap(object):
     """
     SimPEG NonLinearMap
 
     """
 
+    __metaclass__ = Utils.SimPEGMetaClass
+
     counter = None   #: A SimPEG.Utils.Counter object
     mesh = None      #: A SimPEG Mesh
 
@@ -37,7 +31,7 @@ class NonLinearMap(with_metaclass(Utils.SimPEGMetaClass, object)):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.sparse.csr_matrix
+            :rtype: scipy.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -50,7 +44,7 @@ class NonLinearMap(with_metaclass(Utils.SimPEGMetaClass, object)):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.sparse.csr_matrix
+            :rtype: scipy.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -192,7 +186,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
@@ -279,7 +273,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
@@ -349,7 +343,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,12 +1,5 @@
-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
 
 
@@ -68,7 +61,7 @@ class RichardsSurvey(Survey.BaseSurvey):
 
     @Utils.requires('prob')
     def eval(self, U, m):
-        Ds = list(range(len(self.rxList)))
+        Ds = range(len(self.rxList))
         for ii, rx in enumerate(self.rxList):
             Ds[ii] = rx.eval(U, m,
                                 self.prob.mapping,
@@ -80,7 +73,7 @@ class RichardsSurvey(Survey.BaseSurvey):
     @Utils.requires('prob')
     def evalDeriv(self, U, m):
         """The Derivative with respect to the fields."""
-        Ds = list(range(len(self.rxList)))
+        Ds = range(len(self.rxList))
         for ii, rx in enumerate(self.rxList):
             Ds[ii] = rx.evalDeriv(U, m,
                                 self.prob.mapping,
@@ -142,12 +135,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
     @Utils.timeIt
     def fields(self, m):
         tic = time.time()
-        u = list(range(self.nT+1))
+        u = 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
@@ -245,7 +238,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
             f = self.fields(m)
 
         nn = len(f)-1
-        Asubs, Adiags, Bs = list(range(nn)), list(range(nn)), list(range(nn))
+        Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
         for ii in range(nn):
             dt = self.timeSteps[ii]
             bc = self.getBoundaryConditions(ii, f[ii])
@@ -270,7 +263,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
         if f is None:
             f = self.fields(m)
 
-        JvC = list(range(len(f)-1)) # Cell to hold each row of the long vector.
+        JvC = 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])
@@ -302,7 +295,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 = list(range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]))
+            vpart = 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,8 +1,2 @@
-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 *
+import Empirical
+from RichardsProblem import *

SimPEG/FLOW/__init__.py

@@ -1,7 +1 @@
-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
+import Richards

SimPEG/Fields.py

@@ -1,13 +1,4 @@
-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
+import Utils, numpy as np, scipy.sparse as sp
 
 class Fields(object):
     """Fancy Field Storage
@@ -253,7 +244,7 @@ class TimeFields(Fields):
                 out = func(pointerFields, srcII, timeII)
             else: #loop over the time steps
                 nT = pointerShape[2]
-                out = list(range(nT))
+                out = range(nT)
                 for i, TIND_i in enumerate(timeII):
                     fieldI = pointerFields[:,:,i]
                     if fieldI.shape[0] == fieldI.size:

SimPEG/InvProblem.py

@@ -1,21 +1,14 @@
-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
+import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
+from Utils.SolverUtils import *
+import DataMisfit
+import Regularization
 
 
-class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
+class BaseInvProblem(object):
     """BaseInvProblem(dmisfit, reg, opt)"""
 
+    __metaclass__ = Utils.SimPEGMetaClass
+
     beta    = 1.0    #: Trade-off parameter
 
     debug   = False  #: Print debugging information
@@ -61,10 +54,10 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, 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
@@ -72,8 +65,8 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, 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
@@ -94,7 +87,7 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, 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,24 +1,18 @@
-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
-from . import Directives
-from future.utils import with_metaclass
+from Optimization import Remember, IterationPrinters, StoppingCriteria
+import Directives
 
 
-class BaseInversion(with_metaclass(Utils.SimPEGMetaClass, object)):
+class BaseInversion(object):
     """
 
         Inversion Class.
 
     """
 
+    __metaclass__ = Utils.SimPEGMetaClass
+
     name = 'BaseInversion'
 
     debug   = False    #: Print debugging information

SimPEG/MT/BaseMT.py

@@ -1,138 +0,0 @@
-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
-
-
-class BaseMTProblem(BaseFDEMProblem):
-    """
-        Base class for all Natural source problems.
-    """
-
-    def __init__(self, mesh, **kwargs):
-        BaseFDEMProblem.__init__(self, mesh, **kwargs)
-        Utils.setKwargs(self, **kwargs)
-    # Set the default pairs of the problem
-    surveyPair = Survey
-    dataPair = Data
-    fieldsPair = BaseMTFields
-
-    # Set the solver
-    Solver = SimpegSolver
-    solverOpts = {}
-
-    verbose = False
-    # Notes:
-    # Use the forward and devs from BaseFDEMProblem
-    # Might need to add more stuff here.
-
-    ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
-    def Jvec(self, m, v, f=None):
-        """
-        Function to calculate the data sensitivities dD/dm times a vector.
-
-            :param numpy.ndarray m (nC, 1) - conductive model
-            :param numpy.ndarray v (nC, 1) - random vector
-            :param MTfields object (optional) - MT fields object, if not given it is calculated
-            :rtype: MTdata object
-            :return: Data sensitivities wrt m
-        """
-
-        # Calculate the fields
-        if f is None:
-           f= self.fields(m)
-        # Set current model
-        self.curModel = m
-        # Initiate the Jv object
-        Jv = self.dataPair(self.survey)
-
-        # Loop all the frequenies
-        for freq in self.survey.freqs:
-            dA_du = self.getA(freq) #
-
-            dA_duI = self.Solver(dA_du, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                # We need fDeriv_m = df/du*du/dm + df/dm
-                # Construct du/dm, it requires a solve
-                # NOTE: need to account for the 2 polarizations in the derivatives.
-                f_src = f[src,:]
-                # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
-                dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
-                if dRHS_dm is None:
-                    du_dm = dA_duI * ( -dA_dm )
-                else:
-                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
-                # Calculate the projection derivatives
-                for rx in src.rxList:
-                    # Get the projection derivative
-                    # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
-                    Jv[src, rx] = PDeriv_u(mkvc(du_dm))
-            dA_duI.clean()
-        # Return the vectorized sensitivities
-        return mkvc(Jv)
-
-    def Jtvec(self, m, v, f=None):
-        """
-        Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
-
-            :param numpy.ndarray m (nC, 1) - conductive model
-            :param numpy.ndarray v (nD, 1) - vector
-            :param MTfields object u (optional) - MT fields object, if not given it is calculated
-            :rtype: MTdata object
-            :return: Data sensitivities wrt m
-        """
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size)
-
-        for freq in self.survey.freqs:
-            AT = self.getA(freq).T
-
-            ATinv = self.Solver(AT, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                ftype = self._fieldType + 'Solution'
-                f_src = f[src, :]
-
-                for rx in src.rxList:
-                    # Get the adjoint evalDeriv
-                    # PTv needs to be nE,
-                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
-                    # Get the
-                    dA_duIT = ATinv * PTv
-                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
-                    dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
-                    # Make du_dmT
-                    if dRHS_dmT is None:
-                        du_dmT = -dA_dmT
-                    else:
-                        du_dmT = -dA_dmT + dRHS_dmT
-                    # Select the correct component
-                    # du_dmT needs to be of size nC,
-                    real_or_imag = rx.projComp
-                    if real_or_imag == 'real':
-                        Jtv +=  du_dmT.real
-                    elif real_or_imag == 'imag':
-                        Jtv +=  -du_dmT.real
-                    else:
-                        raise Exception('Must be real or imag')
-            # Clean the factorization, clear memory.
-            ATinv.clean()
-        return Jtv

SimPEG/MT/Problem1D/Probs.py

@@ -1,297 +0,0 @@
-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
-from SimPEG.MT.BaseMT import BaseMTProblem
-from SimPEG.MT.SurveyMT import Survey, Data
-from SimPEG.MT.FieldsMT import Fields1D_e
-from SimPEG.MT.Utils.MT1Danalytic import getEHfields
-import numpy as np
-import multiprocessing, sys, time
-
-
-class eForm_psField(BaseMTProblem):
-    """
-    A MT problem soving a e formulation and primary/secondary fields decomposion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly half space ).
-
-
-    """
-    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
-    _fieldType = 'e_1d'
-    _eqLocs    = 'EF'
-    _sigmaPrimary = None
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__init__(self, mesh, **kwargs)
-        self.fieldsPair = Fields1D_e
-        # self._sigmaPrimary = sigmaPrimary
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        return A
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        The derivative of A wrt sigma
-        """
-
-        dsig_dm = self.curModel.sigmaDeriv
-        MeMui = self.MeMui
-        #
-        u_src = u['e_1dSolution']
-        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
-        if adjoint:
-            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
-        # Note: output has to be nN/nF, not nC/nE.
-        # v should be nC
-        return 1j * omega(freq) * ( dMfSigma_dm * v )
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
-            :return: RHS for 1 polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS wrt sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = Fields1D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print('Starting work for {:.3e}'.format(freq))
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
-            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
-
-            # Note curl e = -iwb so b = -curl e /iw
-            # 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))
-                sys.stdout.flush()
-        return F
-
-# Note this is not fully functional.
-# Missing:
-# Fields class corresponding to the fields
-# Update Jvec and Jtvec to include all the derivatives components
-# Other things ...
-class eForm_TotalField(BaseMTProblem):
-    """
-    A MT problem solving a e formulation and a Total bondary domain decompostion.
-
-    Solves the equation:
-
-    Math:
-
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'EF'
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__init__(self, mesh, **kwargs)
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    def getA(self, freq, full=False):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :param logic full: Return full A or the inner part
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        # MeMui = self.MfMui
-        # MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        if full:
-            return A
-        else:
-            return A[1:-1,1:-1]
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        raise NotImplementedError('getADeriv is not implemented')
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-        # Get sources for the frequency
-        # NOTE: Need to use the source information, doesn't really apply in 1D
-        src = self.survey.getSrcByFreq(freq)
-        # Get the full A
-        A = self.getA(freq,full=True)
-        # Define the outer part of the solution matrix
-        Aio = A[1:-1,[0,-1]]
-        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
-        Etot = (Ed + Eu)
-        sourceAmp = 1.0
-        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
-
-        return -Aio*eBC, eBC
-
-    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
-        raise NotImplementedError('getRHSDeriv not implemented yet')
-        return None
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        :param np.ndarray (nC,) m_back: Background conductivity model
-        '''
-        self.curModel = m
-        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
-
-        F = Fields1D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print('Starting work for {:.3e}'.format(freq))
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs, e_o = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_i = Ainv * rhs
-            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)
-            # NOTE: only store e fields
-            F[Src, 'e_1dSolution'] = e[:,0]
-            if self.verbose:
-                print('Ran for {:f} seconds'.format(time.time()-startTime))
-                sys.stdout.flush()
-        return F

SimPEG/MT/Problem1D/__init__.py

@@ -1,7 +0,0 @@
-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,7 +0,0 @@
-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,144 +0,0 @@
-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
-from SimPEG.MT.BaseMT import BaseMTProblem
-from SimPEG.MT.SurveyMT import Survey, Data
-from SimPEG.MT.FieldsMT import Fields3D_e
-import multiprocessing, sys, time
-
-
-
-class eForm_ps(BaseMTProblem):
-    """
-    A MT problem solving a e formulation and a primary/secondary fields decompostion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly as a 1D model).
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'FE'
-    fieldsPair = Fields3D_e
-    _sigmaPrimary = None
-
-    def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__init__(self, mesh, **kwargs)
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A system.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-        Mmui = self.MfMui
-        Msig = self.MeSigma
-        C = self.mesh.edgeCurl
-
-        return C.T*Mmui*C + 1j*omega(freq)*Msig
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        Calculate the derivative of A wrt m.
-
-        """
-
-        # This considers both polarizations and returns a nE,2 matrix for each polarization
-        if adjoint:
-            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
-        else:
-            # Need a nE,2 matrix to be returned
-            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
-        return 1j * omega(freq) * dMe_dsigV
-
-
-    def getRHS(self, freq):
-        """
-        Function to return the right hand side for the system.
-
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS with respect to sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = Fields3D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print('Starting work for {:.3e}'.format(freq))
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            # Solve the system
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # Store the fieldss
-            F[Src, 'e_pxSolution'] = e_s[:,0]
-            F[Src, 'e_pySolution'] = e_s[:,1]
-            # Note curl e = -iwb so b = -curl/iw
-
-            if self.verbose:
-                print('Ran for {:f} seconds'.format(time.time()-startTime))
-                sys.stdout.flush()
-            Ainv.clean()
-        return F
-

SimPEG/MT/Problem3D/__init__.py

@@ -1,7 +0,0 @@
-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/Utils/__init__.py

@@ -1,10 +0,0 @@
-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 .MT1Dsolutions import *  # Add the names of the functions
-from .MT1Danalytic import *
-from .dataUtils import *
-from .ediFilesUtils import *

SimPEG/MT/Utils/srcUtils.py

@@ -1,52 +0,0 @@
-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 as simpeg, numpy as np
-
-def homo1DModelSource(mesh,freq,m_back):
-    '''
-        Function that calculates and return background fields for a 3D mesh and model.
-        The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
-        which is assigned at the fields everywhere for the respective polarizations.2
-
-        :param Simpeg mesh object mesh: Holds information on the discretization
-        :param float freq: The frequency to solve at
-        :param np.array m_back: Background model of conductivity to base the calculations on.
-        :rtype: numpy.ndarray (mesh.nE,2)
-        :return: eBG_bp, E fields for the background model at both polarizations.
-
-    '''
-
-    # import
-    from SimPEG.MT.Utils import get1DEfields
-    # Get a 1d solution for a halfspace background
-    mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
-    # Note: Everything is using e^iwt
-    e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
-    # Setup x (east) polarization (_x)
-    ex_px = np.zeros(mesh.vnEx,dtype=complex)
-    ey_px = np.zeros((mesh.nEy,1),dtype=complex)
-    ez_px = np.zeros((mesh.nEz,1),dtype=complex)
-    # Assign the source to ex_x
-    for i in np.arange(mesh.vnEx[0]):
-        for j in np.arange(mesh.vnEx[1]):
-            ex_px[i,j,:] = -e0_1d
-    eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
-    # Setup y (north) polarization (_py)
-    ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
-    ey_py = np.zeros(mesh.vnEy, dtype='complex128')
-    ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
-    # Assign the source to ey_py
-
-    for i in np.arange(mesh.vnEy[0]):
-        for j in np.arange(mesh.vnEy[1]):
-            ey_py[i,j,:] = e0_1d
-    # ey_py[1:-1,1:-1,1:-1] = 0
-    eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
-
-    # Return the electric fields
-    eBG_bp = np.hstack((eBG_px,eBG_py))
-    return eBG_bp

SimPEG/MT/__init__.py

@@ -1,11 +0,0 @@
-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 Utils
-from .SurveyMT import Rx, Survey, Data
-from .FieldsMT import Fields1D_e, Fields3D_e
-from . import Problem1D, Problem2D, Problem3D
-from . import SrcMT

SimPEG/Maps.py

@@ -1,32 +1,23 @@
-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 str
-from builtins import range
-from builtins import object
-from . import Utils
-import numpy as np, scipy.sparse as sp
+import Utils, numpy as np, scipy.sparse as sp
 from scipy.sparse.linalg import LinearOperator
-from .Tests import checkDerivative
+from Tests import checkDerivative
+from PropMaps import PropMap, Property
 from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
 import warnings
-from future.utils import with_metaclass
 
-class IdentityMap(with_metaclass(Utils.SimPEGMetaClass, object)):
+class IdentityMap(object):
     """
     SimPEG Map
 
     """
+    __metaclass__ = Utils.SimPEGMetaClass
 
     def __init__(self, mesh=None, nP=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
         if nP is not None:
-            assert type(nP) in [int, int], ' Number of parameters must be an integer.'
+            assert type(nP) in [int, long], ' Number of parameters must be an integer.'
 
         self.mesh = mesh
         self._nP  = nP
@@ -50,8 +41,8 @@ class IdentityMap(with_metaclass(Utils.SimPEGMetaClass, object)):
             If this is a meshless mapping (i.e. nP is defined independently)
             the shape will be the the shape (nP,nP).
 
-            :rtype: tuple
-            :return: shape of the operator as a tuple (int,int)
+            :rtype: (int,int)
+            :return: shape of the operator as a tuple
         """
         if self._nP is not None:
             return (self.nP, self.nP)
@@ -95,7 +86,7 @@ class IdentityMap(with_metaclass(Utils.SimPEGMetaClass, object)):
             The derivative of the transformation.
 
             :param numpy.array m: model
-            :rtype: scipy.sparse.csr_matrix
+            :rtype: scipy.csr_matrix
             :return: derivative of transformed model
 
         """
@@ -110,7 +101,7 @@ class IdentityMap(with_metaclass(Utils.SimPEGMetaClass, object)):
             :return: passed the test?
 
         """
-        print('Testing %s' % str(self))
+        print 'Testing %s' % str(self)
         if m is None:
             m = abs(np.random.rand(self.nP))
         if 'plotIt' not in kwargs:
@@ -225,7 +216,7 @@ class ExpMap(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.sparse.csr_matrix
+            :rtype: scipy.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -300,7 +291,7 @@ class LogMap(IdentityMap):
         deriv = np.zeros(mod.shape)
         tol = 1e-16 # zero
         ind = np.greater_equal(np.abs(mod),tol)
-        deriv[ind] = 1.0 / mod[ind]
+        deriv[ind] = 1.0/mod[ind]
         return Utils.sdiag(deriv)
 
     def inverse(self, m):
@@ -375,13 +366,13 @@ class SurjectVertical1D(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.sparse.csr_matrix
+            :rtype: scipy.csr_matrix
             :return: derivative of transformed model
         """
         repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
         repVec = sp.csr_matrix(
                     (np.ones(repNum),
-                    (list(range(repNum)), np.zeros(repNum))
+                    (range(repNum), np.zeros(repNum))
                     ), shape=(repNum, 1))
         return sp.kron(sp.identity(self.nP), repVec)
 
@@ -436,14 +427,14 @@ class Surject2Dto3D(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.sparse.csr_matrix
+            :rtype: scipy.csr_matrix
             :return: derivative of transformed model
         """
         inds = self * np.arange(self.nP)
         nC, nP = self.mesh.nC, self.nP
         P = sp.csr_matrix(
                     (np.ones(nC),
-                    (list(range(nC)), inds)
+                    (range(nC), inds)
                 ), shape=(nC, nP))
         return P
 
@@ -511,13 +502,11 @@ class InjectActiveCells(IdentityMap):
         if Utils.isScalar(valInactive):
             self.valInactive = np.ones(self.nC)*float(valInactive)
         else:
-            self.valInactive = np.ones(self.nC)
-            self.valInactive[self.indInactive] = valInactive.copy()
-
+            self.valInactive = valInactive.copy()
         self.valInactive[self.indActive] = 0
 
         inds = np.nonzero(self.indActive)[0]
-        self.P = sp.csr_matrix((np.ones(inds.size),(inds, list(range(inds.size)))), shape=(self.nC, self.nP))
+        self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
 
     @property
     def shape(self):
@@ -604,14 +593,14 @@ class ComplexMap(IdentityMap):
 
     @property
     def shape(self):
-        return (self.nP // 2,self.nP)
+        return (self.nP/2,self.nP)
 
     def _transform(self, m):
         nC = self.mesh.nC
         return m[:nC] + m[nC:]*1j
 
     def deriv(self, m):
-        nC = self.nP//2
+        nC = self.nP/2
         shp = (nC, nC*2)
         def fwd(v):
             return v[:nC] + v[nC:]*1j
@@ -704,7 +693,7 @@ class PolyMap(IdentityMap):
         self.actInd = actInd
 
         if getattr(self, 'actInd', None) is None:
-            self.actInd = list(range(self.mesh.nC))
+            self.actInd = range(self.mesh.nC)
             self.nC = self.mesh.nC
 
         else:
@@ -740,7 +729,7 @@ class PolyMap(IdentityMap):
             elif self.normal =='Y':
                 f = polynomial.polyval(X, c) - Y
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
             X = self.mesh.gridCC[self.actInd,0]
@@ -753,10 +742,10 @@ class PolyMap(IdentityMap):
             elif self.normal =='Z':
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
 
         else:
-            raise Exception
+            raise(Exception("Only supports 2D"))
 
 
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
@@ -778,7 +767,7 @@ class PolyMap(IdentityMap):
                 f = polynomial.polyval(X, c) - Y
                 V = polynomial.polyvander(X, len(c)-1)
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
             X = self.mesh.gridCC[self.actInd,0]
@@ -795,7 +784,7 @@ class PolyMap(IdentityMap):
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
                 V = polynomial.polyvander2d(X, Y, self.order)
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
 
         if self.logSigma:
             g1 = -(np.arctan(alpha*f)/np.pi + 0.5)*sig1 + sig1
@@ -843,7 +832,7 @@ class SplineMap(IdentityMap):
         elif self.mesh.dim == 3:
             return np.size(self.pts)*2+2
         else:
-            raise Exception
+            raise(Exception("Only supports 2D and 3D"))
 
     def _transform(self, m):
         # Set model parameters
@@ -862,7 +851,7 @@ class SplineMap(IdentityMap):
             elif self.normal =='Y':
                 f = self.spl(X) - Y
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
 
         # 3D:
         # Comments:
@@ -877,7 +866,7 @@ class SplineMap(IdentityMap):
 
             npts = np.size(self.pts)
             if np.mod(c.size, 2):
-                raise Exception
+                raise(Exception("Put even points!"))
 
             self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
                         "splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
@@ -890,9 +879,9 @@ class SplineMap(IdentityMap):
             # elif self.normal =='Y':
             # elif self.normal =='Z':
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
         else:
-            raise Exception
+            raise(Exception("Only supports 2D and 3D"))
 
 
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
@@ -912,7 +901,7 @@ class SplineMap(IdentityMap):
             elif self.normal =='Y':
                 f = self.spl(X) - Y
             else:
-                raise Exception
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
@@ -926,7 +915,7 @@ class SplineMap(IdentityMap):
             # elif self.normal =='Y':
             # elif self.normal =='Z':
             else:
-                raise Exception
+                raise(Exception("Not Implemented for Y and Z, your turn :)"))
 
         if self.logSigma:
             g1 = -(np.arctan(alpha*f)/np.pi + 0.5)*sig1 + sig1
@@ -982,5 +971,11 @@ class SplineMap(IdentityMap):
                     fderiv = (flinesa-flinesb)/(2*dy)
                     g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
         else :
-            raise Exception
+            raise(Exception("Not Implemented for Y and Z, your turn :)"))
         return sp.csr_matrix(np.c_[g1,g2,g3])
+
+
+
+
+
+

SimPEG/Mesh/BaseMesh.py

@@ -1,11 +1,3 @@
-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 builtins import object
 import numpy as np
 from SimPEG import Utils
 
@@ -15,8 +7,8 @@ class BaseMesh(object):
     BaseMesh does all the counting you don't want to do.
     BaseMesh should be inherited by meshes with a regular structure.
 
-    :param numpy.array n: (or list) number of cells in each direction (dim, )
-    :param numpy.array x0: (or list) Origin of the mesh (dim, )
+    :param numpy.array,list n: number of cells in each direction (dim, )
+    :param numpy.array,list x0: Origin of the mesh (dim, )
 
     """
 
@@ -42,8 +34,8 @@ class BaseMesh(object):
         """
         Origin of the mesh
 
-        :rtype: numpy.array
-        :return: x0, (dim, )
+        :rtype: numpy.array (dim, )
+        :return: x0
         """
         return self._x0
 
@@ -124,8 +116,8 @@ class BaseMesh(object):
         """
         Total number of edges in each direction
 
-        :rtype: numpy.array
-        :return: [nEx, nEy, nEz], (dim, )
+        :rtype: numpy.array (dim, )
+        :return: [nEx, nEy, nEz]
 
         .. plot::
             :include-source:
@@ -181,8 +173,8 @@ class BaseMesh(object):
         """
         Total number of faces in each direction
 
-        :rtype: numpy.array
-        :return: [nFx, nFy, nFz], (dim, )
+        :rtype: numpy.array (dim, )
+        :return: [nFx, nFy, nFz]
 
         .. plot::
             :include-source:
@@ -208,8 +200,8 @@ class BaseMesh(object):
         """
         Face Normals
 
-        :rtype: numpy.array
-        :return: normals, (sum(nF), dim)
+        :rtype: numpy.array (sum(nF), dim)
+        :return: normals
         """
         if self.dim == 2:
             nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)]
@@ -226,8 +218,8 @@ class BaseMesh(object):
         """
         Edge Tangents
 
-        :rtype: numpy.array
-        :return: normals, (sum(nE), dim)
+        :rtype: numpy.array (sum(nE), dim)
+        :return: normals
         """
         if self.dim == 2:
             tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)]
@@ -244,9 +236,8 @@ class BaseMesh(object):
         Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
 
         :param numpy.array fV: face vector with shape (nF, dim)
-        :rtype: numpy.array
-        :return: projected face vector, (nF, )
-
+        :rtype: numpy.array with shape (nF, )
+        :return: projected face vector
         """
         assert isinstance(fV, np.ndarray), 'fV must be an ndarray'
         assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
@@ -257,9 +248,8 @@ class BaseMesh(object):
         Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
 
         :param numpy.array eV: edge vector with shape (nE, dim)
-        :rtype: numpy.array
-        :return: projected edge vector, (nE, )
-
+        :rtype: numpy.array with shape (nE, )
+        :return: projected edge vector
         """
         assert isinstance(eV, np.ndarray), 'eV must be an ndarray'
         assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
@@ -305,7 +295,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Total number of cells in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: [nCx, nCy, nCz]
         """
         return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
@@ -345,7 +335,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Total number of nodes in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: [nNx, nNy, nNz]
         """
         return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
@@ -355,7 +345,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of x-edges in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: vnEx
         """
         return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
@@ -365,7 +355,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of y-edges in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: vnEy or None if dim < 2
         """
         return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
@@ -375,7 +365,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of z-edges in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: vnEz or None if dim < 3
         """
         return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
@@ -385,7 +375,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of x-faces in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: vnFx
         """
         return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
@@ -395,7 +385,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of y-faces in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: vnFy or None if dim < 2
         """
         return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
@@ -405,7 +395,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of z-faces in each direction
 
-        :rtype: numpy.array
+        :rtype: numpy.array (dim, )
         :return: vnFz or None if dim < 3
         """
         return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])

SimPEG/Mesh/CurvilinearMesh.py

@@ -1,16 +1,8 @@
-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 Utils, np
-from .BaseMesh import BaseRectangularMesh
-from .DiffOperators import DiffOperators
-from .InnerProducts import InnerProducts
-from .View import CurvView
-from future.utils import with_metaclass
+from BaseMesh import BaseRectangularMesh
+from DiffOperators import DiffOperators
+from InnerProducts import InnerProducts
+from View import CurvView
 
 # Some helper functions.
 length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
@@ -19,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
 
 
-class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (BaseRectangularMesh, DiffOperators, InnerProducts, CurvView), {}))):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
     """
     CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
@@ -34,6 +26,8 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
             M.plotGrid(showIt=True)
     """
 
+    __metaclass__ = Utils.SimPEGMetaClass
+
     _meshType = 'Curv'
 
     def __init__(self, nodes):
@@ -297,7 +291,7 @@ class CurvilinearMesh(with_metaclass(Utils.SimPEGMetaClass, type('NewBase', (Bas
     normals = property(**normals())
 
     def edge():
-        doc = "Edge lengths."
+        doc = "Edge legnths."
 
         def fget(self):
             if(self._edge is None or self._tangents is None):
@@ -351,4 +345,4 @@ if __name__ == '__main__':
         X, Y = Utils.ndgrid(h1, h2, vector=False)
         M = CurvilinearMesh([X, Y])
 
-    print(M.r(M.normals, 'F', 'Fx', 'V'))
+    print M.r(M.normals, 'F', 'Fx', 'V')

SimPEG/Mesh/CylMesh.py

@@ -1,16 +1,10 @@
-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()
 import numpy as np
 import scipy.sparse as sp
 from scipy.constants import pi
 from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
-from .TensorMesh import BaseTensorMesh, BaseRectangularMesh
-from .InnerProducts import InnerProducts
-from .View import CylView
+from TensorMesh import BaseTensorMesh, BaseRectangularMesh
+from InnerProducts import InnerProducts
+from View import CylView
 
 
 class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
@@ -37,7 +31,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         BaseTensorMesh.__init__(self, h, x0)
         assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
         if self.dim == 2:
-            print('Warning, a disk mesh has not been tested thoroughly.')
+            print 'Warning, a disk mesh has not been tested thoroughly.'
         cartesianOrigin = np.zeros(self.dim) if cartesianOrigin is None else cartesianOrigin
         assert len(cartesianOrigin) == self.dim, "cartesianOrigin must be the same length as the dimension of the mesh."
         self.cartesianOrigin = np.array(cartesianOrigin, dtype=float)
@@ -74,8 +68,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of x-faces in each direction
 
-        :rtype: numpy.array
-        :return: vnFx, (dim, )
+        :rtype: numpy.array (dim, )
+        :return: vnFx
         """
         return self.vnC
 
@@ -84,8 +78,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of y-edges in each direction
 
-        :rtype: numpy.array
-        :return: vnEy or None if dim < 2, (dim, )
+        :rtype: numpy.array (dim, )
+        :return: vnEy or None if dim < 2
         """
         nNx = self.nNx if self.isSymmetric else self.nNx - 1
         return np.r_[nNx, self.nCy, self.nNz]
@@ -95,8 +89,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of z-edges in each direction
 
-        :rtype: numpy.array
-        :return: vnEz or None if nCy > 1, (dim, )
+        :rtype: numpy.array (dim, )
+        :return: vnEz or None if nCy > 1
         """
         if self.isSymmetric:
             return np.r_[self.nNx, self.nNy, self.nCz]