trial for secondary data

Conflicts:
	SimPEG/Problem.py

.bumpversion.cfg

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

.travis.yml

@@ -16,12 +16,12 @@ addons:
 
 env:
     - TEST_DIR="tests/mesh tests/base tests/utils"
-    - TEST_DIR=tests/examples
-    - TEST_DIR=tests/em/fdem/forward
     - TEST_DIR=tests/em/fdem/inverse/derivs
-    - TEST_DIR=tests/em/fdem/inverse/adjoint
     - TEST_DIR=tests/em/tdem
     - TEST_DIR=tests/flow
+    - TEST_DIR=tests/examples
+    - TEST_DIR=tests/em/fdem/inverse/adjoint
+    - TEST_DIR=tests/em/fdem/forward
 
 # Setup anaconda
 before_install:
@@ -33,7 +33,7 @@ before_install:
 
 # Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git

AUTHORS.rst

@@ -1,9 +1,13 @@
-- Luz Angelica Caudillo-Mata, (`@lacmajedrez <https://github.com/lacmajedrez/>`_)
 - Rowan Cockett, (`@rowanc1 <https://github.com/rowanc1/>`_)
-- Eldad Haber, (`@ehaber99 <https://github.com/ehaber99/>`_)
 - Lindsey Heagy, (`@lheagy <https://github.com/lheagy/>`_)
 - Seogi Kang, (`@sgkang <https://github.com/sgkang/>`_)
-- Dave Marchant, (`@dwfmarchant <https://github.com/dwfmarchant/>`_)
+- Brendan Smithyman, (`@bsmithyman <https://github.com/bsmithyman/>`_)
 - Gudni Rosenkjaer, (`@grosenkj <https://github.com/grosenkj/>`_)
+- Dom Fournier, (`@fourndo <https://github.com/fourndo/>`_)
+- Dave Marchant, (`@dwfmarchant <https://github.com/dwfmarchant/>`_)
 - Lars Ruthotto, (`@lruthotto <https://github.com/lruthotto/>`_)
 - Mike Wathen, (`@mrwathen <https://github.com/mrwathen/>`_)
+- Luz Angelica Caudillo-Mata, (`@lacmajedrez <https://github.com/lacmajedrez/>`_)
+- Eldad Haber, (`@ehaber99 <https://github.com/ehaber99/>`_)
+- Doug Oldenburg, (`@dougoldenburg <https://github.com/dougoldenburg/>`_)
+- Adam Pidlisecky, (`@aPid1 <https://github.com/aPid1/>`_)

CITATION.rst

@@ -1,14 +1,13 @@
 Citing SimPEG
-=============
+-------------
 
-There is a paper about SimPEG!
+There is a `paper about SimPEG <http://dx.doi.org/10.1016/j.cageo.2015.09.015>`_, if you use this code, please help our scientific visibility by citing our work!
 
 
     Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A., & Oldenburg, D. W. (2015). SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers & Geosciences.
 
 
 BibTex:
--------
 
 .. code::
 

GCEtools/PETScIO.py

@@ -1,470 +0,0 @@
-#!/usr/bin/python
-"""
-Input and output functions.
-"""
-
-
-import os as _os
-import errno as _errno
-import sys as _sys
-import numpy as _np
-
-from petsc4py import PETSc as _PETSc
-import fileinput as _fl
-
-def vecToArray(obj):
-    """ Converts a PETSc vector to a numpy array, available on *all* MPI nodes.
-
-        Args:
-            obj (petsc4py.PETSc.Vec): input vector.
-
-        Returns:
-            numpy.array :
-        """
-    # scatter vector 'obj' to all processes
-    comm = obj.getComm()
-    scatter, obj0 = _PETSc.Scatter.toAll(obj)
-    scatter.scatter(obj, obj0, False, _PETSc.Scatter.Mode.FORWARD)
-
-    return _np.asarray(obj0)
-
-    # deallocate
-    comm.barrier()
-    scatter.destroy()
-    obj0.destroy()
-
-
-def vecToArray0(obj):
-    """ Converts a PETSc vector to a numpy array available on MPI node 0.
-
-        Args:
-            obj (petsc4py.PETSc.Vec): input vector.
-
-        Returns:
-            numpy.array :
-        """
-    # scatter vector 'obj' to process 0
-    comm = obj.getComm()
-    rank = comm.getRank()
-    scatter, obj0 = _PETSc.Scatter.toZero(obj)
-    scatter.scatter(obj, obj0, False, _PETSc.Scatter.Mode.FORWARD)
-
-    if rank == 0:   return _np.asarray(obj0)
-
-    # deallocate
-    comm.barrier()
-    scatter.destroy()
-    obj0.destroy()
-
-
-def arrayToVec(vecArray):
-    """ Converts a (global) array to a PETSc vector over :attr:`petsc4py.PETSc.COMM_WORLD`.
-
-        Args:
-            vecArray (array or numpy.array): input vector.
-
-        Returns:
-            petsc4py.PETSc.Vec() :
-        """
-    vec = _PETSc.Vec().create(comm=_PETSc.COMM_WORLD)
-    vec.setSizes(len(vecArray))
-    vec.setUp()
-    (Istart,Iend) = vec.getOwnershipRange()
-    return vec.createWithArray(vecArray[Istart:Iend],
-            comm=_PETSc.COMM_WORLD)
-    vec.destroy()
-
-
-def arrayToMat(matArray):
-    """ Converts a (global) 2D array to a PETSc matrix over :attr:`petsc4py.PETSc.COMM_WORLD`.
-
-        Args:
-            matArray (array or numpy.array): input square array.
-
-        :rtype: petsc4py.PETSc.Mat()
-
-        .. important::
-            Requires `SciPy <http://www.scipy.org>`_.
-
-        """
-    try:
-        import scipy.sparse as sparse
-    except:
-        print '\nERROR: loading matrices from txt files requires Scipy!'
-        return
-
-    matSparse =matArray
-
-    mat = _PETSc.Mat().createAIJ(size=matSparse.shape,comm=_PETSc.COMM_WORLD)
-    (Istart,Iend) = mat.getOwnershipRange()
-
-    ai = matSparse.indptr[Istart:Iend+1] - matSparse.indptr[Istart]
-    aj = matSparse.indices[matSparse.indptr[Istart]:matSparse.indptr[Iend]]
-    av = matSparse.data[matSparse.indptr[Istart]:matSparse.indptr[Iend]]
-
-    mat.setValuesCSR(ai,aj,av)
-    mat.assemble()
-
-    return mat
-    mat.destroy()
-
-def matToSparse(mat):
-    """ Converts a PETSc matrix to a (global) sparse matrix.
-
-        Args:
-            mat (petsc4py.PETSc.Mat): input PETSc matrix.
-
-        :rtype: scipy.sparse.csr_matrix
-
-        .. important::
-            Requires `SciPy <http://www.scipy.org>`_.
-
-        """
-    import scipy.sparse as sparse
-
-    data = mat.getValuesCSR()
-
-    (Istart,Iend) = mat.getOwnershipRange()
-    columns = mat.getSize()[0]
-    sparseSubMat = sparse.csr_matrix(data[::-1],shape=(Iend-Istart,columns))
-
-    comm = _PETSc.COMM_WORLD
-
-    sparseSubMat = comm.tompi4py().allgather(sparseSubMat)
-
-    return sparse.vstack(sparseSubMat)
-
-def adjToH(adj,d=[0],amp=[0.]):
-    """ Creates a 1 particle PETSc-type Hamiltonian matrix from a PETSc adjacency matrix.
-
-        Args:
-            adj (petsc4py.PETSc.Mat): input PETSc-type adjacency matrix.
-
-            d (array of ints): an array containing *integers* indicating the nodes
-                        where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
-
-            amp (array of floats):   an array containing *floats* indicating the diagonal defect
-                        amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
-
-        Returns:
-             : 1 particle Hamiltonian matrix
-        :rtype: petsc4py.PETSc.Mat()
-
-        Warning:
-            * The size of ``a`` and ``d`` must be identical
-
-                >>> amp = [0.5,-1.,4.2]
-                >>> len(d) == len(amp)
-                True
-
-            * Elements of ``d`` can range from :math:`[0,N-1]` where the adjacency matrix is :math:`N\\times N`.
-
-        """
-    (Istart,Iend) = adj.getOwnershipRange()
-    diagSum = []
-    for i in range(Istart,Iend):
-        diagSum.append(_np.sum(adj.getRow(i)[-1]))
-        for j,val in enumerate(d):
-            if i==val:  diagSum[i-Istart] += amp[j]
-
-    mat = _PETSc.Mat().create(comm=_PETSc.COMM_WORLD)
-    mat.setSizes(adj.getSize())
-    mat.setUp()
-
-    for i in range(Istart,Iend):
-        mat.setValue(i,i,diagSum[i-Istart])
-
-    mat.assemble()
-    mat.axpy(-1,adj)
-
-    return mat
-    mat.destroy()
-
-#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#---------------------- Vec I/O functions ---------------------------
-#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-def exportVec(vec,filename,filetype):
-    """ Export a PETSc vector to a file.
-
-        Args:
-            vec (petsc4py.PETSc.Vec): input vector.
-            filename (str): path to desired output file.
-            filetype (str): the filetype of the exported vector.
-
-                            * ``'txt'`` - a column vector in text format.
-                            * ``'bin'`` - a PETSc binary vector.
-        """
-    if _os.path.isabs(filename):
-        outDir = _os.path.dirname(filename)
-    else:
-        outDir = './'+_os.path.dirname(filename)
-
-    # create output directory if it doesn't exist
-    try:
-        _os.mkdir(outDir)
-    except OSError as exception:
-        if exception.errno != _errno.EEXIST:
-            raise
-
-    if filetype == 'txt':
-        # scatter prob to process 0
-        comm = vec.getComm()
-        rank = comm.getRank()
-        scatter, vec0 = _PETSc.Scatter.toZero(vec)
-        scatter.scatter(vec, vec0, False, _PETSc.Scatter.Mode.FORWARD)
-
-        # use process 0 to write to text file
-        if rank == 0:
-            array0 = _np.asarray(vec0)
-            with open(filename,'w') as f:
-                for i in range(len(array0)):
-                    f.write('{0: .12e}\n'.format(array0[i]))
-
-        # deallocate
-        comm.barrier()
-        scatter.destroy()
-        vec0.destroy()
-
-    elif filetype == 'bin':
-        binSave = _PETSc.Viewer().createBinary(filename, 'w')
-        binSave(vec)
-        binSave.destroy()
-
-    vec.comm.barrier()
-
-def loadVec(filename,filetype):
-    """ Import a PETSc vector from a file.
-
-        Args:
-            filename (str): path to input file.
-            filetype (str): the filetype.
-
-                            * ``'txt'`` - a column vector in text format.
-                            * ``'bin'`` - a PETSc binary vector.
-        """
-    if filetype == 'txt':
-        try:
-            vecArray = _np.loadtxt(filename,dtype=_PETSc.ScalarType)
-            return arrayToVec(vecArray)
-        except:
-            print "\nERROR: input state space file " + filename\
-                + " does not exist or is in an incorrect format"
-            _sys.exit()
-
-    elif filetype == 'bin':
-        binLoad = _PETSc.Viewer().createBinary(filename, 'r')
-        try:
-            return _PETSc.Vec().load(binLoad)
-        except:
-            print "\nERROR: input state space file " + filename\
-                + " does not exist or is in an incorrect format"
-            _sys.exit()
-        binLoad.destroy()
-
-#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#---------------------- Mat I/O functions ---------------------------
-#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-def exportMat(mat,filename,filetype,mattype=None):
-    """ Export a PETSc matrix to a file.
-
-        Args:
-            mat (petsc4py.PETSc.Mat): input matrix.
-            filename (str): path to desired output file.
-            filetype (str): the filetype of the exported vector.
-
-                            * ``'txt'`` - a 2D matrix array in text format.
-                            * ``'bin'`` - a PETSc binary matrix.
-            mattype (str): (``None``,``'adj'``) - if set to ``adj``, only
-                            integers ``0`` and ``1`` are written. Note
-                            that this only applied in ``txt`` mode.
-        """
-    rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
-
-    if _os.path.isabs(filename):
-        outDir = _os.path.dirname(filename)
-    else:
-        outDir = './'+_os.path.dirname(filename)
-
-    # create output directory if it doesn't exist
-    try:
-        _os.mkdir(outDir)
-    except OSError as exception:
-        if exception.errno != _errno.EEXIST:
-            raise
-
-    if filetype == 'txt':
-        txtSave = _PETSc.Viewer().createASCII(filename, 'w',
-            format=_PETSc.Viewer.Format.ASCII_DENSE, comm=_PETSc.COMM_WORLD)
-        txtSave(mat)
-        txtSave.destroy()
-
-        if rank == 0:
-            for line in _fl.FileInput(filename,inplace=1):
-                if line[2] != 't':
-                    if mattype == 'adj':
-                        line = line.replace(" i","j")
-                        line = line.replace(" -","-")
-                        line = line.replace("+-","-")
-                        line = line.replace("0000e+01+0.00000e+00j","")
-                        line = line.replace(".00000e+00+0.00000e+00j","")
-                        line = line.replace(".","")
-                        line = line.replace(" -","\t-")
-                        line = line.replace("  ","\t")
-                        line = line.replace(" ","")
-                        line = line.replace("\t"," ")
-                        print line,
-                    else:
-                        line = line.replace(" i","j")
-                        line = line.replace(" -","-")
-                        line = line.replace("+-","-")
-                        print line,
-
-    elif filetype == 'bin':
-        binSave = _PETSc.Viewer().createBinary(filename, 'w', comm=_PETSc.COMM_WORLD)
-        binSave(mat)
-        binSave.destroy()
-
-    mat.comm.barrier()
-
-def loadMat(filename,filetype,delimiter=None):
-    """ Import a PETSc matrix from a file.
-
-        Args:
-            filename (str): path to input file.
-            filetype (str): the filetype.
-
-                            * ``'txt'`` - a 2D matrix array in text format.
-                            * ``'bin'`` - a PETSc matrix vector.
-
-            delimiter (str): this is passed to `numpy.genfromtxt\
-                <http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_
-                in the case of strange delimiters in an imported ``txt`` file.
-        """
-    if filetype == 'txt':
-        try:
-            try:
-                if delimiter is None:
-                    matArray = _np.genfromtxt(filename,dtype=_PETSc.ScalarType)
-                else:
-                    matArray = _np.genfromtxt(filename,dtype=_PETSc.ScalarType,delimiter=delimiter)
-            except:
-                filefix = []
-                for line in _fl.FileInput(filename,inplace=0):
-                    if line[2] != 't':
-                        line = line.replace(" i","j")
-                        line = line.replace(" -","-")
-                        line = line.replace("+-","-")
-                        filefix.append(line)
-
-                matArray = _np.genfromtxt(filefix,dtype=_PETSc.ScalarType)
-
-            return arrayToMat(matArray)
-        except:
-            print "\nERROR: input state space file " + filename\
-                + " does not exist or is in an incorrect format"
-            _sys.exit()
-
-    elif filetype == 'bin':
-        binLoad = _PETSc.Viewer().createBinary(filename, 'r')
-        try:
-            return _PETSc.Mat().load(binLoad)
-        except:
-            print "\nERROR: input state space file " + filename\
-                + " does not exist or is in an incorrect format"
-            _sys.exit()
-        binLoad.destroy()
-
-def exportVecToMat(vec,filename,filetype):
-    """ Export a :math:`N^2` element PETSc vector as a :math:`N\\times N` matrix.
-
-        This is useful when wanting to view the full statespace of a 2 particle
-        quantum walk.
-
-        Args:
-            vec (petsc4py.PETSc.Vec): input :math:`N^2` element vector.
-            filename (str): path to desired output file.
-            filetype (str): the filetype of the exported vector.
-
-                            * ``'txt'`` - an :math:`N\\times N` 2D matrix array in text format.
-                            * ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
-        """
-    rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
-
-    if _os.path.isabs(filename):
-        outDir = _os.path.dirname(filename)
-    else:
-        outDir = './'+_os.path.dirname(filename)
-
-    # create output directory if it doesn't exist
-    try:
-        _os.mkdir(outDir)
-    except OSError as exception:
-        if exception.errno != _errno.EEXIST:
-            raise
-
-    vecArray = vecToArray(vec)
-    matArray = vecArray.reshape([_np.sqrt(vecArray.size),_np.sqrt(vecArray.size)])
-
-    if filetype == 'txt':
-        #if rank == 0:  _np.savetxt(filename,matArray)
-        txtSave = _PETSc.Viewer().createASCII(filename, 'w',
-            format=_PETSc.Viewer.Format.ASCII_DENSE, comm=_PETSc.COMM_WORLD)
-        txtSave(arrayToMat(matArray))
-        txtSave.destroy()
-
-        if rank == 0:
-            for line in _fl.FileInput(filename,inplace=1):
-                if line[2] != 't':
-                    line = line.replace(" i","j")
-                    line = line.replace(" -","-")
-                    line = line.replace("+-","-")
-                    print line,
-
-    elif filetype == 'bin':
-        binSave = _PETSc.Viewer().createBinary(filename, 'w', comm=_PETSc.COMM_WORLD)
-        binSave(arrayToMat(matArray))
-        binSave.destroy()
-    vec.comm.barrier()
-
-
-def loadMatToVec(filename,filetype):
-    """ Load a :math:`N\\times N` matrix as a :math:`N^2` element PETSc vector.
-
-        This is useful when wanting to import the full statespace of a 2 particle
-        quantum walk to use for propagation.
-
-        Args:
-            filename (str): path to the input file.
-            filetype (str): the filetype
-
-                            * ``'txt'`` - an :math:`N\\times N` 2D matrix array in text format.
-                            * ``'bin'`` - **Not yet implemented! Please use a txt \
-                                         format for this type of import**.
-        """
-    if filetype == 'txt':
-        try:
-            try:
-                matArray = _np.loadtxt(filename,dtype=_PETSc.ScalarType)
-            except:
-                filefix = []
-                for line in _fl.FileInput(filename,inplace=0):
-                    if line[2] != 't':
-                        line = line.replace(" i","j")
-                        line = line.replace(" -","-")
-                        line = line.replace("+-","-")
-                        filefix.append(line)
-
-                matArray = _np.loadtxt(filefix,dtype=_PETSc.ScalarType)
-
-            vecArray = matArray.reshape(matArray.shape[0]**2)
-            return arrayToVec(vecArray)
-        except:
-            print "\nERROR: input state space file " + filename\
-                + " does not exist or is in an incorrect format"
-            _sys.exit()
-
-    elif filetype == 'bin':
-        print '\nERROR: only works for txt storage!'
-        _sys.exit()

