From a595893f884a62933dd0653fd290ce9e30c42a7a Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Wed, 20 Nov 2013 22:46:38 -0800 Subject: [PATCH 01/25] Simple startup script to install necessary software on google compute engine. --- SimPEG/GCEtools/startup.sh | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 SimPEG/GCEtools/startup.sh diff --git a/SimPEG/GCEtools/startup.sh b/SimPEG/GCEtools/startup.sh new file mode 100644 index 00000000..0c593a63 --- /dev/null +++ b/SimPEG/GCEtools/startup.sh @@ -0,0 +1,13 @@ +#! /bin/bash +sudo aptitude -y update +sudo aptitude -y upgrade +sudo aptitude -y install gcc gfortran git libopenmpi-dev python-pip python-dev +sudo aptitude -y install ipython python-scipy python-numpy +sudo aptitude -y install python-matplotlib python-nose python-pip +sudo aptitude -y install libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev +sudo pip install mpi4py +sudo pip install pymumps +git clone https://dwfmarchant@bitbucket.org/rcockett/simpeg.git +cd simpeg/SimPEG/ +python setup.py +cd ~ \ No newline at end of file From 594bf20d5fb5eac7057f088d0bc865b13718b617 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 21 Nov 2013 10:15:30 -0800 Subject: [PATCH 02/25] Updated startup script. --- SimPEG/GCEtools/gceStartup.txt | 11 +++++++++++ SimPEG/GCEtools/startup.sh | 15 ++++++++++++--- 2 files changed, 23 insertions(+), 3 deletions(-) create mode 100644 SimPEG/GCEtools/gceStartup.txt diff --git a/SimPEG/GCEtools/gceStartup.txt b/SimPEG/GCEtools/gceStartup.txt new file mode 100644 index 00000000..76bbcb39 --- /dev/null +++ b/SimPEG/GCEtools/gceStartup.txt @@ -0,0 +1,11 @@ +# Check project status +gcutil getproject --project= --cache_flag_values + +# Start an instance +gcutil addinstance + +# Log in +gcutil ssh + +# Shut down +gcutil deleteinstance \ No newline at end of file diff --git a/SimPEG/GCEtools/startup.sh b/SimPEG/GCEtools/startup.sh index 0c593a63..a4055ea0 100644 --- a/SimPEG/GCEtools/startup.sh +++ b/SimPEG/GCEtools/startup.sh @@ -2,12 +2,21 @@ sudo aptitude -y update sudo aptitude -y upgrade sudo aptitude -y install gcc gfortran git libopenmpi-dev python-pip python-dev -sudo aptitude -y install ipython python-scipy python-numpy -sudo aptitude -y install python-matplotlib python-nose python-pip +sudo aptitude -y install ipython python-scipy python-numpy python-nose python-pip sudo aptitude -y install libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev +sudo aptitude -y install libblas-dev liblapack-dev + sudo pip install mpi4py sudo pip install pymumps + +sudo pip install scipy --upgrade +sudo pip install numpy --upgrade +sudo pip install ipython --upgrade + git clone https://dwfmarchant@bitbucket.org/rcockett/simpeg.git cd simpeg/SimPEG/ python setup.py -cd ~ \ No newline at end of file +cd ~ + +echo export PYTHONPATH=/home/$USER/simpeg/ >> .bashrc +source .profile \ No newline at end of file From 5f64ea8399e114f1909e906d8d53a698941ebb7e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 21 Nov 2013 12:25:56 -0800 Subject: [PATCH 03/25] Documentation Updates. --- docs/api_LogicallyOrthogonalMesh.rst | 2 ++ docs/api_Optimize.rst | 4 ++++ docs/api_Problem.rst | 9 ++++++--- docs/api_TensorMesh.rst | 2 ++ 4 files changed, 14 insertions(+), 3 deletions(-) diff --git a/docs/api_LogicallyOrthogonalMesh.rst b/docs/api_LogicallyOrthogonalMesh.rst index 9d7d516b..c59d7886 100644 --- a/docs/api_LogicallyOrthogonalMesh.rst +++ b/docs/api_LogicallyOrthogonalMesh.rst @@ -4,8 +4,10 @@ Logically Orthogonal Mesh ************************* .. automodule:: SimPEG.mesh.LogicallyOrthogonalMesh + :show-inheritance: :members: :undoc-members: + :inherited-members: LOM View diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst index d726b5f8..283d8ce9 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Optimize.rst @@ -4,16 +4,20 @@ Optimize ******** .. automodule:: SimPEG.inverse.Optimize + :show-inheritance: :members: :undoc-members: + :inherited-members: Inversion ********* .. automodule:: SimPEG.inverse.Inversion + :show-inheritance: :members: :undoc-members: + :inherited-members: Beta Schedule ************* diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst index 3d5a6d2c..b4c96c4a 100644 --- a/docs/api_Problem.rst +++ b/docs/api_Problem.rst @@ -1,21 +1,24 @@ .. _api_Problem: - Problem ******* .. automodule:: SimPEG.forward.Problem + :show-inheritance: :members: :undoc-members: + :inherited-members: DCProblem ********* .. automodule:: SimPEG.forward.DCProblem + :show-inheritance: :members: :undoc-members: + :inherited-members: @@ -23,7 +26,7 @@ Linear Problem ************** .. automodule:: SimPEG.forward.LinearProblem + :show-inheritance: :members: :undoc-members: - - + :inherited-members: diff --git a/docs/api_TensorMesh.rst b/docs/api_TensorMesh.rst index 9de206f1..5dc3b280 100644 --- a/docs/api_TensorMesh.rst +++ b/docs/api_TensorMesh.rst @@ -4,8 +4,10 @@ Tensor Mesh *********** .. automodule:: SimPEG.mesh.TensorMesh + :show-inheritance: :members: :undoc-members: + :inherited-members: Tensor View *********** From fff2375acd1f2b58a334a77416b54ce434b9949e Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 21 Nov 2013 14:34:32 -0800 Subject: [PATCH 04/25] Moved scipy direct code to own function. --- SimPEG/utils/Solver.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 16db872f..85d5f640 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -120,6 +120,7 @@ class Solver(object): """ 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 @@ -129,6 +130,20 @@ class Solver(object): assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' + if backend == 'scipy': + X = self.solveDirect_scipy(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) From 623ff674ea1da5f97cf62453bc3a32411e2bfa8d Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 21 Nov 2013 15:00:51 -0800 Subject: [PATCH 05/25] Added pymumps direct solver backend. --- SimPEG/utils/Solver.py | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 85d5f640..c451a7ca 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -19,6 +19,11 @@ except Exception, e: DEFAULTS['forward'] = 'python' DEFAULTS['backward'] = 'python' +try: + import mumps +except Exception, e: + print 'Warning: mumps solver not available.' + class Solver(object): """ Solver is a light wrapper on the various types of @@ -113,6 +118,9 @@ class Solver(object): 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 @@ -132,6 +140,8 @@ class Solver(object): if backend == 'scipy': X = self.solveDirect_scipy(b, factorize) + elif backend == 'mumps': + X = self.solveDirect_mumps(b, factorize) return X @@ -165,6 +175,40 @@ class Solver(object): 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: + return self.dsolve(b) + else: + return mumps.spsolve(self.A, b) + + else: + raise NotImplementedError('Multiple RHS not yet implemented with mumps solver.') + def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50): if backend is None: backend = DEFAULTS['iter'] From 0bd971c22fcfd2e6c010141103f08604047bc9d5 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 21 Nov 2013 15:12:28 -0800 Subject: [PATCH 06/25] Removed PubSub based communication. --- SimPEG/inverse/Optimize.py | 36 ------------------------------------ 1 file changed, 36 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 75209dc6..e305d923 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -5,12 +5,6 @@ norm = np.linalg.norm import scipy.sparse as sp from SimPEG import Solver -try: - from pubsub import pub - doPub = True -except Exception, e: - print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.' - doPub = False class StoppingCriteria(object): """docstring for StoppingCriteria""" @@ -142,27 +136,6 @@ class Minimize(object): return out if len(out) > 1 else out[0] - Events are fired with the following inputs via pypubsub:: - - Minimize.printInit (minimize) - Minimize.evalFunction (minimize, f, g, H) - Minimize.printIter (minimize) - Minimize.searchDirection (minimize, p) - Minimize.scaleSearchDirection (minimize, p) - Minimize.modifySearchDirection (minimize, xt, passLS) - Minimize.endIteration (minimize, xt) - Minimize.printDone (minimize) - - To hook into one of these events (must have pypubsub installed):: - - from pubsub import pub - def listener(minimize,p): - print 'The search direction is: ', p - pub.subscribe(listener, 'Minimize.searchDirection') - - You can use pubsub communication to debug your code, it is not used internally. - - The algorithm for general minimization is as follows:: startup(x0) @@ -189,20 +162,15 @@ class Minimize(object): while True: self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H) self.printIter() if self.stoppingCriteria(): break p = self.findSearchDirection() - if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) p = self.scaleSearchDirection(p) - if doPub: pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p) xt, passLS = self.modifySearchDirection(p) - if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt, passLS=passLS) if not passLS: xt, caught = self.modifySearchDirectionBreak(p) if not caught: return self.xc self.doEndIteration(xt) - if doPub: pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt) self.printDone() @@ -857,10 +825,6 @@ if __name__ == '__main__': x0 = np.array([2.6, 3.7]) checkDerivative(Rosenbrock, x0, plotIt=False) - # def listener1(minimize,p): - # print 'hi: ', p - # if doPub: pub.subscribe(listener1, 'Minimize.searchDirection') - xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0) From dbaea1fda95945dccd4a5c7470efd4219d979925 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 21 Nov 2013 16:56:01 -0800 Subject: [PATCH 07/25] callHooks generalizes some of the hook calling code in Optimize and Inversion --- SimPEG/inverse/Inversion.py | 16 ++++++---------- SimPEG/inverse/Optimize.py | 20 +++++--------------- SimPEG/mesh/DiffOperators.py | 6 ++---- SimPEG/regularization/Regularization.py | 12 +++++++----- SimPEG/utils/__init__.py | 8 ++++++++ 5 files changed, 28 insertions(+), 34 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 5a71f6ec..254e7651 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,7 +1,7 @@ import numpy as np import scipy.sparse as sp import SimPEG -from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt +from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt, callHooks from Optimize import Remember from BetaSchedule import Cooling @@ -31,9 +31,9 @@ class BaseInversion(object): self.opt.printers.insert(3,SimPEG.inverse.IterationPrinters.phi_m) self.opt.stoppers.append(SimPEG.inverse.StoppingCriteria.phi_d_target_Minimize) - if not hasattr(opt, '_bfgsH0'): # Check if it has been set by the user and the default is not being used. - print 'Setting bfgsH0 to the inverse of the modelObj2Deriv.' - opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv(),doDirect=True,options={'factorize':True}) # False, options={'M':'GS','maxIter':15} + if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used. + print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' + opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv()) @property @@ -94,9 +94,7 @@ class BaseInversion(object): :rtype: None :return: None """ - for method in [posible for posible in dir(self) if '_startup' in posible]: - if self.debug: print 'startup is calling self.'+method - getattr(self,method)(m0) + callHooks(self,'startup',m0) if not hasattr(self.reg, '_mref'): print 'Regularization has not set mref. SimPEG will set it to m0.' @@ -124,9 +122,7 @@ class BaseInversion(object): :rtype: None :return: None """ - for method in [posible for posible in dir(self) if '_doEndIteration' in posible]: - if self.debug: print 'doEndIteration is calling self.'+method - getattr(self,method)() + callHooks(self,'doEndIteration') # store old values self.phi_d_last = self.phi_d diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index e305d923..40ddb310 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt +from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt, callHooks norm = np.linalg.norm import scipy.sparse as sp from SimPEG import Solver @@ -104,7 +104,6 @@ class Minimize(object): counter = None def __init__(self, **kwargs): - self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration] self.stoppersLS = [StoppingCriteria.armijoGoldstein, StoppingCriteria.iterationLS] @@ -208,9 +207,7 @@ class Minimize(object): :rtype: None :return: None """ - for method in [posible for posible in dir(self) if '_startup' in posible]: - if self.debug: print 'startup is calling self.'+method - getattr(self,method)(x0) + callHooks(self,'startup',x0) self._iter = 0 self._iterLS = 0 @@ -230,7 +227,6 @@ class Minimize(object): parent.printInit function and call that. """ - if doPub and not inLS: pub.sendMessage('Minimize.printInit', minimize=self) pad = ' '*10 if inLS else '' name = self.name if not inLS else self.nameLS printTitles(self, self.printers if not inLS else self.printersLS, name, pad) @@ -244,12 +240,8 @@ class Minimize(object): parent.printIter function and call that. """ + callHooks(self,'printIter',inLS) - for method in [posible for posible in dir(self) if '_printIter' in posible]: - if self.debug: print 'printIter is calling self.'+method - getattr(self,method)(inLS) - - if doPub and not inLS: pub.sendMessage('Minimize.printIter', minimize=self) pad = ' '*10 if inLS else '' printLine(self, self.printers if not inLS else self.printersLS, pad=pad) @@ -261,7 +253,6 @@ class Minimize(object): parent.printDone function and call that. """ - if doPub and not inLS: pub.sendMessage('Minimize.printDone', minimize=self) pad = ' '*10 if inLS else '' stop, done = (' STOP! ', ' DONE! ') if not inLS else ('----------------', ' End Linesearch ') stoppers = self.stoppers if not inLS else self.stoppersLS @@ -285,6 +276,7 @@ class Minimize(object): :rtype: numpy.ndarray :return: p, projected search direction """ + callHooks(self,'projection',p) return p @timeIt @@ -415,9 +407,7 @@ class Minimize(object): :rtype: None :return: None """ - for method in [posible for posible in dir(self) if '_doEndIteration' in posible]: - if self.debug: print 'doEndIteration is calling self.'+method - getattr(self,method)(xt) + callHooks(self,'doEndIteration',xt) # store old values self.f_last = self.f diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 19a3ab34..4e017c7c 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -318,8 +318,7 @@ class DiffOperators(object): def cellGrady(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): - if self.dim < 2: - return None + if self.dim < 2: return None if getattr(self, '_cellGrady', None) is None: BC = ['neumann', 'neumann'] n = self.n @@ -338,8 +337,7 @@ class DiffOperators(object): def cellGradz(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): - if self.dim < 3: - return None + if self.dim < 3: return None if getattr(self, '_cellGradz', None) is None: BC = ['neumann', 'neumann'] n = self.n diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 3e429324..3d22a297 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -1,4 +1,4 @@ -from SimPEG.utils import sdiag, count, timeIt +from SimPEG.utils import sdiag, count, timeIt, setKwargs import numpy as np class Regularization(object): @@ -41,15 +41,17 @@ class Regularization(object): return self._Wz alpha_s = 1e-6 - alpha_x = 1 - alpha_y = 1 - alpha_z = 1 + alpha_x = 1.0 + alpha_y = 1.0 + alpha_z = 1.0 counter = None - def __init__(self, mesh): + def __init__(self, mesh, **kwargs): + setKwargs(self, **kwargs) self.mesh = mesh + def pnorm(self, r): return 0.5*r.dot(r) diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 2ba0f44d..f7daf273 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -61,6 +61,12 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'): print pad + stopper['str'] % (l<=r,l,r) print pad + "%s%s%s" % ('-'*25,done,'-'*25) +def callHooks(obj, match, *args, **kwargs): + for method in [posible for posible in dir(obj) if ('_'+match) in posible]: + if getattr(obj,'debug',False): print (match+' is calling self.'+method) + getattr(obj,method)(*args, **kwargs) + + import time import numpy as np @@ -154,6 +160,8 @@ def timeIt(f): 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): From 66050f474f53900adf0053af7c9c4bcf8fa075be Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 22 Nov 2013 01:17:06 +0000 Subject: [PATCH 08/25] Changed to generic bitbucket git repo. --- SimPEG/GCEtools/startup.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/GCEtools/startup.sh b/SimPEG/GCEtools/startup.sh index a4055ea0..899f4b26 100644 --- a/SimPEG/GCEtools/startup.sh +++ b/SimPEG/GCEtools/startup.sh @@ -13,7 +13,7 @@ sudo pip install scipy --upgrade sudo pip install numpy --upgrade sudo pip install ipython --upgrade -git clone https://dwfmarchant@bitbucket.org/rcockett/simpeg.git +git clone https://bitbucket.org/rcockett/simpeg.git cd simpeg/SimPEG/ python setup.py cd ~ From 6c029e5e2c5d0bd379b5e086c630119579215483 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 22 Nov 2013 01:26:22 +0000 Subject: [PATCH 09/25] Added matplotlib into the startup. SimPEG needs it. --- SimPEG/GCEtools/startup.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/GCEtools/startup.sh b/SimPEG/GCEtools/startup.sh index 899f4b26..32164e29 100644 --- a/SimPEG/GCEtools/startup.sh +++ b/SimPEG/GCEtools/startup.sh @@ -2,7 +2,7 @@ sudo aptitude -y update sudo aptitude -y upgrade sudo aptitude -y install gcc gfortran git libopenmpi-dev python-pip python-dev -sudo aptitude -y install ipython python-scipy python-numpy python-nose python-pip +sudo aptitude -y install ipython python-scipy python-numpy python-nose python-pip python-matplotlib sudo aptitude -y install libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev sudo aptitude -y install libblas-dev liblapack-dev From 6216ae977ef52d74a460c0551bf46f8a42b1b060 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 22 Nov 2013 10:51:01 -0800 Subject: [PATCH 10/25] Added hook method to utils so that we can hook new methods dynamically into the classes. --- SimPEG/inverse/Inversion.py | 1 - SimPEG/utils/__init__.py | 17 ++++++++++++++--- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 254e7651..a8e05a1a 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -118,7 +118,6 @@ class BaseInversion(object): Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call. You may also completely overwrite this function. - :param numpy.ndarray xt: tested new iterate that ensures a descent direction. :rtype: None :return: None """ diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index f7daf273..724f6f38 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -12,6 +12,10 @@ import Solver from Solver import Solver import Geophysics +import types +import time +import numpy as np + 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: @@ -66,10 +70,17 @@ def callHooks(obj, match, *args, **kwargs): if getattr(obj,'debug',False): print (match+' is calling self.'+method) getattr(obj,method)(*args, **kwargs) +def hook(obj, method, name=None, overwrite=False): + """ + This dynamically binds a method to the instance of the class. - -import time -import numpy as np + If name is None, the name of the method is used. + """ + if name is None: name = method.__name__ + if not hasattr(obj,name) or overwrite: + setattr(obj, name, types.MethodType( method, obj )) + else: + print 'Method '+name+' was not overwritten.' class Counter(object): From c41b9f12e9b61296df2f039033238b1e9884380c Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 22 Nov 2013 11:58:31 -0800 Subject: [PATCH 11/25] Documentation updates to optimize. Updated the hook method. --- SimPEG/inverse/Optimize.py | 88 +++++++++++++++++++++++--------------- SimPEG/utils/__init__.py | 29 ++++++++----- docs/api_Optimize.rst | 3 +- 3 files changed, 73 insertions(+), 47 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 40ddb310..db708881 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -6,6 +6,9 @@ import scipy.sparse as sp from SimPEG import Solver +__all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters'] + + class StoppingCriteria(object): """docstring for StoppingCriteria""" @@ -79,29 +82,26 @@ class IterationPrinters(object): class Minimize(object): """ - Minimize is a general class for derivative based optimization. - - """ - name = "General Optimization Algorithm" + name = "General Optimization Algorithm" #: The name of the optimization algorithm - maxIter = 20 - maxIterLS = 10 - maxStep = np.inf - LSreduction = 1e-4 - LSshorten = 0.5 - tolF = 1e-1 - tolX = 1e-1 - tolG = 1e-1 - eps = 1e-5 + maxIter = 20 #: Maximum number of iterations + maxIterLS = 10 #: Maximum number of iterations for the line-search + maxStep = np.inf #: Maximum step possible, used in scaling before the line-search. + LSreduction = 1e-4 #: Expected decrease in the line-search + LSshorten = 0.5 #: Line-search step is shortened by this amount each time. + tolF = 1e-1 #: Tolerance on function value decrease + tolX = 1e-1 #: Tolerance on norm(x) movement + tolG = 1e-1 #: Tolerance on gradient norm + eps = 1e-5 #: Small value - debug = False - debugLS = False + debug = False #: Print debugging information + debugLS = False #: Print debugging information for the line-search - comment = '' - counter = None + comment = '' #: Used by some functions to indicate what is going on in the algorithm + counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things def __init__(self, **kwargs): self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration] @@ -114,7 +114,8 @@ class Minimize(object): @timeIt def minimize(self, evalFunction, x0): - """ + """minimize(evalFunction, x0) + Minimizes the function (evalFunction) starting at the location x0. :param def evalFunction: function handle that evaluates: f, g, H = F(x) @@ -231,7 +232,6 @@ class Minimize(object): name = self.name if not inLS else self.nameLS printTitles(self, self.printers if not inLS else self.printersLS, name, pad) - @count def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -258,7 +258,6 @@ class Minimize(object): stoppers = self.stoppers if not inLS else self.stoppersLS printStoppers(self, stoppers, pad='', stop=stop, done=done) - @timeIt def stoppingCriteria(self, inLS=False): if self._iter == 0: self.f0 = self.f @@ -267,7 +266,8 @@ class Minimize(object): @timeIt def projection(self, p): - """ + """projection(p) + projects the search direction. by default, no projection is applied. @@ -281,7 +281,8 @@ class Minimize(object): @timeIt def findSearchDirection(self): - """ + """findSearchDirection() + **findSearchDirection** should return an approximation of: .. math:: @@ -311,7 +312,8 @@ class Minimize(object): @count def scaleSearchDirection(self, p): - """ + """scaleSearchDirection(p) + **scaleSearchDirection** should scale the search direction if appropriate. Set the parameter **maxStep** in the minimize object, to scale back the gradient to a maximum size. @@ -325,11 +327,12 @@ class Minimize(object): p = self.maxStep*p/np.abs(p.max()) return p - nameLS = "Armijo linesearch" + nameLS = "Armijo linesearch" #: The line-search name @timeIt def modifySearchDirection(self, p): - """ + """modifySearchDirection(p) + **modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria. By default, an Armijo backtracking linesearch is preformed with the following parameters: @@ -366,7 +369,8 @@ class Minimize(object): @count def modifySearchDirectionBreak(self, p): - """ + """modifySearchDirectionBreak(p) + Code is called if modifySearchDirection fails to find a descent direction. @@ -387,7 +391,8 @@ class Minimize(object): @count def doEndIteration(self, xt): - """ + """doEndIteration(xt) + **doEndIteration** is called at the end of each minimize iteration. By default, function values and x locations are shuffled to store 1 past iteration in memory. @@ -416,7 +421,6 @@ class Minimize(object): if self.debug: self.printDone() - class Remember(object): """ This mixin remembers all the things you tend to forget. @@ -467,7 +471,6 @@ class Remember(object): self._rememberList[param[0]].append( param[1](self) ) - class ProjectedGradient(Minimize, Remember): name = 'Projected Gradient' @@ -504,22 +507,35 @@ class ProjectedGradient(Minimize, Remember): @count def projection(self, x): - """Make sure we are feasible.""" + """projection(x) + + Make sure we are feasible. + + """ return np.median(np.c_[self.lower,x,self.upper],axis=1) @count def activeSet(self, x): - """If we are on a bound""" + """activeSet(x) + + If we are on a bound + + """ return np.logical_or(x == self.lower, x == self.upper) @count def inactiveSet(self, x): - """The free variables.""" + """inactiveSet(x) + + The free variables. + + """ return np.logical_not(self.activeSet(x)) @count def bindingSet(self, x): - """ + """bindingSet(x) + If we are on a bound and the negative gradient points away from the feasible set. Optimality condition. (Satisfies Kuhn-Tucker) MoreToraldo91 @@ -531,6 +547,10 @@ class ProjectedGradient(Minimize, Remember): @timeIt def findSearchDirection(self): + """findSearchDirection() + + Finds the search direction based on either CG or steepest descent. + """ self.aSet_prev = self.activeSet(self.xc) allBoundsAreActive = sum(self.aSet_prev) == self.xc.size @@ -572,6 +592,7 @@ class ProjectedGradient(Minimize, Remember): @timeIt def _doEndIteration_ProjectedGradient(self, xt): + """_doEndIteration_ProjectedGradient(xt)""" aSet = self.activeSet(xt) bSet = self.bindingSet(xt) @@ -600,7 +621,6 @@ class ProjectedGradient(Minimize, Remember): if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD - class BFGS(Minimize, Remember): name = 'BFGS' nbfgs = 10 diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 724f6f38..0c1e9ad1 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -15,6 +15,19 @@ 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 not hasattr(obj,name) or overwrite: + setattr(obj, name, types.MethodType( method, obj )) + elif not silent: + 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.""" @@ -23,6 +36,9 @@ def setKwargs(obj, **kwargs): setattr(obj, attr, kwargs[attr]) else: raise Exception('%s attr is not recognized' % attr) + hook(obj,callHooks, silent=True) + hook(obj,hook, silent=True) + hook(obj,setKwargs, silent=True) def printTitles(obj, printers, name='Print Titles', pad=''): titles = '' @@ -70,17 +86,6 @@ def callHooks(obj, match, *args, **kwargs): if getattr(obj,'debug',False): print (match+' is calling self.'+method) getattr(obj,method)(*args, **kwargs) -def hook(obj, method, name=None, overwrite=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 not hasattr(obj,name) or overwrite: - setattr(obj, name, types.MethodType( method, obj )) - else: - print 'Method '+name+' was not overwritten.' class Counter(object): @@ -156,6 +161,7 @@ class Counter(object): 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__) @@ -164,6 +170,7 @@ def count(f): 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__) diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst index 283d8ce9..8ca74e27 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Optimize.rst @@ -5,9 +5,9 @@ Optimize .. automodule:: SimPEG.inverse.Optimize :show-inheritance: + :private-members: :members: :undoc-members: - :inherited-members: Inversion @@ -17,7 +17,6 @@ Inversion :show-inheritance: :members: :undoc-members: - :inherited-members: Beta Schedule ************* From aae725ff0535c0d60631803e7cbf73e7188254c3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 22 Nov 2013 13:39:26 -0800 Subject: [PATCH 12/25] TimeSteppingInversion and Estimate Initial Beta based on eigenvalues comparison. --- SimPEG/inverse/Inversion.py | 80 +++++++++++++++++++++++++++++++------ SimPEG/inverse/Optimize.py | 8 ++-- SimPEG/utils/__init__.py | 7 +++- 3 files changed, 76 insertions(+), 19 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index a8e05a1a..933e296c 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -4,16 +4,19 @@ import SimPEG from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt, callHooks from Optimize import Remember from BetaSchedule import Cooling +from SimPEG.inverse import IterationPrinters, StoppingCriteria class BaseInversion(object): """docstring for BaseInversion""" - maxIter = 1 + maxIter = 1 #: Maximum number of iterations name = 'BaseInversion' - debug = False - beta0 = 1e4 - counter = None + debug = False #: Print debugging information + + comment = '' #: Used by some functions to indicate what is going on in the algorithm + counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things + def __init__(self, prob, reg, opt, **kwargs): setKwargs(self, **kwargs) @@ -22,14 +25,13 @@ class BaseInversion(object): self.opt = opt self.opt.parent = self - self.stoppers = [SimPEG.inverse.StoppingCriteria.iteration, SimPEG.inverse.StoppingCriteria.phi_d_target_Inversion] + self.stoppers = [StoppingCriteria.iteration] # Check if we have inserted printers into the optimization - if not np.any([p is SimPEG.inverse.IterationPrinters.phi_d for p in self.opt.printers]): - self.opt.printers.insert(1,SimPEG.inverse.IterationPrinters.beta) - self.opt.printers.insert(2,SimPEG.inverse.IterationPrinters.phi_d) - self.opt.printers.insert(3,SimPEG.inverse.IterationPrinters.phi_m) - self.opt.stoppers.append(SimPEG.inverse.StoppingCriteria.phi_d_target_Minimize) + if IterationPrinters.phi_d not in self.opt.printers: + self.opt.printers.insert(1,IterationPrinters.beta) + self.opt.printers.insert(2,IterationPrinters.phi_d) + self.opt.printers.insert(3,IterationPrinters.phi_m) if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used. print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' @@ -128,9 +130,29 @@ class BaseInversion(object): self.phi_m_last = self.phi_m self._iter += 1 + @property + def beta0(self): + if getattr(self,'_beta0',None) is None: + self._beta0 = self.estimateBeta0() + return self._beta0 + @beta0.setter + def beta0(self, value): + self._beta0 = value + def getBeta(self): return self.beta0 + def estimateBeta0(self): + """ + + + """ + u = self.prob.field(self.m) + v = np.random.rand(*self.m.shape) + t = v.dot(self.dataObj2Deriv(self.m,v,u=u)) + b = v.dot(self.reg.modelObj2Deriv()*v) + return 0.1*(t/b) + def stoppingCriteria(self): if self.debug: print 'checking stoppingCriteria' return checkStoppers(self, self.stoppers) @@ -170,7 +192,7 @@ class BaseInversion(object): return phi_d2Deriv + self._beta * phi_m2Deriv - operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float ) + operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype ) out += (operator,) return out if len(out) > 1 else out[0] @@ -239,8 +261,10 @@ class BaseInversion(object): @timeIt def dataObj2Deriv(self, m, v, u=None): - """ + """dataObj2Deriv(m, v, u=None) + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply :param numpy.array u: fields :rtype: numpy.array :return: data misfit derivative @@ -277,7 +301,7 @@ class BaseInversion(object): R = self.Wd*self.prob.dataResidual(m, u=u) # TODO: abstract to different norms a little cleaner. - # \/ it goes here. in l2 it is the identity. + # \/ it goes here. in l2 it is the identity. dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u) return dmisfit @@ -289,3 +313,33 @@ class Inversion(Cooling, Remember, BaseInversion): def __init__(self, prob, reg, opt, **kwargs): BaseInversion.__init__(self, prob, reg, opt, **kwargs) + + self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) + + if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: + self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) + +class TimeSteppingInversion(Remember, BaseInversion): + """ + A slightly different view on regularization parameters, + let Beta be viewed as 1/dt, and timestep by updating the + reference model every optimization iteration. + """ + maxIter = 1 + name = "Time-Stepping SimPEG Inversion" + + def __init__(self, prob, reg, opt, **kwargs): + BaseInversion.__init__(self, prob, reg, opt, **kwargs) + + self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) + + if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: + self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) + + def _startup_TimeSteppingInversion(self, m0): + + def _doEndIteration_updateMref(self, xt): + if self.debug: 'Updating the reference model.' + self.parent.reg.mref = self.xc + + self.opt.hook(_doEndIteration_updateMref, overwrite=True) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index db708881..5bf3ef4d 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -474,8 +474,8 @@ class Remember(object): class ProjectedGradient(Minimize, Remember): name = 'Projected Gradient' - maxIterCG = 10 - tolCG = 1e-3 + maxIterCG = 5 + tolCG = 1e-1 lower = -np.inf upper = np.inf @@ -716,8 +716,8 @@ class InexactGaussNewton(BFGS, Minimize, Remember): name = 'Inexact Gauss Newton' - maxIterCG = 10 - tolCG = 1e-3 + maxIterCG = 5 + tolCG = 1e-1 @property def approxHinv(self): diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 0c1e9ad1..b1cb47e8 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -26,9 +26,12 @@ def hook(obj, method, name=None, overwrite=False, silent=False): if name is None: name = method.__name__ if not hasattr(obj,name) or overwrite: setattr(obj, name, types.MethodType( method, obj )) - elif not silent: + 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: @@ -97,7 +100,7 @@ class Counter(object): 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): From 82f5a5660fd02bc0d27b1e2aa04d066004848d3c Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 22 Nov 2013 14:06:22 -0800 Subject: [PATCH 13/25] Documentation on how I am picking the initial beta. --- SimPEG/inverse/Inversion.py | 61 +++++++++++++++++++++++++++++++------ 1 file changed, 51 insertions(+), 10 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 933e296c..b5e4906f 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -70,6 +70,11 @@ class BaseInversion(object): @timeIt def run(self, m0): + """run(m0) + + Runs the inversion! + + """ self.startup(m0) while True: self._beta = self.getBeta() @@ -114,7 +119,7 @@ class BaseInversion(object): If you have things that also need to run at the end of every iteration, you can create a method:: - def _doEndIteration*(self, xt): + def _doEndIteration*(self): pass Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call. @@ -142,16 +147,46 @@ class BaseInversion(object): def getBeta(self): return self.beta0 - def estimateBeta0(self): - """ + def estimateBeta0(self, u=None, ratio=0.1): + """estimateBeta0(u=None, ratio=0.1) + + The initial beta is calculated by comparing the estimated + eigenvalues of JtJ and WtW. + + To estimate the eigenvector of **A**, we will use one iteration + of the *Power Method*: + + .. math:: + + \mathbf{x_1 = A x_0} + + Given this (very course) approximation of the eigenvector, + we can use the *Rayleigh quotient* to approximate the largest eigenvalue. + + .. math:: + + \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}} + + We will approximate the largest eigenvalue for both JtJ and WtW, and + use some ratio of the quotient to estimate beta0. + + .. math:: + + \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}} + :param numpy.array u: fields + :param float ratio: desired ratio of the eigenvalues, default is 0.1 + :rtype: float + :return: beta0 """ - u = self.prob.field(self.m) - v = np.random.rand(*self.m.shape) - t = v.dot(self.dataObj2Deriv(self.m,v,u=u)) - b = v.dot(self.reg.modelObj2Deriv()*v) - return 0.1*(t/b) + if u is None: + u = self.prob.field(self.m) + + x0 = np.random.rand(*self.m.shape) + t = x0.dot(self.dataObj2Deriv(self.m,x0,u=u)) + b = x0.dot(self.reg.modelObj2Deriv()*x0) + return ratio*(t/b) def stoppingCriteria(self): if self.debug: print 'checking stoppingCriteria' @@ -167,6 +202,10 @@ class BaseInversion(object): @timeIt def evalFunction(self, m, return_g=True, return_H=True): + """evalFunction(m, return_g=True, return_H=True) + + + """ u = self.prob.field(m) phi_d = self.dataObj(m, u) @@ -198,7 +237,8 @@ class BaseInversion(object): @timeIt def dataObj(self, m, u=None): - """ + """dataObj(m, u=None) + :param numpy.array m: geophysical model :param numpy.array u: fields :rtype: float @@ -220,7 +260,8 @@ class BaseInversion(object): @timeIt def dataObjDeriv(self, m, u=None): - """ + """dataObjDeriv(m, u=None) + :param numpy.array m: geophysical model :param numpy.array u: fields :rtype: numpy.array From 1439052172dc2048494feadf5a63d4c8cb032a3a Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 22 Nov 2013 15:39:06 -0800 Subject: [PATCH 14/25] Multiple RHS using mumps backend. --- SimPEG/utils/Solver.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index c451a7ca..5861c146 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -187,7 +187,7 @@ class Solver(object): if factorize and self.dsolve is None: self.mctx = mumps.DMumpsContext() self.mctx.set_icntl(14, 60) - self.mctx.set_silent() + # self.mctx.set_silent() self.mctx.set_centralized_sparse(self.A) self.mctx.run(job=4) @@ -202,12 +202,20 @@ class Solver(object): if len(b.shape) == 1 or b.shape[1] == 1: # Just one RHS if factorize: - return self.dsolve(b) + X = self.dsolve(b) else: - return mumps.spsolve(self.A, b) + X = mumps.spsolve(self.A, b) else: - raise NotImplementedError('Multiple RHS not yet implemented with mumps solver.') + # 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'] From f5e3a4171de770e92ffccf264203f189e74dd956 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 22 Nov 2013 16:11:43 -0800 Subject: [PATCH 15/25] Raise an exception if name is not provided when hooking a lambda function. --- SimPEG/utils/__init__.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index b1cb47e8..622557d6 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -23,7 +23,10 @@ def hook(obj, method, name=None, overwrite=False, silent=False): If name is None, the name of the method is used. """ - if name is None: name = method.__name__ + if name is None: + name = method.__name__ + if name == '': + 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): From 4bd4c3cd8184a688a45ce1912e8bebf99c342ea3 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 23 Nov 2013 00:45:19 +0000 Subject: [PATCH 16/25] bug fix, we were sourcing the wrong file. --- SimPEG/GCEtools/startup.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/GCEtools/startup.sh b/SimPEG/GCEtools/startup.sh index 32164e29..2f23fab9 100644 --- a/SimPEG/GCEtools/startup.sh +++ b/SimPEG/GCEtools/startup.sh @@ -19,4 +19,4 @@ python setup.py cd ~ echo export PYTHONPATH=/home/$USER/simpeg/ >> .bashrc -source .profile \ No newline at end of file +source .bashrc \ No newline at end of file From 4ca6f090f6c026fdac6cbe537e66a7834f79b003 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 24 Nov 2013 16:07:19 -0800 Subject: [PATCH 17/25] no message --- docs/index.rst | 6 ++++-- docs/simpeg-logo.png | Bin 0 -> 23545 bytes 2 files changed, 4 insertions(+), 2 deletions(-) create mode 100644 docs/simpeg-logo.png diff --git a/docs/index.rst b/docs/index.rst index fb75f312..792f3be9 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,5 +1,7 @@ -SimPEG -====== +.. image:: simpeg-logo.png + :width: 300 px + :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 geophysical applications. diff --git a/docs/simpeg-logo.png b/docs/simpeg-logo.png new file mode 100644 index 0000000000000000000000000000000000000000..a0ec740f135721f5cf7d1c430ae2908c20ed4547 GIT binary patch literal 23545 zcmd>lg;$i{_q9qllG2UR3epWjcMUM4LkL56hjfTENXv-SFvLiwARW>rT_Ox2Aky-C z{Cwa4;$4fyz?yoVd(J&)pS|~u*VR@b!l%W5@ZbTFnkq>D!GlL8z~2kF*uZa5wE{li zfyP_O#M{8l-rMh$r`-d2TQ{g3vzp7R*LM1LuWSQ6hU}ysJTS{s1IZiue?KU|38t{S z9(Wz3b<7?$k){j6RC>HYkx{B&GCu(t_gr1nKAthq>G;Y|C5Dwj#{7G=BvM}q5*|gS zT9bERACz~0%#nB>j@!nSx3YZ|>MxsP;p=O^niH~i+m* zA+Gr^f%i}!nt8no~5M&4Yu5$$2x!bEKMxdvyjE(mjX>DwnP zOA^HmsuE>BrWq^u$7~bg5)sKtxqPqpPF4NQY3n3vo7hT%2P_Tz-3O{FY<6a5ISI$9 zk~S9If|kmxinYzsly{%ADo#&NKaOMz$o|;#vvTSl9!*-VdJ$zufc0(7Hg^}Z-OyG~ zgrx2EgMMk3F>8HZH0#L?DKe83i7z4aU&=-72ncehtEo9UiR9JvX`TDt*{yi9?`s}5XYsd||DUaaoHZ-j>d1o_7n0)e~=Ifx93FuGG%7jI}6M477o zu~Dx$(JoePFDoy1=+6~)@|~pN0Ah?w##{4Izrr^Csm{?BtN<$5w^`Pc+Q?c!fs6Un zTVkaqr9#u~Fug3ohp3QR>=QJRTSf|ztPWl@c(yLgn)>sur>7^=X}0f64-eaDWb;c;rQDbrbYE-Li#-&0ylAWU=Bgkf zbB9Eo!?=w^HqwDF#$j?g_J0d!WRpXtStupm53SSe84)X~R$WOE6azjj3f5#!w``v1^Ye2@cMIJNVO0SjM+q(p`qwC8XNI)jO|gv^|Mddb0N)iaA6F1LtZu6 zznyR?{ZjL*%XUwuit}>C0JDKaWcPE|7r!DyvKI?Be|p-eHiB#kDwCol1Bwkkl62{V zXj4`|5o#VYYE5CP5ec*;4H!~Exj_?;MDZZ zF7u(CGGf0T;&^X1Uv1IvbNJn;C?g}o2Qy6LX>uZMR1wi z7d(3HR4tph!HJ)27x-B1YP_|Y#P}H@?&H!cy>)eUyF!A3++wusD zP%!Vya}vCa&LyN-llF@jW5a?%LITsXioi~LX;hv#Ayl{i-elXj_YP9-;>#`Ztji{@ z@<0HVd1|F2bG1#9nVETS(du1M%4LU#mf@X_wVM3-^+$#6^+fr4L}99)9$lGMo-de} zDl9cSyUG3T`rP~E@&$vWui^R8isYoj%NW$WXN$St_F(E;#aQCKXMVs<{av7mGy7o3 zJ)6#{b>;U)xxoEX$rz+~f~A8W1FN0cM30ea7E1^;f;GVT{dpT0v6c# zbRBcgAL|f;J$X$C|5UcjiW||)2FGXGk$;mwkdjM~s zgrC!ZoUCjfbH@*IC*$)v*JSLa+d{oiCbr3SlD#cdZJW+5EY4=1xWvkOo|u2gpn;Ov z@Sl_7`jcC*gmM`B$yLjAMVsGYv!iuVs{H&5b;N@&pKCqQ1|l3`WQ8jhX@nHFT5zpX zfGt0vG<0|Skh>%na-B}y09lb`hUDF$9sD1HxY#uGwAMphOl(2YEoXx=_}V8;UB3cRdr zB@Nqg0k_b4)ujVW5W=9m@}Td-mvBF3$Qd()Nuhoi4OL`eVOfqk*@_shSAKndQPZ!3 z=Yksi(qkJv^F$!~eHx1)o_5)Ub=}DL&#^aM;3tY%j#}MHvu&4B1p_JnuN}vbA#JcI zJ-T9^37@>^Zc&Vkcw`yn5JEV4-AoM9_^}=B;^OO3s$NM#K=5quLz*$?NZ$|kxZ^qlzlpn|CLxg@K(5VK>3R(3og-^Q9^LmfHD^4ZxTYnU`BXD zeS7=u$j(T1(dqV$f`S6()V@e)5V&?HE2$=rS1o{i?2bF&P^Riy_LvCETHd zU{bwnF|fP4dlQeGlU=cV!30FYdfqEYN9swZ|EU`;;Z^r+&j3e&%JQk$-K?{PrYH8|Khzy&L`3v$?@SJSpFtjz2cRHt-Jii7rDZCu(G4$9-2RD zSWe~iGIxasT-6K`SDvBV+cT!q*uB`BE<;Q@(C`0%jigW+905T{vA)8p;u$|8&8x|1 z{QE-WWfy2u?4!C_BaC;4P;WdfGf_6|jjcyhI&X`)%*9|T+c=|4@PdN43MuBv{lt0- zru?!Y%J*R5g=K}I_Q1%9=&dLnv!CR4CDS(Sd$D4yYBf+(UBe-@K-Xg%5@ZdD?K_t6 zJ$SO(=6C!~yUf7r_s>Zo91_NQ1LFU#<}s!ofp(FMMM8a^C$OeH(&dT&r5|wDiT~2@ z-{UWtH&RdA2aZzg>#p+Ypg8`wHf${#phow1p^xywcI`UiAk3kp{tv*hZHp$XdZqJ`cY}h^0$4tB9>>8KK(fl(~VNiB@akP@B z@o#SZjh0an{xM! zFtM1MX%|Yh1gv$2a1eRW1D$R=@Zu<>3gD#H72Ioh7I;^Ave>()sTlk&7%o{_VpfL1# z^;z|QYRPv`ZNB?c?9}Esfs#jiXJ^OsPF%Loe-kgj-~W<9yDZ&;gjcPGDzB?1Z<-Qc>8dB1 z-~w!D_~h=~e;}Eu=M?h#6{Rzf}_rTd;BqW|O?DixyrK0KADdR12?8-WG`U&p*`;roN z;yp22G>>XG?I1XZx}=gkox7brzRLx1*VMvt+6W-bgq!28E}W=@!|5&6Ns~X@JdaXM z2v{Yn)mB=4wo>0ahjalivMbXrs7(C=Ey@63$rX7Ntl6oHwe~)ux zstPl;mGf$QfR?>|>9}z{DK|^hzC5mxZ2Li^H%`_k;mbaK$rTh`M^y64t%53{n_ zWIIheg9DU9Lp!Z9FaV@-$VfvYKDy!lyKlrdg{NuX(pQ><=^Ot?#MB)*$WA88Nj7Sv zjKgcW+2ipg!ySgcFtz~O$$SBgQlO>Tm-8O2e6du%(zaeFffA>uC1}=?kv$ck7alwr z>PjZVB~L{yN0V1AE0YJw@bfJTmXaH=)s(Wnh)hxkni003v-qs8x|WdSw6x2x-|x9J z1ijdr90*?-*N#Koj)k1M6G+Cp2H=wR=54NL(^`gA*KmFRHJ^EK(GzBf6AMI$8R8?7 zp%P2X&@2;jeamT7Jr{#c8QED*WmA!n25Owf>af|MhJG@mbfMyo&?AA5`;J0^hNMOE z^)itTE55faM=Nav%)u7S`U17e`0Pf~2L>{|o8)nG=5bkb+0_>+i~*p9Qmwqg@+Yx% zc@O5tu%z7ggrQSk%DCFEjdg< z-WQ5_K*#^y-8iv<#}%q%6!b3@3Kc3%jF%KsxWT4%U8um-654^kKLC7iGpa};{rmb` zf#hz!Oh|55yvYzVSdj%>!|YSB8v=`Q+E2Ia2uynN-1>_KDYA$KLKJAp7<##p<+9Ys zfKeAn%+)`V#}s}0L;~6yXPnJv_QyhA{>^y(k$iQlTyjgk3LXEPu}&~i=c`KfAPtoG=<%pEBmZSxw{0PK@6;>F?KASs|)z0`*^ zIu8-izI33+zgtw@$G(=nh;^2pYVq17G6}^<^J(Su>^AZlzl<_p@tL^1ynWI_F~1aX#1$`n31I7np&Sg;CMQux$}=A>hHVZgT%Al@t}@2 zT&Rx@6uP;F1f%5)orwYpR8y&olIg`yU#QN#R0T)&2P;&G#m$DXfX7)X_5-{WZaFO@ ztYs>7pOual$=;`Emz@PB68EZ?>sMrLLg*vdPc}<)oTXiwmWs61)WmTobMYo2@6eO4 zHyc%JapA7p>W^9iA4cD)$ka-*8eT`TnDR0^YHFFGg5|0xq=WthV?~_Na;ITuN*{bZ zOQ{qq*=?@;vQ0C%bwt#x@yOhXiI`qYx7BgFl!V1pL?GTsC~h{J6j7`yLk^8-_rIJ} zVOohx2*n4XpMwxA+uV;l8=gS>Hovv!hz?zb+!78wc{KDLcjz?NTCN(xJZ`K|)#^?} z=jqJ^Ld$G-KP5$mjpqpMO18VYPr^Nb>-89#HQM&P$^pR7_#c~=ef*a11a$$iXQ^2) z!Lh6iDv-pP6>4AH`Q@OLD&wU_7qc2Z8~|Ew6bAgr$tEKsn_>a$xOdWYEj_%da{=ZY z5@i*Ypo^X-5-;gJ8yA&C*kFAn(=PW#!w%}yoI`G2dhkoU1-cT!=EP!5%>JijN2Wwbbsn0y zutJ^axc8>F^+Sa8t>v6j^T68aXP{0e0GVy!(op&8)|hIF;|HVaZt6mM_@PC>;G<2` z51}sjE4U=GH7XL}-h@LYhoy5H5@if6m9GdQoLUw6hc@>S90Z=ToLpRcW^uEM4Y%BX zN}L+e5P^=N0PCB=FGkVOe~jiB%#ad6*iHVv z(bYX?XHMeU33~qUx(XXzkTzs%Ra#b7v)U2lw}+H@|Ni|mdZ1X-8J1LS(uB|Cwe=am zE?1HN+SY~}(l`adpMqE<<$GSFOE zOuHW$vQS2|)a1723tWMj`ZxxO6m#?0H_0|rUaCf&gQ39*TIXY;>|U2@AR$Ovu4z4I zE3pyAtS&ZqnJ38zFtRm2J$(zz?0Nt9?xXve;G+HKj-M4Xn|@@Hu$l7Xsi~=cYAPzr zDz_GZ{`vm1*t`z?>S>U=&qhyl+wjOp7`38k@TI$>gTv3S!0!6)!+jx#GSL82tIpM= zR9&a=5bt`?Y*!i)bw!-7RjYjG$W-?m%@eJeaTc-iv(SA1yclEn1?QskPi)RV7Eb;) zFY%$m4*N3|1}4E&SqwmvR>8!M-%43odt2@<=RMVW;o&>v7kG2B&N@9$R~R2_Xa{wT&Ycw-hD!hwcS_Im04nDD&mmUc%VVf_`a1M!t$zp=ufOa5q6kBC1ajC_0&Bk6?=)qOwPkwptOkA-F@8RZm+sfqtS`eROUuEb!mRky$)RruFFLGRBUj)E|Of?dlj4oji)El-vRYM%V5 z_S^o}aT||Zjb`b^#z(sjLnUhIf7!|6S0fxU{Sm?l`J1PA592}bgXI=<^Qrnd2QU(E zgy=o~1ueEH<}89#D7L^aP6N1xbZK2uou5;mMhFg$i;-TT(<8d!^^ou6tVvl({Hf$VsHH;U&g#>jVw)~5g(n%yQ*g%+qqWO6Zc*z}kY z?=$Ak=%-1L&fu#GfYF@yU$V2aYXksIK~LlE`m$S(MW|Obtk+YUmug&qn>&z+(k!V< zR!~s?d!-Q~^m4VjsY8)GvE}ym79$`aUV zdSxB$S4E~fUBKAe9lkcM*uYCtV_j6Z0;Ca|PQ-bTi>qRh`9{nQyT!cog@}jNE~abNu_~*AfC4Y3b?rNHAjHRO%Q!N?x6OuW(s#y7X93!L8;BayAbT%OVfkcvXzTA zy}v|J<}T6Z0m}*Ao4t5aP0?tEhHk}@ zFxiL6bY83mna&U6925>uAItTLz8-kTaIdMOBT$_f0^KEXd^KJuowvog&yd)+IpmgU zkuI4MS`lqd6ls39gwZk#`Xm`LmMg4jzrVw}*dA~uih)hI4ykuI;jcb>P6EEmK+s+a zf8@V-J($8$z#Ou*%#iSA$Os1857PzDmkSf^NTNmPv4_m@hP<6fb6~zp z%=}rk93jO{!@Mf@EZP7)D9DlUC*NY&ouNq)QpiT=x2CnwwrKzn5Tp66 z;SI3a$|3Gg@Q?a1sSyxQ%VI@Y?cKAd0`=32?nPOAW(hb%v@_9TNk3xDEy8a{Bg|>I zc40(_UbLO72B#@7fj}O`OKxmb9vq36$oz?Kuq}`Z@Z{vKP&}&sT0Nw~6jHB~C*q>e z?yYTN{y}u>&fdMcnV#g5V0T}khE17`qm9cLSdU~D;~k&6R~dKq8M`$fWZV(N;R~A4 zQS6wQm@}|8oBp>n4|0yaowYBm{A_l;7|RoNQ?M8@^7%1d+AP{!b_*JA6>zZnL_vM8q9-|`}p~ItKuv$a{~Vwx@I>C z5iIcy&9UE>2;pSn?mPTd%yzIa=?FY8L4$Z=KqHDKh@PV5+XJrZeKqc$fj71c;kL*( zwgD(wx<{4)OEn&oX-A*tQk@1NjUMpB$W&Du!nx^NOhDBqY=#LAbCJ+q17?c3+gmSe~&8 z{HB+2=GQPNo)k%_ItQzU`O_7>X^q_noBbZv3NKzkN)7&oY%??dG+d5I#7ZKT8Y=%zxzCS9|kv3v50m3 z(qGoEnFz+ceDP3r2B`L0@oF*VS`yFqB1vx zQ&QD1`bJvuEVOdZVIZC;=Jf7Zo?l54a@Lk!3eqa>WaA<3bQVLGK7}B{ zaXemd-a-B%|2P-&?Q&!WM|-jGzH-Yc=RAlYkt+q+CD&*P9_0T^3HNynZy<(e>cD-x z{1iOP#tM?$eBPi&Y8FO}`UI}htnUVmh;kRxw8!0ocgWsKu3SBL>Z6qo^e+3PkzL5e z&0RU>TjptmG1vNaz}%?Gx8O6hFX4<$`K=foFa?b6d9UgRkJ@{AtE;OoU7hXj8iT=e z?B#h_>1}r}5eO%Dy6wuLi#XjKQw(#Y0p3GT7q}rbK0U5#Oi4k7-KhYzK$$}sQ;Xb* znIOPUyl#EP5aly7YpjF%6z$xBH+lU8ex9c9yOC+nvY0zHT%L(VSyyKIOS<4{&H7E@ z(Qmdo&R}V~?0a_rc!}WJhBmpbnDlLv-_*OU$%xqNf&qvnme{B8H2wxjiM9<6zAa#U z>)6{`6{;{?F(d-egUf6Z5)zXe&3Df*0Q?LB>i19}tfRHxxw?{Lr2vPF(Q8G_pMsQrdf5yBMyS*-2r>tLixfVfs+-^p0AD&u2-Jz3riNrv3s}0;qjtjt>f=H5QS<8N0|S zYf+jAF>=J?oyepHnTGO{YPwfxw48;=Qgr2O;1-#@GWDWH0WhTD1&FPJQxJ(eR|p5K{y$>$}kBHU?Y>*alX&nKeuauWDX4tb<|$6 zV71;9Sl5+(*bH$fGWn!*7dg97?`W{|c7p}Z)IQy=tfozlNqT!b*(Y+7uC z(;_`L7ot195;?zcEKzFKuUdEHDq{Zi_Uf$0pEUNo`Kdl<`j9tJ@RpS7jJ0%Yw>?`( z1!u7{EvegWGa1t>Z|JcVQ+w=3Qq7R!@{q&*@~Npr?KkVWVkaVXVhp>*`aK{tQQX|XEC@je zZ#cbxbF2gOBnuUUQtIht#bZlt=Vhip{%<8cOfI~PbWVa5`-tRc#GT*Ps)%3Ue6i&d zXAvgl)oI09^w(TXkcX1Uz8AlUm@rKqmH&3vEj-{+7%f({yC1OrO`_Mk;R?_QysW1p z06sZJzI4GHM$!^h`Ed;Tc(7`Y2mkX``>=x)y(PHO#Qt-gV|y5SAHiun-z~$|xWpGv z&0qDR+%V|#7W8JMofpz6Cp?+^!uQ}k7b&ug1@e&@;;2=QuAMTHTk{nOFpp#bV&rxr zWz_-z?`T6NGEjl%Hx8jAe6ZKIuakoo-n*4Q0)1jJEGc+J_R7b`#s$06TG-^P>uQcp zZLRfcCph^3e(@mVy@*cmIX`b)be?a$mMDfgaeg-Wrcmxl!D$fapwzejprbQ_8W&SE zVVbXtzQ>L_(I$h_0B?>GMS>v0)@pk@+DD4yn+Nf**8b9r=eJj^kbzgkfp@5_MsHaG zUb+nc>d6rf9gGYaqDCJ{tO~4FZ(qBaw-c$dOug{K9N{q410?9e>1vC%wr_>fH)Np1 zM{;=ZBg~P*Zr(-Ji)S@QX06_PlW;g|6U5B?#2D{LO8RHSSl520SP8ImUcY ztVg_{iJ`asqekuWMktBQ$iDRvz#GfD3g*yzt6kSwt9{a+1Sa7&!*nhx(b+;FT7l2a zTrG?F4Di`dYKUdwK-5o(oL2d2*Vr!{){hFavB2Og`Fg|I`{TgL1hhaGHv~&dy*?o^ zF$sVUw~X-2evk$`@3l&t`k7ds$#>IpQahWBc&7;A<#Y{d__=7@oDC_Gngv3)@8qfw z;lSM-sIeAwq?khCX9^3uM?{&la@&UI2ccwWE278yiStAwQm(pVP^c-|*yB$PD zp!_1@FS`j@Jm@(fl+LA8?0XnoNRIu0M$hJcg_K2jx1Q^9a1^b8O4|iNJ@<*={xOh~ zG>E6IB97_^zl%x8essK9BUn`K72S|zOWMLhb7!*b-^(WSe2v@7f!xeYJ1d|Utjb4X zQ*N2CNt<`@{rC3IAvb@7rl+TO4^}(!@Hr#O#S480INC4~y<_7+Q__Q*LLGM$3}{a0 zk(JOA9})lRwllf9Balvyw6^DwpHZ51K!{2&l+&QNnc$V+IkT~_f}23wEO(Tu9na_w z%}+z#ytOD#hP2W^r;H}sN^JP~ML#bjH8u4V$P+(N+N2J1*}>|-k22rjwfXLG^s-z; z!+*qz+*Fj5I!=F2>vR=lXJ>zp1EGye#8vOhH?HHU5_b>)u4rEbo!`KFX}!O8ZS=AZUEij@Yjs@At*c6O~#zY5jTkXCtfUB_9=mxD+vV|5?; z`3}uQn!IbH66Fxh^THsE(fRr^w|%WoD_~)E{tH}qJ=cL!v3~g0u%4Bb6%&A$j>I`o z5~A9y-dK-X25+!Mr;%favwqGu^dp}$;m5f=RR+39<0fr&9Z!7!ZlU(@*IqHwB$L+< z;OAgh*Q0$RI+4qB)h<(0(|#ZH93TBl2VH^4?q70zo!t2Fb8Zk?siXh^^!KaHn2C&%dMuIdOq_@)EcrnmMX=j0XbdxOdc6=HHw#uP+TaZbD0b&QSjTrLBmYH6qHIvTBV&iqKv(b$agwRE*ho$d^VX@kP1;HL3I2Kj0eSvDM#a z1l)Y2_W%*nzzWtQY#gP%;WMx4)6z<*q@?stpz{(c075^Eo5SFr3|uNw;^E;*Gz0d= zHLy35zMhnL`iKtudd1BV#?1w(-spvgcmbSR5+KqE>!m7@c704cWBxmE&i?8Kz=Mq5 zUZ9IFoLNPM*v-4{6k)UO5j*Z;`TvA6A)kH2#$~e_!_bt%{9EvQ<#C-rgR`qh=c!5IQIhj7@2Q zsmJrk;(Uz{i=fL%1QfXt$wHgU>O6<>)J_ot4=g;nxjUNkS%sJCB45M>vIZX$R!dgnp!%wPTOyf()nZf=bQBEWjbo<1Za z9&oxH1C$^=F6Ea#`$Dz(TG=!2Ih~XJR+<-G&uDVX7S=G)w7}>f9`xscj}&=GihRlf zxsP%n@{ONM<>Z7AoU8AY=A%If8W7KIr2`&jyDp~)kH6?ZIzAqrohe`AtN z8n%e!wmq%=;v-`FrX_wSQ}V>95TB}T*-m;ioWpi#eCJ_nz*1@Wwt zDDl{bM87C`u3JpB$5U z+bNM*7mjW+8z&+8m5;lW;pLj}N?!gUdB*1yL@JQmO{SmipWf-#qvmU@WdEKFFj8tvs2g z%D$#mAmP)BP(lRGzxv2jekL&q)nVxGCr*jYvrI${2{N?v&fD*3EKbu-m+6F!&Pw?B zs55_DeDcQ(o$M!akk0m+Z;4JP$VS|GO zIg-b?y1gsl><6!vY<4!k1$3n-&#VE~Km?y98TduD+F8a(U`|DG&&TdWk!2Nb`^jKZ znRHd9!jK;rBJbc;Rikn(@|T9Ph%tK`_&E^>y##{`i?x>6s(p@#_{}*D|1_W7I6i6~ zWplr$Zf>5SXuHRr665$(^Zj~Ny%Ckk?40Jup7%(i*`rF_>SwO*n;Zv*NQXNHxX*+R$cY@}Bn0Al3>Ri_D`SqH9r(j+8j_Zl)~HZ#DGSh6 zJ$E-326*%S)@*j4%No+2@dW3_JaNeQrPt2ItDUNvAq|*NeUB360IY0=OObHlX#-U* zJeb&U>)Y7CBTZ4hs((axJo3R1+sbkEwo^z=jXbPl>)i!vbL)Pe0za(ph1(WZJFS{q zCgkjN>$J(-(Z{ZaU)Ol*0n&ho#mtU0)0o+a&rWdhaok)oDN-6>?3%}!{xHG&Z0>~s z$>S?$;5*|ALtYv517q2L(ZkMp2UYONk!i}z5bZZ|CXDp-UN2ZSiVY5zWJ9V9%92^Y_c5TDyqSF?+Kx^E zZriKVosn2>iUf|p2AM7MST8p`u{RAq%IK!=g-$n7%k0K-t+}IHbpX};i8M1CkvLa$ zf4jQ^EMC{R)w_+&zC>kuMcV($3+rm&8$|B227WsYs7T%#^AKR;a}c}T)YC-F!ffDx zL}EI>fylN)&EvDNqTm|pb1GRG=BKgMijB#X^74AQ>+7~3^>b`9RUzWh3REe-_y7SnKI$dDia2`DsYE=qvXn z+(Ew9o{&lAR=~rc| zH+rNNKpr7%vrzkLi{%+Xt@F(p5C@&OInTH(EBb+G0&0_X8ak=peT^v139Zi;%Q?%UW$-K0*XsNagYnb6Y^(SM)$Hv!yM znSM?v(cc>(MUvDcfP-%dqEk~%iMSyA#p>a@9KZk-jW}N~mt(#x(J!-6@ag;}5(jG% ztVKQVFt{%3&RwS(^URZeIC*VIO*om*_Wi#SDk^GBBA-Ll0>PZjeIR~?J92!J=kVfN zk?eRYX|`N_mR!BWhKXhW5?eE^j6N_W6QrHLq$<;!kd*iD2b{eqfch|h%)l@vTs*nl zz=%(UK!dy+UOwu+#+bYwT!Ld2lTm||8W4tv^ie@mxEAXxYMMN6Z=5JtYA4EbZR(Ro z;=a{EemkeRR|{+X=4VnZRRvAg(Fp_}M9a9|UcQI zl~ZLAAZcD+hd;zv?LE=#cXA}^tnSj7$}&h*d8v@=VH6n?XlU+!KCLu}{rvnQzSy*7 zYsf$F4Smw>YBMGpdy+e2$W}lrM{oy{IA^8a)1%!qIA&2S^2GZ9uj9Mx01H|C85Anf z8pTbc&=;y2oBD18-4&^Ocpjs>qFz7TGlCjc9FdK6Xb5+ZPR%b|u92z=wsQ?W;%{3x zawhyL+Kfl=T`K32b0ZUQXrKXZ)fGC;YzhgJd_OumDpdlEs({r5@IJ-_p*k0x z-yJTuoPCrU7q|xwCX~;x(b84%f?IE1WvcfTSAFB-_H1jm*STHYD-{!V!;I)-RQ3ug z$4l(%(E{vL(-&)DR?aVJ9jE&H*-d7(PkrEax)pyzMP8h^G)(%J1h^h4>wpqNzOiTH zWW+>spp7a4SI*D|zr~OBu6N9by6^nKYa+_@s6K+%|A)|{ltGDG}nEs z^>*PF!#93qaskA{&(H4wB)YSlj>zTy*Y4>Z<1Z-`C%$@+{)?DQzf6QxwYWe3R+Hv z2pAVU%4cE3qIO5s{qOR8$kVKamAhN@>H&*EqGmleR)n$JQ%G|ZCTH|!t}X|PUsjcp zj?j~C%8_NDLFh?RaU!+ApS<95PRKl@J|-qBhB(3j8STLE+jue|DTxv&JndMR9|epl zbpr*WxLcONp>q__;R&3d^Lr8A0i2)I#d zL_Dx4887_t)VHs#qtg%kMAqvVxqkXGy`)V7j7(dQ($dmGULKyaT`4<*&vmtr{EFBQ zwy$}$W@q1Hp=s?u=iL9ubHitbR4F+TgU|r(NoFyT+&n#NlKZf2$HH>>R^1M6#2Wt< zjmv3R-Y)WC^E+0K1l*8Y$)qV+`QzS(j;GH@hfL)ROQ6eM3UbPpYB;F<{78x%V1c0T zV6mwkCaS$F&D;{-akm3~an? zY+qw!t>jUgsR9No=1E4xdWPGNli~pZ8M}920O$K5BT?ZusGYQ9H`F4kYL;hiZ2U7_ zmqq?Ax$X@aFR{~e6-D6`>lNktd{U%$+?@NMeo(eUNYp(vbH*2U7Ju32pod{#Bd_XP z#MfEXY+wc7`-}Awn`tNJM@Em~CoAnL`Mt5d>vHJZVfvN%1o{5`lRm({c^oTdp4Z+Q z^r=E0nrSBR-mugwrQ&$_>#$ceaq(rp&F$SAFDaR)>5vlj!mWm+z1A38$6v4c*b|9# zSiH`fR?`xs;Hi`Eygu6PVJ}1t>>P#qAAWDxiNqoKmi+FW>+z5DsJj&Vs}}3St;$xV zmBXX<$g}$eA+jn+!klB=T*f~}3fAzQnwk&tTf;zNQB#3gzcC*(4t#J9H=?p^>}0O~ zf`&&c@VL*E>vgc;S^P?C(D#Yfr5WY+YL&MeEXPG$>Qa~UAZ0ivl045<^uuMS2VZSAG2Mm?`Z*D(&)utQp|rjg zlD0=HIDjH^9F#r;T8Yxqs*N*C0MFM>)dEXjjDZiur&{muPE}+_q7Rjol<-uN85uBg zTB!ws+VmxnMtL^G&x;U*%pYb|EY7EE>Fg#T@|aAH!nd1+ z?+nc&51n@^?M<7J)-9X0{hvxpE%||D9syTzXk@tq8}rHt3066gI(EIfd199NQLGC4 zG4CI-M|A$5Kaae_TYHSeNU5}X@6fkB@HUzIEs9Esnr|SeOqz=C;XvXcG`?M-`d!?+ z381Sdop0!Pt9>c@tsLT&H{<)Q`iomI{G{M^E7}9*hNQ{5)L=*wN0|r9Jx_y zF07w>g4dk~K3l0fr-josDYhK?EKEz6J)?{3{{k@0Eti*<4}p>I@JK#WJrfx8r)(Gr zp9iuzb4ERtCVdy?WmIPtYC6wriUP^nkG{8oHMDl&tVG>|t_v{{@%UMC%%r)RoacbW zV|efPeX#=`i)tys5cQlY-xL!7H?$8vKUR7i(8|2RAFY-eenHNxjVBU6=3q!iFw8ol z-p~!9dSCtZ22~zWo~0g{w$hTSIm_fc^HbY3v@rrOD zd6IK~kyF5w0S21U(%UT7-Ym=%;g^i2U>v9RUe~YL}zQ_rx#a@ z2K{8K&04@&x#dkqjAYnFb#1{qa$!wND-7BD-=`x95)%?&K{j_Z#+3V;R4G;C%%%lj zdwOmKYqEjmGOd{}dW?(dO0q{hbdAH#qg|}~^2#ZbTd{|45;e%z8nF;zwHRKwoBr+I zkk&MpU=nrd@5hu#U;xBPSKGY-t86kQBV(p2?EddpBJ8Gta~m(AzY4}eleqp271^=q z+U5@szu7u)6>wH;j5T#|WsZhQ|Ab}tFY~9)!5jKry(^lxzXhscglaV@X=(5jRrZ3A zF-PF^VLu_s!wgjI^FqXpM=jupyc*=GX-C!R6v-U6^BU7{!(MK~G9-(#M|NTWyTx%? zd|`Ym_WZh+Vhkb8JY0#OLbOrT`b+pKATv%M`UZu!p=;+5_4V~iJ`;gHL~{w5hexqV z$~H9sH*_b&%R4H-^sDj}&s-F6dhM;=Gm7(23%P1#x#|TezM3CCrXbT;L^fA6nVy?O@th`7(lB%DCfq)|>MZGaheyNue>zgPf- zSnaQhG+&84F2ZV2|K&bw+#(VDC9(l3Ki5ef4{{wRemHb3#Ezbwn3#CrR+tai0HQUE zWbQui76!Q*$;ruq_&O7Iid8xuq*N`YR9QB6BeI>Ko4xTXfYqt zh{G9M=;n98abQ{yLM1pjxM0<}Hj8J_h<2NMGW{Glj+FzfoDaa~{E(2b6j}L>rvnbR z?GOpQy>w$X-8c^sRDMUl_!ydrUBT@xQ4l7aH>PNyuP?FglcXS##cO&h^vx|-j`
>hg zP&!aji(r8E%*4PUrjK&DV;GuNnKDJsz5$XstR0Aoghy1<-aDs8Nf26t_c+RV%mQ#+ z{K?O!Z;H7*TI&j<2z5$#5m`zkMbP@VyvnO7{q*UOiFE>`I0MMo_ zh|!gwoQ0X#++l&|FS3CG!BQ*71if0qolrIwb@=Obkg9Rz0>)i+>DIEG%RyH@@*br%*p4t83a?mXT{I*-j<)5CvGq}q>3V>d z7St>u`R*~`Z4?`%3QCq<9d%t^B;rKVkZuOi^4C8SwAp63j5CAIF00BT|HR`1la zBdv$svmUjE&^D8#o#pSLD}d;hIFB^DV8zeU2|iC@G0p~@pW+4E2@VPCO;;--R>TaQ3^B0VB@6DqurYb0ML6cTVey z{~8}NqSb-Lc3_FFLW?gys^M08<9dl5YlNCJ$bRgt^Oq*8if1rd)i0kCv+03!LZ%-B zmMHF6+m48c{FP|R%lGDo;pRl&fr6kcaTb-m5dNELN{J~2G0tr4MomHCe;2PG%j6#uhAFG#@emdaWN@XKo>asY(I9rIh@IZM~dw7 zu*@}~U$XoRIBF$NBsHVWO+#2@nN5S3O;?yry-uck z(d_IzcxmSh#+0()HToXmKR*Hf?@-TzFbPA=exJtN-_{DpSF9Aw=Ott9ZaLG>_ zJ%_|8J`KWl^V{10dlD$}ZB8~zz^YR~R(5T}>FU9d;O2Nep!kK7(a>B~ece#q*!Oeo z92xat27dhuV==pHsVrN7NysuJ2ADsr zyQYO^KCC;v5_|~n`ym%X3_qs=0Xj?~me*$*W|*WMF@fQ0jOD`^!$Xr=Yxx*3!sXZH z$C(+!|10FX&ZuiOtZdiFc5$z9&FrjDWMpKI>?k8;Wn9<1mvE(U z6+%~CGx^wkT#|G%<9B?%zsL9AuRlHRYu>HQW}c?M*%UpAOH?^(5ZQ za^D?Q)2uSqV?Eg1rkixO7-;xC68U#8{vU;a)2>u$;*SBe>iyo{=(HfNuyqF#y6t+X zTUo=4Xw*8U(V4VFB&u{hj0B*xa}TU-^vUsCIw&M0#H#Z1B&2?0v1WOc^YBscN}Hdw z^*{phZ@d3o_k2nU%ZqiMSkFf3`Bobq^++Pm{#TAB%KIhWwY~~8VFB7!ml0@X9}Y1O zI%AHu%#cEd+XX2D$IVWokBy*hVnayP-u~+{Zq+{%F3fctb9|YXdjRozV|#5`ZfxklesgaQ zrTbAV{TNsQ_W@vSiafV_ni#zTjxM}P`0`sp*~OiyLwjNf`NXEz^?1o+OQR*HRtZw? zl~MjLrnr=#xdu(7KtsV0@uOE~scimm_1Oe`%FCC-M9kl<98SB5!4Fkn;ZVkb#K%)Upc$N25s^w1@t9KrPp~BpxQdL%n6^qmZ#Fb{w(p? z7nfEo55zDRE~K7?Q(#u9V`y?sZS6yke@80-uLW-M+;oFDfvWnIR2V}#mym|dK~BqL z8&q@ybfhK5sPWdfd}QOT$d3Onw7|*c17%WAOjUMvD#f;pS2*iQ{jYLN02+UWaVV%1 z4J{*_Dy2y$|F_@}=<>-{z6p#7BmfACZ~pLCX9oDHOEpuqB$=+AiEk73LLl-G(n6Hv0Sn8Q}xQK#<}CL zl9;Nfhjpn);VdRW3!IRvq#pnv`wm-Hj9v}Gv2s>^@+)zkgxPq@Q>Cna^C}LN1)dcM zyu1}BZb<3*j~UR#=MZ+nuE;F9)S7#zM{jrj6DZFcGa;Q9mzN-KL8Wj6=BOCCXUDTD zK0;Po>Hg#$;qHa;1}E{Etf?~#gL`KtaMjV@!%JSC|EU$@9be=-NGd4?yg%X*udD~Z z#M+_PP%?=_Voz*vleIzcs-{4#Ov` zr>Ac*deN70iRo00zI+Uu)G;sY&2!7LZ6at+&9uYXePXAO{Pye-w?YC!x_)N=6sq8; z@FrftSZ!0J&9!z;;!pa@=cEu??$9)QDPY6t*`y1TgR3rd`KQ~bHPjXF^#N^=sa;+0JN(0@-)Rl1vzrwzsv8CuzC+SvS3V#iTa0E~k_x<+@nzd0>Z~+IWkeghrbK2rqaq7j zcWF>4bYrZ{voUKzl)cdM#8h<4^;Ag*=Ufq%=ZNJa1Q@!K2oP4l0zR`t?tBm(3E z<<-WNC`Qyw^n%-*8;~j_e;WNX$Bv`z2JCixWNFUVJ~qAoj280S`-l-qYz?l_jQJ{r z8J=wR5Qy_sr6cp$a9!9n(?KFfzR%AS6w>V}6FP>Dq8NI5DiRM)_UOwq=cHJTvkMnD z_4*-&l|hGJ5b-#MVDUQxwgSA3)($XtSOwT7<(cx{gWxai7H319<9aeRAq5XFw zMWxB`baCY0*nPD4?c}PTAXK^3?mG|!YL z%FY)1QP}%L{s6^Owo}OB7Pv^ewnY1molS30Da|_6b9XQ1b{ zd=dy8+We4deh;MyQ*Mzh5B`LmJ#|b=bu_+4p!6FrcvPUl*ae zbQ;jybPMG-7^wFP3s;tZ0k4T1sw5~rBacM#2^A}jlQ7R=~J;v&@jzQeNiG*8@M#^9{4nYHA9Z}c?t zS}x#2bJ|1?r$ztqauqaohLhR%rT!#-jCTzH6AxUQF(bi~_L<@Ikzu7}*9Ch!;%l4R zx7%eNT>OPbK;t7!uETTtEtAU$7MEq`x>_`=)Dn>EkC68MD;*#w25>;J>^OBze8SuY z6{uN`QgzEhK{}5MClm+jh35Y(!lbs;Y(qVj^ski?68Y(9Pa8o@s*p8R2ou#a3-q>S;9x}Aw*|WCa&qF1H1H66DNO7@k07RTw<8T+ zmY*UHx`m^YL6{qDcHeSVcDaZT%M*&b7i-uSS9upx?pC>Q>@P5buLTtSK2*f}f*09J z2Jc0CKyZjErz{Oh*v}+5q>cnKqeYH+;#`IjFV+%g<7a~mC*LV3Jx^z*?Cwp-Yhxjx z_hIlBlOAZ*LLegJJ*8qjRawc`WrfH}wy<+#IIB0G4g0R=zQy|foo=~{zr7+r6?~3z z1$AXd41ir$SGZWjP;B#}h7VlYHjUJEddm|qS{oKq`W%)q`VLRN+wI*bmTGH3w6?QYMv;X z5Sk;+JqW!1r}eK^~gHnG>8RoUgF%G zbdFNLBTZJil+&+knt#2bOH5>O$mPMe-4~#f)@h87_b4Lps#9x_{YxjT;teX=#Zb3* z>Tx^uL1CSdGJp0-Zokb9*BOR$l85q7Q$@Y-sRP%rUAW%0GK(@(|JWsi4$ZU z?y7x8#3fpVcQnCVU&OWW;efWh8Jgb8^K< zx@7dop=#{nnS_~R3VJf23K*67*l}8!*hlQRNOqiQcTdmJW-q6y3UmJy`twHTl!1u9 z67D)y-#O>QLyf8@>^0}4OT12(u)N>*@gogEk{2np#99|1ool>zyu{edf`Tw+8j>v@ zVfG@Bc5QL>(g$?oko#8k0=Gg7Cks%!0HW|XphdbK_9EzYkmVT5c~~L5GN{3_a0tyQ zEe}I$qEHdA=SF-n7I*vuW*c&<dIg)ZtgBtveL@JaiKUqmWku!9hsJ{>8w^q`~I## zC*Z0OK{1$Ze{ZV^h`k(~RmE!ktdBlO#JbkukWqm$t64kT%r4}p=C5g?*ugG!4PXf- zeU%}TJjp?Pd_l|?PnKtR(xyr3RBG*^pj8hKmwHId>Lf6Wa1f!{k zZ|QNo#0BGr&Q*FMfEWtnEUoC8(Jk#R4O3D3yZJnEn!lKClDST%#GLv?pCt!q8XASf zBjb?k!^6XB4Pa&ngkzaX55b6S8Ums`T%Y8d&bt+=!^;5QkSFIorgCp-AXjogZr=T$ zzC%-%cxS$mtw(YPEPu2Q81*PTdK4DBoJh_+Ax{Hwp3;f}rJ&V83@*42aC~3+#q}j2 zXo9Z2jJZNi)pX}Ldko1@2BT8Rg4}#CJ1u1<0*U%KV6L|;vE)`ze~wS_Eoep_O*+rr zWBo-nZtxRE^@s2+K;UBtLTX%!;6vS(AJIN%Z;|9#WH_-;t$>aGVhl;M!9dV@j=we- zc1NMyOQzgQv`)29oxda6 zw|S4yZ_xzg-)A{F&2HdaW7_wax;=CZnaHcj1AgZn(sfe=-wA(xQY#>2!#3mkP+}~+ z+ZjJW%Spp?o6`n?MR&luDx*t4fjMsl5v1JzOS@&&w!yFU%t)XTJw1Ky_Mg{L=V6X8 zrg&Z7&|U25KkKtXx&6I$BZ`97Y#I&08VzEI^ccT;=uJcilGl!dmlMVRq*fNL42V?i zr$}EH`^dJows3i$iA*gk0nHr&HuZystJkh=a-b6W*=s56wLSyTc~}fZAO=eR&UrkU zQ{iWBX4QgwbDNp^zcb~hum#r4enMu1EXlig=TNw{!rws=mePth_3!_hDuo_ygiyte z^Otjb)CCkJ@{(>cS;{btWxR((>1Rz~ne+Zx-6MXw!lm`fIc6|+`pRP~+bYOUepbZn zKpb2~=I+KGjatnv271Q=24j;3zs#D|?)wt&exaI(lPR=YF-t9WZHGX78t#O2`VgTS zMx6IoEbb{ZJj&=)p71SWqE|e`bq4!FvUtv+w7m>iVpz^G!9|9NbQy0A{pT)(a1M8x zqxEu=P<5x0Q|P6xpe^zt92b4qn*AWYEgj9XAJZ+|{vPRDuIsOTNtgUf3%#-byFaI1 z;NXDTq0n`(_Q!N%X1_*Zqw8C~j5x}wrO#I2%;5#%>{0!ZiFDPEiI?qK$P;u-wJnM0 z%1a%}~Fu(=$hPDfzFJ+KMlnNeo;lW2*=9n!6&(>R%9% zMfL>|*S_&YwT8Y9Z9cji^5<5Zj_~bP!&uEg+5Bg>GPXk>9A8O} zO=5P!YmcLp8b&EFbwcAPeU$u9O)^t`%kJ9?rs}(28OK{&I{YRErS=IMY{mRv+gn?= zU(Rl92X3o&x`z96U_d%Bwlg4i@tlJx*KH`oDQG;xd5|!+7=o>lv+|L;J6(k#wltC> z!o$ON2LJ9a$QmoVnW2gGX_EJ3Z#HGoo;~29i+4!w36CXvrq9x8I1yy(@LA|>RpHG_TVwh`@jARFn%i`^VIa!UFI&g; zzn_(cJukF;(N5?Q#nl;lsTlw0`2S-jt;VKe`--b+gHJhKj66}h4oc}XlONQAX%el; z&^$nW{HQ`v)X=MYO1!=EWal`w-~Vc67aLpX$_YA*lwf35d0yznHI$Zx{VZ{_2c=!?1T(#a>=~rnwxUlMRLT~6(z;a{{Q`qo*n7&i}H675l16fMMkXG zYQ-+=J25hP@BI64=l(MF$q+8OPky*viSD(X;?1qD#UXOlT0IlBckeus)2X1dwnjQN zg*xTSrG5W)#ka_J;>0H0vG{AQp6LXg@A#h#CY}4@Mt;AVA@BG4D|`%~-U5h@G4)(^(Jt7S2n;0l5bmX+RG-VatVfy7iPO4(PPy35*vD>tOPZdJ-P)TUQ;@`x z26p23hUUHKHzQnzuB6XW-7H{@EVZ9Mtp>YouMX>b-iaM}V1VXE>-WxT&9tx?UowPR z%cz|zV!C~b?kx{v-lZ02eygReET^@Z)$A5@;Ae4PdNcZ) z^u=3kocOJg-n*^Xxa^hu%SiE;Aej6A{v( Date: Mon, 25 Nov 2013 14:59:52 -0800 Subject: [PATCH 18/25] Make tensor mesh copy the h and list so that the outside h is not affected. --- SimPEG/mesh/TensorMesh.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 30eebc13..4c0841b8 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -33,14 +33,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): """ _meshType = 'TENSOR' - def __init__(self, h, x0=None): - for i, h_i in enumerate(h): + def __init__(self, h_in, x0=None): + assert type(h_in) is list, 'h_in must be a list' + h = range(len(h_in)) + for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: # This gives you something over the unit cube. h_i = np.ones(int(h_i))/int(h_i) - h[i] = h_i assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) + h[i] = h_i[:] # make a copy. BaseMesh.__init__(self, np.array([x.size for x in h]), x0) assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" From d81582cc504c78fdfd46e1b928937a514b205032 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 25 Nov 2013 15:01:12 -0800 Subject: [PATCH 19/25] Slice through a tensor mesh using a video. --- SimPEG/mesh/TensorView.py | 43 +- notebooks/SliceThroughModel.ipynb | 1139 +++++++++++++++++++++++++++++ 2 files changed, 1178 insertions(+), 4 deletions(-) create mode 100644 notebooks/SliceThroughModel.ipynb diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 3ae91259..0374d1dd 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -2,7 +2,7 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib from mpl_toolkits.mplot3d import Axes3D -from SimPEG.utils import mkvc +from SimPEG.utils import mkvc, animate class TensorView(object): @@ -14,7 +14,7 @@ class TensorView(object): def __init__(self): pass - def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): + def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False,clim=None): """ Mesh.plotImage(I) @@ -141,7 +141,9 @@ class TensorView(object): C = I[:].reshape(self.nEy, order='F') C = 0.5*(C[:-1,:] + C[1:,:] ) - ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) + if clim is None: + clim = [C.min(),C.max()] + ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T, vmin=clim[0], vmax=clim[1]) ax.axis('tight') ax.set_xlabel("x") ax.set_ylabel("y") @@ -196,7 +198,10 @@ class TensorView(object): xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] # Plot the mesh - ph = ax.pcolormesh(xx, yy, C.T) + + if clim is None: + clim = [C.min(),C.max()] + ph = ax.pcolormesh(xx, yy, C.T, vmin=clim[0], vmax=clim[1]) # Plot the lines gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) @@ -336,3 +341,33 @@ class TensorView(object): ax.set_ylabel('x2') ax.set_zlabel('x3') if showIt: plt.show() + + def Slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None): + assert normal in 'xyz', 'normal must be x, y, or z' + if ax is None: ax = plt.subplot(111) + I = mesh.r(var,'CC','CC','M') + axes = [p for p in 'xyz' if p not in normal.lower()] + if normal is 'x': I = I[index,:,:] + if normal is 'y': I = I[:,index,:] + if normal is 'z': I = I[:,:,index] + if clim is None: clim = [I.min(),I.max()] + p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1]) + ax.axis('tight') + ax.set_xlabel(axes[0]) + ax.set_ylabel(axes[1]) + return p + + def SliceVideo(mesh,var,imageType='CC',normal='z',figsize=(10,8)): + # First set up the figure, the axis, and the plot element we want to animate + fig = plt.figure(figsize=figsize) + ax = plt.axes() + clim = [var.min(),var.max()] + plt.colorbar(mesh.Slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim)) + tlt = plt.title(normal) + + def animateFrame(i): + mesh.Slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim) + tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i]))) + + return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)]) + diff --git a/notebooks/SliceThroughModel.ipynb b/notebooks/SliceThroughModel.ipynb new file mode 100644 index 00000000..bc8548fb --- /dev/null +++ b/notebooks/SliceThroughModel.ipynb @@ -0,0 +1,1139 @@ +{ + "metadata": { + "name": "" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from SimPEG import mesh, utils\n", + "sz = [10,20,30]\n", + "M = mesh.TensorMesh(sz)\n", + "mtrue = utils.ModelBuilder.randomModel(sz,its=20)\n", + "M.SliceVideo(utils.mkvc(mtrue), normal='x')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Warning: mumps solver not available.\n", + "Using a seed of: " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + " 155\n" + ] + }, + { + "html": [ + "" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 1, + "text": [ + "" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file From 623b14c05c76fb5940fd5266a4496f233f4c0696 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 25 Nov 2013 15:02:11 -0800 Subject: [PATCH 20/25] face divergence operators for each direction. Bug fixes for boundary conditions. --- SimPEG/mesh/DiffOperators.py | 71 +++++++++++++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 2 deletions(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 4e017c7c..d384ca17 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -152,6 +152,73 @@ class DiffOperators(object): _faceDiv = None faceDiv = property(**faceDiv()) + def faceDivx(): + doc = "Construct divergence operator in the x component (face-stg to cell-centres)." + + def fget(self): + if(self._faceDivx is None): + # The number of cell centers in each direction + n = self.n + # Compute faceDivergence operator on faces + if(self.dim == 1): + D1 = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]), ddx(n[0])) + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fx', 'V') + V = self.vol + self._faceDivx = sdiag(1/V)*D1*sdiag(S) + + return self._faceDivx + return locals() + _faceDivx = None + faceDivx = property(**faceDivx()) + + def faceDivy(): + doc = "Construct divergence operator in the y component (face-stg to cell-centres)." + + def fget(self): + if(self.dim < 2): return None + if(self._faceDivy is None): + # The number of cell centers in each direction + n = self.n + # Compute faceDivergence operator on faces + if(self.dim == 2): + D2 = sp.kron(ddx(n[1]), speye(n[0])) + elif(self.dim == 3): + D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fy', 'V') + V = self.vol + self._faceDivy = sdiag(1/V)*D2*sdiag(S) + + return self._faceDivy + return locals() + _faceDivy = None + faceDivy = property(**faceDivy()) + + def faceDivz(): + doc = "Construct divergence operator in the z component (face-stg to cell-centres)." + + def fget(self): + if(self.dim < 3): return None + if(self._faceDivz is None): + # The number of cell centers in each direction + n = self.n + # Compute faceDivergence operator on faces + D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fz', 'V') + V = self.vol + self._faceDivz = sdiag(1/V)*D3*sdiag(S) + + return self._faceDivz + return locals() + _faceDivz = None + faceDivz = property(**faceDivz()) + def nodalGrad(): doc = "Construct gradient operator (nodes to edges)." @@ -279,12 +346,12 @@ class DiffOperators(object): elif(self.dim == 2): G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0])) G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0])) - G = sp.vstack((G1, G2), format="csr") + G = sp.block_diag((G1, G2), format="csr") elif(self.dim == 3): G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0])) G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0])) G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0])) - G = sp.vstack((G1, G2, G3), format="csr") + G = sp.block_diag((G1, G2, G3), format="csr") # Compute areas of cell faces & volumes S = self.area V = self.aveCC2F*self.vol # Average volume between adjacent cells From 5b4270560c3db885f99829467dbd1758be1ef167 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 25 Nov 2013 16:51:48 -0800 Subject: [PATCH 21/25] renamed video commands and added a more generic one. --- SimPEG/forward/Problem.py | 5 +- SimPEG/mesh/TensorView.py | 35 +- notebooks/SliceThroughModel.ipynb | 4059 +++++++++++++++++++++-------- 3 files changed, 3017 insertions(+), 1082 deletions(-) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 304f3912..54c623ee 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -226,7 +226,7 @@ class Problem(object): """ return sp.eye(m.size) - def createSyntheticData(self, m, std=0.05): + def createSyntheticData(self, m, std=0.05, u=None): """ Create synthetic data given a model, and a standard deviation. @@ -238,8 +238,7 @@ class Problem(object): Returns the observed data with random Gaussian noise and Wd which is the same size as data, and can be used to weight the inversion. """ - dobs = self.dpred(m) - dobs = dobs + dobs = self.dpred(m,u=u) noise = std*abs(dobs)*np.random.randn(*dobs.shape) dobs = dobs+noise eps = np.linalg.norm(mkvc(dobs),2)*1e-5 diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 0374d1dd..4b1c4fc3 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -342,7 +342,7 @@ class TensorView(object): ax.set_zlabel('x3') if showIt: plt.show() - def Slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None): + def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None): assert normal in 'xyz', 'normal must be x, y, or z' if ax is None: ax = plt.subplot(111) I = mesh.r(var,'CC','CC','M') @@ -357,17 +357,44 @@ class TensorView(object): ax.set_ylabel(axes[1]) return p - def SliceVideo(mesh,var,imageType='CC',normal='z',figsize=(10,8)): + def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)): + assert mesh.dim > 2, 'This is for 3D meshes only.' # First set up the figure, the axis, and the plot element we want to animate fig = plt.figure(figsize=figsize) ax = plt.axes() clim = [var.min(),var.max()] - plt.colorbar(mesh.Slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim)) + plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim)) tlt = plt.title(normal) def animateFrame(i): - mesh.Slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim) + mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim) tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i]))) return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)]) + def video(mesh,var,function,figsize=(10,8)): + """ + Call a function for a list of models to create a video. + + :: + + def function(var, ax, clim, tlt, i): + tlt.set_text('%%d'%%i) + return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim) + + mesh.video([model1, model2, ..., modeln],function) + """ + # First set up the figure, the axis, and the plot element we want to animate + fig = plt.figure(figsize=figsize) + ax = plt.axes() + VAR = np.concatenate(var) + clim = [VAR.min(),VAR.max()] + tlt = plt.title('') + plt.colorbar(function(var[0],ax,clim,tlt,0)) + + def animateFrame(i): + function(var[i],ax,clim,tlt,i) + + return animate(fig, animateFrame, frames=len(var)) + + diff --git a/notebooks/SliceThroughModel.ipynb b/notebooks/SliceThroughModel.ipynb index bc8548fb..1259ebf6 100644 --- a/notebooks/SliceThroughModel.ipynb +++ b/notebooks/SliceThroughModel.ipynb @@ -14,31 +14,16 @@ "from SimPEG import mesh, utils\n", "sz = [10,20,30]\n", "M = mesh.TensorMesh(sz)\n", - "mtrue = utils.ModelBuilder.randomModel(sz,its=20)\n", - "M.SliceVideo(utils.mkvc(mtrue), normal='x')" + "mtrue = utils.ModelBuilder.randomModel(sz,seed=786,its=20)\n", + "M.videoSlicer(utils.mkvc(mtrue), normal='x')" ], "language": "python", "metadata": {}, "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "Warning: mumps solver not available.\n", - "Using a seed of: " - ] - }, - { - "output_type": "stream", - "stream": "stdout", - "text": [ - " 155\n" - ] - }, { "html": [ "" ], "metadata": {}, "output_type": "pyout", - "prompt_number": 1, + "prompt_number": 2, "text": [ - "" + "" ] } ], - "prompt_number": 1 + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "models = [utils.ModelBuilder.randomModel(sz,seed=786,its=i) for i in range(40)]" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def function(var, ax, clim, tlt, i):\n", + " p = M.slicer(var, normal='y', index=5, imageType='CC', ax=ax, clim=clim)\n", + " tlt.set_text('%d smoothing iterations'%(i))\n", + " return p\n", + " \n", + "M.video(models,function)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 5, + "text": [ + "" + ] + } + ], + "prompt_number": 5 }, { "cell_type": "code", From 1088c09864a22d8be6a7c817b4d9648040a0e87a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 26 Nov 2013 14:35:08 -0800 Subject: [PATCH 22/25] Bug fix for integrating over the model, should have a volume term! --- SimPEG/regularization/Regularization.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 3d22a297..41938152 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -22,22 +22,19 @@ class Regularization(object): @property def Wx(self): if getattr(self, '_Wx', None) is None: - a = self.mesh.r(self.mesh.area,'F','Fx','V') - self._Wx = sdiag(a)*self.mesh.cellGradx + self._Wx = self.mesh.cellGradx*sdiag(self.mesh.vol) return self._Wx @property def Wy(self): if getattr(self, '_Wy', None) is None: - a = self.mesh.r(self.mesh.area,'F','Fy','V') - self._Wy = sdiag(a)*self.mesh.cellGrady + self._Wy = self.mesh.cellGrady*sdiag(self.mesh.vol) return self._Wy @property def Wz(self): if getattr(self, '_Wz', None) is None: - a = self.mesh.r(self.mesh.area,'F','Fz','V') - self._Wz = sdiag(a)*self.mesh.cellGradz + self._Wz = self.mesh.cellGradz*sdiag(self.mesh.vol) return self._Wz alpha_s = 1e-6 From a6e2686f09e6deb5b364e1bab1e61f97117b73f3 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 26 Nov 2013 17:30:43 -0800 Subject: [PATCH 23/25] Added doStartIteration to the code, set beta0 to None, and generate it in the evalFunction when fields are available. --- SimPEG/inverse/BetaSchedule.py | 4 ++-- SimPEG/inverse/Inversion.py | 38 +++++++++++++++++++++++++--------- SimPEG/inverse/Optimize.py | 16 ++++++++++++-- 3 files changed, 44 insertions(+), 14 deletions(-) diff --git a/SimPEG/inverse/BetaSchedule.py b/SimPEG/inverse/BetaSchedule.py index af4d883d..c225c0c0 100644 --- a/SimPEG/inverse/BetaSchedule.py +++ b/SimPEG/inverse/BetaSchedule.py @@ -3,8 +3,8 @@ class Cooling(object): """Simple Beta Schedule""" - beta0 = 1.e6 - beta_coolingFactor = 5. + beta0 = None #: The initial beta value, set to none means that it will be approximated in the first iteration. + beta_coolingFactor = 2. def getBeta(self): if self._beta is None: diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index b5e4906f..160678aa 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -17,6 +17,8 @@ class BaseInversion(object): comment = '' #: Used by some functions to indicate what is going on in the algorithm counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things + beta0 = None #: The initial Beta (regularization parameter) + def __init__(self, prob, reg, opt, **kwargs): setKwargs(self, **kwargs) @@ -77,7 +79,7 @@ class BaseInversion(object): """ self.startup(m0) while True: - self._beta = self.getBeta() + self.doStartIteration() self.m = self.opt.minimize(self.evalFunction, self.m) self.doEndIteration() if self.stoppingCriteria(): break @@ -113,6 +115,26 @@ class BaseInversion(object): self.phi_d_last = np.nan self.phi_m_last = np.nan + def doStartIteration(self): + """ + **doStartIteration** is called at the end of each run iteration. + + If you have things that also need to run at the end of every iteration, you can create a method:: + + def _doStartIteration*(self): + pass + + Where the * can be any string. If present, _doStartIteration* will be called at the start of the default doStartIteration call. + You may also completely overwrite this function. + + :rtype: None + :return: None + """ + callHooks(self,'doStartIteration') + + self._beta = self.getBeta() + + def doEndIteration(self): """ **doEndIteration** is called at the end of each run iteration. @@ -135,15 +157,6 @@ class BaseInversion(object): self.phi_m_last = self.phi_m self._iter += 1 - @property - def beta0(self): - if getattr(self,'_beta0',None) is None: - self._beta0 = self.estimateBeta0() - return self._beta0 - @beta0.setter - def beta0(self, value): - self._beta0 = value - def getBeta(self): return self.beta0 @@ -208,9 +221,14 @@ class BaseInversion(object): """ u = self.prob.field(m) + + if self._iter is 0 and self._beta is None: + self._beta = self.beta0 = self.estimateBeta0(u=u) + phi_d = self.dataObj(m, u) phi_m = self.reg.modelObj(m) + self.dpred = self.prob.dpred(m, u=u) # This is a cheap matrix vector calculation. self.phi_d = phi_d self.phi_m = phi_m diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 5bf3ef4d..7e934e13 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -142,6 +142,7 @@ class Minimize(object): printInit() while True: + doStartIteration() f, g, H = evalFunction(xc) printIter() if stoppingCriteria(): break @@ -161,11 +162,12 @@ class Minimize(object): self.printInit() while True: + self.doStartIteration() self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) self.printIter() if self.stoppingCriteria(): break - p = self.findSearchDirection() - p = self.scaleSearchDirection(p) + self.searchDirection = self.findSearchDirection() + p = self.scaleSearchDirection(self.searchDirection) xt, passLS = self.modifySearchDirection(p) if not passLS: xt, caught = self.modifySearchDirectionBreak(p) @@ -219,6 +221,16 @@ class Minimize(object): self.f_last = np.nan self.x_last = x0 + @count + def doStartIteration(self): + """doStartIteration() + + **doStartIteration** is called at the start of each minimize iteration. + + :rtype: None + :return: None + """ + callHooks(self,'doStartIteration') def printInit(self, inLS=False): """ From e2943ba56e1d244d440791200ae460cd7b72aa0e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 26 Nov 2013 17:31:42 -0800 Subject: [PATCH 24/25] Bug fix if no colorer is used in the video. --- SimPEG/mesh/TensorView.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 4b1c4fc3..9956565d 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -372,7 +372,7 @@ class TensorView(object): return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)]) - def video(mesh,var,function,figsize=(10,8)): + def video(mesh,var,function,figsize=(10,8),colorbar=True): """ Call a function for a list of models to create a video. @@ -390,7 +390,8 @@ class TensorView(object): VAR = np.concatenate(var) clim = [VAR.min(),VAR.max()] tlt = plt.title('') - plt.colorbar(function(var[0],ax,clim,tlt,0)) + if colorbar: + plt.colorbar(function(var[0],ax,clim,tlt,0)) def animateFrame(i): function(var[i],ax,clim,tlt,i) From 83bdf73a63a2bf4853eb6c3bd2e541480ed07845 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 26 Nov 2013 17:50:03 -0800 Subject: [PATCH 25/25] Added a travis.yaml file and moved requirements.txt for auto unit testing --- .travis.yaml | 7 +++++++ docs/requirements.txt | 2 +- requirements.txt | 3 +++ 3 files changed, 11 insertions(+), 1 deletion(-) create mode 100644 .travis.yaml create mode 100644 requirements.txt diff --git a/.travis.yaml b/.travis.yaml new file mode 100644 index 00000000..8a4f7cda --- /dev/null +++ b/.travis.yaml @@ -0,0 +1,7 @@ +language: python +python: + - "2.7" +# command to install dependencies +install: "pip install -r requirements.txt" +# command to run tests +script: nosetests diff --git a/docs/requirements.txt b/docs/requirements.txt index 45fa35b8..27ba63e9 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,3 @@ numpy -pypubsub +scipy ipython diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 00000000..45fa35b8 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +numpy +pypubsub +ipython