From df51919d8135ab0d8a94da0ecca3e615246b2f62 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 11:49:55 -0800 Subject: [PATCH] Documentation Updates. Issue #11 --- SimPEG/inverse/Optimize.py | 166 +++++++++++++++++++++++++++++++++-- SimPEG/mesh/InnerProducts.py | 8 +- 2 files changed, 163 insertions(+), 11 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 3e50e3df..d4829b7b 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -3,6 +3,7 @@ import matplotlib.pyplot as plt from SimPEG.utils import mkvc, sdiag norm = np.linalg.norm import scipy.sparse as sp +from SimPEG import Solver try: from pubsub import pub @@ -16,7 +17,7 @@ except Exception, e: class Minimize(object): """ - Minimize is a general class for derivative based optimization. + Minimize is a general class for derivative based optimization. """ @@ -38,6 +39,7 @@ class Minimize(object): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): + """Sets key word arguments (kwargs) that are present in the object.""" # Set the variables, throw an error if they don't exist. for attr in kwargs: if hasattr(self, attr): @@ -47,11 +49,58 @@ class Minimize(object): def minimize(self, evalFunction, x0): """ + Minimizes the function (evalFunction) starting at the location x0. + + :param def evalFunction: function handle that evaluates: f, g, H = F(x) + :param numpy.ndarray x0: starting location + :rtype: numpy.ndarray + :return: x, the last iterate of the optimization algorithm evalFunction is a function handle:: - evalFunction(x, return_g=True, return_H=True ) + (f[, g][, H]) = evalFunction(x, return_g=False, return_H=False ) + + 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: + f, g, H = evalFunction(xc) + printIter() + if stoppingCriteria(): break + p = findSearchDirection() + p = scaleSearchDirection(p) + xt, passLS = modifySearchDirection(p) + if not passLS: + xt, caught = modifySearchDirectionBreak(p) + if not caught: return xc + doEndIteration(xt) + + printDone() + return xc """ self.evalFunction = evalFunction self.startup(x0) @@ -67,7 +116,7 @@ class Minimize(object): 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) + 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 @@ -89,6 +138,19 @@ class Minimize(object): self._parent = value def startup(self, x0): + """ + **startup** is called at the start of any new minimize call. + + This will set:: + + x0 = x0 + xc = x0 + _iter = _iterLS = 0 + + :param numpy.ndarray x0: initial x + :rtype: None + :return: None + """ self._iter = 0 self._iterLS = 0 self._STOP = np.zeros((5,1),dtype=bool) @@ -99,7 +161,10 @@ class Minimize(object): def printInit(self): """ - printIter is called at the beginning of the optimization routine. + printInit is called at the beginning of the optimization routine. + + If there is a parent object, printInit will check for a + parent.printInit function and call that. """ if doPub: pub.sendMessage('Minimize.printInit', minimize=self) @@ -114,6 +179,9 @@ class Minimize(object): """ printIter is called directly after function evaluations. + If there is a parent object, printIter will check for a + parent.printIter function and call that. + """ if doPub: pub.sendMessage('Minimize.printIter', minimize=self) if self.parent is not None and hasattr(self.parent, 'printIter'): @@ -122,6 +190,13 @@ class Minimize(object): print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) def printDone(self): + """ + printDone is called at the end of the optimization routine. + + If there is a parent object, printDone will check for a + parent.printDone function and call that. + + """ if doPub: pub.sendMessage('Minimize.printDone', minimize=self) if self.parent is not None and hasattr(self.parent, 'printDone'): self.parent.printDone() @@ -149,17 +224,79 @@ class Minimize(object): return all(self._STOP[0:3]) | any(self._STOP[3:]) def projection(self, p): + """ + projects the search direction. + + by default, no projection is applied. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray + :return: p, projected search direction + """ return p def findSearchDirection(self): + """ + **findSearchDirection** should return an approximation of: + + .. math:: + + H p = - g + + Where you are solving for the search direction, p + + The default is: + + .. math:: + + H = I + + p = - g + + And corresponds to SteepestDescent. + + The latest function evaluations are present in:: + + self.f, self.g, self.H + + :rtype: numpy.ndarray + :return: p, Search Direction + """ return -self.g def scaleSearchDirection(self, 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. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray + :return: p, Scaled Search Direction + """ + if self.maxStep < np.abs(p.max()): p = self.maxStep*p/np.abs(p.max()) return p def modifySearchDirection(self, 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: + + * maxIterLS, the maximum number of linesearch iterations + * LSreduction, the expected reduction expected, default: 1e-4 + * LSshorten, how much the step is reduced, default: 0.5 + + If the linesearch is completed, and a descent direction is found, passLS is returned as True. + + Else, a modifySearchDirectionBreak call is preformed. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray,bool + :return: (xt, passLS) + """ # Armijo linesearch descent = np.inner(self.g, p) t = 1 @@ -185,7 +322,7 @@ class Minimize(object): and if the searchDirection break has been caught. By default, no additional work is done, and the - function returns a False indicating the break was not caught. + evalFunction returns a False indicating the break was not caught. :param numpy.ndarray p: searchDirection :rtype: numpy.ndarray,bool @@ -195,6 +332,17 @@ class Minimize(object): return p, False def doEndIteration(self, 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. + + self.xc must be updated in this code. + + :param numpy.ndarray xt: tested new iterate that ensures a descent direction. + :rtype: None + :return: None + """ # store old values self.fOld = self.f self.xOld, self.xc = self.xc, xt @@ -204,14 +352,18 @@ class Minimize(object): class GaussNewton(Minimize): name = 'GaussNewton' def findSearchDirection(self): - return np.linalg.solve(self.H,-self.g) + return Solver(self.H).solve(-self.g) class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' + + maxIterCG = 10 + tolCG = 1e-5 + 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=1e-05, maxiter=10) + p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG) return p diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index 9fe84ac3..f0ac4ab0 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -81,9 +81,9 @@ class InnerProducts(object): def getFaceInnerProduct(self, mu=None, returnP=False): """Wrapper function, - :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct` - :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct2D` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D` """ if self.dim == 2: return getFaceInnerProduct2D(self, mu, returnP) @@ -93,9 +93,9 @@ class InnerProducts(object): def getEdgeInnerProduct(self, sigma=None, returnP=False): """Wrapper function, - :py:func:`SimPEG.InnerProducts.getFaceInnerProduct` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct` - :py:func:`SimPEG.InnerProducts.getFaceInnerProduct2D` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D` """ if self.dim == 2: return getEdgeInnerProduct2D(self, sigma, returnP)