Utils

Matrix Utilities

SimPEG.Utils.matutils.mkvc(x, numDims=1)[source]

Creates a vector with the number of dimension specified

e.g.:

a = np.array([1, 2, 3])

mkvc(a, 1).shape
    > (3, )

mkvc(a, 2).shape
    > (3, 1)

mkvc(a, 3).shape
    > (3, 1, 1)
SimPEG.Utils.matutils.sdiag(h)[source]

Sparse diagonal matrix

SimPEG.Utils.matutils.sdInv(M)[source]

Inverse of a sparse diagonal matrix

SimPEG.Utils.matutils.speye(n)[source]

Sparse identity

SimPEG.Utils.matutils.kron3(A, B, C)[source]

Three kron prods

SimPEG.Utils.matutils.spzeros(n1, n2)[source]
SimPEG.Utils.matutils.ddx(n)[source]

Define 1D derivatives, inner, this means we go from n+1 to n

SimPEG.Utils.matutils.av(n)[source]

Define 1D averaging operator from nodes to cell-centers.

SimPEG.Utils.matutils.avExtrap(n)[source]

Define 1D averaging operator from cell-centers to nodes.

SimPEG.Utils.matutils.ndgrid(*args, **kwargs)[source]

Form tensorial grid for 1, 2, or 3 dimensions.

Returns as column vectors by default.

To return as matrix input:

ndgrid(..., vector=False)

The inputs can be a list or separate arguments.

e.g.:

a = np.array([1, 2, 3])
b = np.array([1, 2])

XY = ndgrid(a, b)
    > [[1 1]
       [2 1]
       [3 1]
       [1 2]
       [2 2]
       [3 2]]

X, Y = ndgrid(a, b, vector=False)
    > X = [[1 1]
           [2 2]
           [3 3]]
    > Y = [[1 2]
           [1 2]
           [1 2]]
SimPEG.Utils.matutils.ind2sub(shape, inds)[source]

From the given shape, returns the subscripts of the given index

SimPEG.Utils.matutils.sub2ind(shape, subs)[source]

From the given shape, returns the index of the given subscript

SimPEG.Utils.matutils.getSubArray(A, ind)[source]

subArray

SimPEG.Utils.matutils.inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33, returnMatrix=True)[source]

B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)

inverts a stack of 3x3 matrices

Input:
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
Output:
B - inverse
SimPEG.Utils.matutils.inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True)[source]

B = inv2X2BlockDiagonal(a11, a12, a21, a22)

Inverts a stack of 2x2 matrices by using the inversion formula

inv(A) = (1/det(A)) * cof(A)^T

Input: A - a11, a12, a21, a22

Output: B - inverse

class SimPEG.Utils.matutils.TensorType(M, tensor)[source]
SimPEG.Utils.matutils.makePropertyTensor(M, tensor)[source]
SimPEG.Utils.matutils.invPropertyTensor(M, tensor, returnMatrix=False)[source]
SimPEG.Utils.matutils.diagEst(matFun, n, k=None, approach='Probing')[source]

Estimate the diagonal of a matrix, A. Note that the matrix may be a function which returns A times a vector.

Three different approaches have been implemented:

  1. Probing: cyclic permutations of vectors with 1’s and 0’s (default)
  2. Ones: random +/- 1 entries
  3. Random: random vectors
Parameters:
  • matFun (callable) – takes a (numpy.array) and multiplies it by a matrix to estimate the diagonal
  • n (int) – size of the vector that should be used to compute matFun(v)
  • k (int) – number of vectors to be used to estimate the diagonal
  • approach (str) – approach to be used for getting vectors
Return type:

numpy.array

Returns:

est_diag(A)

Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf

Solver Utilities

SimPEG.Utils.SolverUtils.SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-06, name=None)[source]

Wraps a direct Solver.

Solver   = SolverUtils.SolverWrapD(sp.linalg.spsolve, factorize=False)
SolverLU = SolverUtils.SolverWrapD(sp.linalg.splu, factorize=True)
SimPEG.Utils.SolverUtils.SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-05, name=None)[source]

Wraps an iterative Solver.

SolverCG = SolverUtils.SolverWrapI(sp.linalg.cg)
class SimPEG.Utils.SolverUtils.Solver(A, **kwargs)
clean()
class SimPEG.Utils.SolverUtils.SolverLU(A, **kwargs)
clean()
class SimPEG.Utils.SolverUtils.SolverCG(A, **kwargs)
clean()
class SimPEG.Utils.SolverUtils.SolverBiCG(A, **kwargs)
clean()
class SimPEG.Utils.SolverUtils.SolverDiag(A)[source]