GCEtools/gceStartup.txt

@@ -1,11 +0,0 @@
-# Check project status
-gcutil getproject --project=<ProjectName> --cache_flag_values
-
-# Start an instance
-gcutil addinstance <instanceName>
-
-# Log in 
-gcutil ssh <instanceName>
-
-# Shut down
-gcutil deleteinstance <instanceName>

GCEtools/petscStartup.sh

@@ -1,145 +0,0 @@
-#! /bin/bash
-
-locale-gen en_US en_US.UTF-8 hu_HU hu_HU.UTF-8 > output.t
-dpkg-reconfigure locales >> output.t
-
-sudo apt-get update  >> output.t
-echo " "
-echo " "
-echo "   ============================================"
-echo "   | Installing packages form package manager |"
-echo "   ============================================"
-echo " "
-echo " "
-
-sudo apt-get -y install aptitude >> output.t
-
-packages=(gcc gfortran git libopenmpi-dev python-pip python-dev git flex bison cmake vim cython ipython python-scipy python-numpy python-nose python-pip python-matplotlib python-vtk python-h5py libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev libblas-dev liblapack-dev )
-
-
-for item in ${packages[*]}
-do
-    printf "     %-30s\n" $item
-
-done
-
-for item in ${packages[*]}
-do
-    tput cuu1
-done
-
-for item in ${packages[*]}
-do
-    sudo aptitude -y install $item >> output.t
-    printf "     %-30s              %-4s\n" $item done
-done
-
-
-echo " "
-echo " "
-echo "   ====================================="
-echo "   | Installing extra Python libraries |"
-echo "   ====================================="
-echo " "
-echo " "
-
-
-pipPackages=(mpi4py pymumps)
-
-for item in ${pipPackages[*]}
-do
-    printf "     %-30s\n" $item
-done
-
-for item in ${pipPackages[*]}
-do
-    tput cuu1
-done
-
-for item in ${pipPackages[*]}
-do
-        sudo pip install $item >> output.t
-        printf "     %-30s              %-4s\n" $item done
-done
-
-Upgrade=(scipy numpy ipython)
-
-for item in ${Upgrade[*]}
-do
-    printf "     %-8s%-7s\n" $item upgrade
-
-done
-
-for item in ${Upgrade[*]}
-do
-    tput cuu1
-done
-
-for item in ${Upgrade[*]}
-do
-    sudo pip install $item --upgrade >> output.t
-    printf "     %-8s%-7s                             %-4s\n" $item upgrade done
-done
-
-
-
-echo " "
-echo " "
-echo "   ====================="
-echo "   | Installing SimPEG |"
-echo "   ====================="
-echo " "
-echo " "
-cd ~
-
-
-git clone https://github.com/simpeg/simpeg.git >> output.t
-cd simpeg/SimPEG/
-python setup.py >> output.t
-cd ~
-
-mkdir petsc
-cd petsc
-
-echo " "
-echo " "
-echo "   ===================="
-echo "   | Installing PETSc |"
-echo "   ===================="
-echo " "
-echo " "
-wget http://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.4.3.tar.gz
-
-tar -zxf petsc-3.4.3.tar.gz
-
-cd petsc-3.4.3
-
-./configure --with-debugging=no  --dowload-mpich=yes --download-blacs=yes --download-f-blas-lapack=yes --download-scalapack=yes --download-mumps=yes --download-ml=yes --download-spooles=yes --download-hypre=yes --dowload-trilinos=yes --download-metis=yes --download-parmetis=yes --download-umfpack=yes --download-ptscotch=yes  --download-superlu=yes --download-superlu_dist=yes --download-essl=yes --download-eucild=yes  --download-spai=yes  --download-mpi4py=yes --download-petsc4py=yes --download-scientificpython=yes
-
-
-echo "export PETSC_DIR=/home/${USER}/petsc/petsc-3.4.3" >> ~/.bashrc
-echo "export PETSC_ARCH=arch-linux2-c-opt" >> ~/.bashrc
-export PETSC_DIR=/home/${USER}/petsc/petsc-3.4.3
-export PETSC_ARCH=arch-linux2-c-opt
-. ~/.bashrc
-
-make PETSC_DIR=/home/${USER}/petsc/petsc-3.4.3 PETSC_ARCH=arch-linux2-c-opt all
-make PETSC_DIR=/home/${USER}/petsc/petsc-3.4.3 PETSC_ARCH=arch-linux2-c-opt test
-
-cd ~/petsc
-echo " "
-echo " "
-echo "   ======================="
-echo "   | Installing PETSc4PY |"
-echo "   ======================="
-echo " "
-echo " "
-git clone https://bitbucket.org/petsc/petsc4py.git
-cd petsc4py/
-python setup.py build >> output.t
-python setup.py install --prefix=~/petsc >> output.t
-
-echo "export PYTHONPATH=~/petsc/lib/python2.7/site-packages:/home/$USER/simpeg:${PYTHONPATH}" >> ~/.bashrc
-
-cd ~
-source ~/.bashrc

GCEtools/startup.sh

@@ -1,22 +0,0 @@
-#! /bin/bash
-sudo aptitude -y update
-sudo aptitude -y upgrade
-sudo aptitude -y install gcc gfortran git libopenmpi-dev python-pip python-dev
-sudo aptitude -y install ipython python-scipy python-numpy python-nose python-pip python-matplotlib
-sudo aptitude -y install libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev
-sudo aptitude -y install libblas-dev liblapack-dev
-
-sudo pip install mpi4py
-sudo pip install pymumps
-
-sudo pip install scipy --upgrade
-sudo pip install numpy --upgrade
-sudo pip install ipython --upgrade
-
-git clone https://github.com/simpeg/simpeg.git
-cd simpeg/SimPEG/
-python setup.py
-cd ~
-
-echo export PYTHONPATH=/home/$USER/simpeg/ >> .bashrc
-source .bashrc

LICENSE

@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2013-2015 SimPEG Developers
+Copyright (c) 2013-2016 SimPEG Developers
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in

PROJECTS.rst

@@ -1,36 +0,0 @@
-- Electromagnetics (`simpegEM <http://simpegem.rtfd.org/>`_)
-    .. image:: https://travis-ci.org/simpeg/simpegem.svg?branch=master
-        :target: https://travis-ci.org/simpeg/simpegem
-        :alt: Master Branch
-    .. image:: https://coveralls.io/repos/simpeg/simpegem/badge.png?branch=master
-        :target: https://coveralls.io/r/simpeg/simpegem?branch=master
-- Potential Fields (`simpegPF <http://simpegpf.rtfd.org/>`_)
-    .. image:: https://travis-ci.org/simpeg/simpegpf.svg?branch=master
-        :target: https://travis-ci.org/simpeg/simpegpf
-        :alt: Master Branch
-    .. image:: https://coveralls.io/repos/simpeg/simpegpf/badge.png?branch=master
-        :target: https://coveralls.io/r/simpeg/simpegpf?branch=master
-- Ground Water Flow (`simpegFLOW <http://simpegflow.rtfd.org/>`_)
-    .. image:: https://travis-ci.org/simpeg/simpegflow.svg?branch=master
-        :target: https://travis-ci.org/simpeg/simpegflow
-        :alt: Master Branch
-    .. image:: https://coveralls.io/repos/simpeg/simpegflow/badge.png?branch=master
-        :target: https://coveralls.io/r/simpeg/simpegflow?branch=master
-- Direct Current Resistivity (`simpegDC <http://simpeg-dc.rtfd.org/>`_)
-    .. image:: https://travis-ci.org/simpeg/simpegdc.svg?branch=master
-        :target: https://travis-ci.org/simpeg/simpegdc
-        :alt: Master Branch
-    .. image:: https://coveralls.io/repos/simpeg/simpegdc/badge.png?branch=master
-        :target: https://coveralls.io/r/simpeg/simpegdc?branch=master
-- Electromagnetics 1D (`simpegEM1D <http://simpegem1d.rtfd.org/>`_)
-    .. image:: https://travis-ci.org/simpeg/simpegEM1D.svg?branch=master
-        :target: https://travis-ci.org/simpeg/simpegEM1D
-        :alt: Master Branch
-    .. image:: https://coveralls.io/repos/simpeg/simpegEM1D/badge.png?branch=master
-        :target: https://coveralls.io/r/simpeg/simpegEM1D?branch=master
-- Magnetotellurics (`simpegMT <http://simpegmt.rtfd.org/>`_)
-    .. image:: https://travis-ci.org/simpeg/simpegmt.svg?branch=master
-        :target: https://travis-ci.org/simpeg/simpegmt
-        :alt: Master Branch
-    .. image:: https://coveralls.io/repos/simpeg/simpegmt/badge.png?branch=master
-        :target: https://coveralls.io/r/simpeg/simpegmt?branch=master

README.rst

@@ -17,7 +17,7 @@ SimPEG
     :target: https://github.com/simpeg/simpeg/blob/master/LICENSE
     :alt: BSD 3 clause license.
 
-.. image:: https://img.shields.io/travis/simpeg/simpeg.svg
+.. image:: https://api.travis-ci.org/simpeg/simpeg.svg?branch=master
     :target: https://travis-ci.org/simpeg/simpeg
     :alt: Travis CI build status
 

SimPEG/Directives.py

@@ -206,17 +206,84 @@ 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"""
+
+    def initialize(self):
+        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.smoothModel == 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'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
 
 
-# class UpdateReferenceModel(Parameter):
+class update_IRLS(InversionDirective):
 
-#     mref0 = None
+     m = None
+     eps_min = None
+     factor = None
+     gamma = None
+     phi_m_last = None
+     
+     def initialize(self):
+        
+         # Scale the regularization for changes in norm
+         if getattr(self, 'phi_m_last', None) is not None:
+             self.reg.gamma = 1.
+             phim_new = self.reg.eval(self.invProb.curModel)
+             self.gamma = self.phi_m_last / phim_new
+         
+         self.reg.gamma = self.gamma
+         
+     def endIter(self):
+         # Cool the threshold parameter 
+         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
+                 
+         
+         # Update the model used for the IRLS weights
+         if getattr(self, 'm', None) is None:
+             self.reg.m = self.invProb.curModel
+         
+         # Update the pre-conditioner
+         diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.prob.mesh.nC))**2.
+         PC     = Utils.sdiag(diagA**-1.)
 
-#     def nextIter(self):
-#         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
+         self.opt.approxHinv = PC
+         
+         phim_new = self.reg.eval(self.invProb.curModel)
+         self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new
+
+#==============================================================================
+#          import pylab as plt
+#          plt.figure()
+#          ax = plt.subplot(221)
+#          self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
+#==============================================================================

SimPEG/EM/Analytics/FDEM.py

@@ -2,7 +2,6 @@ from __future__ import division
 import numpy as np
 from scipy.constants import mu_0, pi
 from scipy.special import erf
-import matplotlib.pyplot as plt
 from SimPEG import Utils
 
 
@@ -59,8 +58,9 @@ def MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu
 
         from SimPEG import EM
         import matplotlib.pyplot as plt
+        from scipy.constants import mu_0
         freqs = np.logspace(-2,5,100)
-        Bx, By, Bz = EM.Analytics.FDEM.AnalyticMagDipoleWholeSpace([0,100,0], [0,0,0], 1e-2, freqs, m=1, orientation='Z')
+        Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace([0,100,0], [0,0,0], 1e-2, freqs, moment=1, orientation='Z')
         plt.loglog(freqs, np.abs(Bz.real)/mu_0, 'b')
         plt.loglog(freqs, np.abs(Bz.imag)/mu_0, 'r')
         plt.legend(('real','imag'))

SimPEG/EM/FDEM/FDEM.py

@@ -50,7 +50,7 @@ class BaseFDEMProblem(BaseEMProblem):
             Srcs = self.survey.getSrcByFreq(freq)
             ftype = self._fieldType + 'Solution'
             F[Srcs, ftype] = sol
-
+            Ainv.clean()
         return F
 
     def Jvec(self, m, v, f=None):
@@ -89,6 +89,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
                     Jv[src, rx] = P(Df_Dm)
 
+            Ainv.clean()
         return Utils.mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
@@ -139,7 +140,8 @@ class BaseFDEMProblem(BaseEMProblem):
                         Jtv += - np.array(du_dmT,dtype=complex).real
                     else:
                         raise Exception('Must be real or imag')
-
+            
+            ATinv.clean()
         return Jtv
 
     def getSourceTerm(self, freq):

SimPEG/EM/FDEM/SrcFDEM.py

@@ -277,10 +277,10 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticDipoleVectorPotential
+            srcfct = MagneticLoopVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -26,6 +26,13 @@ class Rx(SimPEG.Survey.BaseRx):
                     'byi':['b', 'Fy', 'imag'],
                     'bzi':['b', 'Fz', 'imag'],
 
+                    'bxr_sec':['bSecondary', 'Fx', 'real'],
+                    'byr_sec':['bSecondary', 'Fy', 'real'],
+                    'bzr_sec':['bSecondary', 'Fz', 'real'],
+                    'bxi_sec':['bSecondary', 'Fx', 'imag'],
+                    'byi_sec':['bSecondary', 'Fy', 'imag'],
+                    'bzi_sec':['bSecondary', 'Fz', 'imag'],
+
                     'jxr':['j', 'Fx', 'real'],
                     'jyr':['j', 'Fy', 'real'],
                     'jzr':['j', 'Fz', 'real'],
@@ -106,7 +113,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
         SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
         _freqDict = {}
-        for src in srcList:
+        for src in self.srcList:
             if src.freq not in _freqDict:
                 _freqDict[src.freq] = []
             _freqDict[src.freq] += [src]

SimPEG/EM/TDEM/BaseTDEM.py

@@ -31,12 +31,20 @@ class FieldsTDEM(Problem.TimeFields):
 
 
 class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
-    """docstring for ProblemTDEM1D"""
+    """docstring for BaseTDEMProblem"""
     def __init__(self, mesh, mapping=None, **kwargs):
         BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
 
     _FieldsForward_pair = FieldsTDEM  #: used for the forward calculation only
 
+    waveformType = "STEPOFF"
+    current = None
+
+    def currentwaveform(self, wave):
+        self._timeSteps = np.diff(wave[:,0])
+        self.current = wave[:,1]
+        self.waveformType = "GENERAL"
+
     def fields(self, m):
         if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
         self.curModel = m

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -79,12 +79,32 @@ class SrcTDEM(Survey.BaseSrc):
 
 class SrcTDEM_VMD_MVP(SrcTDEM):
 
-    def __init__(self,rxList,loc):
+    def __init__(self,rxList,loc,waveformType="STEPOFF"):
         self.loc = loc
+        self.waveformType = waveformType
         SrcTDEM.__init__(self,rxList)
 
     def getInitialFields(self, mesh):
         """Vertical magnetic dipole, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"        
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+            else:
+                raise Exception('Unknown mesh for VMD') 
+            return {"b": mesh.edgeCurl*MVP}               
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
                 MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -93,20 +113,38 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
         else:
-            raise Exception('Unknown mesh for VMD')
-
-        return {"b": mesh.edgeCurl*MVP}
+            raise Exception('Unknown mesh for VMD')            
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 
 class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-
-    def __init__(self,rxList,loc,radius):
+    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
         self.loc = loc
         self.radius = radius
-        SrcTDEM.__init__(self,rxList)
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)        
 
     def getInitialFields(self, mesh):
         """Circular Loop, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            else:
+                raise Exception('Unknown mesh for CircularLoop')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
                 MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -115,9 +153,8 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
         else:
-            raise Exception('Unknown mesh for CircularLoop')
-
-        return {"b": mesh.edgeCurl*MVP}
+            raise Exception('Unknown mesh for CircularLoop')            
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 
 class SurveyTDEM(Survey.BaseSurvey):

SimPEG/Examples/DC_PseudoSection_Simulation.py

@@ -0,0 +1,179 @@
+from SimPEG import *
+import simpegDCIP as DC
+import scipy.interpolate as interpolation
+import matplotlib.pyplot as plt
+import time
+import re
+
+def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp',  plotIt=True):
+    """
+    DC Forward Simulation
+    
+    Forward model conductive spheres in a half-space and plot a pseudo-section
+    
+    Created on Mon Feb 01 19:28:06 2016
+    
+    @fourndo
+    """
+    
+    # First we need to create a mesh and a model.
+    
+    # This is our mesh
+    dx    = 5.
+    
+    hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)]
+    hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)]
+    hzind = [(dx,15,-1.3),(dx, 15)]
+    
+    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
+    
+    
+    # Set background conductivity
+    model = np.ones(mesh.nC) * sig[0]
+    
+    # First anomaly
+    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC)
+    model[ind] = sig[1]
+    
+    # Second anomaly
+    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC)
+    model[ind] = sig[2]
+    
+    # Get index of the center
+    indy = int(mesh.nCy/2)
+
+    
+    # Plot the model for reference
+    # Define core mesh extent
+    xlim = 200
+    zlim = 125
+    
+    # Specify the survey type: "pdp" | "dpdp"
+    
+    
+    # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
+    ends = [(-175,0),(175,0)]
+    ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]]
+    
+    # Snap the endpoints to the grid. Easier to create 2D section.
+    indx = Utils.closestPoints(mesh, ends )
+    locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
+    
+    # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
+    [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
+    
+    # Define some global geometry
+    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) 
+    dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
+    dl_y = ( Tx[-1][1,1] - Tx[0][1,0]  ) / dl_len
+    azm =  np.arctan(dl_y/dl_x)
+
+    #Set boundary conditions
+    mesh.setCellGradBC('neumann')
+    
+    # Define the differential operators needed for the DC problem
+    Div = mesh.faceDiv
+    Grad = mesh.cellGrad
+    Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
+    
+    A = Div*Msig*Grad
+    
+    # Change one corner to deal with nullspace
+    A[0,0] = 1
+    A = sp.csc_matrix(A)
+    
+    # We will solve the system iteratively, so a pre-conditioner is helpful
+    # This is simply a Jacobi preconditioner (inverse of the main diagonal)
+    dA = A.diagonal()
+    P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
+
+    # Now we can solve the system for all the transmitters
+    # We want to store the data
+    data = []
+    
+    # There is probably a more elegant way to do this, but we can just for-loop through the transmitters
+    for ii in range(len(Tx)):
+        
+        start_time = time.time() # Let's time the calculations
+        
+        #print("Transmitter %i / %i\r" % (ii+1,len(Tx)))
+        
+        # Select dipole locations for receiver
+        rxloc_M = np.asarray(Rx[ii][:,0:3])
+        rxloc_N = np.asarray(Rx[ii][:,3:])
+
+      
+        # For usual cases "dpdp" or "gradient"
+        if not re.match(stype,'pdp'): 
+            inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
+            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )   
+            
+        else: 
+            
+            # Create an "inifinity" pole
+            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] )  
+    
+    
+        # Iterative Solve
+        Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
+    
+        # We now have the potential everywhere
+        phi = mkvc(Ainvb[0])
+        
+        # Solve for phi on pole locations
+        P1 = mesh.getInterpolationMat(rxloc_M, 'CC')
+        P2 = mesh.getInterpolationMat(rxloc_N, 'CC')
+            
+        # Compute the potential difference
+        dtemp = (P1*phi - P2*phi)*np.pi
+        
+        data.append( dtemp )     
+        print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
+        
+    print '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
+    [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(Tx,Rx)
+    
+
+    # Here is an example for the first tx-rx array
+    if plotIt:
+        fig = plt.figure()
+        ax = plt.subplot(2,1,1, aspect='equal')
+        mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
+        ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
+        plt.gca().set_aspect('equal', adjustable='box')
+
+        plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
+        plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y')
+        plt.xlim([-xlim,xlim])
+        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+        
+        
+        ax = plt.subplot(2,1,2, aspect='equal')   
+        
+        # Plot the location of the spheres for reference
+        circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
+        ax.add_artist(circle1)
+        ax.add_artist(circle2)  
+        
+        # Add the speudo section
+        DC.plot_pseudoSection(Tx2d,Rx2d,data,mesh.vectorNz[-1],stype)
+        
+        plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
+        plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
+        plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
+        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+        
+        plt.show()
+        
+        return fig, ax
+        
+if __name__ == '__main__':
+    run()

SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py

@@ -0,0 +1,43 @@
+from SimPEG import *
+import SimPEG.EM as EM
+
+def run(XYZ=None, loc=np.r_[0.,0.,0.], sig=1.0, freq=1.0, orientation='Z', plotIt=True):
+    """
+        EM: Magnetic Dipole in a Whole-Space
+        ====================================
+
+        Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
+
+    """
+
+    if XYZ is None:
+        x = np.arange(-100.5,100.5,step = 1.) #(avoid putting measurement points where source is located)
+        y = np.r_[0]
+        z = x
+        XYZ = Utils.ndgrid(x,y,z)
+
+
+    Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, loc, sig, freq, orientation=orientation)
+    absB = np.sqrt(Bx*Bx.conj()+By*By.conj()+Bz*Bz.conj()).real
+
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        from matplotlib.colors import LogNorm
+        fig, ax = plt.subplots(1,1,figsize=(6,5))
+        bxplt = Bx.reshape(x.size,z.size)
+        bzplt = Bz.reshape(x.size,z.size)
+        pc = ax.pcolor(x,z,absB.reshape(x.size,z.size),norm=LogNorm())
+        ax.streamplot(x,z,bxplt.real,bzplt.real,color='k',density=1)
+        ax.set_xlim([x.min(),x.max()])
+        ax.set_ylim([z.min(),z.max()])
+        ax.set_xlabel('x')
+        ax.set_ylabel('z')
+        cb = plt.colorbar(pc,ax = ax)
+        cb.set_label('|B| (T)')
+        plt.show()
+
+        return fig, ax
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -1,14 +1,14 @@
 from SimPEG import *
 import SimPEG.EM as EM
 from scipy.constants import mu_0
-import matplotlib.pyplot as plt
+
 
 def run(plotIt=True):
     """
-        EM: FDEM: 1D: Inversion
+        EM: TDEM: 1D: Inversion
         =======================
 
-        Here we will create and run a FDEM 1D inversion.
+        Here we will create and run a TDEM 1D inversion.
 
     """
 
@@ -31,6 +31,7 @@ def run(plotIt=True):
 
 
     if plotIt:
+        import matplotlib.pyplot as plt
         fig, ax = plt.subplots(1,1, figsize = (3, 6))
         plt.semilogx(sigma[active], mesh.vectorCCz[active])
         ax.set_ylim(-600, 0)
@@ -60,6 +61,7 @@ def run(plotIt=True):
     survey.Wd = 1/(abs(survey.dobs)*std)
 
     if plotIt:
+        import matplotlib.pyplot as plt
         fig, ax = plt.subplots(1,1, figsize = (10, 6))
         ax.loglog(rx.times, dtrue, 'b.-')
         ax.loglog(rx.times, survey.dobs, 'r.-')
@@ -88,6 +90,7 @@ def run(plotIt=True):
     mopt = inv.run(m0)
 
     if plotIt:
+        import matplotlib.pyplot as plt
         fig, ax = plt.subplots(1,1, figsize = (3, 6))
         plt.semilogx(sigma[active], mesh.vectorCCz[active])
         plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])

SimPEG/Examples/FLOW_Richards_1D_Celia1990.py

@@ -1,6 +1,5 @@
 from SimPEG import *
 from SimPEG.FLOW import Richards
-import matplotlib.pyplot as plt
 
 def run(plotIt=True):
     """
@@ -61,6 +60,7 @@ def run(plotIt=True):
     Hs_H120= getFields(120.,'head')
 
     if not plotIt:return
+    import matplotlib.pyplot as plt
     plt.figure(figsize=(13,5))
     plt.subplot(121)
     plt.plot(40-M.gridCC, Hs_M10[-1],'b-')

SimPEG/Examples/Forward_BasicDirectCurrent.py

@@ -1,7 +1,4 @@
 from SimPEG import Mesh, Utils, np, SolverLU
-import matplotlib.pyplot as plt
-import matplotlib
-from matplotlib.mlab import griddata
 
 ## 2D DC forward modeling example with Tensor and Curvilinear Meshes
 
@@ -39,6 +36,10 @@ def run(plotIt=True):
 
     if not plotIt: return
 
+    import matplotlib.pyplot as plt
+    import matplotlib
+    from matplotlib.mlab import griddata
+
     #Step4: Making Figure
     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
     label = ["(a)", "(b)"]

SimPEG/Examples/__init__.py

@@ -1,7 +1,8 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import EM_FDEM_1D_Inversion
+import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
 import Inversion_Linear
@@ -13,7 +14,7 @@ import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
 
-__examples__ = ["EM_FDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
+__examples__ = ["EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
 
 ##### AUTOIMPORTS #####
 
@@ -28,16 +29,17 @@ if __name__ == '__main__':
     from SimPEG import Examples
 
     # Create the examples dir in the docs folder.
-    docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples'])
+    fName = os.path.realpath(__file__)
+    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
     shutil.rmtree(docExamplesDir)
     os.makedirs(docExamplesDir)
 
     # Get all the python examples in this folder
-    thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1])
+    thispath = os.path.sep.join(fName.split(os.path.sep)[:-1])
     exfiles  = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')]
 
     # Add the imports to the top in the AUTOIMPORTS section
-    f = file(__file__, 'r')
+    f = file(fName, 'r')
     inimports = False
     out = ''
     for line in f:
@@ -52,7 +54,7 @@ if __name__ == '__main__':
                 out += '\n##### AUTOIMPORTS #####\n'
     f.close()
 
-    f = file(__file__, 'w')
+    f = file(fName, 'w')
     f.write(out)
     f.close()
 

SimPEG/Exceptions.py

@@ -0,0 +1,13 @@
+
+
+class SimPEGException(Exception):
+
+    def __init__(self, reason=''):
+        self.reason = reason
+
+    def __str__(self):
+        return '%s: %s' %(self.__class__.__name__,  self.reason)
+
+
+class PairingException(SimPEGException):
+    pass

SimPEG/InvProblem.py

@@ -66,8 +66,8 @@ class BaseInvProblem(object):
         self.curModel = m0
 
         print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
-                    ***Done using same solver as the problem***"""
-        self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel))
+                    ***Done using same Solver and solverOpts as the problem***"""
+        self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts)
 
     @property
     def warmstart(self):

SimPEG/Maps.py

@@ -4,27 +4,32 @@ from Tests import checkDerivative
 from PropMaps import PropMap, Property
 from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
+from scipy.spatial import cKDTree
 
 class IdentityMap(object):
     """
     SimPEG Map
 
     """
-
     __metaclass__ = Utils.SimPEGMetaClass
 
-    mesh = None      #: A SimPEG Mesh
-
-    def __init__(self, mesh, **kwargs):
+    def __init__(self, mesh=None, nP=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
+
+        if nP is not None:
+            assert type(nP) in [int, long], ' Number of parameters must be an integer.'
+
         self.mesh = mesh
+        self._nP  = nP
 
     @property
     def nP(self):
         """
             :rtype: int
-            :return: number of parameters in the model
+            :return: number of parameters that the mapping accepts
         """
+        if self._nP is not None:
+            return self._nP
         if self.mesh is None:
             return '*'
         return self.mesh.nC
@@ -32,11 +37,15 @@ class IdentityMap(object):
     @property
     def shape(self):
         """
-            The default shape is (mesh.nC, nP).
+            The default shape is (mesh.nC, nP) if the mesh is defined.
+            If this is a meshless mapping (i.e. nP is defined independently)
+            the shape will be the the shape (nP,nP).
 
             :rtype: (int,int)
             :return: shape of the operator as a tuple
         """
+        if self._nP is not None:
+            return (self.nP, self.nP)
         if self.mesh is None:
             return ('*', self.nP)
         return (self.mesh.nC, self.nP)
@@ -113,11 +122,12 @@ class IdentityMap(object):
             if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
             return self._transform(val)
-        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
+        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray! Not a %s'%type(val))
 
     def __str__(self):
         return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
 
+
 class ComboMap(IdentityMap):
     """Combination of various maps."""
 
@@ -449,6 +459,89 @@ class Mesh2Mesh(IdentityMap):
         return self.P
 
 
+class Mesh2MeshTopo(IdentityMap):
+    """
+        Takes a model on one mesh are translates it to another mesh
+        with consideration of topography
+
+    """
+    tree = None
+    nIterpPts = 6
+    P = None #: The CSR projection matrix.
+
+    def __init__(self, meshes, actinds, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+
+        assert type(meshes) is list, "meshes must be a list of two meshes"
+        assert len(meshes) == 2, "meshes must be a list of two meshes"
+        assert type(actinds) is list, "actinds must be a list of two meshes"
+        assert len(actinds) == 2, "actinds must be a list of two meshes"
+        assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
+
+        self.mesh  = meshes[0]
+        self.mesh2 = meshes[1]
+        self.actind = actinds[0]
+        self.actind2 = actinds[1]
+        self._createProjection()
+
+        # Old version using SimPEG interpolation
+        # self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
+
+    def genActiveindfromTopo(mesh, xyztopo):
+        #TODO: This possibly needs to be improved use vtk(?)
+        if mesh.dim==3:
+            nCxy = mesh.nCx*mesh.nCy
+            Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F')
+            Ftopo = NearestNDInterpolator(xyztopo[:,:2], xyztopo[:,2])
+            XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
+            XY.shape
+            topo = Ftopo(XY)
+            actind = []
+            for ixy in range(nCxy):
+                actind.append(topo[ixy] <= Zcc[ixy,:])
+        else:
+            raise NotImplementedError("Only 3D is working")
+
+        return Utils.mkvc(np.vstack(actind))
+
+    #Question .. is it only generated once?
+    def _createProjection(self):
+        """
+            KD Tree interpolation onto the active cells.
+        """
+        if self.tree==None:
+            self.tree = cKDTree(zip(self.mesh.gridCC[self.actind,0], self.mesh.gridCC[self.actind,1], self.mesh.gridCC[self.actind,2]))
+        d, inds = self.tree.query(zip(self.mesh2.gridCC[self.actind2,0],self.mesh2.gridCC[self.actind2,1],self.mesh2.gridCC[self.actind2,2]), k=self.nIterpPts)
+        # Not sure consideration of the volume ...
+        # vol = np.zeros((self.actind2.sum(), self.nIterpPts))
+        # for i in range(self.nIterpPts):
+        #     vol[:,i] = self.mesh.vol[inds[:,i]]
+        w = 1. / d**2
+        w = Utils.sdiag(1./np.sum(w, axis=1)) * (w)
+        I = Utils.mkvc(np.arange(inds.shape[0]).reshape([-1,1]).repeat(self.nIterpPts, axis=1))
+        J = Utils.mkvc(inds)
+        P = sp.coo_matrix( (Utils.mkvc(w),(I, J)), shape=(inds.shape[0], (self.actind).sum()) )
+        # self.P = Utils.sdiag(self.mesh2.vol[self.actind2])*P.tocsc()
+        self.P = P.tocsr()
+
+    @property
+    def shape(self):
+        """Number of parameters in the model."""
+        # return (self.mesh.nC, self.mesh2.nC)
+        return (self.actind2.sum(), self.actind.sum())
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        # return self.mesh2.nC
+        return self.actind2.sum()
+
+    def _transform(self, m):
+        return self.P*m
+
+    def deriv(self, m):
+        return self.P
+
 class ActiveCells(IdentityMap):
     """
         Active model parameters.
@@ -475,7 +568,7 @@ class ActiveCells(IdentityMap):
         else:
             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, range(inds.size))), shape=(self.nC, self.nP))
 
@@ -708,7 +801,7 @@ class PolyMap(IdentityMap):
         Parameterize the model space using a polynomials in a wholespace.
 
         ..math::
-            
+
             y = \mathbf{V} c
 
         Define the model as:
@@ -752,10 +845,10 @@ class PolyMap(IdentityMap):
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]            
-            Z = self.mesh.gridCC[:,2]            
+            Y = self.mesh.gridCC[:,1]
+            Z = self.mesh.gridCC[:,2]
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
             elif self.normal =='Y':
@@ -766,43 +859,43 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         else:
             raise(Exception("Only supports 2D"))
-                    
+
 
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
-    
+
     def deriv(self, m):
         alpha = self.slope
         sig1,sig2, c = m[0],m[1],m[2:]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
-        if self.mesh.dim == 2:            
+        if self.mesh.dim == 2:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
 
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
-                V = polynomial.polyvander(Y, len(c)-1)        
+                V = polynomial.polyvander(Y, len(c)-1)
             elif self.normal =='Y':
                 f = polynomial.polyval(X, c) - Y
-                V = polynomial.polyvander(X, len(c)-1)        
+                V = polynomial.polyvander(X, len(c)-1)
             else:
-                raise(Exception("Input for normal = X or Y or Z"))            
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
             Z = self.mesh.gridCC[:,2]
 
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
-                V = polynomial.polyvander2d(Y, Z, self.order)        
+                V = polynomial.polyvander2d(Y, Z, self.order)
             elif self.normal =='Y':
                 f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
-                V = polynomial.polyvander2d(X, Z, self.order)        
+                V = polynomial.polyvander2d(X, Z, self.order)
             elif self.normal =='Z':
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
-                V = polynomial.polyvander2d(X, Y, self.order)        
+                V = polynomial.polyvander2d(X, Y, self.order)
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
 
@@ -815,16 +908,16 @@ class PolyMap(IdentityMap):
 
         g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V
 
-        return sp.csr_matrix(np.c_[g1,g2,g3])  
+        return sp.csr_matrix(np.c_[g1,g2,g3])
 
 class SplineMap(IdentityMap):
 
     """SplineMap
 
-        Parameterize the boundary of two geological units using a spline interpolation 
+        Parameterize the boundary of two geological units using a spline interpolation
 
         ..math::
-            
+
             g = f(x)-y
 
         Define the model as:
@@ -849,7 +942,7 @@ class SplineMap(IdentityMap):
     def nP(self):
         if self.mesh.dim == 2:
             return np.size(self.pts)+2
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             return np.size(self.pts)*2+2
         else:
             raise(Exception("Only supports 2D and 3D"))
@@ -866,28 +959,28 @@ class SplineMap(IdentityMap):
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
             self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0)
-            if self.normal =='X':                
+            if self.normal =='X':
                 f = self.spl(Y) - X
             elif self.normal =='Y':
                 f = self.spl(X) - Y
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
 
-        # 3D: 
-        # Comments: 
+        # 3D:
+        # Comments:
         # Make two spline functions and link them using linear interpolation.
         # This is not quite direct extension of 2D to 3D case
         # Using 2D interpolation  is possible
 
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]            
+            Y = self.mesh.gridCC[:,1]
             Z = self.mesh.gridCC[:,2]
 
-            npts = np.size(self.pts)            
+            npts = np.size(self.pts)
             if np.mod(c.size, 2):
                 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)}
 
@@ -902,7 +995,7 @@ class SplineMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         else:
             raise(Exception("Only supports 2D and 3D"))
-                    
+
 
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
 
@@ -912,7 +1005,7 @@ class SplineMap(IdentityMap):
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
-        if self.mesh.dim == 2:            
+        if self.mesh.dim == 2:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
 
@@ -921,9 +1014,9 @@ class SplineMap(IdentityMap):
             elif self.normal =='Y':
                 f = self.spl(X) - Y
             else:
-                raise(Exception("Input for normal = X or Y or Z"))            
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
             Z = self.mesh.gridCC[:,2]
@@ -931,7 +1024,7 @@ class SplineMap(IdentityMap):
                 zb = self.ptsv[0]
                 zt = self.ptsv[1]
                 flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
-                f = flines - X            
+                f = flines - X
             # elif self.normal =='Y':
             # elif self.normal =='Z':
             else:
@@ -944,7 +1037,7 @@ class SplineMap(IdentityMap):
             g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
             g2 = (np.arctan(alpha*f)/np.pi + 0.5)
 
