mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 11:43:07 +08:00
Documentation Updates. Issue #11
This commit is contained in:
Rowan Cockett
+159
-7
@@ -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
|
||||
|
||||
|
||||
|
||||
@@ -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)
|
||||
|
||||
Reference in New Issue
Block a user