docstring for SolverDiag

clean()[source]

Curv Utilities

SimPEG.Utils.curvutils.volTetra(xyz, A, B, C, D)[source]

Returns the volume for tetrahedras volume specified by the indexes A to D.

Parameters:
Return type:

numpy.array

Returns:

V, volume of the tetrahedra

Algorithm https://en.wikipedia.org/wiki/Tetrahedron#Volume

\[ \begin{align}\begin{aligned}V = {1 \over 3} A h\\V = {1 \over 6} | ( a - d ) \cdot ( ( b - d ) ( c - d ) ) |\end{aligned}\end{align} \]
SimPEG.Utils.curvutils.indexCube(nodes, gridSize, n=None)[source]

Returns the index of nodes on the mesh.

Input:
nodes - string of which nodes to return. e.g. ‘ABCD’ gridSize - size of the nodal grid n - number of nodes each i,j,k direction: [ni,nj,nk]
Output:
index - index in the order asked e.g. ‘ABCD’ –> (A,B,C,D)

TWO DIMENSIONS:

node(i,j)          node(i,j+1)
     A -------------- B
     |                |
     |    cell(i,j)   |
     |        I       |
     |                |
    D -------------- C
node(i+1,j)        node(i+1,j+1)

THREE DIMENSIONS:

      node(i,j,k+1)       node(i,j+1,k+1)
          E --------------- F
         /|               / |
        / |              /  |
       /  |             /   |
node(i,j,k)         node(i,j+1,k)
     A -------------- B     |
     |    H ----------|---- G
     |   /cell(i,j)   |   /
     |  /     I       |  /
     | /              | /
     D -------------- C
node(i+1,j,k)      node(i+1,j+1,k)
SimPEG.Utils.curvutils.faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True)[source]

function [N] = faceInfo(y,A,B,C,D)

Returns the averaged normal, area, and edge lengths for a given set of faces.

If average option is FALSE then N is a cell array {nA,nB,nC,nD}

Input:
xyz - X,Y,Z vertex vector A,B,C,D - vert index of the face (counter clockwize)
Options:
average - [true]/false, toggles returning all normals or the average
Output:
N - average face normal or {nA,nB,nC,nD} if average = false area - average face area edgeLengths - exact edge Lengths, 4 column vector [AB, BC, CD, DA]

see also testFaceNormal testFaceArea

@author Rowan Cockett

Last modified on: 2013/07/26

Mesh Utilities

SimPEG.Utils.meshutils.exampleLrmGrid(nC, exType)[source]
SimPEG.Utils.meshutils.meshTensor(value)[source]

meshTensor takes a list of numbers and tuples that have the form:

mT = [ float, (cellSize, numCell), (cellSize, numCell, factor) ]

For example, a time domain mesh code needs many time steps at one time:

[(1e-5, 30), (1e-4, 30), 1e-3]

Means take 30 steps at 1e-5 and then 30 more at 1e-4, and then one step of 1e-3.

Tensor meshes can also be created by increase factors:

[(10.0, 5, -1.3), (10.0, 50), (10.0, 5, 1.3)]

When there is a third number in the tuple, it refers to the increase factor, if this number is negative this section of the tensor is flipped right-to-left.

(Source code, png, hires.png, pdf)

../../_images/api_Utils-1.png
SimPEG.Utils.meshutils.closestPoints(mesh, pts, gridLoc='CC')[source]

Move a list of points to the closest points on a grid.

Parameters:
  • mesh (BaseMesh) – The mesh
  • pts (numpy.ndarray) – Points to move
  • gridLoc (string) – [‘CC’, ‘N’, ‘Fx’, ‘Fy’, ‘Fz’, ‘Ex’, ‘Ex’, ‘Ey’, ‘Ez’]
Return type:

numpy.ndarray

Returns:

nodeInds

SimPEG.Utils.meshutils.ExtractCoreMesh(xyzlim, mesh, meshType='tensor')[source]

Extracts Core Mesh from Global mesh

Parameters:

This function ouputs:

- actind: corresponding boolean index from global to core
- meshcore: core SimPEG mesh

Model Builder Utilities

SimPEG.Utils.ModelBuilder.addBlock(gridCC, modelCC, p0, p1, blockProp)[source]

Add a block to an exsisting cell centered model, modelCC

Parameters:
BlockProp float blockProp:
 

property to assign to the model

Return numpy.array, modelBlock:
 

model with block

SimPEG.Utils.ModelBuilder.getIndicesBlock(p0, p1, ccMesh)[source]

