wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-09-12 12:51:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of https://bitbucket.org/rcockett/simpeg into seismicExample

Browse Source 
Conflicts:
	docs/api_TestResults.rst
This commit is contained in:
Rowan Cockett
2013-11-12 13:43:06 -08:00
parent 291cabe931 cea1141dcd
commit 5248ac81ce
20 changed files with 2004 additions and 913 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/utils/Solver.py
+71 -19
View File
@@ -1,7 +1,22 @@
import numpy as np
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
from SimPEG.utils import mkvc
DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
try:
import TriSolve
except Exception, e:
try:
import os
# Note: this may not work from SublimeText, if that is the case, just run the command in your shell.
os.system('f2py -c TriSolve.f -m TriSolve')
import TriSolve
except Exception, e:
print 'Warning: Python backend is being used for solver. Run setup.py from the command line.'
DEFAULTS['forward'] = 'python'
DEFAULTS['backward'] = 'python'
class Solver(object):
"""
@@ -55,11 +70,11 @@ class Solver(object):
elif self.flag is None and not self.doDirect:
return self.solveIter(b, **self.options)
elif self.flag == 'U':
return self.solveBackward(b)
return self.solveBackward(b, **self.options)
elif self.flag == 'L':
return self.solveForward(b)
return self.solveForward(b, **self.options)
elif self.flag == 'D':
return self.solveDiagonal(b)
return self.solveDiagonal(b, **self.options)
else:
raise Exception('Unknown flag.')
pass
@@ -69,7 +84,7 @@ class Solver(object):
del self.dsolve
self.dsolve = None
def solveDirect(self, b, factorize=False, backend='scipy'):
def solveDirect(self, b, factorize=False, backend=None):
"""
Use solve instead of this interface.
@@ -78,6 +93,8 @@ class Solver(object):
: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 factorize and self.dsolve is None:
@@ -104,7 +121,7 @@ class Solver(object):
def solveIter(self, b, M=None, iterSolver='CG'):
pass
def solveBackward(self, b, backend='python'):
def solveBackward(self, b, backend=None):
"""
Use solve instead of this interface.
@@ -114,21 +131,29 @@ class Solver(object):
:rtype: numpy.ndarray
:return: x
"""
if backend is None: backend = DEFAULTS['backward']
if type(self.A) is not sparse.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
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]
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='python'):
def solveForward(self, b, backend=None):
"""
Use solve instead of this interface.
@@ -138,21 +163,29 @@ class Solver(object):
:rtype: numpy.ndarray
:return: x
"""
if backend is None: backend = DEFAULTS['forward']
if type(self.A) is not sparse.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
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]
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='python'):
def solveDiagonal(self, b, backend=None):
"""
Use solve instead of this interface.
@@ -162,6 +195,8 @@ class Solver(object):
: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
@@ -205,3 +240,20 @@ if __name__ == '__main__':
print np.linalg.norm(e-x,np.inf)
n = 6000
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 = sparse.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
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
+52 -2
View File
@@ -3,9 +3,59 @@ import sputils
import lomutils
import interputils
import ModelBuilder
import Solver
from Solver import Solver
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag, ddx, av
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
from interputils import interpmat
from ipythonUtils import easyAnimate as animate
import Solver
from Solver import Solver
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)
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)
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