-        
+
         if self.mesh.dim ==2:
             g3 = np.zeros((self.mesh.nC, self.npts))
             if self.normal =='Y':
@@ -958,7 +1051,7 @@ class SplineMap(IdentityMap):
                     cb = c.copy()
                     dy = self.mesh.hy[ind]*1.5
                     ca[i] = ctemp+dy
-                    cb[i] = ctemp-dy                
+                    cb[i] = ctemp-dy
                     spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
                     splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
                     fderiv = (spla(X)-splb(X))/(2*dy)
@@ -968,7 +1061,7 @@ class SplineMap(IdentityMap):
             g3 = np.zeros((self.mesh.nC, self.npts*2))
             if self.normal =='X':
                 # Here we use perturbation to compute sensitivity
-                for i in range(self.npts*2):                    
+                for i in range(self.npts*2):
                     ctemp = c[i]
                     ind = np.argmin(abs(self.mesh.vectorCCy-ctemp))
                     ca = c.copy()
@@ -982,20 +1075,20 @@ class SplineMap(IdentityMap):
                         splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0)
                         flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X
                         flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X
-                    #treat top boundary                        
+                    #treat top boundary
                     else:
                         splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
                         spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
                         flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X
-                        flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X                        
-                    fderiv = (flinesa-flinesb)/(2*dy)                                
+                        flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
+                    fderiv = (flinesa-flinesb)/(2*dy)
                     g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
         else :
             raise(Exception("Not Implemented for Y and Z, your turn :)"))
-        return sp.csr_matrix(np.c_[g1,g2,g3])  
+        return sp.csr_matrix(np.c_[g1,g2,g3])
+
 
 
 
 
 
-            

SimPEG/Mesh/DiffOperators.py

@@ -565,7 +565,58 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
-
+    
+    def unitCellGradx():
+        doc = """Cell centered Gradient in the x dimension used for 
+                regularization. The gradient operator is square (nC-by-nC)"""
+        def fget(self):
+            if self.dim < 3: return None
+            if getattr(self, '_unitCellGradx', None) is None:
+                
+                n = self.vnC
+                gx = ddx(n[0]-1)
+                gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr")
+                
+                self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square)
+        
+            return self._unitCellGradx
+        return locals()
+    unitCellGradx = property(**unitCellGradx())
+    
+    def unitCellGrady():
+        doc = """Cell centered Gradient in they dimension used for 
+                regularization. The gradient operator is square (nC-by-nC)"""
+        def fget(self):
+            if self.dim < 3: return None
+            if getattr(self, '_unitCellGrady', None) is None:
+        
+                n = self.vnC
+                gy = ddx(n[1]-1)
+                gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr")
+                
+                self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0]))
+            
+            return self._unitCellGrady
+        return locals()
+    unitCellGrady = property(**unitCellGrady())
+        
+    def unitCellGradz():
+        doc = """Cell centered Gradient in they dimension used for 
+                regularization. The gradient operator is square (nC-by-nC)"""
+        def fget(self):
+            if self.dim < 3: return None
+            if getattr(self, '_unitCellGradz', None) is None:
+        
+                n = self.vnC
+                gz = ddx(n[2]-1)
+                gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr")
+                
+                self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0]))
+            
+            return self._unitCellGradz
+        return locals()
+    unitCellGradz = property(**unitCellGradz())        
+        
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Mesh/MeshIO.py

@@ -0,0 +1,416 @@
+import numpy as np, os
+from SimPEG import Utils
+
+class TensorMeshIO(object):
+
+    @classmethod
+    def readUBC(TensorMesh, fileName):
+        """
+            Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
+
+            Input:
+            :param fileName, path to the UBC GIF mesh file
+
+            Output:
+            :param SimPEG TensorMesh object
+        """
+
+        # Interal function to read cell size lines for the UBC mesh files.
+        def readCellLine(line):
+            for seg in line.split():
+                if '*' in seg:
+                    st = seg
+                    sp = seg.split('*')
+                    re = np.array(sp[0],dtype=int)*(' ' + sp[1])
+                    line = line.replace(st,re.strip())
+            return np.array(line.split(),dtype=float)
+
+        # Read the file as line strings, remove lines with comment = !
+        msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
+
+        # Fist line is the size of the model
+        sizeM = np.array(msh[0].split(),dtype=float)
+        # Second line is the South-West-Top corner coordinates.
+        x0 = np.array(msh[1].split(),dtype=float)
+        # Read the cell sizes
+        h1 = readCellLine(msh[2])
+        h2 = readCellLine(msh[3])
+        h3temp = readCellLine(msh[4])
+        h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom.
+        # Adjust the reference point to the bottom south west corner
+        x0[2] = x0[2] - np.sum(h3)
+        # Make the mesh
+        tensMsh = TensorMesh([h1,h2,h3],x0)
+        return tensMsh
+
+    @classmethod
+    def readVTK(TensorMesh, fileName):
+        """
+            Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
+
+            Input:
+            :param vtrFileName, path to the vtr model file to write to
+
+            Output:
+            :return SimPEG TensorMesh object
+            :return SimPEG model dictionary
+
+        """
+        # Import
+        from vtk import vtkXMLRectilinearGridReader as vtrFileReader
+        from vtk.util.numpy_support import vtk_to_numpy
+
+        # Read the file
+        vtrReader = vtrFileReader()
+        vtrReader.SetFileName(fileName)
+        vtrReader.Update()
+        vtrGrid = vtrReader.GetOutput()
+        # Sort information
+        hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates())))
+        xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0]
+        hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates())))
+        yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0]
+        zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates()))
+        # Check the direction of hz
+        if np.all(zD < 0):
+            hz = np.abs(zD[::-1])
+            zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1]
+        else:
+            hz = np.abs(zD)
+            zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0]
+        x0 = np.array([xR,yR,zR])
+
+        # Make the SimPEG object
+        tensMsh = TensorMesh([hx,hy,hz],x0)
+
+        # Grap the models
+        models = {}
+        for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()):
+            modelName = vtrGrid.GetCellData().GetArrayName(i)
+            if np.all(zD < 0):
+                modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
+                tM = tensMsh.r(modFlip,'CC','CC','M')
+                modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V')
+            else:
+                modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
+            models[modelName] = modArr
+
+        # Return the data
+        return tensMsh, models
+
+    def writeVTK(mesh, fileName, models=None):
+        """
+        Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
+
+        Input:
+        :param str, path to the output vtk file
+        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
+
+        """
+        # Import
+        from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION
+        from vtk.util.numpy_support import numpy_to_vtk
+
+        # Deal with dimensionalities
+        if mesh.dim >= 1:
+            vX = mesh.vectorNx
+            xD = mesh.nNx
+            yD,zD = 1,1
+            vY, vZ = np.array([0,0])
+        if mesh.dim >= 2:
+            vY = mesh.vectorNy
+            yD = mesh.nNy
+        if mesh.dim == 3:
+            vZ = mesh.vectorNz
+            zD = mesh.nNz
+        # Use rectilinear VTK grid.
+        # Assign the spatial information.
+        vtkObj = rectGrid()
+        vtkObj.SetDimensions(xD,yD,zD)
+        vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
+        vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
+        vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
+
+        # Assign the model('s) to the object
+        if models is not None:
+            for item in models.iteritems():
+                # Convert numpy array
+                vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
+                vtkDoubleArr.SetName(item[0])
+                vtkObj.GetCellData().AddArray(vtkDoubleArr)
+            # Set the active scalar
+            vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
+        # vtkObj.Update()
+
+        # Check the extension of the fileName
+        ext = os.path.splitext(fileName)[1]
+        if ext is '':
+            fileName = fileName + '.vtr'
+        elif ext not in '.vtr':
+            raise IOError('{:s} is an incorrect extension, has to be .vtr')
+        # Write the file.
+        vtrWriteFilter = rectWriter()
+        if float(VTK_VERSION.split('.')[0]) >=6:
+            vtrWriteFilter.SetInputData(vtkObj)
+        else:
+            vtuWriteFilter.SetInput(vtuObj)
+        vtrWriteFilter.SetFileName(fileName)
+        vtrWriteFilter.Update()
+
+
+    def readModelUBC(mesh, fileName):
+        """
+            Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
+
+            Input:
+            :param fileName, path to the UBC GIF mesh file to read
+            :param mesh, TensorMesh object, mesh that coresponds to the model
+
+            Output:
+            :return numpy array, model with TensorMesh ordered
+        """
+        f = open(fileName, 'r')
+        model = np.array(map(float, f.readlines()))
+        f.close()
+        model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F')
+        model = model[::-1,:,:]
+        model = np.transpose(model, (1, 2, 0))
+        model = Utils.mkvc(model)
+        return model
+
+    def writeModelUBC(mesh, fileName, model):
+        """
+            Writes a model associated with a SimPEG TensorMesh
+            to a UBC-GIF format model file.
+
+            :param str fileName: File to write to
+            :param simpeg.Mesh.TensorMesh mesh: The mesh
+            :param numpy.ndarray model: The model
+        """
+
+        # Reshape model to a matrix
+        modelMat = mesh.r(model,'CC','CC','M')
+        # Transpose the axes
+        modelMatT = modelMat.transpose((2,0,1))
+        # Flip z to positive down
+        modelMatTR = Utils.mkvc(modelMatT[::-1,:,:])
+
+        np.savetxt(fileName, modelMatTR.ravel())
+
+    def writeUBC(mesh, fileName, models=None):
+        """
+            Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
+
+            :param str fileName: File to write to
+            :param simpeg.Mesh.TensorMesh mesh: The mesh
+
+        """
+        assert mesh.dim == 3
+        s = ''
+        s += '%i %i %i\n' %tuple(mesh.vnC)
+        origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
+        origin.dtype = float
+
+        s += '%.2f %.2f %.2f\n' %tuple(origin)
+        s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
+        s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
+        s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
+        f = open(fileName, 'w')
+        f.write(s)
+        f.close()
+
+        if models is None: return
+        assert type(models) is dict, 'models must be a dict'
+        for key in models:
+            assert type(key) is str, 'The dict key is a file name'
+            mesh.writeModelUBC(key, models[key])
+
+class TreeMeshIO(object):
+
+    def writeUBC(mesh, fileName, models=None):
+        """
+            Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model.
+
+            :param str fileName: File to write to
+            :param simpeg.Mesh.TreeMesh mesh: The mesh
+            :param dictionary models: The models in a dictionary, where the keys is the name of the of the model file
+        """
+
+        # Calculate information to write in the file.
+        # Number of cells in the underlying mesh
+        nCunderMesh = np.array([h.size for h in mesh.h],dtype=np.int64)
+        # The top-south-west most corner of the mesh
+        tswCorn = mesh.x0 + np.array([0,0,np.sum(mesh.h[2])])
+        # Smallest cell size
+        smallCell = np.array([h.min() for h in mesh.h])
+        # Number of cells
+        nrCells = mesh.nC
+
+        ## Extract iformation about the cells.
+        # cell pointers
+        cellPointers = np.array([c._pointer for c in mesh])
+        # cell with
+        cellW = np.array([ mesh._levelWidth(i) for i in cellPointers[:,-1] ])
+        # Need to shift the pointers to work with UBC indexing
+        # UBC Octree indexes always the top-left-close (top-south-west) corner first and orders the cells in z(top-down),x,y vs x,y,z(bottom-up).
+        # Shift index up by 1
+        ubcCellPt = cellPointers[:,0:-1].copy() + np.array([1.,1.,1.])
+        # Need reindex the z index to be from the top-left-close corner and to be from the global top.
+        ubcCellPt[:,2] = ( nCunderMesh[-1] + 2) - (ubcCellPt[:,2] + cellW)
+
+        # Reorder the ubcCellPt
+        ubcReorder = np.argsort(ubcCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0]
+        # Make a array with the pointers and the withs, that are order in the ubc ordering
+        indArr = np.concatenate((ubcCellPt[ubcReorder,:],cellW[ubcReorder].reshape((-1,1)) ),axis=1)
+
+        ## Write the UBC octree mesh file
+        with open(fileName,'w') as mshOut:
+            mshOut.write('{:.0f} {:.0f} {:.0f}\n'.format(nCunderMesh[0],nCunderMesh[1],nCunderMesh[2]))
+            mshOut.write('{:.4f} {:.4f} {:.4f}\n'.format(tswCorn[0],tswCorn[1],tswCorn[2]))
+            mshOut.write('{:.3f} {:.3f} {:.3f}\n'.format(smallCell[0],smallCell[1],smallCell[2]))
+            mshOut.write('{:.0f} \n'.format(nrCells))
+            np.savetxt(mshOut,indArr,fmt='%i')
+
+        ## Print the models
+        # Assign the model('s) to the object
+        if models is not None:
+            # indUBCvector = np.argsort(cX0[np.argsort(np.concatenate((cX0[:,0:2],cX0[:,2:3].max() - cX0[:,2:3]),axis=1).view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]].view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]
+            for item in models.iteritems():
+                # Save the data
+                np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e')
+
+    @classmethod
+    def readUBC(TreeMesh, meshFile):
+        """
+            Read UBC 3D OcTree mesh and/or modelFiles
+
+            Input:
+            :param str meshFile: path to the UBC GIF OcTree mesh file to read
+
+            Output:
+            :return SimPEG.Mesh.TreeMesh mesh: The octree mesh
+            :return list of ndarray's: models as a list of numpy array's
+        """
+
+        ## Read the file lines
+        fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n')
+        # Extract the data
+        nCunderMesh = np.array(fileLines[0].split(),dtype=float)
+        # I think this is the case?
+        if np.unique(nCunderMesh).size >1:
+            raise Exception('SimPEG TreeMeshes have the same number of cell in all directions')
+        tswCorn = np.array(fileLines[1].split(),dtype=float)
+        smallCell = np.array(fileLines[2].split(),dtype=float)
+        nrCells = np.array(fileLines[3].split(),dtype=float)
+        # Read the index array
+        indArr = np.genfromtxt(fileLines[4::],dtype=np.int)
+
+        ## Calculate simpeg parameters
+        h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)]
+        x0 = tswCorn - np.array([0,0,np.sum(h3)])
+        # Need to convert the index array to a points list that complies with SimPEG TreeMesh.
+        # Shift to start at 0
+        simpegCellPt = indArr[:,0:-1].copy()
+        simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] + indArr[:,3])
+        # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom.
+        simpegCellPt = simpegCellPt - np.array([1.,1.,1.])
+
+        # Calculate the cell level
+        simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3])
+        # Make a pointer matrix
+        simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1)
+
+        ## Make the tree mesh
+        mesh = TreeMesh([h1,h2,h3],x0)
+        mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()])
+
+        # Figure out the reordering
+        mesh._simpegReorderUBC = np.argsort(np.array([mesh._index(i) for i in simpegPointers.tolist()]))
+        # mesh._simpegReorderUBC = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0]
+
+        return mesh
+
+
+    def readModelUBC(mesh, fileName):
+        """
+            Read UBC OcTree model and get vector
+
+            Input:
+            :param fileName, path to the UBC GIF model file to read
+
+            Output:
+            :return numpy array, OcTree model
+        """
+
+        if type(fileName) is list:
+            out = {}
+            for f in fileName:
+                out[f] = mesh.readModelUBC(f)
+            return out
+
+        assert hasattr(mesh, '_simpegReorderUBC'), 'The file must have been loaded from a UBC format.'
+        assert mesh.dim == 3
+
+        modList = []
+        modArr = np.loadtxt(fileName)
+        if len(modArr.shape) == 1:
+            modList.append(modArr[mesh._simpegReorderUBC])
+        else:
+            modList.append(modArr[mesh._simpegReorderUBC,:])
+        return modList
+
+    def writeVTK(mesh, fileName, models=None):
+        """
+        Function to write a VTU file from a SimPEG TreeMesh and model.
+        """
+        import vtk
+        from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION
+        from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray
+
+        if str(type(mesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh':
+            raise IOError('mesh is not a SimPEG TreeMesh.')
+
+        # Make the data parts for the vtu object
+        # Points
+        mesh.number()
+        ptsMat = mesh._gridN + mesh.x0
+
+        vtkPts = vtk.vtkPoints()
+        vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True))
+        # Cells
+        cellConn = np.array([c.nodes for c in mesh],dtype=np.int64)
+
+        cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel()
+        cellsArr = vtk.vtkCellArray()
+        cellsArr.SetNumberOfCells(cellConn.shape[0])
+        cellsArr.SetCells(cellConn.shape[0],numpy_to_vtkIdTypeArray(cellsMat,deep=True))
+
+        # Make the object
+        vtuObj = vtk.vtkUnstructuredGrid()
+        vtuObj.SetPoints(vtkPts)
+        vtuObj.SetCells(vtk.VTK_VOXEL,cellsArr)
+        # Add the level of refinement as a cell array
+        cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())])
+        uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True)
+        refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel,deep=1)
+        refineLevelArr.SetName('octreeLevel')
+        vtuObj.GetCellData().AddArray(refineLevelArr)
+        # Assign the model('s) to the object
+        if models is not None:
+            for item in models.iteritems():
+                # Convert numpy array
+                vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
+                vtkDoubleArr.SetName(item[0])
+                vtuObj.GetCellData().AddArray(vtkDoubleArr)
+
+        # Make the writer
+        vtuWriteFilter = Writer()
+        if float(VTK_VERSION.split('.')[0]) >=6:
+            vtuWriteFilter.SetInputData(vtuObj)
+        else:
+            vtuWriteFilter.SetInput(vtuObj)
+        vtuWriteFilter.SetFileName(fileName)
+        # Write the file
+        vtuWriteFilter.Update()
+

SimPEG/Mesh/TreeMesh.py

@@ -90,11 +90,6 @@
 #
 
 from SimPEG import np, sp, Utils, Solver
-import matplotlib.pyplot as plt
-import matplotlib
-from mpl_toolkits.mplot3d import Axes3D
-import matplotlib.colors as colors
-import matplotlib.cm as cmx
 
 try:
     import TreeUtils
@@ -105,11 +100,12 @@ except Exception, e:
 
 from InnerProducts import InnerProducts
 from TensorMesh import TensorMesh, BaseTensorMesh