Creates a vector containing the block indices in the cell centers mesh. Returns a tuple

The block is defined by the points

p0, describe the position of the left upper front corner, and

p1, describe the position of the right bottom back corner.

ccMesh represents the cell-centered mesh

The points p0 and p1 must live in the the same dimensional space as the mesh.

SimPEG.Utils.ModelBuilder.defineBlock(ccMesh, p0, p1, vals=None)[source]

Build a block with the conductivity specified by condVal. Returns an array. vals[0] conductivity of the block vals[1] conductivity of the ground

SimPEG.Utils.ModelBuilder.defineElipse(ccMesh, center=None, anisotropy=None, slope=10.0, theta=0.0)[source]
SimPEG.Utils.ModelBuilder.getIndicesSphere(center, radius, ccMesh)[source]

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.

SimPEG.Utils.ModelBuilder.defineTwoLayers(ccMesh, depth, vals=None)[source]

Define a two layered model. Depth of the first layer must be specified. CondVals vector with the conductivity values of the layers. Eg:

Convention to number the layers:

<----------------------------|------------------------------------>
0                          depth                                 zf
     1st layer                       2nd layer
SimPEG.Utils.ModelBuilder.scalarConductivity(ccMesh, pFunction)[source]

Define the distribution conductivity in the mesh according to the analytical expression given in pFunction

SimPEG.Utils.ModelBuilder.layeredModel(ccMesh, layerTops, layerValues)[source]

Define a layered model from layerTops (z-positive up)

Parameters:
  • ccMesh (numpy.array) – cell-centered mesh
  • layerTops (numpy.array) – z-locations of the tops of each layer
  • layerValue (numpy.array) – values of the property to assign for each layer (starting at the top)
Return type:

numpy.array

Returns:

M, layered model on the mesh

SimPEG.Utils.ModelBuilder.randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None)[source]

Create a random model by convolving a kernel with a uniformly distributed model.

Parameters:
  • shape (tuple) – shape of the model.
  • seed (int) – pick which model to produce, prints the seed if you don’t choose.
  • anisotropy (numpy.ndarray) – this is the (3 x n) blurring kernel that is used.
  • its (int) – number of smoothing iterations
  • bounds (list) – bounds on the model, len(list) == 2
Return type:

numpy.ndarray

Returns:

M, the model

(Source code, png, hires.png, pdf)

../../_images/api_Utils-2.png

Interpolation Utilities

SimPEG.Utils.interputils.interpmat(locs, x, y=None, z=None)[source]

Local interpolation computed for each receiver point in turn

Parameters:
  • loc (numpy.ndarray) – Location of points to interpolate to
  • x (numpy.ndarray) – Tensor vector of 1st dimension of grid.
  • y (numpy.ndarray) – Tensor vector of 2nd dimension of grid. None by default.
  • z (numpy.ndarray) – Tensor vector of 3rd dimension of grid. None by default.
Return type:

scipy.sparse.csr_matrix

Returns:

Interpolation matrix

(Source code, png, hires.png, pdf)

../../_images/api_Utils-3.png

Counter Utilities

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
class MyClass(object):
    def __init__(self, url):
        self.counter = Counter()

    @count
    def MyMethod(self):
        pass

    @timeIt
    def MySecondMethod(self):
        pass

c = MyClass('blah')
for i in range(100): c.MyMethod()
for i in range(300): c.MySecondMethod()
c.counter.summary()
1
2
3
4
5
Counters:
  MyClass.MyMethod                        :      100

Times:                                        mean      sum
  MyClass.MySecondMethod                  : 1.70e-06, 5.10e-04,  300x

The API

class SimPEG.Utils.CounterUtils.Counter[source]

Counter allows anything that calls it to record iterations and timings in a simple way.

Also has plotting functions that allow quick recalls of data.

If you want to use this, import count or timeIt and use them as decorators on class methods.

class MyClass(object):
    def __init__(self, url):
        self.counter = Counter()

    @count
    def MyMethod(self):
        pass

    @timeIt
    def MySecondMethod(self):
        pass

c = MyClass('blah')
for i in range(100): c.MyMethod()
for i in range(300): c.MySecondMethod()
c.counter.summary()
count(prop)[source]

Increases the count of the property.

countTic(prop)[source]

Times a property call, this is the init call.

countToc(prop)[source]

Times a property call, this is the end call.

summary()[source]

Provides a text summary of the current counters and timers.

SimPEG.Utils.CounterUtils.count(f)[source]
SimPEG.Utils.CounterUtils.timeIt(f)[source]