From 025cb5c58355d07586a4cb44945092a43b83ee2a Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 14 Mar 2014 15:19:38 -0700 Subject: [PATCH 1/8] Removed hold off for consistent behaviour with the rest of matplotlib. --- SimPEG/Mesh/TensorView.py | 1 - 1 file changed, 1 deletion(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index a67d4609..074065dd 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -471,7 +471,6 @@ class TensorView(object): ax.set_zlabel('x3') ax.grid(True) - ax.hold(False) if showIt: plt.show() def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None): From e30dfb447d69ae91b62f3f06f0dca8ba150f542d Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 17 Mar 2014 12:47:36 -0700 Subject: [PATCH 2/8] Added points2nodes method to meshutils --- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/meshutils.py | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 209d1cff..f96ca68a 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,5 +1,5 @@ from matutils import * -from meshutils import exampleLrmGrid, meshTensors +from meshutils import exampleLrmGrid, meshTensors, points2nodes from lrmutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 67d747c3..6b47d690 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -53,6 +53,27 @@ def meshTensors(*args): return list(tensors) if len(tensors) > 1 else tensors[0] +def points2nodes(mesh, pts): + """ + Move a list of the nearest nodes to a set of points + + :param simpeg.Mesh.TensorMesh mesh: The mesh + :param numpy.ndarray pts: Points to move} + :rtype: numpy.ndarray + :return: nodeInds + """ + + pts = np.atleast_2d(pts) + + assert mesh._meshType == 'TENSOR' + assert pts.shape[1] == mesh.dim + + nodeInds = np.empty(pts.shape[0], dtype=int) + + for i, pt in enumerate(pts): + nodeInds[i] = ((np.tile(pt, (mesh.nN,1)) - mesh.gridN)**2).sum(axis=1).argmin() + + return nodeInds if __name__ == '__main__': from SimPEG import mesh import matplotlib.pyplot as plt From 44ab7de6bcfb304fc72f6b35108452dfc53c18c0 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 19 Mar 2014 10:34:37 -0700 Subject: [PATCH 3/8] Wrapper utils for Solvers --- SimPEG/Solver.py | 26 +++++++++++-- SimPEG/Utils/SolverUtils.py | 75 +++++++++++++++++++++++++++++++++++++ SimPEG/Utils/__init__.py | 1 + docs/api_Utils.rst | 8 ++++ 4 files changed, 107 insertions(+), 3 deletions(-) create mode 100644 SimPEG/Utils/SolverUtils.py diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 68f36798..cc6d6338 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -1,7 +1,7 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg -from Utils.matutils import mkvc, sdiag +from Utils import mkvc, sdiag import warnings DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'} @@ -214,12 +214,18 @@ class Solver(object): 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} + algorithms = {'CG':sp.linalg.cg, 'QMR':sp.linalg.qmr} assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!" alg = algorithms[iterSolver] + if iterSolver == 'CG': + opts = {'M':M} + elif iterSolver == 'QMR': + #TODO: make preconditioner better. + opts = {'M1':np.sqrt(M), 'M2':np.sqrt(M)} + if len(b.shape) == 1 or b.shape[1] == 1: - x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter) + x, self.info = alg(self.A, b, tol=tol, maxiter=maxIter) else: x = np.empty_like(b) for i in range(b.shape[1]): @@ -381,3 +387,17 @@ if __name__ == '__main__': toc = time() - tic print x print 'Error CG = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info + + + + 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={'iterSolver': 'QMR', 'M':'J'}) + tic = time() + x = iSolve.solve(b) + toc = time() - tic + print x + print 'Error QMR = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py new file mode 100644 index 00000000..69ca6508 --- /dev/null +++ b/SimPEG/Utils/SolverUtils.py @@ -0,0 +1,75 @@ +import numpy as np +from matutils import mkvc +import warnings + +def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-8): + + def __init__(self, A, **kwargs): + self.A = A.tocsc() + self.kwargs = kwargs + if factorize: + self.solver = fun(self.A, **kwargs) + + def solve(self, b): + if len(b.shape) == 1 or b.shape[1] == 1: + # Just one RHS + if factorize: + X = self.solver.solve(b.flatten(), **self.kwargs) + else: + X = fun(self.A, b.flatten(), **self.kwargs) + else: # Multiple RHSs + X = np.empty_like(b) + for i in range(b.shape[1]): + if factorize: + X[:,i] = self.solver.solve(b[:,i]) + else: + X[:,i] = fun(self.A, b[:,i], **self.kwargs) + + if checkAccuracy: + nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(Utils.mkvc(b)) + if nrm > accuracyTol: + msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol) + print msg + warnings.warn(msg, RuntimeWarning) + return X + + return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve}) + + + +def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5): + + def __init__(self, A, **kwargs): + self.A = A.tocsc() + self.kwargs = kwargs + + def solve(self, b): + if len(b.shape) == 1 or b.shape[1] == 1: + # Just one RHS + out = fun(self.A, b.flatten(), **self.kwargs) + if type(out) is tuple and len(out) == 2: + # We are dealing with scipy output with an info! + X = out[0] + self.info = out[1] + else: + X = out + else: # Multiple RHSs + X = np.empty_like(b) + for i in range(b.shape[1]): + out = fun(self.A, b[:,i], **self.kwargs) + if type(out) is tuple and len(out) == 2: + # We are dealing with scipy output with an info! + X[:,i] = out[0] + self.info = out[1] + else: + X[:,i] = out + + if checkAccuracy: + nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(Utils.mkvc(b)) + if nrm > accuracyTol: + msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol) + print msg + warnings.warn(msg, RuntimeWarning) + return X + + return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve}) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 209d1cff..56a30bc2 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -4,6 +4,7 @@ from lrmutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate import ModelBuilder +import SolverUtils import types import time diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index f9e24e04..0ea3e22e 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -15,6 +15,14 @@ Matrix Utilities :members: :undoc-members: + +Solver Utilities +================ + +.. automodule:: SimPEG.Utils.SolverUtils + :members: + :undoc-members: + LRM Utilities ============= From 8ff74f12d9c1b1f48c15a46ea93a6ffb961cf4ad Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 19 Mar 2014 11:47:12 -0700 Subject: [PATCH 4/8] minor documentation updates --- docs/index.rst | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index e8da6b25..d4f25039 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,9 +3,11 @@ :alt: SimPEG :align: center -SimPEG (Simulation and Parameter Estimation in Geophysics) is a python -package for simulation and gradient based parameter estimation in the -context of geoscience applications. +Simulation and Parameter Estimation in Geophysics (SimPEG) +********************************************************** + +SimPEG is a python package for simulation and gradient +based parameter estimation in the context of geoscience applications. The vision is to create a package for finite volume simulation and parameter estimation with applications to geophysical imaging and subsurface flow. To enable From ecc7aaf9b0c11f1e48e97811ce7d5eec40e03319 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 19 Mar 2014 13:33:45 -0700 Subject: [PATCH 5/8] Updates to Solver Wrappers. bug fix in dependentProperty --- SimPEG/Utils/SolverUtils.py | 14 ++++++++------ SimPEG/Utils/__init__.py | 3 ++- docs/index.rst | 4 ++-- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index 69ca6508..53ffceba 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -2,7 +2,7 @@ import numpy as np from matutils import mkvc import warnings -def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-8): +def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): def __init__(self, A, **kwargs): self.A = A.tocsc() @@ -12,11 +12,12 @@ def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-8): def solve(self, b): if len(b.shape) == 1 or b.shape[1] == 1: + b = b.flatten() # Just one RHS if factorize: - X = self.solver.solve(b.flatten(), **self.kwargs) + X = self.solver.solve(b, **self.kwargs) else: - X = fun(self.A, b.flatten(), **self.kwargs) + X = fun(self.A, b, **self.kwargs) else: # Multiple RHSs X = np.empty_like(b) for i in range(b.shape[1]): @@ -26,7 +27,7 @@ def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-8): X[:,i] = fun(self.A, b[:,i], **self.kwargs) if checkAccuracy: - nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(Utils.mkvc(b)) + nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(mkvc(b)) if nrm > accuracyTol: msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol) print msg @@ -45,8 +46,9 @@ def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5): def solve(self, b): if len(b.shape) == 1 or b.shape[1] == 1: + b = b.flatten() # Just one RHS - out = fun(self.A, b.flatten(), **self.kwargs) + out = fun(self.A, b, **self.kwargs) if type(out) is tuple and len(out) == 2: # We are dealing with scipy output with an info! X = out[0] @@ -65,7 +67,7 @@ def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5): X[:,i] = out if checkAccuracy: - nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(Utils.mkvc(b)) + nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(mkvc(b)) if nrm > accuracyTol: msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol) print msg diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index ce0a6b10..94d016f6 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -136,7 +136,8 @@ def callHooks(match, mainFirst=False): def dependentProperty(name, value, children, doc): def fget(self): return getattr(self,name,value) def fset(self, val): - if getattr(self,name,value) == val: return # it is the same! + if isScalar(val) and getattr(self,name,value) == val: + return # it is the same! for child in children: if hasattr(self, child): delattr(self, child) diff --git a/docs/index.rst b/docs/index.rst index d4f25039..71aaa282 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,8 +3,8 @@ :alt: SimPEG :align: center -Simulation and Parameter Estimation in Geophysics (SimPEG) -********************************************************** +Simulation and Parameter Estimation in Geophysics +************************************************* SimPEG is a python package for simulation and gradient based parameter estimation in the context of geoscience applications. From 61a08408d016fbc47b7dc75369e192d1037f906d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Mar 2014 14:06:18 -0700 Subject: [PATCH 6/8] Changed BaseRxList --> BaseRx --- SimPEG/Survey.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 42ee50be..936560d6 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -159,7 +159,7 @@ class BaseSurvey(object): return self.mtrue is not None - #TODO: Move this to the data class? + #TODO: Move this to the survey class? # @property # def phi_d_target(self): # """ @@ -178,8 +178,8 @@ class BaseSurvey(object): # self._phi_d_target = value -class BaseRxList(object): - """SimPEG Receiver List Object""" +class BaseRx(object): + """SimPEG Receiver Object""" locs = None #: Locations (nRx x 3) @@ -207,12 +207,15 @@ class BaseTx(object): loc = None #: Location [x,y,z] rxList = None #: SimPEG Receiver List - rxListPair = BaseRxList + rxPair = BaseRx knownTxTypes = None #: Set this to a list of strings to ensure that txType is known def __init__(self, loc, txType, rxList, **kwargs): - assert isinstance(rxList, self.rxListPair), 'rxList must be a %s'%self.rxListPair.__name__ + assert type(rxList) is list, 'rxList must be a list' + for rx in rxList: + assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxListPair.__name__ + assert len(set(rxList)) == len(rxList), 'The rxList must be unique' self.loc = loc self.txType = txType From a3e7b96d2e4d418fb965829c537953c08cf6d66b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 23 Mar 2014 10:22:39 -0700 Subject: [PATCH 7/8] Key word args for the model.test() function that go to checkDeriv --- SimPEG/Model.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index ce2349d0..db87f460 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -58,11 +58,13 @@ class BaseModel(object): def example(self): return np.random.rand(self.nP) - def test(self, m=None): + def test(self, m=None, **kwargs): print 'Testing the %s Class!' % self.__class__.__name__ if m is None: m = self.example() - return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, plotIt=False) + if 'plotIt' not in kwargs: + kwargs['plotIt'] = False + return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, **kwargs) class BaseNonLinearModel(object): """ From 9ec3e20ad98a3b5d73a9a0e5972723ace5258ca0 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 23 Mar 2014 15:00:09 -0700 Subject: [PATCH 8/8] bug fix in ActiveModel --- SimPEG/Model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index db87f460..00f8f302 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -310,7 +310,7 @@ class ActiveModel(BaseModel): indActive = z self.indActive = indActive self.indInactive = np.logical_not(indActive) - if type(valInactive) in [float, int, long]: + if Utils.isScalar(valInactive): valInactive = np.ones(self.nC)*float(valInactive) valInactive[self.indActive] = 0