+from MeshIO import TreeMeshIO
 import time
 
 MAX_BITS = 20
 
-class TreeMesh(BaseTensorMesh, InnerProducts):
+class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
 
     _meshType = 'TREE'
 
@@ -569,15 +565,18 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         return [p - (p % mod) for p in pointer[:-1]] + [pointer[-1]-1]
 
     def _cellN(self, p):
+        """Node location [x,y(,z)] of a single cell, closest to origin, given a pointer."""
         p = self._asPointer(p)
         return [hi[:p[ii]].sum() for ii, hi in enumerate(self.h)]
 
     def _cellH(self, p):
+        """Widths of a single cell given a pointer."""
         p = self._asPointer(p)
         w = self._levelWidth(p[-1])
         return [hi[p[ii]:p[ii]+w].sum() for ii, hi in enumerate(self.h)]
 
     def _cellC(self, p):
+        """Cell center of a single cell (without origin correction), given a pointer."""
         return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist()
 
     def _levelWidth(self, level):
@@ -832,8 +831,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
     def _numberCells(self, force=False):
         if not self.__dirtyCells__ and not force: return
         self._cc2i = dict()
+        self._i2cc = dict()
         for ii, c in enumerate(sorted(self._cells)):
             self._cc2i[c] = ii
+            self._i2cc[ii] = c
         self.__dirtyCells__ = False
 
     def _numberNodes(self, force=False):
@@ -1709,9 +1710,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         "Construct the averaging operator on cell faces to cell centers."
         if getattr(self, '_aveF2CC', None) is None:
             if self.dim == 2:
-                self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC])
+                self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]).tocsr()
             elif self.dim == 3:
-                self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC])
+                self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr()
         return self._aveF2CC
 
     @property
@@ -1719,9 +1720,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         "Construct the averaging operator on cell faces to cell centers."
         if getattr(self, '_aveF2CCV', None) is None:
             if self.dim == 2:
-                self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC])
+                self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]).tocsr()
             elif self.dim == 3:
-                self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC])
+                self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr()
         return self._aveF2CCV
 
     @property
@@ -1973,6 +1974,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         facesX=False, facesY=False, facesZ=False,
         edgesX=False, edgesY=False, edgesZ=False):
 
+
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        import matplotlib.colors as colors
+        import matplotlib.cm as cmx
+
         # self.number()
 
         axOpts = {'projection':'3d'} if self.dim == 3 else {}
@@ -2094,6 +2102,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
     def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
         if self.dim == 3: raise Exception('Use plot slice?')
 
+
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        import matplotlib.colors as colors
+        import matplotlib.cm as cmx
+
         if ax is None: ax = plt.subplot(111)
         jet = cm = plt.get_cmap('jet')
         cNorm  = colors.Normalize(
@@ -2123,6 +2138,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         assert vType in ['CC','F','E']
         assert self.dim == 3
 
+
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        import matplotlib.colors as colors
+        import matplotlib.cm as cmx
+
         szSliceDim = len(getattr(self, 'h'+normal.lower())) #: Size of the sliced dimension
         if ind is None: ind = int(szSliceDim/2)
         assert type(ind) in [int, long], 'ind must be an integer'
@@ -2202,6 +2224,25 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         if showIt: plt.show()
         return tuple(out)
 
+    def __len__(self): return self.nC
+
+    def __getitem__(self, key):
+        if isinstance( key, slice ) :
+            #Get the start, stop, and step from the slice
+            return [self[ii] for ii in xrange(*key.indices(len(self)))]
+        elif isinstance( key, int ) :
+            if key < 0 : #Handle negative indices
+                key += len( self )
+            if key >= len( self ) :
+                raise IndexError, "The index (%d) is out of range."%key
+
+            self._numberCells() # no-op if numbered
+            index   = self._i2cc[key]
+            pointer = self._asPointer(index)
+            return Cell(self, index, pointer)
+        else:
+            raise TypeError, "Invalid argument type."
+
 
 class Cell(object):
     def __init__(self, mesh, index, pointer):
@@ -2209,6 +2250,35 @@ class Cell(object):
         self._index   = index
         self._pointer = pointer
 
+    @property
+    def nodes(self):
+        """The node index in _gridN (this may include hanging nodes)."""
+        M = self.mesh
+        M._numberNodes()
+        p = self._pointer
+        i = self._index
+        w = M._levelWidth(p[-1])
+
+        if M.dim == 2:
+            n = [
+                    i,
+                    M._index([ p[0] + w, p[1]   , p[2]]),
+                    M._index([ p[0]    , p[1]+ w, p[2]]),
+                    M._index([ p[0] + w, p[1]+ w, p[2]]),
+                ]
+        elif self.dim == 3:
+            n = [
+                    i,
+                    M._index([ p[0] + w, p[1]    , p[2]    ,p[3]]),
+                    M._index([ p[0]    , p[1] + w, p[2]    ,p[3]]),
+                    M._index([ p[0] + w, p[1] + w, p[2]    ,p[3]]),
+                    M._index([ p[0]    , p[1]    , p[2] + w,p[3]]),
+                    M._index([ p[0] + w, p[1]    , p[2] + w,p[3]]),
+                    M._index([ p[0]    , p[1] + w, p[2] + w,p[3]]),
+                    M._index([ p[0] + w, p[1] + w, p[2] + w,p[3]]),
+                ]
+        return [M._n2i[_] for _ in n]
+
     @property
     def center(self):
         if getattr(self, '_center', None) is None:
@@ -2266,114 +2336,3 @@ class NotBalancedException(TreeException):
     pass
 class CellLookUpException(TreeException):
     pass
-
-if __name__ == '__main__':
-
-    def topo(x):
-        return np.sin(x*(2.*np.pi))*0.3 + 0.5
-
-    def function(cell):
-        r = cell.center - np.array([0.5]*len(cell.center))
-        dist = np.sqrt(r.dot(r))
-        # dist2 = np.abs(cell.center[-1] - topo(cell.center[0]))
-
-        # dist = min([dist1,dist2])
-        # if dist < 0.05:
-        #     return 5
-        if dist < 0.1:
-            return 5
-        if dist < 0.2:
-            return 4
-        if dist < 0.4:
-            return 3
-        return 2
-
-    # T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7)
-    # T = TreeMesh([128,128,128])
-    # T = TreeMesh([64,64],levels=6)
-    T = TreeMesh([4,4,4])
-    # T = TreeMesh([[(1,128)],[(1,128)]],levels=7)
-    # T.refine(lambda xc:2, balance=False)
-    # T._index([0,0,0])
-    # T._pointer(0)
-
-
-    # tic = time.time()
-    T.refine(function)#, balance=False)
-    # print time.time() - tic
-    # print T.nC
-    T.plotSlice(np.log(T.vol))#np.random.rand(T.nC))
-
-    plt.show()
-    blah
-
-    # T.plotImage(np.arange(len(T.vol)),showIt=True)
-
-    # print T.getFaceInnerProduct()
-    # print T.gridFz
-
-
-    # T._refineCell([8,0,1])
-    # T._refineCell([8,0,2])
-    # T._refineCell([12,0,2])
-    # T._refineCell([8,4,2])
-    # T._refineCell([6,0,3])
-    # T._refineCell([8,8,1])
-    # T._refineCell([0,0,0,1])
-    # T.__dirty__ = True
-
-
-    # print T.gridFx.shape[0], T.nFx
-
-
-
-    ax = plt.subplot(211)
-    ax.spy(T.edgeCurl)
-
-    # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense()
-    # print T.edgeCurl.todense()
-    # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - T.edgeCurl.todense()
-    # print T.gridEy - Mesh.TensorMesh([2,2,2]).gridEy
-
-    # print T.edge
-    # T.plotGrid(ax=ax)
-
-    # R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i)
-    # print R
-
-    ax = plt.subplot(212)#, projection='3d')
-    ax.spy(Mesh.TensorMesh([2,2,2]).edgeCurl)
-
-    # ax = plt.subplot(313)
-    # ax.spy(T.faceDiv[:,:T.nFx] * R)
-
-
-    # T.balance()
-    # T.plotGrid(ax=ax)
-
-    # cx = T._getNextCell([0,0,1],direction=0,positive=True)
-    # print cx
-    # # print [T._asPointer(_) for _ in cx]
-    # cx = T._getNextCell([8,0,3],direction=0,positive=False)
-    # print T._asPointer(cx)
-    # cx = T._getNextCell([8,8,1],direction=1,positive=False)
-    # print cx, #[T._asPointer(_) for _ in cx]
-    # cm = T._getNextCell([64,80,4],direction=0,positive=False)
-    # cy = T._getNextCell([64,80,4],direction=1,positive=True)
-    # cp = T._getNextCell([64,80,4],direction=1,positive=False)
-
-    # ax.plot( T._cellN([4,0,1])[0],T._cellN([4,0,1])[1], 'yd')
-    # ax.plot( T._cellN(cx)[0],T._cellN(cx)[1], 'ys')
-    # ax.plot( T._cellN(cm)[0],T._cellN(cm)[1], 'ys')
-    # ax.plot( T._cellN(cy)[0],T._cellN(cy)[1], 'ys')
-    # ax.plot( T._cellN(cp[0])[0],T._cellN(cp[0])[1], 'ys')
-    # ax.plot( T._cellN(cp[1])[0],T._cellN(cp[1])[1], 'ys')
-
-
-
-
-
-    # print T.nN
-
-    plt.show()
-

SimPEG/Mesh/View.py

@@ -1,8 +1,11 @@
 import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib
-from mpl_toolkits.mplot3d import Axes3D
-from SimPEG.Utils import mkvc, animate
+from SimPEG.Utils import mkvc
+try:
+    import matplotlib.pyplot as plt
+    import matplotlib
+    from mpl_toolkits.mplot3d import Axes3D
+except ImportError, e:
+    print 'Trouble importing matplotlib.'
 
 
 class TensorView(object):
@@ -479,63 +482,6 @@ class TensorView(object):
         ax.grid(True)
         if showIt: plt.show()
 
-    def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None):
-        assert normal in 'xyz', 'normal must be x, y, or z'
-        if ax is None: ax = plt.subplot(111)
-        I = mesh.r(var,'CC','CC','M')
-        axes = [p for p in 'xyz' if p not in normal.lower()]
-        if normal is 'x': I = I[index,:,:]
-        if normal is 'y': I = I[:,index,:]
-        if normal is 'z': I = I[:,:,index]
-        if clim is None: clim = [I.min(),I.max()]
-        p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1])
-        ax.axis('tight')
-        ax.set_xlabel(axes[0])
-        ax.set_ylabel(axes[1])
-        return p
-
-    def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)):
-        assert mesh.dim > 2, 'This is for 3D meshes only.'
-        # First set up the figure, the axis, and the plot element we want to animate
-        fig = plt.figure(figsize=figsize)
-        ax = plt.axes()
-        clim = [var.min(),var.max()]
-        plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim))
-        tlt = plt.title(normal)
-
-        def animateFrame(i):
-            mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
-            tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
-
-        return animate(fig, animateFrame, frames=mesh.vnC['xyz'.index(normal)])
-
-    def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1):
-        """
-        Call a function for a list of models to create a video.
-
-        ::
-
-            def function(var, ax, clim, tlt, i):
-                tlt.set_text('%d'%i)
-                return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim)
-
-            mesh.video([model1, model2, ..., modeln],function)
-        """
-        # First set up the figure, the axis, and the plot element we want to animate
-        fig = plt.figure(figsize=figsize)
-        ax = plt.axes()
-        VAR = np.concatenate(var)
-        clim = [VAR.min(),VAR.max()]
-        tlt = plt.title('')
-        if colorbar:
-            plt.colorbar(function(var[0],ax,clim,tlt,0))
-
-        frames = np.arange(0,len(var),skip)
-        def animateFrame(j):
-            i = frames[j]
-            function(var[i],ax,clim,tlt,i)
-
-        return animate(fig, animateFrame, frames=len(frames))
 
 class CylView(object):
 

SimPEG/Optimization.py

@@ -989,5 +989,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                 if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG):
                     cgFlag = 1
                 # End CG Iterations
+                    
+            # Take a gradient step on the active cells if exist
+            if temp != self.xc.size:
+        
+                rhs_a = (Active) * -self.g
+                
+                dm_i = max( abs( delx ) )
+                dm_a = max( abs(rhs_a) )           
+    
+                delx = delx + rhs_a * dm_i / dm_a /10.
+                    
+            # Only keep gradients going in the right direction on the active set
+            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
+            delx[indx] = 0.
 
-            return delx
+            return delx

SimPEG/Problem.py

@@ -1,6 +1,6 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
-import Maps, Mesh
+import Maps, Mesh, Exceptions
 from Fields import Fields, TimeFields
 
 class BaseProblem(object):
@@ -18,10 +18,14 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    mesh    = None    #: A SimPEG.Mesh instance.
-
     PropMap = None    #: A SimPEG PropertyMap class.
 
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        self.mesh = mesh
+        self.mapping = mapping or Maps.IdentityMap(mesh)
+
     @property
     def mapping(self):
         "A SimPEG.Map instance or a property map is PropMap is not None"
@@ -32,13 +36,8 @@ class BaseProblem(object):
             val._assertMatchesPair(self.mapPair)
             self._mapping = val
         else:
-            self._mapping = self.PropMap(val) 
-    
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Utils.setKwargs(self, **kwargs)
-        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        self.mesh = mesh
-        self.mapping = mapping or Maps.IdentityMap(mesh)
+            self._propMapMapping = val
+            self._mapping = self.PropMap(val)
 
     @property
     def survey(self):
