wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-13 17:45:30 +08:00
Code Issues Packages Projects Releases Wiki Activity

renaming to ensure capitals

Browse Source
This commit is contained in:
rowanc1
2014-01-16 13:22:46 -08:00
parent 7432591450
commit fa8a5cd7cb
48 changed files with 0 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Utils/Geophysics/__init__.py
+1
View File
@@ -0,0 +1 @@
import emSources
SimPEG/Utils/Geophysics/emSources/__init__.py
+1
View File
@@ -0,0 +1 @@
from emSources import MagneticDipoleVectorPotential
SimPEG/Utils/Geophysics/emSources/emSources.py
+40
View File
@@ -0,0 +1,40 @@
import numpy as np
from scipy.constants import mu_0, pi
def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
"""
Calculate the vector potential of a set of magnetic dipoles
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
:param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z)
:param str component: The component to calculate - 'x', 'y', or 'z'
:param numpy.ndarray dipoleMoment: The vector dipole moment
:rtype: numpy.ndarray
:return: The vector potential each dipole at each observation location
"""
if component=='x':
dimInd = 0
elif component=='y':
dimInd = 1
elif component=='z':
dimInd = 2
else:
raise ValueError('Invalid component')
txLoc = np.atleast_2d(txLoc)
obsLoc = np.atleast_2d(obsLoc)
dipoleMoment = np.atleast_2d(dipoleMoment)
nEdges = obsLoc.shape[0]
nTx = txLoc.shape[0]
m = np.array(dipoleMoment).repeat(nEdges, axis=0)
A = np.empty((nEdges, nTx))
for i in range(nTx):
dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0)
mCr = np.cross(m, dR)
r = np.sqrt((dR**2).sum(axis=1))
A[:, i] = -(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
return A
SimPEG/Utils/ModelBuilder.py
+270
View File
@@ -0,0 +1,270 @@
import numpy as np
import scipy.ndimage as ndi
import scipy.sparse as sp
def getIndecesBlock(p0,p1,ccMesh):
"""
Creates a vector containing the block indexes in the cell centerd 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.
"""
# Validation: p0 and p1 live in the same dimensional space
assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)"
# Validation: mesh and points live in the same dimensional space
dimMesh = np.size(ccMesh[0,:])
assert len(p0) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
if dimMesh == 1:
# Define the reference points
x1 = p0[0]
x2 = p1[0]
indX = (x1 <= ccMesh[:,0]) & (ccMesh[:,0] <= x2)
ind = np.where(indX)
elif dimMesh == 2:
# Define the reference points
x1 = p0[0]
y1 = p0[1]
x2 = p1[0]
y2 = p1[1]
indX = (x1 <= ccMesh[:,0]) & (ccMesh[:,0] <= x2)
indY = (y1 <= ccMesh[:,1]) & (ccMesh[:,1] <= y2)
ind = np.where(indX & indY)
elif dimMesh == 3:
# Define the points
x1 = p0[0]
y1 = p0[1]
z1 = p0[2]
x2 = p1[0]
y2 = p1[1]
z2 = p1[2]
indX = (x1 <= ccMesh[:,0]) & (ccMesh[:,0] <= x2)
indY = (y1 <= ccMesh[:,1]) & (ccMesh[:,1] <= y2)
indZ = (z1 <= ccMesh[:,2]) & (ccMesh[:,2] <= z2)
ind = np.where(indX & indY & indZ)
# Return a tuple
return ind
def defineBlockConductivity(p0,p1,ccMesh,condVals):
"""
Build a block with the conductivity specified by condVal. Returns an array.
condVals[0] conductivity of the block
condVals[1] conductivity of the ground
"""
sigma = np.zeros(ccMesh.shape[0]) + condVals[1]
ind = getIndecesBlock(p0,p1,ccMesh)
sigma[ind] = condVals[0]
return sigma
def defineTwoLayeredConductivity(depth,ccMesh,condVals):
"""
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
"""
sigma = np.zeros(ccMesh.shape[0]) + condVals[1]
dim = np.size(ccMesh[0,:])
p0 = np.zeros(dim)
p1 = np.zeros(dim)
# Identify 1st cell centered reference point
p0[0] = ccMesh[0,0]
if dim>1: p0[1] = ccMesh[0,1]
if dim>2: p0[2] = ccMesh[0,2]
# Identify the last cell-centered reference point
p1[0] = ccMesh[-1,0]
if dim>1: p1[1] = ccMesh[-1,1]
if dim>2: p1[2] = ccMesh[-1,2]
# The depth is always defined on the last one.
p1[len(p1)-1] -= depth
ind = getIndecesBlock(p0,p1,ccMesh)
sigma[ind] = condVals[0];
return sigma
def scalarConductivity(ccMesh,pFunction):
"""
Define the distribution conductivity in the mesh according to the
analytical expression given in pFunction
"""
dim = np.size(ccMesh[0,:])
CC = [ccMesh[:,0]]
if dim>1: CC.append(ccMesh[:,1])
if dim>2: CC.append(ccMesh[:,2])
sigma = pFunction(*CC)
return sigma
def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
"""
Create a random model by convolving a kernal with a
uniformly distributed model.
:param int,tuple shape: shape of the model.
:param int seed: pick which model to produce, prints the seed if you don't choose.
:param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
:param int its: number of smoothing iterations
:param list bounds: bounds on the model, len(list) == 2
:rtype: numpy.ndarray
:return: M, the model
.. plot::
import matplotlib.pyplot as plt
import SimPEG.Utils.ModelBuilder as MB
plt.colorbar(plt.imshow(MB.randomModel((50,50),bounds=[-4,0])))
plt.title('A very cool, yet completely random model.')
plt.show()
"""
if seed is None:
seed = np.random.randint(1e3)
print 'Using a seed of: ', seed
if type(shape) in [int, long, float]:
shape = (shape,) # make it a tuple for consistency
np.random.seed(seed)
mr = np.random.rand(*shape)
if anisotropy is None:
if len(shape) is 1:
smth = np.array([1,10.,1],dtype=float)
elif len(shape) is 2:
smth = np.array([[1,7,1],[2,10,2],[1,7,1]],dtype=float)
elif len(shape) is 3:
kernal = np.array([1,4,1], dtype=float).reshape((1,3))
smth = np.array(sp.kron(sp.kron(kernal,kernal.T).todense()[:],kernal).todense()).reshape((3,3,3))
else:
assert len(anisotropy.shape) is len(shape), 'Anisotropy must be the same shape.'
smth = np.array(anisotropy,dtype=float)
smth = smth/smth.sum() # normalize
mi = mr
for i in range(its):
mi = ndi.convolve(mi, smth)
# scale the model to live between the bounds.
mi = (mi - mi.min())/(mi.max()-mi.min()) # scaled between 0 and 1
mi = mi*(bounds[1]-bounds[0])+bounds[0]
return mi
if __name__ == '__main__':
from SimPEG.mesh import TensorMesh
from matplotlib import pyplot as plt
# Define the mesh
testDim = 2
h1 = 0.3*np.ones(7)
h1[0] = 0.5
h1[-1] = 0.6
h2 = .5 * np.ones(4)
h3 = .4 * np.ones(6)
x0 = np.zeros(3)
if testDim == 1:
h = [h1]
x0 = x0[0]
elif testDim == 2:
h = [h1, h2]
x0 = x0[0:2]
else:
h = [h1, h2, h3]
M = TensorMesh(h, x0)
ccMesh = M.gridCC
# ------------------- Test conductivities! --------------------------
print('Testing 1 block conductivity')
p0 = np.array([0.5,0.5,0.5])[:testDim]
p1 = np.array([1.0,1.0,1.0])[:testDim]
condVals = np.array([100,1e-6])
sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
# Plot sigma model
print sigma.shape
M.plotImage(sigma)
print 'Done with block! :)'
plt.show()
# -----------------------------------------
print('Testing the two layered model')
condVals = np.array([100,1e-5]);
depth = 1.0;
sigma = defineTwoLayeredConductivity(depth,ccMesh,condVals)
M.plotImage(sigma)
print sigma
print 'layer model!'
plt.show()
# -----------------------------------------
print('Testing scalar conductivity')
if testDim == 1:
pFunction = lambda x: np.exp(x)
elif testDim == 2:
pFunction = lambda x,y: np.exp(x+y)
elif testDim == 3:
pFunction = lambda x,y,z: np.exp(x+y+z)
sigma = scalarConductivity(ccMesh,pFunction)
# Plot sigma model
M.plotImage(sigma)
print sigma
print 'Scalar conductivity defined!'
plt.show()
# -----------------------------------------
SimPEG/Utils/Save.py
+352
View File
@@ -0,0 +1,352 @@
import numpy as np
import time
import re
try:
import h5py
except Exception, e:
print 'Warning: SimPEG.Utils.Save needs h5py to be installed.'
SAVEABLES = {}
def natural_keys(text):
'''
alist.sort(key=natural_keys) sorts in human order
http://nedbatchelder.com/blog/200712/human_sorting.html
(See Toothy's implementation in the comments)
'''
atoi = lambda text: int(text) if text.isdigit() else text
return [ atoi(c) for c in re.split('(\d+)', text) ]
def preIteration(group):
group.attrs['complete'] = False
group.attrs['time'] = time.time()
def postIteration(group):
group.attrs['time'] = time.time() - group.attrs['time']
group.attrs['date'] = time.ctime()
group.attrs['complete'] = True
class SimPEGTable:
"""
This is a wrapper class on the HDF5 file.
"""
def __init__(self, name, mode='a'):
if '.hdf5' not in name:
name += '.hdf5'
self.f = h5py.File(name, mode)
self.root = hdf5Group(self,self.f)
self.inversions = hdf5InversionGroup(self,self.root.addGroup('inversions',soft=True))
def show(self): self.root.show()
def saveInversion(self, invObj):
# Create a new inversion anytime this is run.
def _startup_hdf5_inv(invObj, m0):
node = self.inversions.addGroup('%d'%self.inversions.numChildren)
saveSavable(invObj,node.addGroup('rebuild'))
results = node.addGroup('results')
preIteration(results)
invObj._invNode = results
self.f.flush()
invObj.hook(_startup_hdf5_inv, overwrite=True)
# At the start of every iteration we will create a inversion iteration node.
def _doStartIteration_hdf5_inv(invObj):
invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1))
preIteration(invObj._invNodeIt)
invObj.hook(_doStartIteration_hdf5_inv, overwrite=True)
def _doEndIteration_hdf5_inv(invObj):
invObj.save(invObj._invNodeIt)
postIteration(invObj._invNodeIt)
self.f.flush()
invObj.hook(_doEndIteration_hdf5_inv, overwrite=True)
# Delete all iterates that did not finish.
def _finish_hdf5_inv(invObj):
postIteration(invObj._invNode)
for it in invObj._invNode:
if not it.attrs['complete']:
del self.f[it.path]
del invObj._invNode
self.f.flush()
invObj.hook(_finish_hdf5_inv, overwrite=True)
def _doStartIteration_hdf5_opt(optObj):
optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter))
preIteration(optObj._optNodeIt)
invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True)
def _doEndIteration_hdf5_opt(optObj, xt):
optObj.save(optObj._optNodeIt)
postIteration(optObj._optNodeIt)
self.f.flush()
invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True)
class hdf5Group(object):
"""Has some low level support for wrapping the native HDF5-Group class"""
def __init__(self, T, groupNode):
self.T = T
# check if you are inputing a hdf5Group rather than a raw node, and act accordingly
if issubclass(groupNode.__class__, hdf5Group):
self.node = groupNode.node
else:
self.node = groupNode
self.childClass = hdf5Group
self.parentClass = hdf5Group
@property
def children(self):
"""Children names in a list
Use obj[name] to return the actual node.
"""
myChildren = [c for c in self.node]
myChildren.sort(key=natural_keys)
return myChildren
@property
def numChildren(self):
"""Returns the len(children)"""
return len(self.children)
@property
def parent(self):
"""Returns parent node"""
return self.parentClass(self.T, self.node.parent)
@property
def name(self):
return self.path.split('/')[-1]
@property
def path(self):
"""Returns the root path of the group"""
return self.node.name
@property
def attrs(self):
"""Returns a list of attributes in the group"""
return self.node.attrs
def addGroup(self, name, soft=False):
"""Adds a child group to the current node."""
if name in self.children and soft:
return self[name]
assert name not in self.children, 'Already a child called: '+self.path+'/'+name
return self.childClass(self.T, self.node.create_group(name))
def setArray(self, name, data):
a = self.node.create_dataset(name, data.shape)
a[...] = data
return a
def __getitem__(self, val):
if type(val) is int:
val = self.children[val]
child = self.node[val]
if type(child) is h5py.Group:
child = self.childClass(self.T, child)
return child
def __contains__(self, key):
return key in self.children
def show(self, pad='', maxDepth=1, depth=0):
"""
Recursively show the structure of the database.
For example::
<hdf5InversionGroup group "/inversions" (1 member)>
- <hdf5Inversion group "/inversions/0" (4 members)>
- <hdf5InversionIteration group "/inversions/0/0.0" (3 members)>
- <hdf5InversionIteration group "/inversions/0/0.1" (3 members)>
- <hdf5InversionIteration group "/inversions/0/0.2" (3 members)>
- <hdf5InversionIteration group "/inversions/0/0.3" (3 members)>
"""
s = self.__str__()
pad += ' '*4
if maxDepth <= 0: print s
if depth >= maxDepth: return s
for c in self.children:
if issubclass(self[c].__class__, hdf5Group):
s += '\n%s- %s' % (pad, self[c].show(pad=pad,depth=depth+1,maxDepth=maxDepth))
else:
s += '\n%s- %s' % (pad, self[c].__str__())
if depth is 0:
print s
else:
return s
def __str__(self):
return '<%s "%s" (%d member%s, attrs=[%s])>' % (self.__class__.__name__, self.path, self.numChildren, '' if self.numChildren == 1 else 's',', '.join([a for a in self.attrs]))
class hdf5InversionGroup(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.childClass = hdf5Inversion
class hdf5Inversion(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5InversionGroup
self.childClass = hdf5InversionResults
def rebuild(self):
return loadSavable(self['rebuild'])
@property
def results(self): return self['results']
class hdf5InversionResults(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5Inversion
self.childClass = hdf5InversionIteration
class hdf5InversionIteration(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5InversionResults
class Savable(type):
def __new__(cls, name, bases, attrs):
__init__ = attrs['__init__']
def init_wrapper(self, *args, **kwargs):
self._args_init = args
self._kwargs_init = kwargs
return __init__(self, *args,**kwargs)
attrs['__init__'] = init_wrapper
newClass = super(Savable, cls).__new__(cls, name, bases, attrs)
SAVEABLES[name] = newClass
return newClass
def saveSavable(obj, group, debug=False):
"""
This creates softlinks if _savable exists in children object.
The first object is always created.
"""
assert type(obj.__class__) is Savable, 'Can only save objects that are Savable objects.'
def doSave(grp, name, val):
if debug: print name, val
if type(val.__class__) is Savable:
link = getattr(val,'_savable',None)
if link is not None:
group.node[name] = h5py.SoftLink(link.path)
if debug: 'Created a softlink path to %s' % link.path
else:
subgrp = grp.addGroup(name)
saveSavable(val, subgrp, debug=debug)
elif type(val) in [list, tuple]:
# Split up, and save each element
for i, v in enumerate(val):
doSave(grp, name + '[%d]'%i, v)
elif type(val) is np.ndarray:
grp.setArray(name, val)
elif val is None:
grp.attrs[name] = 'None'
else:
# just try saving it as an attr
try:
grp.attrs[name] = val
except Exception, e:
print 'Warning: Could not save %s, problems may arise in loading.' % name
group.attrs['__class__'] = obj.__class__.__name__
for arg in obj._kwargs_init:
doSave(group, '_kwarg_'+arg, obj._kwargs_init[arg])
for i, arg in enumerate(obj._args_init):
doSave(group, '_arg%d'%i, arg)
obj._savable = group
def loadSavable(node, pointers=None):
"""
pointers allow things that point to the same node in the h5py file to
be returned as the same object, if they have already been created.
"""
if pointers is None: pointers = []
for pointer in pointers:
if pointer._savable.node == node.node: return pointer
args = ([a for a in node.attrs if '_arg' in a] + [a for a in node.children if '_arg' in a])
kwargs = ([a for a in node.attrs if '_kwarg' in a] + [a for a in node.children if '_kwarg' in a])
args.sort(key=natural_keys)
kwargs.sort(key=natural_keys)
def get(node,key):
if key in node.children: return node[key]
elif key in node.attrs: return node.attrs[key]
ARGS = []
for name in args:
val = get(node, name)
if val.__class__ is h5py.Dataset: val = val[:]
if val is 'None': val = None
if '[' in name: # We are reloading a list
ind = int(name[4:name.index('[')])
if len(ARGS) is ind: # Create the list
ARGS.append([val])
else:
ARGS[ind].append(val)
elif issubclass(val.__class__,hdf5Group):
ARGS.append(loadSavable(val,pointers=pointers))
else:
ind = int(name[4:])
ARGS.append(val)
KWARGS = {}
for name in kwargs:
val = get(node, name)
if val.__class__ is h5py.Dataset: val = val[:]
if val is 'None': val = None
if '[' in name: # We are reloading a list
key = name[7:name.index('[')]
if key not in KWARGS: # Create the list
KWARGS[key] = [val]
else:
KWARGS[key].append(val)
elif issubclass(val.__class__,hdf5Group):
key = name[7:]
KWARGS[key] = loadSavable(val,pointers=pointers)
else:
key = name[7:]
KWARGS[key] = val
cls = get(node, '__class__')
if cls in SAVEABLES:
try:
out = SAVEABLES[cls](*ARGS, **KWARGS)
out._savable = node
pointers.append(out) # Because this is recursive.
return out
except Exception, e:
print 'Warning: %s Class could not be initiated.' % cls
print 'ARGS: ', ARGS
print 'KWARGS: ', KWARGS
return (cls, ARGS, KWARGS, node)
else:
print 'Warning: %s Class not found in SimPEG.Utils.Save.SAVABLES' % cls
return (cls, ARGS, KWARGS, node)
SimPEG/Utils/Solver.py
+384
View File
@@ -0,0 +1,384 @@
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
from matutils import mkvc
from sputils import sdiag
import warnings
DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'}
OPTIONS = {'direct':['scipy'], 'iter':['scipy'], 'triangular':['python'], 'diagonal':['python']}
try:
import TriSolve
OPTIONS['triangular'].append('fortran')
except Exception, e:
print 'Warning: Python backend is being used for solver. Run setup.py from the command line.'
DEFAULTS['triangular'] = 'python'
try:
import mumps
OPTIONS['direct'].append('mumps')
except Exception, e:
print 'Warning: mumps solver not available.'
class Solver(object):
"""
Solver is a light wrapper on the various types of
linear solvers available in python.
:param scipy.sparse A: Matrix
:param bool doDirect: if you want a direct solver
:param string flag: Matrix type flag for special solves: [None, 'L', 'U', 'D']
:param dict options: options which are passed to each sub solver, see each for details.
:rtype: Solver
:return: Solver
To use for direct solvers::
solve = Solver(A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'})
x = solve.solve(rhs)
Or in one line::
x = Solver(A).solve(rhs)
The flag can be set to None, 'L', 'U', or 'D', for general, lower, upper, and diagonal matrices, respectively.
"""
def __init__(self, A, doDirect=True, flag=None, options={}):
assert type(doDirect) is bool, 'doDirect must be a boolean'
assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'"
assert type(options) is dict, 'options must be a dictionary object'
self.A = A
self.dsolve = None
self.doDirect = doDirect
self.flag = flag
self.options = options
if doDirect: return
# Now deal with iterative stuff only
if 'M' not in options:
warnings.warn("You should provide a preconditioner, M.", UserWarning)
return
M = options['M']
if type(M) is sp.linalg.LinearOperator:
return
PreconditionerList = ['J','GS']
if type(M) is str:
assert M in PreconditionerList, "M must be in the known preconditioner list. ['J','GS']"
M = (M,A) # use A as the base for the preconditioner.
if type(M) is tuple:
assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix."
if M[0] is 'J':
Jacobi = sdiag(1.0/M[1].diagonal())
options['M'] = Jacobi
elif M[0] is 'GS':
DD = sdiag(M[1].diagonal())
Uinv = Solver(M[1], flag='U')
Linv = Solver(M[1], flag='L')
def GS(f):
return Uinv.solve(DD*Linv.solve(f))
options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype )
else:
raise Exception('M must be a LinearOperator or a tuple')
def solve(self, b):
"""
Solves the linear system.
.. math::
Ax=b
:param numpy.ndarray b: the right hand side
:rtype: numpy.ndarray
:return: x
"""
if self.flag is None and self.doDirect:
return self.solveDirect(b, **self.options)
elif self.flag is None and not self.doDirect:
return self.solveIter(b, **self.options)
elif self.flag == 'U':
return self.solveBackward(b, **self.options)
elif self.flag == 'L':
return self.solveForward(b, **self.options)
elif self.flag == 'D':
return self.solveDiagonal(b, **self.options)
else:
raise Exception('Unknown flag.')
pass
def clean(self):
"""Cleans up the memory"""
if self.options.has_key('backend'):
if self.options['backend'] == 'mumps':
self.mctx.destroy()
del self.dsolve
self.dsolve = None
def solveDirect(self, b, factorize=False, backend=None):
"""
Use solve instead of this interface.
:param numpy.ndarray b: the right hand side
:param bool factorize: if you want to factorize and store factors
:param str backend: which backend to use. Default is scipy
:rtype: numpy.ndarray
:return: x
"""
if backend is None: backend = DEFAULTS['direct']
assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
if backend == 'scipy':
X = self.solveDirect_scipy(b, factorize)
elif backend == 'mumps':
X = self.solveDirect_mumps(b, factorize)
return X
def solveDirect_scipy(self, b, factorize):
"""
Use solve instead of this interface.
:param numpy.ndarray b: the right hand side
:param bool factorize: if you want to factorize and store factors
:rtype: numpy.ndarray
:return: x
"""
if factorize and self.dsolve is None:
self.A = self.A.tocsc() # for efficiency
self.dsolve = linalg.factorized(self.A)
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
if factorize:
return self.dsolve(b)
else:
return linalg.dsolve.spsolve(self.A, b)
# Multiple RHSs
X = np.empty_like(b)
for i in range(b.shape[1]):
if factorize:
X[:,i] = self.dsolve(b[:,i])
else:
X[:,i] = linalg.dsolve.spsolve(self.A,b[:,i])
return X
def solveDirect_mumps(self, b, factorize):
"""
Use solve instead of this interface.
:param numpy.ndarray b: the right hand side
:param bool factorize: if you want to factorize and store factors
:rtype: numpy.ndarray
:return: x
"""
if factorize and self.dsolve is None:
self.mctx = mumps.DMumpsContext()
self.mctx.set_icntl(14, 60)
# self.mctx.set_silent()
self.mctx.set_centralized_sparse(self.A)
self.mctx.run(job=4)
def mdsolve(rhs):
x = rhs.copy()
self.mctx.set_rhs(x)
self.mctx.run(job=3)
return x
self.dsolve = mdsolve
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
if factorize:
X = self.dsolve(b)
else:
X = mumps.spsolve(self.A, b)
else:
# Multiple RHSs
X = np.empty_like(b)
for i in range(b.shape[1]):
if factorize:
X[:,i] = self.dsolve(b[:,i])
else:
X[:,i] = mumps.spsolve(self.A,b[:,i])
return X
def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50):
if backend is None: backend = DEFAULTS['iter']
algorithms = {'CG':sp.linalg.cg}
assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!"
alg = algorithms[iterSolver]
if len(b.shape) == 1 or b.shape[1] == 1:
x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter)
else:
x = np.empty_like(b)
for i in range(b.shape[1]):
x[:,i], self.info = alg(self.A, b[:,i], M=M, tol=tol, maxiter=maxIter)
return x
def solveBackward(self, b, backend=None):
"""
Use solve instead of this interface.
Perform a backwards solve with upper triangular A in CSR format (best, if not, it will be converted).
:param str backend: which backend to use. Default is python.
:rtype: numpy.ndarray
:return: x
"""
if backend is None: backend = DEFAULTS['triangular']
if backend not in OPTIONS['triangular']:
print 'Warning: %s-backend not being used, %s-default will be used instead.'%(backend,DEFAULTS['triangular'])
backend = DEFAULTS['triangular']
if type(self.A) is not sp.csr.csr_matrix:
self.A = sp.csr_matrix(self.A)
vals = self.A.data
rowptr = self.A.indptr
colind = self.A.indices
if backend == 'fortran':
if len(b.shape) == 1 or b.shape[1] == 1:
x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size, 1)
x = mkvc(x)
else:
x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1])
elif backend == 'python':
x = np.empty_like(b) # empty() is faster than zeros().
for i in reversed(xrange(self.A.shape[0])):
ith_row = vals[rowptr[i] : rowptr[i+1]]
cols = colind[rowptr[i] : rowptr[i+1]]
x_vals = x[cols]
x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0]
return x
def solveForward(self, b, backend=None):
"""
Use solve instead of this interface.
Perform a forward solve with lower triangular A in CSR format (best, if not, it will be converted).
:param str backend: which backend to use. Default is python.
:rtype: numpy.ndarray
:return: x
"""
if backend is None: backend = DEFAULTS['triangular']
if backend not in OPTIONS['triangular']:
print 'Warning: %s-backend not being used, %s-default will be used instead.'%(backend,DEFAULTS['triangular'])
backend = DEFAULTS['triangular']
if type(self.A) is not sp.csr.csr_matrix:
from scipy.sparse import csr_matrix
self.A = csr_matrix(self.A)
vals = self.A.data
rowptr = self.A.indptr
colind = self.A.indices
if backend == 'fortran':
if len(b.shape) == 1 or b.shape[1] == 1:
x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size, 1)
x = mkvc(x)
else:
x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1])
elif backend == 'python':
x = np.empty_like(b) # empty() is faster than zeros().
for i in xrange(self.A.shape[0]):
ith_row = vals[rowptr[i] : rowptr[i+1]]
cols = colind[rowptr[i] : rowptr[i+1]]
x_vals = x[cols]
x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1]
return x
def solveDiagonal(self, b, backend=None):
"""
Use solve instead of this interface.
Perform a diagonal solve with diagonal matrix A.
:param str backend: which backend to use. Default is python.
:rtype: numpy.ndarray
:return: x
"""
if backend is None: backend = DEFAULTS['diagonal']
diagA = self.A.diagonal()
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
return b/diagA
# Multiple RHSs
X = np.empty_like(b)
for i in range(b.shape[1]):
X[:,i] = b[:,i]/diagA
return X
if __name__ == '__main__':
from SimPEG.mesh import TensorMesh
from time import time
h1 = np.ones(20)*100.
h2 = np.ones(20)*100.
h3 = np.ones(20)*100.
h = [h1,h2,h3]
M = TensorMesh(h)
D = M.faceDiv
G = M.cellGrad
Msig = M.getFaceMass()
A = D*Msig*G
A[0,0] *= 10 # remove the constant null space from the matrix
e = np.ones(M.nC)
rhs = A.dot(e)
tic = time()
solve = Solver(A, options={'factorize':True})
x = solve.solve(rhs)
print 'Factorized', time() - tic
print np.linalg.norm(e-x,np.inf)
tic = time()
solve = Solver(A, options={'factorize':False})
x = solve.solve(rhs)
print 'spsolve', time() - tic
print np.linalg.norm(e-x,np.inf)
n = 600
A_dense = np.random.random((n,n))
L = np.tril(np.dot(A_dense, A_dense)) # Positive definite is better conditioned.
e = np.ones(n)
b = np.dot(L, e)
A = sp.csr_matrix(L)
pSolve = Solver(A,flag='L',options={'backend':'python'});
fSolve = Solver(A,flag='L',options={'backend':'fortran'})
tic = time()
x = pSolve.solve(b)
toc = time() - tic
print 'Error Forward Python = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc
tic = time()
x = fSolve.solve(b)
toc = time() - tic
print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc
A = -D*D.T
A[0,0] *= 10 # remove the constant null space from the matrix
e = np.ones(M.nC)
b = A.dot(e)
iSolve = Solver(A, doDirect=False,options={'M':('GS',A)})
tic = time()
x = iSolve.solve(b)
toc = time() - tic
print x
print 'Error CG = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info
SimPEG/Utils/TriSolve.f
+64
View File
@@ -0,0 +1,64 @@
c File TriSolve.f
subroutine forward(al, ial, jal, b, nv, n, nRHS, x)
double precision al(nv)
integer ial(n+1)
integer jal(nv)
double precision b(n,nRHS)
double precision x(n,nRHS)
integer nv
integer n
integer nRHS
integer rhs
cf2py intent(in) :: al
cf2py intent(in) :: ial
cf2py intent(in) :: jal
cf2py intent(in) :: b
cf2py intent(in) :: nv
cf2py intent(in) :: n
cf2py intent(in) :: nRHS
cf2py intent(out) :: x
real ( kind = 8 ) t
do rhs = 1, nRHS
do k = 1, n
t = b(k,rhs)
do j = ial(k)+1, ial(k+1)
t = t - al(j) * x(jal(j)+1,rhs)
end do
x(k,rhs) = t/al(ial(k+1))
end do
end do
end subroutine forward
subroutine backward(au,iau, jau, b, nv, n, nRHS, x)
double precision au(nv)
integer iau(n+1)
integer jau(nv)
double precision b(n,nRHS)
double precision x(n,nRHS)
integer nv
integer n
integer nRHS
integer rhs
cf2py intent(in) :: au
cf2py intent(in) :: iau
cf2py intent(in) :: jau
cf2py intent(in) :: b
cf2py intent(in) :: nv
cf2py intent(in) :: n
cf2py intent(in) :: nRHS
cf2py intent(out) :: x
real ( kind = 8 ) t
do rhs = 1, nRHS
do k = n, 1, -1
t = b(k,rhs)
do j = iau(k)+1, iau(k+1)
t = t - au(j) * x(jau(j)+1,rhs)
end do
x(k,rhs) = t/au(iau(k)+1)
end do
end do
end subroutine backward
SimPEG/Utils/__init__.py
+250
View File
@@ -0,0 +1,250 @@
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap
from meshutils import exampleLomGird, meshTensors
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube
from interputils import interpmat
from ipythonutils import easyAnimate as animate
from Solver import Solver
import Save
import Geophysics
import types
import time
import numpy as np
from functools import wraps
def hook(obj, method, name=None, overwrite=False, silent=False):
"""
This dynamically binds a method to the instance of the class.
If name is None, the name of the method is used.
"""
if name is None:
name = method.__name__
if name == '<lambda>':
raise Exception('Must provide name to hook lambda functions.')
if not hasattr(obj,name) or overwrite:
setattr(obj, name, types.MethodType( method, obj ))
if getattr(obj,'debug',False):
print 'Method '+name+' was added to class.'
elif not silent or getattr(obj,'debug',False):
print 'Method '+name+' was not overwritten.'
def setKwargs(obj, **kwargs):
"""Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
for attr in kwargs:
if hasattr(obj, attr):
setattr(obj, attr, kwargs[attr])
else:
raise Exception('%s attr is not recognized' % attr)
hook(obj,hook, silent=True)
hook(obj,setKwargs, silent=True)
def printTitles(obj, printers, name='Print Titles', pad=''):
titles = ''
widths = 0
for printer in printers:
titles += ('{:^%i}'%printer['width']).format(printer['title']) + ''
widths += printer['width']
print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name)
print pad + titles
print pad + "%s" % '-'*widths
def printLine(obj, printers, pad=''):
values = ''
for printer in printers:
values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj))
print pad + values
def checkStoppers(obj, stoppers):
# check stopping rules
optimal = []
critical = []
for stopper in stoppers:
l = stopper['left'](obj)
r = stopper['right'](obj)
if stopper['stopType'] == 'optimal':
optimal.append(l <= r)
if stopper['stopType'] == 'critical':
critical.append(l <= r)
if obj.debug: print 'checkStoppers.optimal: ', optimal
if obj.debug: print 'checkStoppers.critical: ', critical
return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical))
def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
print pad + "%s%s%s" % ('-'*25,stop,'-'*25)
for stopper in stoppers:
l = stopper['left'](obj)
r = stopper['right'](obj)
print pad + stopper['str'] % (l<=r,l,r)
print pad + "%s%s%s" % ('-'*25,done,'-'*25)
def callHooks(match, mainFirst=False):
"""
Use this to wrap a funciton::
@callHooks('doEndIteration')
def doEndIteration(self):
pass
This will call everything named _doEndIteration* at the beginning of the function call.
By default the master method (doEndIteration) is run after all of the sub methods (_doEndIteration*).
This can be reversed by adding the mainFirst=True kwarg.
"""
def callHooksWrap(f):
@wraps(f)
def wrapper(self,*args,**kwargs):
if not mainFirst:
for method in [posible for posible in dir(self) if ('_'+match) in posible]:
if getattr(self,'debug',False): print (match+' is calling self.'+method)
getattr(self,method)(*args, **kwargs)
return f(self,*args,**kwargs)
else:
out = f(self,*args,**kwargs)
for method in [posible for posible in dir(self) if ('_'+match) in posible]:
if getattr(self,'debug',False): print (match+' is calling self.'+method)
getattr(self,method)(*args, **kwargs)
return out
extra = """
If you have things that also need to run in the method %s, you can create a method::
def _%s*(self, ... ):
pass
Where the * can be any string. If present, _%s* will be called at the start of the default %s call.
You may also completely overwrite this function.
""" % (match, match, match, match)
doc = wrapper.__doc__
wrapper.__doc__ = ('' if doc is None else doc) + extra
return wrapper
return callHooksWrap
def dependentProperty(name, value, children, doc):
def fget(self): return getattr(self,name,value)
def fset(self, val):
for child in children:
if hasattr(self, child):
delattr(self, child)
setattr(self, name, val)
return property(fget=fget, fset=fset, doc=doc)
class Counter(object):
"""
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()
"""
def __init__(self):
self._countList = {}
self._timeList = {}
def count(self, prop):
"""
Increases the count of the property.
"""
assert type(prop) is str, 'The property must be a string.'
if prop not in self._countList:
self._countList[prop] = 0
self._countList[prop] += 1
def countTic(self, prop):
"""
Times a property call, this is the init call.
"""
assert type(prop) is str, 'The property must be a string.'
if prop not in self._timeList:
self._timeList[prop] = []
self._timeList[prop].append(-time.time())
def countToc(self, prop):
"""
Times a property call, this is the end call.
"""
assert type(prop) is str, 'The property must be a string.'
assert prop in self._timeList, 'The property must already be in the dictionary.'
self._timeList[prop][-1] += time.time()
def summary(self):
"""
Provides a text summary of the current counters and timers.
"""
print 'Counters:'
for prop in sorted(self._countList):
print " {0:<40}: {1:8d}".format(prop,self._countList[prop])
print '\nTimes:'+' '*40+'mean sum'
for prop in sorted(self._timeList):
l = len(self._timeList[prop])
a = np.array(self._timeList[prop])
print " {0:<40}: {1:4.2e}, {2:4.2e}, {3:4d}x".format(prop,a.mean(),a.sum(),l)
def count(f):
@wraps(f)
def wrapper(self,*args,**kwargs):
counter = getattr(self,'counter',None)
if type(counter) is Counter: counter.count(self.__class__.__name__+'.'+f.__name__)
out = f(self,*args,**kwargs)
return out
return wrapper
def timeIt(f):
@wraps(f)
def wrapper(self,*args,**kwargs):
counter = getattr(self,'counter',None)
if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
out = f(self,*args,**kwargs)
if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
return out
return wrapper
if __name__ == '__main__':
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()
SimPEG/Utils/interputils.py
+181
View File
@@ -0,0 +1,181 @@
import numpy as np
import scipy.sparse as sp
from sputils import spzeros
from matutils import mkvc, sub2ind
def _interp_point_1D(x, xr_i):
"""
given a point, xr_i, this will find which two integers it lies between.
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param float xr_i: Location of a point
:rtype: int,int,float,float
:return: index1, index2, portion1, portion2
"""
# TODO: This fails if the point is on the outside of the mesh. We may want to replace this by extrapolation?
im = np.argmin(abs(x-xr_i))
if xr_i - x[im] >= 0: # Point on the left
ind_x1 = im
ind_x2 = im+1
elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
dx1 = xr_i - x[ind_x1]
dx2 = x[ind_x2] - xr_i
return ind_x1, ind_x2, dx1, dx2
def interpmat(locs, x, y=None, z=None):
"""
Local interpolation computed for each receiver point in turn
:param numpy.ndarray loc: Location of points to interpolate to
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param numpy.ndarray y: Tensor vector of 2nd dimension of grid. None by default.
:param numpy.ndarray z: Tensor vector of 3rd dimension of grid. None by default.
:rtype: scipy.sparse.csr.csr_matrix
:return: Interpolation matrix
.. plot::
import SimPEG
import numpy as np
import matplotlib.pyplot as plt
locs = np.random.rand(50)*0.8+0.1
x = np.linspace(0,1,7)
dense = np.linspace(0,1,200)
fun = lambda x: np.cos(2*np.pi*x)
Q = SimPEG.Utils.interpmat(locs, x)
plt.plot(x, fun(x), 'bs-')
plt.plot(dense, fun(dense), 'y:')
plt.plot(locs, Q*fun(x), 'mo')
plt.plot(locs, fun(locs), 'rx')
plt.show()
"""
if y is None and z is None:
return _interpmat1D(locs, x)
elif z is None:
return _interpmat2D(locs, x, y)
else:
return _interpmat3D(locs, x, y, z)
def _interpmat1D(locs, x):
"""Use interpmat with only x component provided."""
nx = x.size
locs = mkvc(locs)
npts = locs.shape[0]
Q = sp.lil_matrix((npts, nx))
for i in range(npts):
ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i])
dv = (x[ind_x2] - x[ind_x1])
Dx = x[ind_x2] - x[ind_x1]
# Get the row in the matrix
inds = [ind_x1, ind_x2]
vals = [(1-dx1/Dx),(1-dx2/Dx)]
Q[i, inds] = vals
return Q.tocsr()
def _interpmat2D(locs, x, y):
"""Use interpmat with only x and y components provided."""
nx = x.size
ny = y.size
npts = locs.shape[0]
Q = sp.lil_matrix((npts, nx*ny))
for i in range(npts):
ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1])
dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1])
Dx = x[ind_x2] - x[ind_x1]
Dy = y[ind_y2] - y[ind_y1]
# Get the row in the matrix
inds = sub2ind((nx,ny),[
( ind_x1, ind_y2),
( ind_x1, ind_y1),
( ind_x2, ind_y1),
( ind_x2, ind_y2)])
vals = [(1-dx1/Dx)*(1-dy2/Dy),
(1-dx1/Dx)*(1-dy1/Dy),
(1-dx2/Dx)*(1-dy1/Dy),
(1-dx2/Dx)*(1-dy2/Dy)]
Q[i, mkvc(inds)] = vals
return Q.tocsr()
def _interpmat3D(locs, x, y, z):
"""Use interpmat."""
nx = x.size
ny = y.size
nz = z.size
npts = locs.shape[0]
Q = sp.lil_matrix((npts, nx*ny*nz))
for i in range(npts):
ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, dz1, dz2 = _interp_point_1D(z, locs[i, 2])
dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
Dx = x[ind_x2] - x[ind_x1]
Dy = y[ind_y2] - y[ind_y1]
Dz = z[ind_z2] - z[ind_z1]
# Get the row in the matrix
inds = sub2ind((nx,ny,nz),[
( ind_x1, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
( ind_x1, ind_y2, ind_z2),
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)])
vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz),
(1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz),
(1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz),
(1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz),
(1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz),
(1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz),
(1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz),
(1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)]
Q[i, mkvc(inds)] = vals
return Q.tocsr()
if __name__ == '__main__':
import SimPEG
import numpy as np
import matplotlib.pyplot as plt
locs = np.random.rand(50)*0.8+0.1
x = np.linspace(0,1,7)
dense = np.linspace(0,1,200)
fun = lambda x: np.cos(2*np.pi*x)
Q = SimPEG.Utils.interpmat(locs, x)
plt.plot(x, fun(x), 'bs-')
plt.plot(dense, fun(dense), 'y:')
plt.plot(locs, Q*fun(x), 'mo')
plt.plot(locs, fun(locs), 'rx')
plt.show()
SimPEG/Utils/ipythonutils.py
+28
View File
@@ -0,0 +1,28 @@
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/lomutils.py
+265
View File
@@ -0,0 +1,265 @@
import numpy as np
from scipy import sparse as sp
from matutils import mkvc, ndgrid, sub2ind
from sputils import sdiag
def volTetra(xyz, A, B, C, D):
"""
Returns the volume for tetrahedras volume specified by the indexes A to D.
:param numpy.array xyz: X,Y,Z vertex vector
:param numpy.array A,B,C,D: vert index of the tetrahedra
:rtype: numpy.array
:return: V, volume of the tetrahedra
Algorithm http://en.wikipedia.org/wiki/Tetrahedron#Volume
.. math::
V = {1 \over 3} A h
V = {1 \over 6} | ( a - d ) \cdot ( ( b - d ) ( c - d ) ) |
"""
AD = xyz[A, :] - xyz[D, :]
BD = xyz[B, :] - xyz[D, :]
CD = xyz[C, :] - xyz[D, :]
V = (BD[:, 0]*CD[:, 1] - BD[:, 1]*CD[:, 0])*AD[:, 2] - (BD[:, 0]*CD[:, 2] - BD[:, 2]*CD[:, 0])*AD[:, 1] + (BD[:, 1]*CD[:, 2] - BD[:, 2]*CD[:, 1])*AD[:, 0]
return V/6
def indexCube(nodes, gridSize, n=None):
"""
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)
"""
assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'"
assert type(gridSize) == np.ndarray, "Number of nodes must be an ndarray"
nodes = nodes.upper()
# Make sure that we choose from the possible nodes.
possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH'
for node in nodes:
assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes
dim = gridSize.size
if n is None:
n = gridSize - 1
if dim == 2:
ij = ndgrid(np.arange(n[0]), np.arange(n[1]))
i, j = ij[:, 0], ij[:, 1]
elif dim == 3:
ijk = ndgrid(np.arange(n[0]), np.arange(n[1]), np.arange(n[2]))
i, j, k = ijk[:, 0], ijk[:, 1], ijk[:, 2]
else:
raise Exception('Only 2 and 3 dimensions supported.')
nodeMap = {'A': [0, 0, 0], 'B': [0, 1, 0], 'C': [1, 1, 0], 'D': [1, 0, 0],
'E': [0, 0, 1], 'F': [0, 1, 1], 'G': [1, 1, 1], 'H': [1, 0, 1]}
out = ()
for node in nodes:
shift = nodeMap[node]
if dim == 2:
out += (sub2ind(gridSize, np.c_[i+shift[0], j+shift[1]]).flatten(), )
elif dim == 3:
out += (sub2ind(gridSize, np.c_[i+shift[0], j+shift[1], k+shift[2]]).flatten(), )
return out
def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True):
"""
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
"""
assert type(average) is bool, 'average must be a boolean'
assert type(normalizeNormals) is bool, 'normalizeNormals must be a boolean'
# compute normal that is pointing away from you.
#
# A -------A-B------- B
# | |
# | |
# D-A (X) B-C
# | |
# | |
# D -------C-D------- C
AB = xyz[B, :] - xyz[A, :]
BC = xyz[C, :] - xyz[B, :]
CD = xyz[D, :] - xyz[C, :]
DA = xyz[A, :] - xyz[D, :]
def cross(X, Y):
return np.c_[X[:, 1]*Y[:, 2] - X[:, 2]*Y[:, 1],
X[:, 2]*Y[:, 0] - X[:, 0]*Y[:, 2],
X[:, 0]*Y[:, 1] - X[:, 1]*Y[:, 0]]
nA = cross(AB, DA)
nB = cross(BC, AB)
nC = cross(CD, BC)
nD = cross(DA, CD)
length = lambda x: np.sqrt(x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)
normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(x), 2))
if average:
# average the normals at each vertex.
N = (nA + nB + nC + nD)/4 # this is intrinsically weighted by area
# normalize
N = normalize(N)
else:
if normalizeNormals:
N = [normalize(nA), normalize(nB), normalize(nC), normalize(nD)]
else:
N = [nA, nB, nC, nD]
# Area calculation
#
# Approximate by 4 different triangles, and divide by 2.
# Each triangle is one half of the length of the cross product
#
# So also could be viewed as the average parallelogram.
#
# TODO: This does not compute correctly for concave quadrilaterals
area = (length(nA)+length(nB)+length(nC)+length(nD))/4
return N, area
def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
""" 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
"""
a11 = mkvc(a11)
a12 = mkvc(a12)
a13 = mkvc(a13)
a21 = mkvc(a21)
a22 = mkvc(a22)
a23 = mkvc(a23)
a31 = mkvc(a31)
a32 = mkvc(a32)
a33 = mkvc(a33)
detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
b11 = +(a22*a33 - a23*a32)/detA
b12 = -(a12*a33 - a13*a32)/detA
b13 = +(a12*a23 - a13*a22)/detA
b21 = +(a31*a23 - a21*a33)/detA
b22 = -(a31*a13 - a11*a33)/detA
b23 = +(a21*a13 - a11*a23)/detA
b31 = -(a31*a22 - a21*a32)/detA
b32 = +(a31*a12 - a11*a32)/detA
b33 = -(a21*a12 - a11*a22)/detA
B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12), sdiag(b13))),
sp.hstack((sdiag(b21), sdiag(b22), sdiag(b23))),
sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33)))))
return B
def inv2X2BlockDiagonal(a11, a12, a21, a22):
""" 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, a13, a21, a22, a23, a31, a32, a33
Output:
B - inverse
"""
a11 = mkvc(a11)
a12 = mkvc(a12)
a21 = mkvc(a21)
a22 = mkvc(a22)
# compute inverse of the determinant.
detAinv = 1./(a11*a22 - a21*a12)
b11 = +detAinv*a22
b12 = -detAinv*a12
b21 = -detAinv*a21
b22 = +detAinv*a11
B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
sp.hstack((sdiag(b21), sdiag(b22)))))
return B
SimPEG/Utils/matutils.py
+140
View File
@@ -0,0 +1,140 @@
import numpy as np
def mkvc(x, numDims=1):
"""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)
"""
if type(x) == np.matrix:
x = np.array(x)
assert type(x) == np.ndarray, "Vector must be a numpy array"
if numDims == 1:
return x.flatten(order='F')
elif numDims == 2:
return x.flatten(order='F')[:, np.newaxis]
elif numDims == 3:
return x.flatten(order='F')[:, np.newaxis, np.newaxis]
def ndgrid(*args, **kwargs):
"""
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]]
"""
# Read the keyword arguments, and only accept a vector=True/False
vector = kwargs.pop('vector', True)
assert type(vector) == bool, "'vector' keyword must be a bool"
assert len(kwargs) == 0, "Only 'vector' keyword accepted"
# you can either pass a list [x1, x2, x3] or each seperately
if type(args[0]) == list:
xin = args[0]
else:
xin = args
# Each vector needs to be a numpy array
assert np.all([type(x) == np.ndarray for x in xin]), "All vectors must be numpy arrays."
if len(xin) == 1:
return xin[0]
elif len(xin) == 2:
XY = np.broadcast_arrays(mkvc(xin[1], 1), mkvc(xin[0], 2))
if vector:
X2, X1 = [mkvc(x) for x in XY]
return np.c_[X1, X2]
else:
return XY[1], XY[0]
elif len(xin) == 3:
XYZ = np.broadcast_arrays(mkvc(xin[2], 1), mkvc(xin[1], 2), mkvc(xin[0], 3))
if vector:
X3, X2, X1 = [mkvc(x) for x in XYZ]
return np.c_[X1, X2, X3]
else:
return XYZ[2], XYZ[1], XYZ[0]
def ind2sub(shape, ind):
"""From the given shape, returns the subscrips of the given index"""
revshp = []
revshp.extend(shape)
mult = [1]
for i in range(0, len(revshp)-1):
mult.extend([mult[i]*revshp[i]])
mult = np.array(mult).reshape(len(mult))
sub = []
for i in range(0, len(shape)):
sub.extend([np.math.floor(ind / mult[i])])
ind = ind - (np.math.floor(ind/mult[i]) * mult[i])
return sub
def sub2ind(shape, subs):
"""From the given shape, returns the index of the given subscript"""
revshp = list(shape)
mult = [1]
for i in range(0, len(revshp)-1):
mult.extend([mult[i]*revshp[i]])
mult = np.array(mult).reshape(len(mult), 1)
idx = np.dot((subs), (mult))
return idx
def getSubArray(A, ind):
"""subArray"""
assert type(ind) == list, "ind must be a list of vectors"
assert len(A.shape) == len(ind), "ind must have the same length as the dimension of A"
if len(A.shape) == 2:
return A[ind[0], :][:, ind[1]]
elif len(A.shape) == 3:
return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]
else:
raise Exception("getSubArray does not support dimension asked.")
SimPEG/Utils/meshutils.py
+58
View File
@@ -0,0 +1,58 @@
import numpy as np
from scipy import sparse as sp
from matutils import mkvc, ndgrid, sub2ind
from sputils import sdiag
def exampleLomGird(nC, exType):
assert type(nC) == list, "nC must be a list containing the number of nodes"
assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions"
exType = exType.lower()
possibleTypes = ['rect', 'rotate']
assert exType in possibleTypes, "Not a possible example type."
if exType == 'rect':
return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
elif exType == 'rotate':
if len(nC) == 2:
X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2)
amt[amt < 0] = 0
return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt
elif len(nC) == 3:
X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2)
amt[amt < 0] = 0
return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt
def meshTensors(*args):
"""
**meshTensors** takes any number of tuples that have the form::
h1 = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCore), (numPad, sizeStart [, increaseFactor]) )
.. plot::
from SimPEG import mesh, Utils
M = mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
M.plotGrid()
"""
def padding(num, start, factor=1.3, reverse=False):
pad = ((np.ones(num)*factor)**np.arange(num))*start
if reverse: pad = pad[::-1]
return pad
tensors = tuple()
for i, arg in enumerate(args):
tensors += (np.r_[padding(*arg[0],reverse=True),np.ones(arg[1][0])*arg[1][1],padding(*arg[2])],)
return list(tensors) if len(tensors) > 1 else tensors[0]
if __name__ == '__main__':
from SimPEG import mesh
import matplotlib.pyplot as plt
M = mesh.TensorMesh(meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
M.plotGrid()
plt.gca().axis('tight')
plt.show()
SimPEG/Utils/sputils.py
+38
View File
@@ -0,0 +1,38 @@
from scipy import sparse as sp
from matutils import mkvc
import numpy as np
def sdiag(h):
"""Sparse diagonal matrix"""
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
def speye(n):
"""Sparse identity"""
return sp.identity(n, format="csr")
def kron3(A, B, C):
"""Three kron prods"""
return sp.kron(sp.kron(A, B), C, format="csr")
def spzeros(n1, n2):
"""spzeros"""
return sp.coo_matrix((n1, n2)).tocsr()
def ddx(n):
"""Define 1D derivatives, inner, this means we go from n+1 to n"""
return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
def av(n):
"""Define 1D averaging operator from nodes to cell-centers."""
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
def avExtrap(n):
"""Define 1D averaging operator from cell-centers to nodes."""
Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
return Av