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/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 new file mode 100644 index 00000000..2f23fab9 --- /dev/null +++ b/SimPEG/GCEtools/startup.sh @@ -0,0 +1,22 @@ +#! /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 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 + +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://bitbucket.org/rcockett/simpeg.git +cd simpeg/SimPEG/ +python setup.py +cd ~ + +echo export PYTHONPATH=/home/$USER/simpeg/ >> .bashrc +source .bashrc \ No newline at end of file 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/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 5a71f6ec..160678aa 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,19 +1,24 @@ 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 +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 + + beta0 = None #: The initial Beta (regularization parameter) + def __init__(self, prob, reg, opt, **kwargs): setKwargs(self, **kwargs) @@ -22,18 +27,17 @@ 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'): # 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 @@ -68,9 +72,14 @@ class BaseInversion(object): @timeIt def run(self, m0): + """run(m0) + + Runs the inversion! + + """ 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 @@ -94,9 +103,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.' @@ -108,25 +115,42 @@ 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. 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. You may also completely overwrite this function. - :param numpy.ndarray xt: tested new iterate that ensures a descent direction. :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 @@ -136,6 +160,47 @@ class BaseInversion(object): def getBeta(self): return self.beta0 + 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 + """ + 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' return checkStoppers(self, self.stoppers) @@ -150,11 +215,20 @@ 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) + + 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 @@ -175,13 +249,14 @@ 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] @timeIt def dataObj(self, m, u=None): - """ + """dataObj(m, u=None) + :param numpy.array m: geophysical model :param numpy.array u: fields :rtype: float @@ -203,7 +278,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 @@ -244,8 +320,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 @@ -282,7 +360,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 @@ -294,3 +372,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 75209dc6..7e934e13 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -1,16 +1,13 @@ 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 -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 + +__all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters'] + class StoppingCriteria(object): """docstring for StoppingCriteria""" @@ -85,32 +82,28 @@ 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._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] @@ -121,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) @@ -142,33 +136,13 @@ 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) printInit() while True: + doStartIteration() f, g, H = evalFunction(xc) printIter() if stoppingCriteria(): break @@ -188,21 +162,17 @@ class Minimize(object): self.printInit() while True: + self.doStartIteration() 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) + self.searchDirection = self.findSearchDirection() + p = self.scaleSearchDirection(self.searchDirection) 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() @@ -240,9 +210,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 @@ -253,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): """ @@ -262,12 +240,10 @@ 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) - @count def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -276,12 +252,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) @@ -293,13 +265,11 @@ 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 printStoppers(self, stoppers, pad='', stop=stop, done=done) - @timeIt def stoppingCriteria(self, inLS=False): if self._iter == 0: self.f0 = self.f @@ -308,7 +278,8 @@ class Minimize(object): @timeIt def projection(self, p): - """ + """projection(p) + projects the search direction. by default, no projection is applied. @@ -317,11 +288,13 @@ class Minimize(object): :rtype: numpy.ndarray :return: p, projected search direction """ + callHooks(self,'projection',p) return p @timeIt def findSearchDirection(self): - """ + """findSearchDirection() + **findSearchDirection** should return an approximation of: .. math:: @@ -351,7 +324,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. @@ -365,11 +339,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: @@ -406,7 +381,8 @@ class Minimize(object): @count def modifySearchDirectionBreak(self, p): - """ + """modifySearchDirectionBreak(p) + Code is called if modifySearchDirection fails to find a descent direction. @@ -427,7 +403,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. @@ -447,9 +424,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 @@ -458,7 +433,6 @@ class Minimize(object): if self.debug: self.printDone() - class Remember(object): """ This mixin remembers all the things you tend to forget. @@ -509,12 +483,11 @@ class Remember(object): self._rememberList[param[0]].append( param[1](self) ) - class ProjectedGradient(Minimize, Remember): name = 'Projected Gradient' - maxIterCG = 10 - tolCG = 1e-3 + maxIterCG = 5 + tolCG = 1e-1 lower = -np.inf upper = np.inf @@ -546,22 +519,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 @@ -573,6 +559,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 @@ -614,6 +604,7 @@ class ProjectedGradient(Minimize, Remember): @timeIt def _doEndIteration_ProjectedGradient(self, xt): + """_doEndIteration_ProjectedGradient(xt)""" aSet = self.activeSet(xt) bSet = self.bindingSet(xt) @@ -642,7 +633,6 @@ class ProjectedGradient(Minimize, Remember): if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD - class BFGS(Minimize, Remember): name = 'BFGS' nbfgs = 10 @@ -738,8 +728,8 @@ class InexactGaussNewton(BFGS, Minimize, Remember): name = 'Inexact Gauss Newton' - maxIterCG = 10 - tolCG = 1e-3 + maxIterCG = 5 + tolCG = 1e-1 @property def approxHinv(self): @@ -857,10 +847,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) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 19a3ab34..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 @@ -318,8 +385,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 +404,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/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)" diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 3ae91259..9956565d 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,61 @@ 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 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)) + 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)]) + + def video(mesh,var,function,figsize=(10,8),colorbar=True): + """ + 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('') + if colorbar: + 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/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 3e429324..41938152 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): @@ -22,34 +22,33 @@ 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 - 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/Solver.py b/SimPEG/utils/Solver.py index 16db872f..5861c146 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 @@ -120,6 +128,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 +138,22 @@ class Solver(object): assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' + if backend == 'scipy': + X = self.solveDirect_scipy(b, factorize) + elif backend == 'mumps': + X = self.solveDirect_mumps(b, factorize) + + return X + + def solveDirect_scipy(self, b, factorize): + """ + Use solve instead of this interface. + + :param numpy.ndarray b: the right hand side + :param bool factorize: if you want to factorize and store factors + :rtype: numpy.ndarray + :return: x + """ if factorize and self.dsolve is None: self.A = self.A.tocsc() # for efficiency self.dsolve = linalg.factorized(self.A) @@ -150,6 +175,48 @@ 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: + X = self.dsolve(b) + else: + X = mumps.spsolve(self.A, b) + + else: + # Multiple RHSs + X = np.empty_like(b) + for i in range(b.shape[1]): + if factorize: + X[:,i] = self.dsolve(b[:,i]) + else: + X[:,i] = mumps.spsolve(self.A,b[:,i]) + + return X + def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50): if backend is None: backend = DEFAULTS['iter'] diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 2ba0f44d..622557d6 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -12,6 +12,29 @@ import Solver from Solver import Solver import Geophysics +import types +import time +import numpy as np +from functools import wraps + +def hook(obj, method, name=None, overwrite=False, silent=False): + """ + This dynamically binds a method to the instance of the class. + + If name is None, the name of the method is used. + """ + if name is None: + name = method.__name__ + if name == '': + raise Exception('Must provide name to hook lambda functions.') + if not hasattr(obj,name) or overwrite: + setattr(obj, name, types.MethodType( method, obj )) + if getattr(obj,'debug',False): + print 'Method '+name+' was added to class.' + elif not silent or getattr(obj,'debug',False): + print 'Method '+name+' was not overwritten.' + + def setKwargs(obj, **kwargs): """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" for attr in kwargs: @@ -19,6 +42,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 = '' @@ -61,9 +87,11 @@ 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 class Counter(object): @@ -75,7 +103,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): @@ -139,6 +167,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__) @@ -147,6 +176,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__) @@ -154,6 +184,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): 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..8ca74e27 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Optimize.rst @@ -4,6 +4,8 @@ Optimize ******** .. automodule:: SimPEG.inverse.Optimize + :show-inheritance: + :private-members: :members: :undoc-members: @@ -12,6 +14,7 @@ Inversion ********* .. automodule:: SimPEG.inverse.Inversion + :show-inheritance: :members: :undoc-members: 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 *********** 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/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/docs/simpeg-logo.png b/docs/simpeg-logo.png new file mode 100644 index 00000000..a0ec740f Binary files /dev/null and b/docs/simpeg-logo.png differ diff --git a/notebooks/SliceThroughModel.ipynb b/notebooks/SliceThroughModel.ipynb new file mode 100644 index 00000000..1259ebf6 --- /dev/null +++ b/notebooks/SliceThroughModel.ipynb @@ -0,0 +1,3048 @@ +{ + "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,seed=786,its=20)\n", + "M.videoSlicer(utils.mkvc(mtrue), normal='x')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 2, + "text": [ + "" + ] + } + ], + "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", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file 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