@@ -47,13 +46,22 @@ class BaseProblem(object):
         """
         return getattr(self, '_survey', None)
 
-    def pair(self, d):
+    def pair(self, survey):
         """Bind a survey to this problem instance using pointers."""
-        assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
-        if d.ispaired:
+        assert isinstance(survey, self.surveyPair), "Survey must be an instance of a %s class."%(self.surveyPair.__name__)
+        if survey.ispaired:
             raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
-        self._survey = d
-        d._prob = self
+        try:
+            self._survey = survey
+            self._validatePairing()
+        except Exceptions.PairingException, e:
+            self._survey = None
+            raise e
+        survey._prob = self
+
+    def _validatePairing(self):
+        """Called when the pair is done, raise a SimPEG.Exceptions.PairingException if unsuccessful"""
+        pass
 
     def unpair(self):
         """Unbind a survey from this problem instance."""
@@ -158,6 +166,9 @@ class BaseProblem(object):
 
 class BaseTimeProblem(BaseProblem):
     """Sets up that basic needs of a time domain problem."""
+    
+    waveformType = "STEPOFF"    
+    current = None
 
     @property
     def timeSteps(self):
@@ -184,6 +195,11 @@ class BaseTimeProblem(BaseProblem):
         self._timeSteps = Utils.meshTensor(value)
         del self.timeMesh
 
+    def currentwaveform(self, wave):
+        self._timeSteps = np.diff(wave[:,0])
+        self.current = wave[:,1]
+        self.waveformType = "GENERAL"
+
     @property
     def nT(self):
         "Number of time steps."
@@ -214,4 +230,191 @@ class BaseTimeProblem(BaseProblem):
             del self._timeMesh
 
 
+class GlobalProblem(BaseProblem):
+    """
+
+        The GlobalProblem allows you to run a whole bunch of SubProblems,
+        potentially in parallel, potentially of different meshes.
+
+        This is handy for working with lots of sources,
+
+    """
+
+    surveyKwargs = {}
+    probKwargs   = {}
+
+    def __init__(self, SubProblem, globalMesh, mapping=None, **kwargs):
+
+        # assert isclass??(SubProblem, BaseProblem), "SubProblem must be a SimPEG.Problem.BaseProblem object."
+        self.surveyPair = SubProblem.surveyPair
+        self.PropMap = SubProblem.PropMap
+        self.mapPair = SubProblem.mapPair
+        self.SubProblem = SubProblem
+
+        Utils.setKwargs(self, **kwargs)
+        assert isinstance(globalMesh, Mesh.BaseMesh), "globalMesh must be a SimPEG.Mesh object."
+        self.globalMesh = globalMesh
+        self.mapping = mapping or Maps.IdentityMap()
+
+    @property
+    def groups(self):
+        """
+            List of lists/integers to say how the sources are grouped.
+
+            e.g.
+
+            survey.srcList = [s0,s1,s2,s3,s4]
+            groups = [ [0,4], [1,3], 2 ]
+        """
+        if getattr(self, '_groups', None) is None:
+            if not self.ispaired: return None
+            self._groups = range(self.survey.nSrc)
+        return self._groups
+    @groups.setter
+    def groups(self, val):
+        assert type(val) is list, 'This should be an list of groups'
+        if self.ispaired:
+            for g in val:
+                assert type(g) in [int, list], 'Must be an integer or a list'
+                if type(g) is int:
+                    assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
+                if type(g) is list:
+                    for sg in g:
+                        assert type(g) is int, 'Must be an integer or a list'
+                        assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
+            assert len(val) == len(self.survey.srcList), 'The groups must be the same length as the srcList in the survey'
+        self._groups = val
+        self._nGroups = None
+
+    @property
+    def meshes(self):
+        if getattr(self, '_meshes', None) is None:
+            if not self.ispaired: return None
+            self._meshes = [self.globalMesh]*self.nGroups
+        return self._meshes
+    @meshes.setter
+    def meshes(self, val):
+        assert type(val) is list
+        if self.ispaired:
+            assert len(val) == self.nGroups
+        self._meshes = val
+
+    @property
+    def nGroups(self):
+        if getattr(self, '_groups', None) is None:
+            return None
+        return len(self.groups)
+
+    def _validatePairing(self):
+        try:
+            self.groups = self.groups # check the assumptions for the grouping
+        except Exception, e:
+            raise Exceptions.PairingException(reason='The grouping does not match the survey')
+        if self.nGroups is not len(self.meshes):
+            raise Exceptions.PairingException(reason='The meshes are not the the same length as the number of groups')
+
+    def getSubProblemandSubSurvey(self, subMap, ind):
+        #This is a core place that we can proceed parallelization
+        assert self.ispaired, 'You must be paired to a survey'
+        assert type(ind) in [int,long] and ind >= 0 and ind < self.nGroups, 'ind must be an index into the group list'
+
+        subMesh = self.meshes[ind]
+        # subMap  = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
+        # subMap  = self.getSubMap(subMesh, ind)
+
+        if self.PropMap is None:
+            prob = self.SubProblem(subMesh, mapping=subMap * self.mapping, **self.probKwargs)
+        else:
+            # This will not work with a fancier propmap...
+            prob = self.SubProblem(subMesh, mapping=subMap * self._propMapMapping, **self.probKwargs)
+
+        survey = self.survey.__class__(srcList=self.survey.srcList[self.groups[ind]], **self.surveyKwargs)
+        prob.pair(survey)
+
+        return prob, survey
+
+    # Not sure we need this here ... 
+
+    def getSubMap(self, subMesh, ind):
+        """The sub"""
+        mesh2mesh = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
+
+        if self.PropMap is None:
+            subMap = mesh2mesh * self.mapping
+        else:
+            subMap = mesh2mesh * self._propMapMapping
+
+        return subMap
+
+if __name__ == '__main__':
+
+
+    from SimPEG import *
+    from SimPEG import EM
+    from scipy.constants import mu_0
+    from pymatsolver import MumpsSolver
+
+    cs = 10.
+    ncx, ncy, ncz = 10, 10, 10
+    npad = 4
+    freq = 1e2
+
+    hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
+
+    mapping = Maps.ExpMap(mesh)
+
+    x = np.linspace(-10,10,5)
+    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
+    rxList = EM.FDEM.Rx(XYZ, 'exi')
+    Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
+    Src1 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
+
+
+    prb0 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
+    survey = EM.FDEM.Survey([Src0])
+    prb0.pair(survey)
+    prb1 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
+    survey = EM.FDEM.Survey([Src1])
+    prb1.pair(survey)
+
+
+
+    sig = 1e-1
+    sigma = np.ones(mesh.nC)*sig
+    sigma[mesh.gridCC[:,2] > 0] = 1e-8
+    m = np.log(sigma)
+
+    GP = GlobalProblem(EM.FDEM.Problem_b, mesh, mapping=mapping, meshes=[mesh,mesh])
+    survey = EM.FDEM.Survey([Src0, Src1])
+    GP.pair(survey)
+
+    gp1 = GP.getSubProblem(0)
+    gp1.Solver = MumpsSolver
+
+    pu = prb0.fields(m)
+    gpu = gp1.fields(m)
+
+    bfz = mesh.r(pu[Src0, 'b'],'F','Fz','M')
+    bfz = mesh.r(gpu[Src0, 'b'],'F','Fz','M')
+    x = np.linspace(-55,55,12)
+    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
+    P = mesh.getInterpolationMat(XYZ, 'Fz')
+
+    # an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, Src0.freq, sig)
+
+    # diff = np.log10(np.abs(P*np.imag(pu[Src0, 'b']) - mu_0*np.imag(an)))
+    # diff = np.log10(np.abs(P*np.imag(gpu[Src0, 'b']) - mu_0*np.imag(an)))
+
+    import matplotlib.pyplot as plt
+    plt.plot(x,np.log10(np.abs(P*np.imag(pu[Src0, 'b']))), 'r-s')
+    plt.plot(x,np.log10(np.abs(P*np.imag(gpu[Src0, 'b']))), 'b')
+    # plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
+    # plt.plot(x,diff,'g')
+    plt.show()
+
+
+
 

SimPEG/Regularization.py

@@ -20,12 +20,13 @@ class BaseRegularization(object):
     mesh    = None    #: A SimPEG.Mesh instance.
     mref = None       #: Reference model.
 
-    def __init__(self, mesh, mapping=None, **kwargs):
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        self.mapping = mapping or Maps.IdentityMap(mesh)
+        self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
+        self.indActive = indActive
 
     @property
     def parent(self):
@@ -112,89 +113,8 @@ class BaseRegularization(object):
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
 
-
-
 class Tikhonov(BaseRegularization):
-    """**Tikhonov Regularization**
-
-        Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same:
-
-        .. math::
-
-            R(m) = \int_\Omega \\frac{\\alpha_x}{2}\left(\\frac{\partial m}{\partial x}\\right)^2 + \\frac{\\alpha_y}{2}\left(\\frac{\partial m}{\partial y}\\right)^2 \partial v
-
-        Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this:
-
-        .. math::
-
-            R(m) \\approx \sum_{ij} \left[\\frac{\\alpha_x}{2}\left[\left(\\frac{m_{i+1,j} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i-1,j}}{h}\\right)^2\\right]
-            + \\frac{\\alpha_y}{2}\left[\left(\\frac{m_{i,j+1} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i,j-1}}{h}\\right)^2\\right]
-            \\right]h^2
-
-        If we let D_1 be the derivative matrix in the x direction
-
-        .. math::
-
-            \mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1
-
-        .. math::
-
-            \mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1
-
-        Where d_1 is the one dimensional derivative:
-
-        .. math::
-
-            \mathbf{d}_1 = \\frac{1}{h} \left[ \\begin{array}{cccc}
-            -1 & 1 & & \\\\
-             & \ddots & \ddots&\\\\
-             &  & -1 & 1\end{array} \\right]
-
-        .. math::
-
-            R(m) \\approx \mathbf{v}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \\frac{\\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \\right]
-
-        Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise)
-
-        .. math::
-
-            \mathbf{a\odot b} = \\text{diag}(\mathbf{a})\mathbf{b} = \\text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a}
-
-        and the transpose also is true (but the sizes have to make sense...):
-
-        .. math::
-
-            \mathbf{a}^\\top\\text{diag}(\mathbf{b}) = \mathbf{b}^\\top\\text{diag}(\mathbf{a})
-
-        So R(m) can simplify to:
-
-        .. math::
-
-            R(m) \\approx  \mathbf{m}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{D}_1^\\top \\text{diag}(\mathbf{A}_1^\\top\mathbf{v}) \mathbf{D}_1 +  \\frac{\\alpha_y}{2}\mathbf{D}_2^\\top \\text{diag}(\mathbf{A}_2^\\top \mathbf{v}) \mathbf{D}_2 \\right] \mathbf{m}
-
-        We will define W_x as:
-
-        .. math::
-
-            \mathbf{W}_x = \sqrt{\\alpha_x}\\text{diag}\left(\sqrt{\mathbf{A}_1^\\top\mathbf{v}}\\right) \mathbf{D}_1
-
-
-        And then W as a tall matrix of all of the different regularization terms:
-
-        .. math::
-
-            \mathbf{W} = \left[ \\begin{array}{c}
-            \mathbf{W}_s\\\\
-            \mathbf{W}_x\\\\
-            \mathbf{W}_y\end{array} \\right]
-
-        Then we can write
-
-        .. math::
-
-            R(m) \\approx \\frac{1}{2}\mathbf{m^\\top W^\\top W m}
-
-
+    """
     """
     smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
     alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
@@ -205,9 +125,182 @@ class Tikhonov(BaseRegularization):
     alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
 
+    def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
+        self.indActive = indActive
+
+    @property
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self,'_Ws', None) is None:
+            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Ws = Pac.T * self._Ws * Pac
+        return self._Ws
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
+            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
+
+            if self.indActive is not None:
+                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
+                self._Wx = Pafx.T*self._Wx*Pac
+
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
+            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
+
+            if self.indActive is not None:
+                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
+                self._Wy = Pafy.T*self._Wy*Pac
+
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
+            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
+
+            if self.indActive is not None:
+                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
+                self._Wz = Pafz.T*self._Wz*Pac
+
+        return self._Wz
+
+    @property
+    def Wxx(self):
+        """Regularization matrix Wxx"""
+        if getattr(self, '_Wxx', None) is None:
+            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
+
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Wxx = Pac.T*self._Wxx*Pac
+
+        return self._Wxx
+
+    @property
+    def Wyy(self):
+        """Regularization matrix Wyy"""
+        if getattr(self, '_Wyy', None) is None:
+            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
+
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Wyy = Pac.T*self._Wyy*Pac
+
+        return self._Wyy
+
+    @property
+    def Wzz(self):
+        """Regularization matrix Wzz"""
+        if getattr(self, '_Wzz', None) is None:
+            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
+
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Wzz = Pac.T*self._Wzz*Pac
+
+        return self._Wzz
+
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx, self.Wxx)
+            if self.mesh.dim > 1:
+                wlist += (self.Wy, self.Wyy)
+            if self.mesh.dim > 2:
+                wlist += (self.Wz, self.Wzz)
+            self._Wsmooth = sp.vstack(wlist)
+        return self._Wsmooth
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Ws, self.Wsmooth)
+            self._W = sp.vstack(wlist)
+        return self._W
+
+    @Utils.timeIt
+    def eval(self, m):
+        if self.smoothModel == True:
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
+        elif self.smoothModel == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        if self.smoothModel == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1+out2
+        elif self.smoothModel == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
+
+class Simple(BaseRegularization):
+    """
+        Only for tensor mesh
+    """     
+    
+    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s  = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
+    alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+    alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+    alpha_z  = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+    alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
+    alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
+    alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
+
     def __init__(self, mesh, mapping=None, **kwargs):
         BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
 
+    
+        
     @property
     def Ws(self):
         """Regularization matrix Ws"""
@@ -219,24 +312,21 @@ class Tikhonov(BaseRegularization):
     def Wx(self):
         """Regularization matrix Wx"""
         if getattr(self, '_Wx', None) is None:
-            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
-            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
+            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
         if getattr(self, '_Wy', None) is None:
-            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
-            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
+            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
         if getattr(self, '_Wz', None) is None:
-            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
-            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
+            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
         return self._Wz
 
     @property
@@ -311,7 +401,7 @@ class Tikhonov(BaseRegularization):
         if self.smoothModel == True:
             mD1 = self.mapping.deriv(m)
             mD2 = self.mapping.deriv(m - self.mref)
-            r1 = self.Wsmooth * ( self.mapping * (m)) 
+            r1 = self.Wsmooth * ( self.mapping * (m))
             r2 = self.Ws * ( self.mapping * (m - self.mref) )
             out1 = mD1.T * ( self.Wsmooth.T * r1 )
             out2 = mD2.T * ( self.Ws.T * r2 )
@@ -322,3 +412,109 @@ class Tikhonov(BaseRegularization):
             out = mD.T * ( self.W.T * r )
         return out
 
+class SparseRegularization(Simple):
+    
+    eps = 1e-1
+
+    m = None
+    gamma = 1.
+    p = 0.
+    qx = 2.
+    qy = 2.
+    qz = 2.
+    
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Simple.__init__(self, mesh, mapping=mapping, **kwargs)
+        
+    
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx, self.Wxx)
+            if self.mesh.dim > 1:
+                wlist += (self.Wy, self.Wyy)
+            if self.mesh.dim > 2:
+                wlist += (self.Wz, self.Wzz)
+            self._Wsmooth = sp.vstack(wlist)
+        return self._Wsmooth
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Ws, self.Wsmooth)
+            self._W = sp.vstack(wlist)
+        return self._W
+    
+    @property
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self, 'm', None) is None:
+            self.Rs = Utils.speye(self.mesh.nC)            
+            
+        else:
+            f_m = self.m
+            self.rs = self.R(f_m , self.p, self.eps)
+            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
+            self.Rs = Utils.sdiag( self.rs )
+            
+        self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
+            
+        return self._Ws
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        
+        if getattr(self, 'm', None) is None:
+            self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])            
+            
+        else:
+            f_m = self.mesh.unitCellGradx * self.m
+            self.rx = self.R( f_m , self.qx, self.eps)
+            self.Rx = Utils.sdiag( self.rx )
+            
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        
+        if getattr(self, 'm', None) is None:
+            self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])            
+            
+        else:
+            f_m = self.mesh.unitCellGrady * self.m
+            self.ry = self.R( f_m , self.qy, self.eps)
+            self.Ry = Utils.sdiag( self.ry )
+            
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+
+        if getattr(self, 'm', None) is None:
+            self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])            
+            
+        else:
+            f_m = self.mesh.unitCellGradz * self.m
+            self.rz = self.R( f_m , self.qz, self.eps)
+            self.Rz = Utils.sdiag( self.rz )
+        
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
+        return self._Wz    
+
+    
+    def R(self, f_m , p, dec):
+
+        eta = (self.eps**(1-p/2.))**0.5
+        r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
+
+        return r

SimPEG/Survey.py

@@ -1,6 +1,5 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
 
-
 class BaseRx(object):
     """SimPEG Receiver Object"""
 
@@ -223,6 +222,8 @@ class BaseSurvey(object):
 
     @srcList.setter
     def srcList(self, value):
+        if isinstance(value, self.srcPair):
+            value = [value]
         assert type(value) is list, 'srcList must be a list'
         assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
         assert len(set(value)) == len(value), 'The srcList must be unique'
@@ -374,3 +375,11 @@ class BaseSurvey(object):
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
         return self.dobs
+
+class LinearSurvey(BaseSurvey):
+    def projectFields(self, u):
+        return u
+    
+    @property
+    def nD(self):
+        return self.prob.G.shape[0]

SimPEG/Tests.py

@@ -1,5 +1,4 @@
 import numpy as np
-import matplotlib.pyplot as plt
 from numpy.linalg import norm
 from SimPEG.Utils import mkvc, sdiag, diagEst
 from SimPEG import Utils
@@ -311,6 +310,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
 
 
     if plotIt:
+        import matplotlib.pyplot as plt
         ax = ax or plt.subplot(111)
         ax.loglog(h, E0, 'b')
         ax.loglog(h, E1, 'g--')

SimPEG/Utils/ModelBuilder.py

@@ -118,6 +118,44 @@ def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.
     D = np.sqrt(np.sum(G**2,axis=1))
     return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5
 
