From 102d5721c56f27374d74817603939c09b43f9bde Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 16 May 2014 19:33:29 -0700 Subject: [PATCH] Solver updates. Remove .solve() replace with a multiplication. --- SimPEG/Mesh/TreeMesh.py | 2 +- SimPEG/Mesh/View.py | 4 +-- SimPEG/Optimization.py | 10 +++--- SimPEG/Solver.py | 6 ---- SimPEG/Tests/test_boundaryPoisson.py | 8 ++--- SimPEG/Tests/test_tensorMesh.py | 2 +- SimPEG/Utils/SolverUtils.py | 50 +++++++++++++++++++--------- SimPEG/__init__.py | 2 +- docs/api_Solver.rst | 42 +++++++++++++++++++++-- 9 files changed, 87 insertions(+), 39 deletions(-) delete mode 100644 SimPEG/Solver.py diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 8ac841d4..0ddd337b 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1127,7 +1127,7 @@ if __name__ == '__main__': q = np.zeros(M.nC) q[208] = -1.0 q[291] = 1.0 - b = Solver(-DIV*Mf*DIV.T).solve(q) + b = Solver(-DIV*Mf*DIV.T) * (q) plt.figure() M.plotImage(b) # plt.gca().invert_yaxis() diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index cc1c7279..eed8d06a 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -186,7 +186,7 @@ class TensorView(object): q[[4,4],[4,4],[2,6]]=[-1,1] q = Utils.mkvc(q) A = M.faceDiv*M.cellGrad - b = Solver(A).solve(q) + b = Solver(A) * (q) M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8}) """ @@ -628,7 +628,7 @@ if __name__ == '__main__': q[[4,4],[4,4],[2,6]]=[-1,1] q = Utils.mkvc(q) A = M.faceDiv*M.cellGrad - b = Solver(A).solve(q) + b = Solver(A) * (q) M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, pcolorOpts={'alpha':0.8}) M2 = Mesh.TensorMesh([10,20],x0=[10,5]) diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 12dc385e..c8314af5 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -1,5 +1,5 @@ import Utils, numpy as np, scipy.sparse as sp -from Solver import Solver, SolverCG +from Utils.SolverUtils import Solver, SolverCG norm = np.linalg.norm @@ -692,7 +692,7 @@ class BFGS(Minimize, Remember): def bfgsrec(self,k,n,nn,S,Y,d): """BFGS recursion""" if k < 0: - d = self.bfgsH0.solve(d) + d = self.bfgsH0 * (d) else: khat = 0 if nn is 0 else np.mod(n-nn+k,nn) gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]) @@ -731,7 +731,7 @@ class GaussNewton(Minimize, Remember): @Utils.timeIt def findSearchDirection(self): - return Solver(self.H).solve(-self.g) + return Solver(self.H) * (-self.g) class InexactGaussNewton(BFGS, Minimize, Remember): @@ -779,7 +779,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember): @Utils.timeIt def findSearchDirection(self): Hinv = SolverCG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG) - p = Hinv.solve(-self.g) + p = Hinv * (-self.g) return p @@ -846,7 +846,7 @@ class NewtonRoot(object): r, J = fun(x, return_g=True) Jinv = self.Solver(J, **self.solverOpts) - dh = - Jinv.solve(r) + dh = - (Jinv * r) muLS = 1. LScnt = 1 diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py deleted file mode 100644 index 8f96923c..00000000 --- a/SimPEG/Solver.py +++ /dev/null @@ -1,6 +0,0 @@ -import scipy.sparse as sp -from SimPEG.Utils import SolverUtils - -Solver = SolverUtils.DSolverWrap(sp.linalg.spsolve, factorize=False) -SolverLU = SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True) -SolverCG = SolverUtils.ISolverWrap(sp.linalg.cg) diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index d4ac504b..e604c35b 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -55,12 +55,12 @@ class Test1D_InhomogeneousDirichlet(OrderTest): elif self.myTest == 'xc': #TODO: fix the null space solver = SolverCG(A, maxiter=1000) - xc = solver.solve(rhs) + xc = solver * (rhs) print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) err = np.linalg.norm((xc-xc_anal), np.inf) elif self.myTest == 'xcJ': #TODO: fix the null space - xc = Solver(A).solve(rhs) + xc = Solver(A) * (rhs) print np.linalg.norm(Utils.mkvc(A*xc) - rhs) j = McI*(G*xc + P*phi_bc) err = np.linalg.norm((j-j_anal), np.inf) @@ -140,10 +140,10 @@ class Test2D_InhomogeneousDirichlet(OrderTest): elif self.myTest == 'q': err = np.linalg.norm((q-q_anal), np.inf) elif self.myTest == 'xc': - xc = Solver(A).solve(rhs) + xc = Solver(A) * (rhs) err = np.linalg.norm((xc-xc_anal), np.inf) elif self.myTest == 'xcJ': - xc = Solver(A).solve(rhs) + xc = Solver(A) * (rhs) j = McI*(G*xc + P*bc) err = np.linalg.norm((j-j_anal), np.inf) diff --git a/SimPEG/Tests/test_tensorMesh.py b/SimPEG/Tests/test_tensorMesh.py index 23658739..3e3286e5 100644 --- a/SimPEG/Tests/test_tensorMesh.py +++ b/SimPEG/Tests/test_tensorMesh.py @@ -79,7 +79,7 @@ class TestPoissonEqn(OrderTest): err = np.linalg.norm((sA - sN), np.inf) else: fA = fun(self.M.gridCC) - fN = Solver(D*G).solve(sol(self.M.gridCC)) + fN = Solver(D*G) * (sol(self.M.gridCC)) err = np.linalg.norm((fA - fN), np.inf) return err diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index c2f17b6d..50317b63 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -1,4 +1,4 @@ -import numpy as np +import numpy as np, scipy.sparse as sp from matutils import mkvc import warnings @@ -13,7 +13,16 @@ def _checkAccuracy(A, b, X, accuracyTol): warnings.warn(msg, RuntimeWarning) -def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): +def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): + """ + Wraps a direct Solver. + + :: + + Solver = SolverUtils.SolverWrapD(sp.linalg.spsolve, factorize=False) + SolverLU = SolverUtils.SolverWrapD(sp.linalg.splu, factorize=True) + + """ def __init__(self, A, **kwargs): self.A = A.tocsc() @@ -21,7 +30,10 @@ def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): if factorize: self.solver = fun(self.A, **kwargs) - def solve(self, b): + def __mul__(self, b): + if type(b) is not np.ndarray: + raise TypeError('Can only multiply by a numpy array.') + if len(b.shape) == 1 or b.shape[1] == 1: b = b.flatten() # Just one RHS @@ -45,22 +57,28 @@ def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): if factorize and hasattr(self.solver, 'clean'): return self.solver.clean() - def __mul__(self, val): - if type(val) is np.ndarray: - return self.solve(val) - raise TypeError('Can only multiply by a numpy array.') - - return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "clean": clean, "__mul__": __mul__}) + return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__}) -def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5): +def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5): + """ + Wraps an iterative Solver. + + :: + + SolverCG = SolverUtils.SolverWrapI(sp.linalg.cg) + + """ def __init__(self, A, **kwargs): self.A = A self.kwargs = kwargs - def solve(self, b): + def __mul__(self, b): + if type(b) is not np.ndarray: + raise TypeError('Can only multiply by a numpy array.') + if len(b.shape) == 1 or b.shape[1] == 1: b = b.flatten() # Just one RHS @@ -90,9 +108,9 @@ def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5): if hasattr(self.solver, 'clean'): return self.solver.clean() - def __mul__(self, val): - if type(val) is np.ndarray: - return self.solve(val) - raise TypeError('Can only multiply by a numpy array.') + return type(fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__}) - return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "clean": clean, "__mul__": __mul__}) + +Solver = SolverWrapD(sp.linalg.spsolve, factorize=False) +SolverLU = SolverWrapD(sp.linalg.splu, factorize=True) +SolverCG = SolverWrapI(sp.linalg.cg) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 299bd01e..7d894e65 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,7 +1,7 @@ import numpy as np import scipy.sparse as sp import Utils -from Solver import * +from Utils.SolverUtils import * import Mesh import Maps import Problem diff --git a/docs/api_Solver.rst b/docs/api_Solver.rst index be009f97..9909b65e 100644 --- a/docs/api_Solver.rst +++ b/docs/api_Solver.rst @@ -4,7 +4,43 @@ Solver ****** -.. automodule:: SimPEG.Solver - :members: - :undoc-members: +The numerical linear algebra solver that you use will ultimately be the +bottleneck of your large scale inversion. To be the most flexible, SimPEG +provides wrappers rather than a comprehensive set of solvers (i.e. BYOS). + +The interface is as follows:: + + A # Where A is a sparse matrix (or linear operator) + Ainv = Solver(A, **solverOpts) # Create a solver object with key word arguments + x = Ainv * b # Where b is a numpy array of shape (n,) or (n,*) + Ainv.clean() # This cleans the memory footprint (if any) + +.. note:: + + This is somewhat an abuse of notation for solvers as we never actually + create A inverse. Instead we are creating an object that acts like A + inverse, whether that be a Krylov subspace solver or an LU decomposition. + +To wrap up solvers in scipy.sparse.linalg it takes one line of code:: + + + Solver = SolverWrapD(sp.linalg.spsolve, factorize=False) + SolverLU = SolverWrapD(sp.linalg.splu, factorize=True) + SolverCG = SolverWrapI(sp.linalg.cg) + +.. note:: + + The above solvers are loaded into the base name space of SimPEG. + +.. seealso:: + + https://github.com/rowanc1/pymatsolver + + +The API +======= + +.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD + +.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI