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/__init__.py b/SimPEG/__init__.py index 7f059a74..1ea0601c 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -2,5 +2,11 @@ import utils from utils import Solver import mesh import inverse +import visualize import forward import regularization +import examples + +import scipy.version as _v +if _v.version < '0.13.0': + print 'Warning: upgrade your scipy to 0.13.0' diff --git a/SimPEG/examples/__init__.py b/SimPEG/examples/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 9ddb5332..074cdb8b 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -1,11 +1,10 @@ from SimPEG.mesh import TensorMesh -from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms +from SimPEG.forward import Problem, ModelTransforms from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag, mkvc from SimPEG import Solver import numpy as np import scipy.sparse as sp -import scipy.sparse.linalg as linalg class DCProblem(ModelTransforms.LogModel, Problem): @@ -202,23 +201,17 @@ if __name__ == '__main__': P = Q.T # Create some data - class syntheticDCProblem(DCProblem, SyntheticProblem): - pass + problem = DCProblem(mesh) + problem.P = P + problem.RHS = q + dobs, Wd = problem.createSyntheticData(mSynth, std=0.05) - synthetic = syntheticDCProblem(mesh); - synthetic.P = P - synthetic.RHS = q - dobs, Wd = synthetic.createData(mSynth, std=0.05) - - u = synthetic.field(mSynth) - u = synthetic.reshapeFields(u) + u = problem.field(mSynth) + u = problem.reshapeFields(u) mesh.plotImage(u[:,10]) # plt.show() # Now set up the problem to do some minimization - problem = DCProblem(mesh) - problem.P = P - problem.RHS = q problem.dobs = dobs problem.std = dobs*0 + 0.05 m0 = mesh.gridCC[:,0]*0+sig2 diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index cf22baae..54c623ee 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -1,5 +1,5 @@ import numpy as np -from SimPEG.utils import mkvc, sdiag +from SimPEG.utils import mkvc, sdiag, count, timeIt import scipy.sparse as sp norm = np.linalg.norm @@ -37,6 +37,8 @@ class Problem(object): to (locally) find how model parameters change the data, and optimize! """ + counter = None + def __init__(self, mesh): self.mesh = mesh @@ -83,6 +85,7 @@ class Problem(object): def dobs(self, value): self._dobs = value + @count def dpred(self, m, u=None): """ Predicted data. @@ -94,6 +97,7 @@ class Problem(object): u = self.field(m) return self.P*u + @count def dataResidual(self, m, u=None): """ :param numpy.array m: geophysical model @@ -113,6 +117,7 @@ class Problem(object): return self.dpred(m, u=u) - self.dobs + @timeIt def J(self, m, v, u=None): """ :param numpy.array m: model @@ -142,6 +147,7 @@ class Problem(object): """ raise NotImplementedError('J is not yet implemented.') + @timeIt def Jt(self, m, v, u=None): """ :param numpy.array m: model @@ -155,6 +161,7 @@ class Problem(object): raise NotImplementedError('Jt is not yet implemented.') + @timeIt def J_approx(self, m, v, u=None): """ @@ -169,6 +176,7 @@ class Problem(object): """ return self.J(m, v, u) + @timeIt def Jt_approx(self, m, v, u=None): """ :param numpy.array m: model @@ -218,32 +226,19 @@ class Problem(object): """ return sp.eye(m.size) - - - -class SyntheticProblem(object): - """ - Has helpful functions when dealing with synthetic problems - - To use this class, inherit to your problem:: - - class mySyntheticExample(Problem, SyntheticProblem): - pass - """ - def createData(self, m, std=0.05): + def createSyntheticData(self, m, std=0.05, u=None): """ + Create synthetic data given a model, and a standard deviation. + :param numpy.array m: geophysical model :param numpy.array std: standard deviation :rtype: numpy.array, numpy.array :return: dobs, Wd - Create synthetic data given a model, and a standard deviation. - 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/Inversion.py b/SimPEG/inverse/Inversion.py index e3d500b1..b5e4906f 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,17 +1,22 @@ import numpy as np import scipy.sparse as sp import SimPEG -from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers +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 + + 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) @@ -20,14 +25,18 @@ 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.' + opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv()) + @property def Wd(self): @@ -38,6 +47,9 @@ class BaseInversion(object): eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5 self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) return self._Wd + @Wd.setter + def Wd(self, value): + self._Wd = value @property def phi_d_target(self): @@ -56,7 +68,13 @@ class BaseInversion(object): def phi_d_target(self, value): self._phi_d_target = value + @timeIt def run(self, m0): + """run(m0) + + Runs the inversion! + + """ self.startup(m0) while True: self._beta = self.getBeta() @@ -83,13 +101,17 @@ 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.' + self.reg.mref = m0 self.m = m0 self._iter = 0 self._beta = None + self.phi_d_last = np.nan + self.phi_m_last = np.nan def doEndIteration(self): """ @@ -97,28 +119,75 @@ 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. 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 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=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) @@ -131,8 +200,12 @@ class BaseInversion(object): """ printStoppers(self, self.stoppers) - + @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) @@ -154,17 +227,18 @@ class BaseInversion(object): if return_H: def H_fun(v): phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) - phi_m2Deriv = self.reg.modelObj2Deriv(m)*v + phi_m2Deriv = self.reg.modelObj2Deriv()*v 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 @@ -184,8 +258,10 @@ class BaseInversion(object): R = mkvc(R) return 0.5*np.vdot(R, R) + @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 @@ -224,9 +300,12 @@ class BaseInversion(object): return dmisfit + @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 @@ -263,7 +342,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 @@ -275,3 +354,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 37c8b296..5bf3ef4d 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 +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""" @@ -76,7 +73,7 @@ class IterationPrinters(object): itType = {"title": "itType", "value": lambda M: M._itType, "width": 8, "format": "%s"} aSet = {"title": "aSet", "value": lambda M: np.sum(M.activeSet(M.xc)), "width": 8, "format": "%d"} bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"} - comment = {"title": "Comment", "value": lambda M: M.projComment, "width": 7, "format": "%s"} + comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"} beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format": "%1.2e"} phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"} @@ -85,29 +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 = '' #: 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] @@ -116,8 +112,10 @@ class Minimize(object): setKwargs(self, **kwargs) + @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) @@ -138,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) @@ -185,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() @@ -236,9 +208,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 @@ -258,7 +228,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) @@ -271,7 +240,8 @@ class Minimize(object): parent.printIter function and call that. """ - if doPub and not inLS: pub.sendMessage('Minimize.printIter', minimize=self) + callHooks(self,'printIter',inLS) + pad = ' '*10 if inLS else '' printLine(self, self.printers if not inLS else self.printersLS, pad=pad) @@ -283,22 +253,21 @@ 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) - def stoppingCriteria(self, inLS=False): if self._iter == 0: self.f0 = self.f self.g0 = self.g return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) - + @timeIt def projection(self, p): - """ + """projection(p) + projects the search direction. by default, no projection is applied. @@ -307,10 +276,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:: @@ -338,8 +310,10 @@ class Minimize(object): """ return -self.g + @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. @@ -353,10 +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: @@ -391,8 +367,10 @@ class Minimize(object): return self._LS_xt, self._iterLS < self.maxIterLS + @count def modifySearchDirectionBreak(self, p): - """ + """modifySearchDirectionBreak(p) + Code is called if modifySearchDirection fails to find a descent direction. @@ -411,8 +389,10 @@ class Minimize(object): print 'The linesearch got broken. Boo.' return p, False + @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. @@ -432,9 +412,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 @@ -443,7 +421,6 @@ class Minimize(object): if self.debug: self.printDone() - class Remember(object): """ This mixin remembers all the things you tend to forget. @@ -494,12 +471,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 @@ -525,24 +501,41 @@ class ProjectedGradient(Minimize, Remember): self.stopDoingPG = False self._itType = 'SD' - self.projComment = '' + self.comment = '' self.aSet_prev = self.activeSet(x0) + @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 @@ -552,7 +545,12 @@ class ProjectedGradient(Minimize, Remember): bind_low = np.logical_and(x == self.upper, self.g <= 0) return np.logical_or(bind_up, bind_low) + @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 @@ -586,13 +584,15 @@ class ProjectedGradient(Minimize, Remember): def reduceHess(v): # Z is tall and skinny return Z.T*(self.H*(Z*v)) - operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=float ) + operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=self.xc.dtype ) p, info = sp.linalg.cg(operator, -Z.T*self.g, tol=self.tolCG, maxiter=self.maxIterCG) p = Z*p # bring up to full size # aSet_after = self.activeSet(self.xc+p) return p + @timeIt def _doEndIteration_ProjectedGradient(self, xt): + """_doEndIteration_ProjectedGradient(xt)""" aSet = self.activeSet(xt) bSet = self.bindingSet(xt) @@ -600,7 +600,7 @@ class ProjectedGradient(Minimize, Remember): self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient f_current_decrease = self.f_last - self.f - self.projComment = '' + self.comment = '' if self._iter < 1: # Note that this is reset on every CG iteration. self.f_decrease_max = -np.inf @@ -608,7 +608,7 @@ class ProjectedGradient(Minimize, Remember): self.f_decrease_max = max(self.f_decrease_max, f_current_decrease) self.stopDoingPG = f_current_decrease < 0.25 * self.f_decrease_max if self.stopDoingPG: - self.projComment = 'Stop SD' + self.comment = 'Stop SD' self.explorePG = False self.exploreCG = True # implement 3.8, MoreToraldo91 @@ -620,40 +620,229 @@ class ProjectedGradient(Minimize, Remember): if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD + +class BFGS(Minimize, Remember): + name = 'BFGS' + nbfgs = 10 + + @property + def bfgsH0(self): + """ + Approximate Hessian used in preconditioning the problem. + + Must be a SimPEG.Solver + """ + _bfgsH0 = getattr(self,'_bfgsH0',None) + if _bfgsH0 is None: + return Solver(sp.identity(self.xc.size).tocsc(), flag='D') + return _bfgsH0 + @bfgsH0.setter + def bfgsH0(self, value): + assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver' + self._bfgsH0 = value + + def _startup_BFGS(self,x0): + self._bfgscnt = -1 + self._bfgsY = np.zeros((x0.size, self.nbfgs)) + self._bfgsS = np.zeros((x0.size, self.nbfgs)) + if not np.any([p is IterationPrinters.comment for p in self.printers]): + self.printers.append(IterationPrinters.comment) + + def bfgs(self, d): + n = self._bfgscnt + nn = ktop = min(self._bfgsS.shape[1],n) + return self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d) + + def bfgsrec(self,k,n,nn,S,Y,d): + """BFGS recursion""" + if k < 0: + d = self.bfgsH0.solve(d) + else: + khat = 0 if nn is 0 else np.mod(n-nn+k,nn) + gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]) + d = d - gamma*Y[:,khat] + d = self.bfgsrec(k-1,n,nn,S,Y,d) + d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat] + return d + + def findSearchDirection(self): + return self.bfgs(-self.g) + + def _doEndIteration_BFGS(self, xt): + if self._iter is 0: + self.g_last = self.g + return + + yy = self.g - self.g_last; + ss = self.xc - xt; + self.g_last = self.g + + if yy.dot(ss) > 0: + self._bfgscnt += 1 + ktop = np.mod(self._bfgscnt,self.nbfgs) + self._bfgsY[:,ktop] = yy + self._bfgsS[:,ktop] = ss + self.comment = '' + else: + self.comment = 'Skip BFGS' + + class GaussNewton(Minimize, Remember): name = 'Gauss Newton' + + @timeIt def findSearchDirection(self): return Solver(self.H).solve(-self.g) -class InexactGaussNewton(Minimize, Remember): +class InexactGaussNewton(BFGS, Minimize, Remember): + """ + Minimizes using CG as the inexact solver of + + .. math:: + + \mathbf{H p = -g} + + By default BFGS is used as the preconditioner. + + Use *nbfgs* to set the memory limitation of BFGS. + + To set the initial H0 to be used in BFGS, set *bfgsH0* to be a SimPEG.Solver + + """ + + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + name = 'Inexact Gauss Newton' - maxIterCG = 10 - tolCG = 1e-5 + maxIterCG = 5 + tolCG = 1e-1 + @property + def approxHinv(self): + """ + The approximate Hessian inverse is used to precondition CG. + + Default uses BFGS, with an initial H0 of *bfgsH0*. + + Must be a scipy.sparse.linalg.LinearOperator + """ + _approxHinv = getattr(self,'_approxHinv',None) + if _approxHinv is None: + M = sp.linalg.LinearOperator( (self.xc.size, self.xc.size), self.bfgs, dtype=self.xc.dtype ) + return M + return _approxHinv + @approxHinv.setter + def approxHinv(self, value): + self._approxHinv = value + + @timeIt def findSearchDirection(self): - # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly - p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG) + Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG}) + p = Hinv.solve(-self.g) return p class SteepestDescent(Minimize, Remember): name = 'Steepest Descent' + + @timeIt def findSearchDirection(self): return -self.g + +class NewtonRoot(object): + """ + Newton Method - Root Finding + + root = newtonRoot(fun,x); + + Where fun is the function that returns the function value as well as the + gradient. + + For iterative solving of dh = -J\\r, use O.solveTol = TOL. For direct + solves, use SOLVETOL = 0 (default) + + Rowan Cockett + 16-May-2013 16:29:51 + University of British Columbia + rcockett@eos.ubc.ca + + """ + + tol = 1.000e-06 + solveTol = 0 # Default direct solve. + maxIter = 20 + stepDcr = 0.5 + maxLS = 30 + comments = False + doLS = True + + def __init__(self, **kwargs): + setKwargs(self, **kwargs) + + def root(self, fun, x): + if self.comments: print 'Newton Method:\n' + + self._iter = 0 + while True: + + [r,J] = fun(x); + if self.solveTol == 0: + Jinv = Solver(J) + dh = - Jinv.solve(r) + else: + raise NotImplementedError('Iterative solve on NewtonRoot is not yet implemented.') + # M = @(x) tril(J)\(diag(J).*(triu(J)\x)); + # [dh, ~] = bicgstab(J,-r,O.solveTol,500,M); + + muLS = 1. + LScnt = 1 + xt = x + dh + rt, Jt = fun(xt) # TODO: get rid of Jt + + if self.comments: print '\tLinesearch:\n' + # Enter Linesearch + while True and self.doLS: + if self.comments: + print '\t\tResid: %e\n'%norm(rt) + if norm(rt) <= norm(r) or norm(rt) < self.tol: + break + + muLS = muLS*self.stepDcr + LScnt = LScnt + 1 + print '.' + if LScnt > self.maxLS: + print 'Newton Method: Line search break.' + root = NaN + return + xt = x + muLS*dh + rt, Jt = fun(xt) # TODO: get rid of Jt + + x = xt + self._iter += 1 + if norm(rt) < self.tol or self._iter > self.maxIter: + break + + return x + + + if __name__ == '__main__': from SimPEG.tests import Rosenbrock, checkDerivative import matplotlib.pyplot as plt 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) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) + + + print 'test the newtonRoot finding.' + fun = lambda x: (np.sin(x), sdiag(np.cos(x))) + x = np.array([np.pi-0.3, np.pi+0.1, 0]) + pnt = NewtonRoot(comments=False).root(fun,x) + print pnt diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 6a9a8032..62ec4251 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -102,7 +102,7 @@ class BaseMesh(object): elif xType in ['F', 'E']: # This will only deal with components of fields, not full 'F' or 'E' xx = mkvc(xx) # unwrap it in case it is a matrix - nn = self.nF if xType == 'F' else self.nE + nn = self.nFv if xType == 'F' else self.nEv nn = np.r_[0, nn] nx = [0, 0, 0] @@ -228,6 +228,17 @@ class BaseMesh(object): return locals() nC = property(**nC()) + def nCv(): + doc = """ + Total number of cells in each direction + + :rtype: numpy.array (dim, ) + :return: [nCx, nCy, nCz] + """ + fget = lambda self: np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None]) + return locals() + nCv = property(**nCv()) + def nNx(): doc = """ Number of nodes in the x-direction @@ -288,6 +299,17 @@ class BaseMesh(object): return locals() nN = property(**nN()) + def nNv(): + doc = """ + Total number of nodes in each direction + + :rtype: numpy.array (dim, ) + :return: [nNx, nNy, nNz] + """ + fget = lambda self: np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None]) + return locals() + nNv = property(**nNv()) + def nEx(): doc = """ Number of x-edges in each direction @@ -331,7 +353,7 @@ class BaseMesh(object): return locals() nEz = property(**nEz()) - def nE(): + def nEv(): doc = """ Total number of edges in each direction @@ -346,6 +368,18 @@ class BaseMesh(object): """ fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None]) return locals() + nEv = property(**nEv()) + + def nE(): + doc = """ + Total number of edges. + + :rtype: int + :return: sum([prod(nEx), prod(nEy), prod(nEz)]) + + """ + fget = lambda self: np.sum(self.nEv) + return locals() nE = property(**nE()) def nFx(): @@ -391,7 +425,7 @@ class BaseMesh(object): return locals() nFz = property(**nFz()) - def nF(): + def nFv(): doc = """ Total number of faces in each direction @@ -406,6 +440,19 @@ class BaseMesh(object): """ fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None]) return locals() + nFv = property(**nFv()) + + + def nF(): + doc = """ + Total number of faces. + + :rtype: int + :return: sum([prod(nFx), prod(nFy), prod(nFz)]) + + """ + fget = lambda self: np.sum(self.nFv) + return locals() nF = property(**nF()) def normals(): @@ -418,13 +465,13 @@ class BaseMesh(object): def fget(self): if self.dim == 2: - nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0])] - nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1])] + nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])] + nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])] return np.r_[nX, nY] elif self.dim == 3: - nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0]), np.zeros(self.nF[0])] - nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1]), np.zeros(self.nF[1])] - nZ = np.c_[np.zeros(self.nF[2]), np.zeros(self.nF[2]), np.ones(self.nF[2])] + nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])] + nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])] + nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])] return np.r_[nX, nY, nZ] return locals() normals = property(**normals()) @@ -439,13 +486,13 @@ class BaseMesh(object): def fget(self): if self.dim == 2: - tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0])] - tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1])] + tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])] + tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])] return np.r_[tX, tY] elif self.dim == 3: - tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0]), np.zeros(self.nE[0])] - tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1]), np.zeros(self.nE[1])] - tZ = np.c_[np.zeros(self.nE[2]), np.zeros(self.nE[2]), np.ones(self.nE[2])] + tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])] + tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])] + tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])] return np.r_[tX, tY, tZ] return locals() tangents = property(**tangents()) diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 93b82b25..e22e12b9 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -62,11 +62,11 @@ class Cyl1DMesh(object): # Counting #################################################### - def nCr(): + def nCx(): doc = "Number of cells in the radial direction" fget = lambda self: self.hr.size return locals() - nCr = property(**nCr()) + nCx = property(**nCx()) def nCz(): doc = "Number of cells in the z direction" @@ -76,10 +76,16 @@ class Cyl1DMesh(object): def nC(): doc = "Total number of cells" - fget = lambda self: self.nCr * self.nCz + fget = lambda self: self.nCx * self.nCz return locals() nC = property(**nC()) + def nCv(): + doc = "Total number of cells in each direction" + fget = lambda self: np.array([self.nCx, self.nCz]) + return locals() + nCv = property(**nCv()) + def nNr(): doc = "Number of nodes in the radial direction" fget = lambda self: self.hr.size @@ -106,10 +112,16 @@ class Cyl1DMesh(object): def nFz(): doc = "Number of z faces" - fget = lambda self: self.nNz * self.nCr + fget = lambda self: self.nNz * self.nCx return locals() nFz = property(**nFz()) + def nFv(): + doc = "Total number of faces in each direction" + fget = lambda self: np.array([self.nFr, self.nFz]) + return locals() + nFv = property(**nFv()) + def nF(): doc = "Total number of faces" fget = lambda self: self.nFr + self.nFz @@ -236,12 +248,12 @@ class Cyl1DMesh(object): def fget(self): if self._edgeCurl is None: #1D Difference matricies - dr = sp.spdiags((np.ones((self.nCr+1, 1))*[-1, 1]).T, [-1,0], self.nCr, self.nCr, format="csr") + dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1,0], self.nCx, self.nCx, format="csr") dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0,1], self.nCz, self.nCz+1, format="csr") #2D Difference matricies Dr = sp.kron(sp.eye(self.nNz), dr) - Dz = -sp.kron(dz, sp.eye(self.nCr)) #Not sure about this negative + Dz = -sp.kron(dz, sp.eye(self.nCx)) #Not sure about this negative #Edge curl operator self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0) @@ -255,7 +267,7 @@ class Cyl1DMesh(object): def fget(self): if self._aveE2CC is None: az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') ar[0,0] = 1 self._aveE2CC = sp.kron(az, ar).T return self._aveE2CC @@ -268,10 +280,10 @@ class Cyl1DMesh(object): def fget(self): if self._aveF2CC is None: az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') ar[0,0] = 1 Afr = sp.kron(sp.eye(self.nCz),ar) - Afz = sp.kron(az,sp.eye(self.nCr)) + Afz = sp.kron(az,sp.eye(self.nCx)) self._aveF2CC = sp.vstack((Afr,Afz)).T return self._aveF2CC return locals() @@ -305,7 +317,7 @@ class Cyl1DMesh(object): elif type(materialProp) is float: materialProp = np.ones(self.nC)*materialProp elif materialProp.shape == (self.nCz,): - materialProp = materialProp.repeat(self.nCr) + materialProp = materialProp.repeat(self.nCx) materialProp = mkvc(materialProp) assert materialProp.shape == (self.nC,), "materialProp incorrect shape" @@ -377,7 +389,7 @@ class Cyl1DMesh(object): dFz = np.sum(dFz**2, axis=1) indBL = np.argmin(dFz) # Face below and to the left - indAL = indBL + self.nCr # Face above and to the left + indAL = indBL + self.nCx # Face above and to the left zF_BL = self.gridFz[indBL,:] zF_AL = self.gridFz[indAL,:] diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 598b392a..d384ca17 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -1,15 +1,16 @@ import numpy as np from scipy import sparse as sp -from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros - - -def ddx(n): - """Define 1D derivatives, inner, this means we go from n+1 to n+1""" - return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") +from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap def checkBC(bc): - """ Checks if boundary condition 'bc' is valid. """ + """ + + Checks if boundary condition 'bc' is valid. + + Each bc must be either 'dirichlet' or 'neumann' + + """ if(type(bc) is str): bc = [bc, bc] assert type(bc) is list, 'bc must be a list' @@ -22,7 +23,33 @@ def checkBC(bc): def ddxCellGrad(n, bc): - """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" + """ + Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + + For Cell-Centered **Dirichlet**, use a ghost point:: + + (u_1 - u_g)/hf = grad + + u_g u_1 u_2 + * | * | * ... + ^ + 0 + + u_g = - u_1 + grad = 2*u1/dx + negitive on the other side. + + For Cell-Centered **Neumann**, use a ghost point:: + + (u_1 - u_g)/hf = 0 + + u_g u_1 u_2 + * | * | * ... + + u_g = u_1 + grad = 0; put a zero in. + + """ bc = checkBC(bc) D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") @@ -38,10 +65,55 @@ def ddxCellGrad(n, bc): D[-1, -1] = 0 return D +def ddxCellGradBC(n, bc): + """ -def av(n): - """Define 1D averaging operator from cell-centres to nodes.""" - return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") + Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + + For Cell-Centered **Dirichlet**, use a ghost point:: + + (u_1 - u_g)/hf = grad + + u_g u_1 u_2 + * | * | * ... + ^ + u_b + + We know the value at the boundary (u_b):: + + (u_g+u_1)/2 = u_b (the average) + u_g = 2*u_b - u_1 + + So plug in to gradient: + + (u_1 - (2*u_b - u_1))/hf = grad + 2*(u_1-u_b)/hf = grad + + Separate, because BC are known (and can move to RHS later):: + + ( 2/hf )*u_1 + ( -2/hf )*u_b = grad + + ( ^ ) JUST RETURN THIS + + + """ + bc = checkBC(bc) + + ij = (np.array([0, n]),np.array([0, 1])) + vals = np.zeros(2) + + # Set the first side + if(bc[0] == 'dirichlet'): + vals[0] = -2 + elif(bc[0] == 'neumann'): + vals[0] = 0 + # Set the second side + if(bc[1] == 'dirichlet'): + vals[1] = 2 + elif(bc[1] == 'neumann'): + vals[1] = 0 + D = sp.csr_matrix((vals, ij), shape=(n+1,2)) + return D class DiffOperators(object): @@ -80,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)." @@ -107,6 +246,38 @@ class DiffOperators(object): _nodalGrad = None nodalGrad = property(**nodalGrad()) + def nodalLaplacian(): + doc = "Construct laplacian operator (nodes to edges)." + + def fget(self): + if(self._nodalLaplacian is None): + print 'Warning: Laplacian has not been tested rigorously.' + # The number of cell centers in each direction + n = self.n + # Compute divergence operator on faces + if(self.dim == 1): + D1 = sdiag(1./self.hx) * ddx(mesh.nCx) + L = - D1.T*D1 + elif(self.dim == 2): + D1 = sdiag(1./self.hx) * ddx(n[0]) + D2 = sdiag(1./self.hy) * ddx(n[1]) + L1 = sp.kron(speye(n[1]+1), - D1.T * D1) + L2 = sp.kron(- D2.T * D2, speye(n[0]+1)) + L = L1 + L2 + elif(self.dim == 3): + D1 = sdiag(1./self.hx) * ddx(n[0]) + D2 = sdiag(1./self.hy) * ddx(n[1]) + D3 = sdiag(1./self.hz) * ddx(n[2]) + L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1) + L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1)) + L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1)) + L = L1 + L2 + L3 + self._nodalLaplacian = L + return self._nodalLaplacian + return locals() + _nodalLaplacian = None + nodalLaplacian = property(**nodalLaplacian()) + def setCellGradBC(self, BC): """ Function that sets the boundary conditions for cell-centred derivative operators. @@ -129,17 +300,19 @@ class DiffOperators(object): for i, bc_i in enumerate(BC): BC[i] = checkBC(bc_i) - self._cellGrad = None # ensure we create a new gradient next time we call it - self._cellGradBC = BC + # ensure we create a new gradient next time we call it + self._cellGrad = None + self._cellGradBC = None + self._cellGradBC_list = BC return BC - _cellGradBC = 'neumann' + _cellGradBC_list = 'neumann' def cellGrad(): doc = "The cell centered Gradient, takes you to cell faces." def fget(self): if(self._cellGrad is None): - BC = self.setCellGradBC(self._cellGradBC) + BC = self.setCellGradBC(self._cellGradBC_list) n = self.n if(self.dim == 1): G = ddxCellGrad(n[0], BC[0]) @@ -154,13 +327,40 @@ class DiffOperators(object): G = sp.vstack((G1, G2, G3), format="csr") # Compute areas of cell faces & volumes S = self.area - V = self.vol - self._cellGrad = sdiag(S)*G*sdiag(1/V) + V = self.aveCC2F*self.vol # Average volume between adjacent cells + self._cellGrad = sdiag(S/V)*G return self._cellGrad return locals() _cellGrad = None cellGrad = property(**cellGrad()) + def cellGradBC(): + doc = "The cell centered Gradient boundary condition matrix" + + def fget(self): + if(self._cellGradBC is None): + BC = self.setCellGradBC(self._cellGradBC_list) + n = self.n + if(self.dim == 1): + G = ddxCellGradBC(n[0], BC[0]) + 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.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.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 + self._cellGradBC = sdiag(S/V)*G + return self._cellGradBC + return locals() + _cellGradBC = None + cellGradBC = property(**cellGradBC()) + def cellGradx(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." @@ -175,19 +375,17 @@ class DiffOperators(object): elif(self.dim == 3): G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fx', 'V') - V = self.vol - self._cellGradx = sdiag(S)*G1*sdiag(1/V) + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fx', 'V') + self._cellGradx = sdiag(L)*G1 return self._cellGradx return locals() cellGradx = property(**cellGradx()) - 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 @@ -196,33 +394,29 @@ class DiffOperators(object): elif(self.dim == 3): G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fy', 'V') - V = self.vol - self._cellGrady = sdiag(S)*G2*sdiag(1/V) + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fy', 'V') + self._cellGrady = sdiag(L)*G2 return self._cellGrady return locals() cellGrady = property(**cellGrady()) - - 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 G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fz', 'V') - V = self.vol - self._cellGradz = sdiag(S)*G3*sdiag(1/V) + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fz', 'V') + self._cellGradz = sdiag(L)*G3 return self._cellGradz return locals() cellGradz = property(**cellGradz()) - def edgeCurl(): doc = "Construct the 3D curl operator." @@ -265,6 +459,8 @@ class DiffOperators(object): _edgeCurl = None edgeCurl = property(**edgeCurl()) + # --------------- Averaging --------------------- + def aveF2CC(): doc = "Construct the averaging operator on cell faces to cell centers." @@ -274,17 +470,37 @@ class DiffOperators(object): if(self.dim == 1): self._aveF2CC = av(n[0]) elif(self.dim == 2): - self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") + self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])), + sp.kron(av(n[1]), speye(n[0]))), format="csr") elif(self.dim == 3): - self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") + self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") return self._aveF2CC return locals() _aveF2CC = None aveF2CC = property(**aveF2CC()) + def aveCC2F(): + doc = "Construct the averaging operator on cell cell centers to faces." + + def fget(self): + if(self._aveCC2F is None): + n = self.n + if(self.dim == 1): + self._aveCC2F = avExtrap(n[0]) + elif(self.dim == 2): + self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])), + sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr") + elif(self.dim == 3): + self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])), + kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])), + kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr") + return self._aveCC2F + return locals() + _aveCC2F = None + aveCC2F = property(**aveCC2F()) + def aveE2CC(): doc = "Construct the averaging operator on cell edges to cell centers." @@ -295,10 +511,10 @@ class DiffOperators(object): if(self.dim == 1): raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') elif(self.dim == 2): - self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])), sp.kron(speye(n[1]), av(n[0]))), format="csr") elif(self.dim == 3): - self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), kron3(av(n[2]), speye(n[1]), av(n[0])), kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._aveE2CC @@ -316,37 +532,57 @@ class DiffOperators(object): if(self.dim == 1): self._aveN2CC = av(n[0]) elif(self.dim == 2): - self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") + self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr() elif(self.dim == 3): - self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") + self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr() return self._aveN2CC return locals() _aveN2CC = None aveN2CC = property(**aveN2CC()) - def aveN2CCv(): - doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." + def aveN2E(): + doc = "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate." def fget(self): - if(self._aveN2CCv is None): + if(self._aveN2E is None): # The number of cell centers in each direction n = self.n if(self.dim == 1): - self._aveN2CCv = av(n[0]) + self._aveN2E = av(n[0]) elif(self.dim == 2): - self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") + self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])), + sp.kron(av(n[1]), speye(n[0]+1))), format="csr") elif(self.dim == 3): - self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveN2CCv + self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])), + kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)), + kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr") + return self._aveN2E return locals() - _aveN2CCv = None - aveN2CCv = property(**aveN2CCv()) + _aveN2E = None + aveN2E = property(**aveN2E()) + + def aveN2F(): + doc = "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate." + + def fget(self): + if(self._aveN2F is None): + # The number of cell centers in each direction + n = self.n + if(self.dim == 1): + self._aveN2F = av(n[0]) + elif(self.dim == 2): + self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)), + sp.kron(speye(n[1]+1), av(n[0]))), format="csr") + elif(self.dim == 3): + self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)), + kron3(av(n[2]), speye(n[1]+1), av(n[0])), + kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") + return self._aveN2F + return locals() + _aveN2F = None + aveN2F = property(**aveN2F()) + + # --------------- Methods --------------------- def getMass(self, materialProp=None, loc='e'): """ Produces mass matricies. diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index f0ac4ab0..5c2e7b5f 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -174,8 +174,8 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): def Pxxx(pos): ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0] - ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1] + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nFv[0] + ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nFv[0] + mesh.nFv[1] IND = np.r_[ind1, ind2, ind3].flatten() @@ -286,7 +286,7 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False): def Pxx(pos): ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nF[0] + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nFv[0] IND = np.r_[ind1, ind2].flatten() @@ -387,8 +387,8 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False): def Pxxx(pos): ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0] - ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1] + ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nEv[0] + ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nEv[0] + mesh.nEv[1] IND = np.r_[ind1, ind2, ind3].flatten() @@ -499,7 +499,7 @@ def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): def Pxx(pos): ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nE[0] + ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nEv[0] IND = np.r_[ind1, ind2].flatten() diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index b510a754..5c4a73db 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -45,8 +45,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridCC is None: - ccV = (self.aveN2CCv*mkvc(self.gridN)) - self._gridCC = ccV.reshape((-1, self.dim), order='F') + self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F') return self._gridCC return locals() _gridCC = None # Store grid by default diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 6cc5d8a6..4c0841b8 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -26,17 +26,26 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): .. plot:: examples/mesh/plot_TensorMesh.py + For a quick tensor mesh on a (10x12x15) unit cube:: + + mesh = TensorMesh([10, 12, 15]) + """ _meshType = 'TENSOR' - def __init__(self, h, x0=None): - super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0) - - assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" - - 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) 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)" # Ensure h contains 1D vectors self._h = [mkvc(x.astype(float)) for x in h] @@ -396,7 +405,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: - nF_nE = self.nF if 'F' in locType else self.nE + nF_nE = self.nFv if 'F' in locType else self.nEv components = [spzeros(loc.shape[0], n) for n in nF_nE] components[ind] = interpmat(loc, *self.getTensor(locType)) Q = sp.hstack(components) diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 0b9ff7b3..4b1c4fc3 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) @@ -89,18 +89,18 @@ class TensorView(object): # Determine the subplot number: 131, 121 numPlots = 130 if plotAll else len(imageType)/2*10+100 pltNum = 1 - fx, fy, fz = self.r(I,'F','F','M') + fxyz = self.r(I,'F','F','M') if plotAll or 'Fx' in imageType: ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(fx, imageType='Fx', ax=ax_x, **options) + self.plotImage(fxyz[0], imageType='Fx', ax=ax_x, **options) pltNum +=1 if plotAll or 'Fy' in imageType: ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(fy, imageType='Fy', ax=ax_y, **options) + self.plotImage(fxyz[1], imageType='Fy', ax=ax_y, **options) pltNum +=1 if plotAll or 'Fz' in imageType: ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(fz, imageType='Fz', ax=ax_z, **options) + self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options) pltNum +=1 return else: @@ -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,60 @@ 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)): + """ + 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/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 6f5970e6..3d22a297 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -1,4 +1,4 @@ -from SimPEG.utils import sdiag +from SimPEG.utils import sdiag, count, timeIt, setKwargs import numpy as np class Regularization(object): @@ -7,7 +7,7 @@ class Regularization(object): @property def mref(self): if getattr(self, '_mref', None) is None: - self._mref = np.zeros(self.mesh.nC); + return np.zeros(self.mesh.nC); return self._mref @mref.setter def mref(self, value): @@ -40,21 +40,22 @@ class Regularization(object): self._Wz = sdiag(a)*self.mesh.cellGradz return self._Wz + alpha_s = 1e-6 + 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 - self._Wx = None - self._Wy = None - self._Wz = None - self.alpha_s = 1e-6 - self.alpha_x = 1 - self.alpha_y = 1 - self.alpha_z = 1 + def pnorm(self, r): return 0.5*r.dot(r) + @timeIt def modelObj(self, m): mresid = m - self.mref @@ -69,6 +70,7 @@ class Regularization(object): return mobj + @timeIt def modelObjDeriv(self, m): """ @@ -104,8 +106,8 @@ class Regularization(object): return mobjDeriv - def modelObj2Deriv(self, m): - mresid = m - self.mref + @timeIt + def modelObj2Deriv(self): mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py index af384971..05ae09df 100644 --- a/SimPEG/tests/HTMLTestRunner.py +++ b/SimPEG/tests/HTMLTestRunner.py @@ -241,7 +241,8 @@ function showClassDetail(cid, count) { for (var i = 0; i < count; i++) { tid = id_list[i]; if (toHide) { - document.getElementById('div_'+tid).style.display = 'none' + var divTid = document.getElementById('div_'+tid); + if(divTid !== null){divTid.style.display = 'none';} document.getElementById(tid).className = 'hiddenRow'; } else { diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 1cc2bbae..2723b0b7 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -189,7 +189,7 @@ def Rosenbrock(x, return_g=True, return_H=True): out += (H,) return out if len(out) > 1 else out[0] -def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10): """ Basic derivative check @@ -201,6 +201,9 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): :param int num: number of times to reduce step length, h :param bool plotIt: if you would like to plot :param numpy.array dx: step direction + :param int expectedOrder: The order that you expect the derivative to yield. + :param float tolerance: The tolerance on the expected order. + :param float eps: What is zero? :rtype: bool :return: did you pass the test?! @@ -243,9 +246,6 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): order1 = np.log10(E1[:-1]/E1[1:]) print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) - tolerance = 0.9 - expectedOrder = 2 - eps = 1e-10 order0 = order0[E0[1:] > eps] order1 = order1[E1[1:] > eps] belowTol = order1.size == 0 and order0.size > 0 @@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): -def getQuadratic(A, b): +def getQuadratic(A, b, c=0): """ - Given A and b, this returns a quadratic, Q + Given A, b and c, this returns a quadratic, Q .. math:: - \mathbf{Q( x ) = 0.5 x A x + b x} + \mathbf{Q( x ) = 0.5 x A x + b x} + c """ def Quadratic(x, return_g=True, return_H=True): - f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + c out = (f,) if return_g: g = A.dot(x) + b diff --git a/SimPEG/tests/runTests.py b/SimPEG/tests/runTests.py index c44f9a1e..58d94303 100644 --- a/SimPEG/tests/runTests.py +++ b/SimPEG/tests/runTests.py @@ -4,47 +4,54 @@ import unittest import HTMLTestRunner # This code will run all tests in directory named test_*.py +def main(html=False): + TITLE = 'Test Results' + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) -TITLE = 'Test Results' -test_file_strings = glob.glob('test_*.py') -module_strings = [str[0:len(str)-3] for str in test_file_strings] -suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str - in module_strings] -testSuite = unittest.TestSuite(suites) -unittest.TextTestRunner(verbosity=2).run(testSuite) + if not html: + unittest.TextTestRunner(verbosity=2).run(testSuite) + return -outfile = open("report.html", "w") -runner = HTMLTestRunner.HTMLTestRunner( - stream=outfile, - title=TITLE, - description='SimPEG Test Report was automatically generated.' - ) + outfile = open("report.html", "w") + runner = HTMLTestRunner.HTMLTestRunner( + stream=outfile, + title=TITLE, + description='SimPEG Test Report was automatically generated.', + verbosity=2 + ) -runner.run(testSuite) -outfile.close() + runner.run(testSuite) + outfile.close() -reader = open("report.html", "r") -writer = open("../../docs/api_TestResults.rst", "w") + reader = open("report.html", "r") + writer = open("../../docs/api_TestResults.rst", "w") -writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n') + writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n') -go = False -for line in reader: - skip = False - if line == '