+def getIndicesSphere(center,radius,ccMesh):
+    """
+        Creates a vector containing the sphere indices in the cell centers mesh.
+        Returns a tuple
+
+        The sphere is defined by the points
+
+        p0, describe the position of the center of the cell
+
+        r, describe the radius of the sphere.
+
+        ccMesh represents the cell-centered mesh
+
+        The points p0 must live in the the same dimensional space as the mesh.
+
+    """
+
+    # Validation: mesh and point (p0) live in the same dimensional space
+    dimMesh = np.size(ccMesh[0,:])
+    assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
+
+    if dimMesh == 1:
+       # Define the reference points
+        
+        ind  = np.abs(center[0] - ccMesh[:,0]) < radius
+
+    elif dimMesh == 2:
+       # Define the reference points
+
+        ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius
+
+    elif dimMesh == 3:
+        # Define the points
+        ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius
+
+    # Return a tuple
+    return ind
+
 def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     """
     Define a two layered model.  Depth of the first layer must be specified.

SimPEG/Utils/SolverUtils.py

@@ -26,7 +26,14 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
 
     def __init__(self, A, **kwargs):
         self.A = A.tocsc()
+
+        self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy)
+        if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"]
+        self.accuracyTol = kwargs.get("accuracyTol", accuracyTol)
+        if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"]
+
         self.kwargs = kwargs
+
         if factorize:
             self.solver = fun(self.A, **kwargs)
 
@@ -57,8 +64,8 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
                 else:
                     X[:,i] = fun(self.A, b[:,i], **self.kwargs)
 
-        if checkAccuracy:
-            _checkAccuracy(self.A, b, X, accuracyTol)
+        if self.checkAccuracy:
+            _checkAccuracy(self.A, b, X, self.accuracyTol)
         return X
 
     def clean(self):
@@ -81,6 +88,12 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
 
     def __init__(self, A, **kwargs):
         self.A = A
+
+        self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy)
+        if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"]
+        self.accuracyTol = kwargs.get("accuracyTol", accuracyTol)
+        if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"]
+
         self.kwargs = kwargs
 
     def __mul__(self, b):
@@ -108,8 +121,8 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
                 else:
                     X[:,i] = out
 
-        if checkAccuracy:
-            _checkAccuracy(self.A, b, X, accuracyTol)
+        if self.checkAccuracy:
+            _checkAccuracy(self.A, b, X, self.accuracyTol)
         return X
 
     def clean(self):

SimPEG/Utils/__init__.py

@@ -1,9 +1,8 @@
 from matutils import *
 from codeutils import *
-from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
+from meshutils import *
 from curvutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
-from ipythonutils import easyAnimate as animate
 from CounterUtils import *
 import ModelBuilder
 import SolverUtils

SimPEG/Utils/codeutils.py

@@ -17,7 +17,7 @@ def memProfileWrapper(towrap, *funNames):
 
         For example::
 
-            foo_mem = memProfile(foo,'my_func')
+            foo_mem = memProfileWrapper(foo,['my_func'])
             fooi = foo_mem()
             for i in range(5):
                 fooi.my_func()

SimPEG/Utils/ipythonutils.py

@@ -1,28 +0,0 @@
-from tempfile import NamedTemporaryFile
-import matplotlib.pyplot as plt
-from matplotlib import animation
-
-# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
-# http://www.renevolution.com/how-to-install-ffmpeg-on-mac-os-x/
-
-VIDEO_TAG = """<video controls loop>
- <source src="data:video/x-m4v;base64,{0}" type="video/mp4">
- Your browser does not support the video tag.
-</video>"""
-
-def anim_to_html(anim):
-    if not hasattr(anim, '_encoded_video'):
-        with NamedTemporaryFile(suffix='.mp4') as f:
-            anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
-            video = open(f.name, "rb").read()
-        anim._encoded_video = video.encode("base64")
-
-    return VIDEO_TAG.format(anim._encoded_video)
-
-def display_animation(anim):
-    plt.close(anim._fig)
-    return anim_to_html(anim)
-
-animation.Animation._repr_html_ = display_animation
-
-easyAnimate = animation.FuncAnimation

SimPEG/Utils/meshutils.py

@@ -102,223 +102,6 @@ def closestPoints(mesh, pts, gridLoc='CC'):
 
     return nodeInds
 
-def readUBCTensorMesh(fileName):
-    """
-        Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
-
-        Input:
-        :param fileName, path to the UBC GIF mesh file
-
-        Output:
-        :param SimPEG TensorMesh object
-        :return
-    """
-
-    # Interal function to read cell size lines for the UBC mesh files.
-    def readCellLine(line):
-        for seg in line.split():
-            if '*' in seg:
-                st = seg
-                sp = seg.split('*')
-                re = np.array(sp[0],dtype=int)*(' ' + sp[1])
-                line = line.replace(st,re.strip())
-        return np.array(line.split(),dtype=float)
-
-    # Read the file as line strings, remove lines with comment = !
-    msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
-
-    # Fist line is the size of the model
-    sizeM = np.array(msh[0].split(),dtype=float)
-    # Second line is the South-West-Top corner coordinates.
-    x0 = np.array(msh[1].split(),dtype=float)
-    # Read the cell sizes
-    h1 = readCellLine(msh[2])
-    h2 = readCellLine(msh[3])
-    h3temp = readCellLine(msh[4])
-    h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom.
-    # Adjust the reference point to the bottom south west corner
-    x0[2] = x0[2] - np.sum(h3)
-    # Make the mesh
-    from SimPEG import Mesh
-    tensMsh = Mesh.TensorMesh([h1,h2,h3],x0)
-    return tensMsh
-
-def readUBCTensorModel(fileName, mesh):
-    """
-        Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
-
-        Input:
-        :param fileName, path to the UBC GIF mesh file to read
-        :param mesh, TensorMesh object, mesh that coresponds to the model
-
-        Output:
-        :return numpy array, model with TensorMesh ordered
-    """
-    f = open(fileName, 'r')
-    model = np.array(map(float, f.readlines()))
-    f.close()
-    model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F')
-    model = model[::-1,:,:]
-    model = np.transpose(model, (1, 2, 0))
-    model = mkvc(model)
-
-    return model
-
-def writeUBCTensorMesh(fileName, mesh):
-    """
-        Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
-
-        :param str fileName: File to write to
-        :param simpeg.Mesh.TensorMesh mesh: The mesh
-
-    """
-    assert mesh.dim == 3
-    s = ''
-    s += '%i %i %i\n' %tuple(mesh.vnC)
-    origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
-    origin.dtype = float
-
-    s += '%.2f %.2f %.2f\n' %tuple(origin)
-    s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
-    s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
-    s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
-    f = open(fileName, 'w')
-    f.write(s)
-    f.close()
-
-def writeUBCTensorModel(fileName, mesh, model):
-    """
-        Writes a model associated with a SimPEG TensorMesh
-        to a UBC-GIF format model file.
-
-        :param str fileName: File to write to
-        :param simpeg.Mesh.TensorMesh mesh: The mesh
-        :param numpy.ndarray model: The model
-    """
-
-    # Reshape model to a matrix
-    modelMat = mesh.r(model,'CC','CC','M')
-    # Transpose the axes
-    modelMatT = modelMat.transpose((2,0,1))
-    # Flip z to positive down
-    modelMatTR = mkvc(modelMatT[::-1,:,:])
-
-    np.savetxt(fileName, modelMatTR.ravel())
-
-
-def readVTRFile(fileName):
-    """
-        Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
-
-        Input:
-        :param vtrFileName, path to the vtr model file to write to
-
-        Output:
-        :return SimPEG TensorMesh object
-        :return SimPEG model dictionary
-
-    """
-    # Import
-    from vtk import vtkXMLRectilinearGridReader as vtrFileReader
-    from vtk.util.numpy_support import vtk_to_numpy
-
-    # Read the file
-    vtrReader = vtrFileReader()
-    vtrReader.SetFileName(fileName)
-    vtrReader.Update()
-    vtrGrid = vtrReader.GetOutput()
-    # Sort information
-    hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates())))
-    xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0]
-    hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates())))
-    yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0]
-    zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates()))
-    # Check the direction of hz
-    if np.all(zD < 0):
-        hz = np.abs(zD[::-1])
-        zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1]
-    else:
-        hz = np.abs(zD)
-        zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0]
-    x0 = np.array([xR,yR,zR])
-
-    # Make the SimPEG object
-    from SimPEG import Mesh
-    tensMsh = Mesh.TensorMesh([hx,hy,hz],x0)
-
-    # Grap the models
-    modelDict = {}
-    for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()):
-        modelName = vtrGrid.GetCellData().GetArrayName(i)
-        if np.all(zD < 0):
-            modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
-            tM = tensMsh.r(modFlip,'CC','CC','M')
-            modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V')
-        else:
-            modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
-        modelDict[modelName] = modArr
-
-    # Return the data
-    return tensMsh, modelDict
-
-def writeVTRFile(fileName,mesh,model=None):
-    """
-    Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
-
-    Input:
-    :param str, path to the output vtk file
-    :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
-    :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells
-
-    """
-    # Import
-    from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter
-    from vtk.util.numpy_support import numpy_to_vtk
-
-    # Deal with dimensionalities
-    if mesh.dim >= 1:
-        vX = mesh.vectorNx
-        xD = mesh.nNx
-        yD,zD = 1,1
-        vY, vZ = np.array([0,0])
-    if mesh.dim >= 2:
-        vY = mesh.vectorNy
-        yD = mesh.nNy
-    if mesh.dim == 3:
-        vZ = mesh.vectorNz
-        zD = mesh.nNz
-    # Use rectilinear VTK grid.
-    # Assign the spatial information.
-    vtkObj = rectGrid()
-    vtkObj.SetDimensions(xD,yD,zD)
-    vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
-    vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
-    vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
-
-    # Assign the model('s) to the object
-    for item in model.iteritems():
-        # Convert numpy array
-        vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
-        vtkDoubleArr.SetName(item[0])
-        vtkObj.GetCellData().AddArray(vtkDoubleArr)
-    # Set the active scalar
-    vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
-    vtkObj.Update()
-
-
-    # Check the extension of the fileName
-    ext = os.path.splitext(fileName)[1]
-    if ext is '':
-        fileName = fileName + '.vtr'
-    elif ext not in '.vtr':
-        raise IOError('{:s} is an incorrect extension, has to be .vtr')
-    # Write the file.
-    vtrWriteFilter = rectWriter()
-    vtrWriteFilter.SetInput(vtkObj)
-    vtrWriteFilter.SetFileName(fileName)
-    vtrWriteFilter.Update()
-
-
 def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
     """
         Extracts Core Mesh from Global mesh

SimPEG/__init__.py

@@ -15,7 +15,7 @@ import Directives
 import Inversion
 import Tests
 
-__version__   = '0.1.4'
+__version__   = '0.1.9'
 __author__    = 'Rowan Cockett'
 __license__   = 'MIT'
 __copyright__ = 'Copyright 2014 Rowan Cockett'

docs/api_FiniteVolume.rst

@@ -0,0 +1,19 @@
+.. _api_FiniteVolume:
+
+Finite Volume
+*************
+
+Any numerical implementation requires the discretization of continuous functions into discrete approximations. These approximations are typically organized in a mesh, which defines boundaries, locations, and connectivity. Of specific interest to geophysical simulations, we require that averaging, interpolation and differential operators be defined for any mesh. In SimPEG, we have implemented a staggered mimetic finite volume approach (`Hyman and Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This approach requires the definitions of variables at either cell-centers, nodes, faces, or edges as seen in the figure below.
+
+.. image:: images/finitevolrealestate.png
+   :width: 400 px
+   :alt: FiniteVolume
+   :align: center
+
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Mesh
+   api_DiffOps
+   api_InnerProducts

docs/api_ForwardProblem.rst

@@ -61,11 +61,6 @@ If the forward problem is invertible, then we can rearrange for \\(\\frac{\\part
 This can often be computed given a vector (i.e. \\(J(v)\\)) rather than stored, as \\(J\\) is a large dense matrix.
 
 
-.. math::
-
-    u(m)
-
-
 
 The API
 =======
@@ -78,7 +73,7 @@ Problem
 
 Survey
 ------
-
 .. automodule:: SimPEG.Survey
     :members:
+    :undoc-members:
 

docs/api_Inversion.rst

@@ -1,21 +1,19 @@
-.. _api_Inverse:
 
 
-Regularization
-**************
+InvProblem
+**********
 
-.. automodule:: SimPEG.Regularization
+.. automodule:: SimPEG.InvProblem
     :show-inheritance:
     :members:
     :undoc-members:
 
 
-Optimize
-********
+Inversion
+*********
 
-.. automodule:: SimPEG.Optimization
+.. automodule:: SimPEG.Inversion
     :show-inheritance:
-    :private-members:
     :members:
     :undoc-members:
 
@@ -27,12 +25,3 @@ Directives
     :members:
     :undoc-members:
 
-Inversion
-*********
-
-.. automodule:: SimPEG.Inversion
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-

docs/api_InversionComponents.rst

@@ -0,0 +1,11 @@
+Inversion Components
+********************
+
+.. toctree::
+   :maxdepth: 3
+
+   api_DataMisfit
+   api_Regularization
+   api_Optimization
+   api_Inversion
+

docs/api_Mesh.rst

@@ -24,8 +24,7 @@ the implementations.
 .. plot::
 
     from SimPEG import Examples
-    Examples.Mesh_ThreeMeshes.run()
-
+    Examples.Mesh_Basic_Types.run()
 
 
 Variable Locations and Terminology

docs/api_MeshCode.rst

@@ -9,6 +9,15 @@ Tensor Mesh
     :undoc-members:
 
 
+Cylindrical Mesh
+================
+
+.. automodule:: SimPEG.Mesh.CylMesh
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+
 Tree Mesh
 =========
 
@@ -21,16 +30,7 @@ Tree Mesh
 Curvilinear Mesh
 ================
 
-.. automodule:: SimPEG.Mesh.Curvilinear
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Cylindrical Mesh
-================
-
-.. automodule:: SimPEG.Mesh.CylMesh
+.. automodule:: SimPEG.Mesh.CurvilinearMesh
     :show-inheritance:
     :members:
     :undoc-members:

docs/api_Optimization.rst

@@ -0,0 +1,9 @@
+
+Optimize
+********
+
+.. automodule:: SimPEG.Optimization
+    :show-inheritance:
+    :private-members:
+    :members:
+    :undoc-members:

docs/api_Regularization.rst

@@ -0,0 +1,100 @@
+
+Regularization
+**************
+
+If there is one model that has a misfit that equals the desired tolerance, then there are infinitely many other models which can fit to the same degree. The challenge is to find that model which has the desired characteristics and is compatible with a priori information. A single model can be selected from an infinite ensemble by measuring the length, or norm, of each model. Then a smallest, or sometimes largest, member can be isolated. Our goal is to design a norm that embodies our prior knowledge and, when minimized, yields a realistic candidate for the solution of our problem. The norm can penalize variation from a reference model, spatial derivatives of the model, or some combination of these.
+
+Tikhonov Regularization
+=======================
+
+Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same:
+
+.. math::
+
+    R(m) = \int_\Omega \frac{\alpha_x}{2}\left(\frac{\partial m}{\partial x}\right)^2 + \frac{\alpha_y}{2}\left(\frac{\partial m}{\partial y}\right)^2 \partial v
+
+Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this:
+
+.. math::
+
+    R(m) \approx \sum_{ij} \left[\frac{\alpha_x}{2}\left[\left(\frac{m_{i+1,j} - m_{i,j}}{h}\right)^2 + \left(\frac{m_{i,j} - m_{i-1,j}}{h}\right)^2\right] \\
+    + \frac{\alpha_y}{2}\left[\left(\frac{m_{i,j+1} - m_{i,j}}{h}\right)^2 + \left(\frac{m_{i,j} - m_{i,j-1}}{h}\right)^2\right]
+    \right]h^2
+
+If we let D_1 be the derivative matrix in the x direction
+
+.. math::
+
+    \mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1
+
+.. math::
+
+    \mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1
+
+Where d_1 is the one dimensional derivative:
+
+.. math::
+
+    \mathbf{d}_1 = \frac{1}{h} \left[ \begin{array}{cccc}
+    -1 & 1 & & \\
+     & \ddots & \ddots&\\
+     &  & -1 & 1\end{array} \right]
+
+.. math::
+
+    R(m) \approx \mathbf{v}^\top \left[\frac{\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \frac{\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \right]
+
+Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise)
+
+.. math::
+
+    \mathbf{a\odot b} = \text{diag}(\mathbf{a})\mathbf{b} = \text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a}
+
+and the transpose also is true (but the sizes have to make sense...):
+
+.. math::
+
+    \mathbf{a}^\top\text{diag}(\mathbf{b}) = \mathbf{b}^\top\text{diag}(\mathbf{a})
+
+So R(m) can simplify to:
+
+.. math::
+
+    R(m) \approx  \mathbf{m}^\top \left[\frac{\alpha_x}{2}\mathbf{D}_1^\top \text{diag}(\mathbf{A}_1^\top\mathbf{v}) \mathbf{D}_1 +  \frac{\alpha_y}{2}\mathbf{D}_2^\top \text{diag}(\mathbf{A}_2^\top \mathbf{v}) \mathbf{D}_2 \right] \mathbf{m}
+
+We will define W_x as:
+
+.. math::
+
+    \mathbf{W}_x = \sqrt{\alpha_x}\text{diag}\left(\sqrt{\mathbf{A}_1^\top\mathbf{v}}\right) \mathbf{D}_1
+
+
+And then W as a tall matrix of all of the different regularization terms:
+
+.. math::
+
+    \mathbf{W} = \left[ \begin{array}{c}
+    \mathbf{W}_s\\
+    \mathbf{W}_x\\
+    \mathbf{W}_y\end{array} \right]
+
+Then we can write
+
+.. math::
+
+    R(m) \approx \frac{1}{2}\mathbf{m^\top W^\top W m}
+
+The API
+-------
+
+.. autoclass:: SimPEG.Regularization.BaseRegularization
+    :members:
+    :undoc-members:
+
+
+.. autoclass:: SimPEG.Regularization.Tikhonov
+    :show-inheritance:
+    :members:
+
+
+

docs/api_Utilities.rst

@@ -0,0 +1,10 @@
+Utilities
+*********
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Solver
+   api_Maps
+   api_Utils
+   api_Tests

docs/api_Utils.rst

@@ -1,8 +1,5 @@
-.. _api_Utils:
-
-
-Utilities
-*********
+Utils
+*****
 
 .. automodule:: SimPEG.Utils
     :members:
@@ -52,7 +49,7 @@ Interpolation Utilities
     :undoc-members:
 
 Counter Utilities
-=======================
+=================
 
 ::
     class MyClass(object):

docs/api_bigPicture.rst

@@ -1,17 +1,69 @@
-.. _api_license:
-
 Why SimPEG?
-***********
+===========
 
+Our essential functions as researchers are the pursuit and dissemination of knowledge through research and education. As scientists we
+seek to find models that reproduce the observations that we make in the world. In geophysics, we use inverse theory to mathematically
+create models of the earth from measured data. It is a difficult problem with many moving pieces: physics, discretization, simulation,
+regularization, optimization, computer science, linear algebra, geology. Exploring each of these disciplines can take a career, if you
+are so inclined, but as geophysicists we care about the combination: how to pull these disciplines together to answer our questions.
+This is the first problem we hope to help solve: to create a toolbox for the geophysicist that allows you to work at a high level and
+keep your geophysical question in focus. However, a toolbox is not enough. The research questions that we are interested in surround
+the integration of information to make better decisions.
+
+We believe that the feedback loops in the geosciences could use some serious work. For example, collect multiple data-sets from the
+same field area (geology, seismic, electromagnetics, hydrogeology), process the data separately, and then reconvene with your
+multidisciplinary team. You may be rather surprised (or not) that the everyone has a (completely!?) different model. Dissonant at best,
+but often conflicting in the details. Therein lies the second problem: how do we integrate these geoscience fields? Not by force or
+even by default, but at least to have the option of quantitative communication and built in feedback loops. What we require is an
+implementation that is inherently and unequivocally modular, with all pieces available to manipulation. Black-box software, where the
+implementations are hidden, obfuscated, or difficult to manipulate, do not promote experimentation and investigation. We are working on
+a framework that exposes the details of the implementation to the geophysicist in a manner that promotes productivity and question
+based interrogation. This framework can be easily extended to encompass many geophysical problems and is built with the inverse problem
+as the fundamental goal.
+
+The future we see is a mix of tools that span our disciplines, and a framework that allows us to integrate many different types of
+geophysical data so that we can communicate effectively and experiment efficiently. A toolbox combined with a framework that allows you
+to solve your own problems, and creates opportunities for us to work together to better image and understand the subsurface. What we
+are building is called SimPEG, simulation and parameter estimation in geophysics. We are building it in the open. We are testing it.
+Breaking it. Building it. Fixing it. Using it. If you believe, like we do, that geophysics can be more innovative and informative in
+the open and that these tools are necessary and invaluable in education as well as research, then you should get in touch. There is a
+lot of work to do!
 
 The Big Picture
-===============
+---------------
+
+Defining a well-posed inverse problem and solving it is a complex task that requires many components that must interact. It is helpful
+to view this task as a workflow in which various elements are explicitly identified and integrated. The figure below outlines the inversion components that consists of inputs, implementation, and evaluation. The inputs are composed of the geophysical data, the equations which are a mathematical description of the governing physics, and prior knowledge or assumptions about the setting. The implementation consists of two broad categories: the forward simulation and the inversion. The **forward simulation** is the means by which we solve the governing equations given a model and the **inversion components** evaluate and update this model. We are considering a gradient based approach, which updates the model through an optimization routine. The output of this implementation is a model, which, prior to interpretation, must be evaluated. This requires considering, and often re-assessing, the choices and assumptions made in both the input and implementation stages.
+
+.. image:: InversionWorkflow-PreSimPEG.png
+   :width: 400 px
+   :alt: Components
+   :align: center
+
+
+A Comprehensive Framework
+-------------------------
+
+There are an overwhelming amount of choices to be made as one works through the forward modeling and inversion process (see figure above). As a result, software implementations of this workflow often become complex and highly interdependent, making it difficult to interact with and to ask other scientists to pick up and change. Our approach to handling this complexity is to propose a framework, (see below), that compartmentalizes the implementation of inversions into various units. We present it in this specific modular style, as each unit contains a targeted subset of choices crucial to the inversion process.
 
 .. image:: InversionWorkflow.png
    :width: 400 px
    :alt: Framework
    :align: center
 
+The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The (:class:`SimPEG.Mesh.BaseMesh`) class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the (:class:`SimPEG.Survey.BaseSurvey`) and the (:class:`SimPEG.Problem.BaseProblem`). The (:class:`SimPEG.Survey.BaseSurvey`) class handles the geometry of a geophysical problem as well as sources. The (:class:`SimPEG.Problem.BaseProblem`) class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The (:class:`SimPEG.Problem.BaseProblem`) creates geophysical fields given a source from the (:class:`SimPEG.Survey.BaseSurvey`). The (:class:`SimPEG.Survey.BaseSurvey`) interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a (:class:`SimPEG.DataMisfit.BaseDataMisfit`) is chosen to capture the goodness of fit of the predicted data and a (:class:`SimPEG.Regularization.BaseRegularization`) is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class (:class:`SimPEG.InvProblem.BaseInvProblem`). (:class:`SimPEG.InvProblem.BaseInvProblem`) is the mathematical statement that will be numerically solved by running an Inversion. The (:class:`SimPEG.Inversion.BaseInversion`) class handles organization and dispatch of directives between all of the various pieces of the framework.
 
-Explaining The Big Picture
-==========================
+The arrows in the figure above indicate what each class takes as a primary argument. For example, both the (:class:`SimPEG.Problem.BaseProblem`) and (:class:`SimPEG.Regularization.BaseRegularization`) classes take a (:class:`SimPEG.Mesh.BaseMesh`) class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The (:class:`SimPEG.Mesh.BaseMesh`) class, for example, could be a regular Cartesian mesh (:class:`SimPEG.Mesh.TensorMesh`) or a cylindrical coordinate mesh (:class:`SimPEG.Mesh.CylMesh`), which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
+
+
+.. include:: ../CITATION.rst
+
+Authors
+-------
+
+.. include:: ../AUTHORS.rst
+
+License
+-------
+
+.. include:: ../LICENSE

docs/api_installing.rst

@@ -1,7 +1,7 @@
 .. _api_installing:
 
-Installation
-************
+Getting Started with SimPEG
+***************************
 
 Dependencies
 ============

docs/api_license.rst

@@ -1,17 +0,0 @@
-.. _api_license:
-
-License
-*******
-
-.. include:: ../LICENSE
-
-Authors
-*******
-
-.. include:: ../AUTHORS.rst
-
-
-Projects Using SimPEG
-*********************
-
-.. include:: ../PROJECTS.rst

docs/conf.py

@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.4'
+version = '0.1.9'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.4'
+release = '0.1.9'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/em/api_Utils.rst

@@ -7,7 +7,7 @@ sources, and analytic functions.
 Analytic Functions - Time
 =========================
 
-.. automodule:: SimPEG.EM.Utils.Ana.TEM
+.. automodule:: SimPEG.EM.Analytics.TDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -17,7 +17,7 @@ Analytic Functions - Time
 Analytic Functions - Frequency
 ==============================
 
-.. automodule:: SimPEG.EM.Utils.Ana.FEM
+.. automodule:: SimPEG.EM.Analytics.FDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -27,8 +27,7 @@ Analytic Functions - Frequency
 Sources
 =======
 
-.. automodule:: SimPEG.EM.Utils.Sources.magneticDipole
+.. autoclass:: SimPEG.EM.FDEM.SrcFDEM.MagDipole
     :show-inheritance:
     :members:
     :undoc-members:
-    :inherited-members:

docs/em/index.rst

@@ -9,7 +9,7 @@ Time Domian Electromagnetics
 ----------------------------
 
 .. toctree::
-   :maxdepth: 2
+    :maxdepth: 2
 
     api_TDEM_derivation
 
@@ -18,7 +18,7 @@ Code for Time Domian Electromagnetics
 -------------------------------------
 
 .. toctree::
-   :maxdepth: 2
+    :maxdepth: 2
 
     api_TDEM
 
@@ -28,7 +28,6 @@ Frequency Domian Electromagnetics
 .. toctree::
    :maxdepth: 2
 
-   api_ForwardProblem
    api_FDEM
 
 

docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst

@@ -0,0 +1,26 @@
+.. _examples_EM_FDEM_Analytic_MagDipoleWholespace:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: Magnetic Dipole in a Whole-Space
+====================================
+
+Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_Analytic_MagDipoleWholespace.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
+    :language: python
+    :linenos:

docs/examples/EM_TDEM_1D_Inversion.rst

@@ -1,4 +1,4 @@
-.. _examples_EM_FDEM_1D_Inversion:
+.. _examples_EM_TDEM_1D_Inversion:
 
 .. --------------------------------- ..
 ..                                   ..
@@ -9,18 +9,18 @@
 .. --------------------------------- ..
 
 
-EM: FDEM: 1D: Inversion
+EM: TDEM: 1D: Inversion
 =======================
 
-Here we will create and run a FDEM 1D inversion.
+Here we will create and run a TDEM 1D inversion.
 
 
 
 .. plot::
 
     from SimPEG import Examples
-    Examples.EM_FDEM_1D_Inversion.run()
+    Examples.EM_TDEM_1D_Inversion.run()
 
-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+.. literalinclude:: ../../SimPEG/Examples/EM_TDEM_1D_Inversion.py
     :language: python
     :linenos:

docs/flow/index.rst

@@ -41,7 +41,7 @@ Here we reproduce the results from Celia et al. (1990):
 Richards
 ========
 
-.. automodule:: simpegFLOW.Richards.Empirical
+.. automodule:: SimPEG.FLOW.Richards.Empirical
     :show-inheritance:
     :members:
     :undoc-members:

docs/index.rst

@@ -1,25 +1,38 @@
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
+    :alt: SimPEG Logo
+
 SimPEG Documentation
 ********************
 
-.. image:: simpeg-logo.png
-   :width: 300 px
-   :alt: SimPEG
-   :align: center
+.. image:: https://img.shields.io/pypi/v/SimPEG.svg
+    :target: https://crate.io/packages/SimPEG/
+    :alt: Latest PyPI version
 
-SimPEG: Simulation and Parameter Estimation in Geophysics
+.. image:: https://img.shields.io/pypi/dm/SimPEG.svg
+    :target: https://crate.io/packages/SimPEG/
+    :alt: Number of PyPI downloads
 
-SimPEG is a framework and a collection of tools that aid in the development of
-large-scale geophysical inversion codes.
-The vision is to create a modular and extensible package for
-finite volume simulation and parameter estimation with
-applications to geophysical imaging and subsurface flow. To enable
-these goals, this package has the following features:
+.. image:: https://img.shields.io/badge/license-MIT-blue.svg
+    :target: https://github.com/simpeg/simpeg/blob/master/LICENSE
+    :alt: BSD 3 clause license.
 
-- is modular with respect to ... everything!
-- is built with the (large-scale) inverse problem in mind
-- provides a framework for geophysical and hydrogeologic problems
-- supports 1D, 2D and 3D problems
-- provides a set of commonly used visualization utilities
+.. image:: https://img.shields.io/travis/simpeg/simpeg.svg
+    :target: https://travis-ci.org/simpeg/simpeg?branch=master
+    :alt: Travis CI build status
+
+.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
+    :target: https://coveralls.io/r/simpeg/simpeg?branch=master
+    :alt: Coverage status
+
+Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+
+Our 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:
+
+* modular with respect to the spacial discretization, optimization routine, and geophysical problem
+* built with the inverse problem in mind
+* provides a framework for geophysical and hydrogeologic problems
+* supports 1D, 2D and 3D problems
+* designed for large-scale inversions
 
 
 About SimPEG
@@ -29,55 +42,16 @@ About SimPEG
    :maxdepth: 2
 
    api_bigPicture
-   api_license
-
-
-Getting Started with SimPEG
-***************************
-
-.. toctree::
-   :maxdepth: 2
-
    api_installing
 
-Discretization
-**************
-
-.. toctree::
-   :maxdepth: 3
-
-   api_Mesh
-   api_DiffOps
-   api_InnerProducts
-
-Forward Problems
-****************
+Examples
+********
 
 .. toctree::
    :maxdepth: 2
 
-   api_ForwardProblem
-
-Inversion
-*********
-
-.. toctree::
-   :maxdepth: 3
-
-   api_DataMisfit
-   api_Inverse
-
-Utility Codes
-*************
-
-.. toctree::
-   :maxdepth: 2
-
-   api_Solver
-   api_Maps
-   api_Utils
-   api_Tests
-
+   api_Examples
+   
 
 Packages
 ********
@@ -88,20 +62,38 @@ Packages
    em/index
    flow/index
 
-Developer's Documentation
-*************************
 
-* Travis-CI Testing
-   .. image:: https://travis-ci.org/simpeg/simpeg.svg?branch=master
-      :target: https://travis-ci.org/simpeg/simpeg
-      :alt: Master Branch
-      :align: center
+Finite Volume
+*************
 
-* Coveralls Testing
-   .. image:: https://coveralls.io/repos/simpeg/simpeg/badge.png?branch=master
-      :target: https://coveralls.io/r/simpeg/simpeg?branch=master
-      :alt: Coveralls
-      :align: center
+.. toctree::
+   :maxdepth: 3
+
+   api_FiniteVolume
+
+Forward Problems
+****************
+
+.. toctree::
+   :maxdepth: 3
+
+   api_ForwardProblem
+
+Inversion Components
+********************
+
+.. toctree::
+   :maxdepth: 3
+
+   api_InversionComponents
+
+Utility Codes
+*************
+
+.. toctree::
+   :maxdepth: 3
+
+   api_Utilities
 
 
 Project Index & Search
@@ -110,11 +102,3 @@ Project Index & Search
 * :ref:`genindex`
 * :ref:`modindex`
 * :ref:`search`
-
-Examples
-********
-
-.. toctree::
-   :maxdepth: 2
-
-   api_Examples

setup.py

@@ -64,6 +64,7 @@ cython_files = [
                     "SimPEG/Mesh/TreeUtils"
                ]
 extensions = [Extension(f, [f+ext]) for f in cython_files]
+scripts = [f+'.pyx' for f in cython_files]
 
 if USE_CYTHON and "cleanall" not in args:
     from Cython.Build import cythonize
@@ -76,7 +77,7 @@ with open("README.rst") as f:
 
 setup(
     name = "SimPEG",
-    version = "0.1.4",
+    version = "0.1.9",
     packages = find_packages(),
     install_requires = ['numpy>=1.7',
                         'scipy>=0.13',
@@ -95,5 +96,6 @@ setup(
     use_2to3 = False,
     include_dirs=[np.get_include()],
     ext_modules = extensions,
+    scripts=scripts,
     **cythonKwargs
 )

tests/base/test_regularization.py

@@ -4,11 +4,17 @@ from SimPEG import *
 from scipy.sparse.linalg import dsolve
 import inspect
 
+TOL = 1e-20
 
 class RegularizationTests(unittest.TestCase):
 
     def setUp(self):
-        self.mesh2 = Mesh.TensorMesh([3, 2])
+        hx, hy, hz = np.random.rand(10), np.random.rand(9), np.random.rand(8)
+        hx, hy, hz = hx/hx.sum(), hy/hy.sum(), hz/hz.sum()
+        mesh1 = Mesh.TensorMesh([hx])
+        mesh2 = Mesh.TensorMesh([hx, hy])
+        mesh3 = Mesh.TensorMesh([hx, hy, hz])
+        self.meshlist = [mesh1,mesh2, mesh3]
 
     def test_regularization(self):
         for R in dir(Regularization):
@@ -16,18 +22,63 @@ class RegularizationTests(unittest.TestCase):
             if not inspect.isclass(r): continue
             if not issubclass(r, Regularization.BaseRegularization):
                 continue
-            # if 'Regularization' not in R: continue
-            mapping = r.mapPair(self.mesh2)
-            reg = r(self.mesh2, mapping=mapping)
-            m = np.random.rand(mapping.nP)
-            reg.mref = m[:]*np.mean(m)
 
-            print 'Check:', R
-            passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-            self.assertTrue(passed)
-            print 'Check 2 Deriv:', R
-            passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-            self.assertTrue(passed)
+            for i, mesh in enumerate(self.meshlist):
+
+                print 'Testing %iD'%mesh.dim
+
+                mapping = r.mapPair(mesh)
+                reg = r(mesh, mapping=mapping)
+                m = np.random.rand(mapping.nP)
+                reg.mref = np.ones_like(m)*np.mean(m)
+
+                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                passed = reg.eval(reg.mref) < TOL
+                self.assertTrue(passed)
+
+                print 'Check:', R
+                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+                print 'Check 2 Deriv:', R
+                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+    def test_regularization_ActiveCells(self):
+        for R in dir(Regularization):
+            r = getattr(Regularization, R)
+            if not inspect.isclass(r): continue
+            if not issubclass(r, Regularization.BaseRegularization):
+                continue
+
+            for i, mesh in enumerate(self.meshlist):
+
+                print 'Testing Active Cells %iD'%(mesh.dim)
+
+                if mesh.dim == 1:
+                    indAct = Utils.mkvc(mesh.gridCC <= 0.8)
+                elif mesh.dim == 2:
+                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                elif mesh.dim == 3:
+                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
+
+                mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
+
+                reg = r(mesh, mapping=mapping, indActive=indAct)
+                m = np.random.rand(mesh.nC)[indAct]
+                reg.mref = np.ones_like(m)*np.mean(m)
+
+                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                passed = reg.eval(reg.mref) < TOL
+                self.assertTrue(passed)
+
+                print 'Check:', R
+                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+                print 'Check 2 Deriv:', R
+                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
 
 
 if __name__ == '__main__':

tests/mesh/test_MeshIO.py

@@ -0,0 +1,100 @@
+import numpy as np
+import unittest, os
+import SimPEG as simpeg
+from SimPEG.Mesh import TensorMesh, TreeMesh
+
+
+class TestTensorMeshIO(unittest.TestCase):
+
+    def setUp(self):
+        h = np.ones(16)
+        mesh = TensorMesh([h,2*h,3*h])
+        self.mesh = mesh
+
+    def test_UBCfiles(self):
+
+        mesh = self.mesh
+        # Make a vector
+        vec = np.arange(mesh.nC)
+        # Write and read
+        mesh.writeUBC('temp.msh', {'arange.txt':vec})
+        meshUBC = TensorMesh.readUBC('temp.msh')
+        vecUBC  = meshUBC.readModelUBC('arange.txt')
+
+        # The mesh
+        assert mesh.__str__() == meshUBC.__str__()
+        assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0
+        assert np.sum(vec - vecUBC) == 0
+        assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0)
+
+
+        vecUBC  = mesh.readModelUBC('arange.txt')
+        assert np.sum(vec - vecUBC) == 0
+
+        mesh.writeModelUBC('arange2.txt', vec + 1)
+        vec2UBC  = mesh.readModelUBC('arange2.txt')
+        assert np.sum(vec + 1 - vec2UBC) == 0
+
+        print 'IO of UBC tensor mesh files is working'
+        os.remove('temp.msh')
+        os.remove('arange.txt')
+        os.remove('arange2.txt')
+
+    def test_VTKfiles(self):
+        mesh = self.mesh
+        vec = np.arange(mesh.nC)
+
+        mesh.writeVTK('temp.vtr', {'arange.txt':vec})
+        meshVTR, models = TensorMesh.readVTK('temp.vtr')
+
+        assert mesh.__str__() == meshVTR.__str__()
+        assert np.all(np.array(mesh.h) - np.array(meshVTR.h) == 0)
+
+        assert 'arange.txt' in models
+        vecVTK = models['arange.txt']
+        assert np.sum(vec - vecVTK) == 0
+
+        print 'IO of VTR tensor mesh files is working'
+        os.remove('temp.vtr')
+
+
+class TestOcTreeMeshIO(unittest.TestCase):
+
+    def setUp(self):
+        h = np.ones(16)
+        mesh = TreeMesh([h,2*h,3*h])
+        mesh.refine(3)
+        mesh._refineCell([0,0,0,3])
+        mesh._refineCell([0,2,0,3])
+        self.mesh = mesh
+
+    def test_UBCfiles(self):
+
+        mesh = self.mesh
+        # Make a vector
+        vec = np.arange(mesh.nC)
+        # Write and read
+        mesh.writeUBC('temp.msh', {'arange.txt':vec})
+        meshUBC = TreeMesh.readUBC('temp.msh')
+        vecUBC  = meshUBC.readModelUBC('arange.txt')
+
+        # The mesh
+        assert mesh.__str__() == meshUBC.__str__()
+        assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0
+        assert np.sum(vec - vecUBC) == 0
+        assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0)
+        print 'IO of UBC octree files is working'
+        os.remove('temp.msh')
+        os.remove('arange.txt')
+
+    def test_VTUfiles(self):
+        mesh = self.mesh
+        vec = np.arange(mesh.nC)
+        mesh.writeVTK('temp.vtu',{'arange':vec})
+        print 'Writing of VTU files is working'
+        os.remove('temp.vtu')
+
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mesh/test_TreeMesh.py

@@ -26,6 +26,27 @@ class TestSimpleQuadTree(unittest.TestCase):
 
         assert np.allclose(np.r_[M._areaFxFull, M._areaFyFull], M._deflationMatrix('F') * M.area)
 
+    def test_getitem(self):
+        M = Mesh.TreeMesh([4,4])
+        M.refine(1)
+        assert M.nC == 4
+        assert len(M) == M.nC
+        assert np.allclose(M[0].center, [0.25,0.25])
+        actual = [[0,0],[0.5,0],[0,0.5],[0.5,0.5]]
+        for i, n in enumerate(M[0].nodes):
+            assert np.allclose(M._gridN[n,:], actual[i])
+
+    def test_getitem3D(self):
+        M = Mesh.TreeMesh([4,4,4])
+        M.refine(1)
+        assert M.nC == 8
+        assert len(M) == M.nC
+        assert np.allclose(M[0].center, [0.25,0.25,0.25])
+        actual = [[0,0,0],[0.5,0,0],[0,0.5,0],[0.5,0.5,0],
+                  [0,0,0.5],[0.5,0,0.5],[0,0.5,0.5],[0.5,0.5,0.5]]
+        for i, n in enumerate(M[0].nodes):
+            assert np.allclose(M._gridN[n,:], actual[i])
+
     def test_refine(self):
         M = Mesh.TreeMesh([4,4,4])
         M.refine(1)