mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-02 09:38:24 +08:00
Merged in BoundConstraint (pull request #21)
Optimization Framework and Bug Fixes.
This commit is contained in:
rowanc1
+3
-2
@@ -1,5 +1,6 @@
|
||||
*.pyc
|
||||
SimPEG.sublime-project
|
||||
SimPEG.sublime-workspace
|
||||
*.so
|
||||
*.sublime-project
|
||||
*.sublime-workspace
|
||||
docs/_build/
|
||||
myNotebooks/*
|
||||
|
||||
@@ -16,7 +16,7 @@ class DCProblem(ModelTransforms.LogModel, Problem):
|
||||
|
||||
"""
|
||||
def __init__(self, mesh):
|
||||
super(DCProblem, self).__init__(mesh)
|
||||
Problem.__init__(self, mesh)
|
||||
self.mesh.setCellGradBC('neumann')
|
||||
|
||||
def reshapeFields(self, u):
|
||||
|
||||
@@ -13,13 +13,13 @@ class LinearProblem(Problem):
|
||||
return self.G.dot(m)
|
||||
|
||||
def J(self, m, v, u=None):
|
||||
return G.dot(v)
|
||||
return self.G.dot(v)
|
||||
|
||||
def Jt(self, m, v, u=None):
|
||||
return G.T.dot(v)
|
||||
return self.G.T.dot(v)
|
||||
|
||||
if __name__ == '__main__':
|
||||
N = 100
|
||||
|
||||
def example(N):
|
||||
h = np.ones(N)/N
|
||||
M = TensorMesh([h])
|
||||
|
||||
@@ -28,8 +28,6 @@ if __name__ == '__main__':
|
||||
p = -0.25
|
||||
q = 0.25
|
||||
|
||||
|
||||
|
||||
g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx)
|
||||
|
||||
G = np.empty((nk, M.nC))
|
||||
@@ -38,12 +36,6 @@ if __name__ == '__main__':
|
||||
G[i,:] = g(i)
|
||||
|
||||
|
||||
|
||||
plt.figure(1)
|
||||
for i in range(nk):
|
||||
plt.plot(G[i,:])
|
||||
|
||||
|
||||
m_true = np.zeros(M.nC)
|
||||
m_true[M.vectorCCx > 0.3] = 1.
|
||||
m_true[M.vectorCCx > 0.45] = -0.5
|
||||
@@ -55,29 +47,29 @@ if __name__ == '__main__':
|
||||
|
||||
d_obs = d_true + noise
|
||||
|
||||
# plt.figure(3)
|
||||
# plt.plot(d_true,'-o')
|
||||
# plt.plot(d_obs,'r-o')
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
prob = LinearProblem(M)
|
||||
prob.G = G
|
||||
prob.dobs = d_obs
|
||||
prob.std = np.ones_like(d_obs)*0.1
|
||||
|
||||
return prob, m_true
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
prob, m_true = example(100)
|
||||
M = prob.mesh
|
||||
|
||||
reg = Regularization(M)
|
||||
|
||||
opt = InexactGaussNewton(maxIter=20)
|
||||
|
||||
inv = Inversion(prob,reg,opt,beta0=1e-4)
|
||||
|
||||
m0 = np.zeros_like(m_true)
|
||||
|
||||
mrec = inv.run(m0)
|
||||
|
||||
plt.figure(1)
|
||||
for i in range(prob.G.shape[0]):
|
||||
plt.plot(prob.G[i,:])
|
||||
|
||||
plt.figure(2)
|
||||
|
||||
|
||||
@@ -140,7 +140,7 @@ class Problem(object):
|
||||
This can often be computed given a vector (i.e. J(v)) rather than stored, as J is a large dense matrix.
|
||||
|
||||
"""
|
||||
pass
|
||||
raise NotImplementedError('J is not yet implemented.')
|
||||
|
||||
def Jt(self, m, v, u=None):
|
||||
"""
|
||||
@@ -152,7 +152,7 @@ class Problem(object):
|
||||
|
||||
Effect of transpose of J on a vector v.
|
||||
"""
|
||||
pass
|
||||
raise NotImplementedError('Jt is not yet implemented.')
|
||||
|
||||
|
||||
def J_approx(self, m, v, u=None):
|
||||
|
||||
@@ -8,5 +8,5 @@ class Cooling(object):
|
||||
|
||||
def getBeta(self):
|
||||
if self._beta is None:
|
||||
return beta0
|
||||
return self._beta / beta_coolingFactor
|
||||
return self.beta0
|
||||
return self._beta / self.beta_coolingFactor
|
||||
|
||||
+94
-38
@@ -1,35 +1,33 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from SimPEG.utils import sdiag, mkvc
|
||||
import SimPEG
|
||||
from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers
|
||||
from Optimize import Remember
|
||||
from BetaSchedule import Cooling
|
||||
|
||||
class Inversion(object):
|
||||
"""docstring for Inversion"""
|
||||
class BaseInversion(object):
|
||||
"""docstring for BaseInversion"""
|
||||
|
||||
maxIter = 10
|
||||
name = 'SimPEG Inversion'
|
||||
maxIter = 1
|
||||
name = 'BaseInversion'
|
||||
debug = False
|
||||
beta0 = 1e4
|
||||
|
||||
def __init__(self, prob, reg, opt, **kwargs):
|
||||
setKwargs(self, **kwargs)
|
||||
self.prob = prob
|
||||
self.reg = reg
|
||||
self.opt = opt
|
||||
self.opt.parent = self
|
||||
self.setKwargs(**kwargs)
|
||||
|
||||
def setKwargs(self, **kwargs):
|
||||
"""Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
|
||||
for attr in kwargs:
|
||||
if hasattr(self, attr):
|
||||
setattr(self, attr, kwargs[attr])
|
||||
else:
|
||||
raise Exception('%s attr is not recognized' % attr)
|
||||
self.stoppers = [SimPEG.inverse.StoppingCriteria.iteration, SimPEG.inverse.StoppingCriteria.phi_d_target_Inversion]
|
||||
|
||||
def printInit(self):
|
||||
print "%s %s %s" % ('='*22, self.name, '='*22)
|
||||
print " # beta phi_d phi_m f norm(dJ) #LS"
|
||||
print "%s" % '-'*62
|
||||
|
||||
def printIter(self):
|
||||
print "%3d %1.2e %1.2e %1.2e %1.2e %1.2e %3d" % (self.opt._iter, self._beta, self._phi_d_last, self._phi_m_last, self.opt.f, np.linalg.norm(self.opt.g), self.opt._iterLS)
|
||||
# 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)
|
||||
|
||||
@property
|
||||
def Wd(self):
|
||||
@@ -53,34 +51,85 @@ class Inversion(object):
|
||||
if getattr(self, '_phi_d_target', None) is None:
|
||||
return self.prob.dobs.size #
|
||||
return self._phi_d_target
|
||||
|
||||
@phi_d_target.setter
|
||||
def phi_d_target(self, value):
|
||||
self._phi_d_target = value
|
||||
|
||||
def run(self, m0):
|
||||
m = m0
|
||||
self._iter = 0
|
||||
self._beta = None
|
||||
self.startup(m0)
|
||||
while True:
|
||||
self._beta = self.getBeta()
|
||||
m = self.opt.minimize(self.evalFunction,m)
|
||||
self.m = self.opt.minimize(self.evalFunction, self.m)
|
||||
self.doEndIteration()
|
||||
if self.stoppingCriteria(): break
|
||||
self._iter += 1
|
||||
return m
|
||||
|
||||
beta0 = 1.e2
|
||||
beta_coolingFactor = 5.
|
||||
self.printDone()
|
||||
return self.m
|
||||
|
||||
def startup(self, m0):
|
||||
"""
|
||||
**startup** is called at the start of any new run call.
|
||||
|
||||
If you have things that also need to run on startup, you can create a method::
|
||||
|
||||
def _startup*(self, x0):
|
||||
pass
|
||||
|
||||
Where the * can be any string. If present, _startup* will be called at the start of the default startup call.
|
||||
You may also completely overwrite this function.
|
||||
|
||||
:param numpy.ndarray x0: initial x
|
||||
: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)
|
||||
|
||||
self.m = m0
|
||||
self._iter = 0
|
||||
self._beta = None
|
||||
|
||||
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):
|
||||
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)()
|
||||
|
||||
# store old values
|
||||
self.phi_d_last = self.phi_d
|
||||
self.phi_m_last = self.phi_m
|
||||
self._iter += 1
|
||||
|
||||
def getBeta(self):
|
||||
if self._beta is None:
|
||||
return self.beta0
|
||||
return self._beta / self.beta_coolingFactor
|
||||
return self.beta0
|
||||
|
||||
def stoppingCriteria(self):
|
||||
self._STOP = np.zeros(2,dtype=bool)
|
||||
self._STOP[0] = self._iter >= self.maxIter
|
||||
self._STOP[1] = self._phi_d_last <= self.phi_d_target
|
||||
return np.any(self._STOP)
|
||||
if self.debug: print 'checking stoppingCriteria'
|
||||
return checkStoppers(self, self.stoppers)
|
||||
|
||||
|
||||
def printDone(self):
|
||||
"""
|
||||
**printDone** is called at the end of the inversion routine.
|
||||
|
||||
"""
|
||||
printStoppers(self, self.stoppers)
|
||||
|
||||
|
||||
def evalFunction(self, m, return_g=True, return_H=True):
|
||||
@@ -89,8 +138,8 @@ class Inversion(object):
|
||||
phi_d = self.dataObj(m, u)
|
||||
phi_m = self.reg.modelObj(m)
|
||||
|
||||
self._phi_d_last = phi_d
|
||||
self._phi_m_last = phi_m
|
||||
self.phi_d = phi_d
|
||||
self.phi_m = phi_m
|
||||
|
||||
f = phi_d + self._beta * phi_m
|
||||
|
||||
@@ -111,7 +160,7 @@ class Inversion(object):
|
||||
|
||||
operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float )
|
||||
out += (operator,)
|
||||
return out
|
||||
return out if len(out) > 1 else out[0]
|
||||
|
||||
|
||||
def dataObj(self, m, u=None):
|
||||
@@ -219,3 +268,10 @@ class Inversion(object):
|
||||
|
||||
return dmisfit
|
||||
|
||||
class Inversion(Cooling, Remember, BaseInversion):
|
||||
|
||||
maxIter = 10
|
||||
name = "SimPEG Inversion"
|
||||
|
||||
def __init__(self, prob, reg, opt, **kwargs):
|
||||
BaseInversion.__init__(self, prob, reg, opt, **kwargs)
|
||||
|
||||
+351
-79
@@ -1,6 +1,6 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG.utils import mkvc, sdiag
|
||||
from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers
|
||||
norm = np.linalg.norm
|
||||
import scipy.sparse as sp
|
||||
from SimPEG import Solver
|
||||
@@ -12,6 +12,75 @@ except Exception, e:
|
||||
print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.'
|
||||
doPub = False
|
||||
|
||||
class StoppingCriteria(object):
|
||||
"""docstring for StoppingCriteria"""
|
||||
|
||||
iteration = { "str": "%d : maxIter = %3d <= iter = %3d",
|
||||
"left": lambda M: M.maxIter, "right": lambda M: M._iter,
|
||||
"stopType": "critical"}
|
||||
|
||||
iterationLS = { "str": "%d : maxIterLS = %3d <= iterLS = %3d",
|
||||
"left": lambda M: M.maxIterLS, "right": lambda M: M._iterLS,
|
||||
"stopType": "critical"}
|
||||
|
||||
armijoGoldstein = { "str": "%d : ft = %1.4e <= alp*descent = %1.4e",
|
||||
"left": lambda M: M._LS_ft, "right": lambda M: M.f + M.LSreduction * M._LS_descent,
|
||||
"stopType": "optimal"}
|
||||
|
||||
tolerance_f = { "str": "%d : |fc-fOld| = %1.4e <= tolF*(1+|f0|) = %1.4e",
|
||||
"left": lambda M: 1 if M._iter==0 else abs(M.f-M.f_last), "right": lambda M: 0 if M._iter==0 else M.tolF*(1+abs(M.f0)),
|
||||
"stopType": "optimal"}
|
||||
|
||||
moving_x = { "str": "%d : |xc-x_last| = %1.4e <= tolX*(1+|x0|) = %1.4e",
|
||||
"left": lambda M: 1 if M._iter==0 else norm(M.xc-M.x_last), "right": lambda M: 0 if M._iter==0 else M.tolX*(1+norm(M.x0)),
|
||||
"stopType": "optimal"}
|
||||
|
||||
tolerance_g = { "str": "%d : |proj(x-g)-x| = %1.4e <= tolG = %1.4e",
|
||||
"left": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "right": lambda M: M.tolG,
|
||||
"stopType": "optimal"}
|
||||
|
||||
norm_g = { "str": "%d : |proj(x-g)-x| = %1.4e <= 1e3*eps = %1.4e",
|
||||
"left": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "right": lambda M: 1e3*M.eps,
|
||||
"stopType": "critical"}
|
||||
|
||||
bindingSet = { "str": "%d : probSize = %3d <= bindingSet = %3d",
|
||||
"left": lambda M: M.xc.size, "right": lambda M: np.sum(M.bindingSet(M.xc)),
|
||||
"stopType": "critical"}
|
||||
|
||||
bindingSet_LS = { "str": "%d : probSize = %3d <= bindingSet = %3d",
|
||||
"left": lambda M: M._LS_xt.size, "right": lambda M: np.sum(M.bindingSet(M._LS_xt)),
|
||||
"stopType": "critical"}
|
||||
|
||||
phi_d_target_Minimize = { "str": "%d : phi_d = %1.4e <= phi_d_target = %1.4e ",
|
||||
"left": lambda M: M.parent.phi_d, "right": lambda M: M.parent.phi_d_target,
|
||||
"stopType": "critical"}
|
||||
|
||||
phi_d_target_Inversion = { "str": "%d : phi_d = %1.4e <= phi_d_target = %1.4e ",
|
||||
"left": lambda I: I.phi_d, "right": lambda I: I.phi_d_target,
|
||||
"stopType": "critical"}
|
||||
|
||||
|
||||
class IterationPrinters(object):
|
||||
"""docstring for IterationPrinters"""
|
||||
|
||||
iteration = {"title": "#", "value": lambda M: M._iter, "width": 5, "format": "%3d"}
|
||||
f = {"title": "f", "value": lambda M: M.f, "width": 10, "format": "%1.2e"}
|
||||
norm_g = {"title": "|proj(x-g)-x|", "value": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "width": 15, "format": "%1.2e"}
|
||||
totalLS = {"title": "LS", "value": lambda M: M._iterLS, "width": 5, "format": "%d"}
|
||||
|
||||
iterationLS = {"title": "#", "value": lambda M: (M._iter, M._iterLS), "width": 5, "format": "%3d.%d"}
|
||||
LS_ft = {"title": "ft", "value": lambda M: M._LS_ft, "width": 10, "format": "%1.2e"}
|
||||
LS_t = {"title": "t", "value": lambda M: M._LS_t, "width": 10, "format": "%0.5f"}
|
||||
LS_armijoGoldstein = {"title": "f + alp*g.T*p", "value": lambda M: M.f + M.LSreduction*M._LS_descent, "width": 16, "format": "%1.2e"}
|
||||
|
||||
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"}
|
||||
|
||||
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"}
|
||||
phi_m = {"title": "phi_m", "value": lambda M: M.parent.phi_m, "width": 10, "format": "%1.2e"}
|
||||
|
||||
|
||||
class Minimize(object):
|
||||
@@ -22,7 +91,7 @@ class Minimize(object):
|
||||
|
||||
"""
|
||||
|
||||
name = "GeneralOptimizationAlgorithm"
|
||||
name = "General Optimization Algorithm"
|
||||
|
||||
maxIter = 20
|
||||
maxIterLS = 10
|
||||
@@ -34,17 +103,18 @@ class Minimize(object):
|
||||
tolG = 1e-1
|
||||
eps = 1e-5
|
||||
|
||||
debug = False
|
||||
debugLS = False
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub
|
||||
self.setKwargs(**kwargs)
|
||||
self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration]
|
||||
self.stoppersLS = [StoppingCriteria.armijoGoldstein, StoppingCriteria.iterationLS]
|
||||
|
||||
def setKwargs(self, **kwargs):
|
||||
"""Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
|
||||
for attr in kwargs:
|
||||
if hasattr(self, attr):
|
||||
setattr(self, attr, kwargs[attr])
|
||||
else:
|
||||
raise Exception('%s attr is not recognized' % attr)
|
||||
self.printers = [IterationPrinters.iteration, IterationPrinters.f, IterationPrinters.norm_g, IterationPrinters.totalLS]
|
||||
self.printersLS = [IterationPrinters.iterationLS, IterationPrinters.LS_ft, IterationPrinters.LS_t, IterationPrinters.LS_armijoGoldstein]
|
||||
|
||||
setKwargs(self, **kwargs)
|
||||
|
||||
def minimize(self, evalFunction, x0):
|
||||
"""
|
||||
@@ -59,6 +129,14 @@ class Minimize(object):
|
||||
|
||||
(f[, g][, H]) = evalFunction(x, return_g=False, return_H=False )
|
||||
|
||||
def evalFunction(x, return_g=False, return_H=False):
|
||||
out = (f,)
|
||||
if return_g:
|
||||
out += (g,)
|
||||
if return_H:
|
||||
out += (H,)
|
||||
return out if len(out) > 1 else out[0]
|
||||
|
||||
|
||||
Events are fired with the following inputs via pypubsub::
|
||||
|
||||
@@ -146,19 +224,33 @@ class Minimize(object):
|
||||
xc = x0
|
||||
_iter = _iterLS = 0
|
||||
|
||||
If you have things that also need to run on startup, you can create a method::
|
||||
|
||||
def _startup*(self, x0):
|
||||
pass
|
||||
|
||||
Where the * can be any string. If present, _startup* will be called at the start of the default startup call.
|
||||
You may also completely overwrite this function.
|
||||
|
||||
:param numpy.ndarray x0: initial x
|
||||
: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)
|
||||
|
||||
self._iter = 0
|
||||
self._iterLS = 0
|
||||
self._STOP = np.zeros((5,1),dtype=bool)
|
||||
|
||||
x0 = self.projection(x0) # ensure that we start of feasible.
|
||||
self.x0 = x0
|
||||
self.xc = x0
|
||||
self.xOld = x0
|
||||
self.f_last = np.nan
|
||||
self.x_last = x0
|
||||
|
||||
def printInit(self):
|
||||
|
||||
def printInit(self, inLS=False):
|
||||
"""
|
||||
**printInit** is called at the beginning of the optimization routine.
|
||||
|
||||
@@ -166,15 +258,12 @@ class Minimize(object):
|
||||
parent.printInit function and call that.
|
||||
|
||||
"""
|
||||
if doPub: pub.sendMessage('Minimize.printInit', minimize=self)
|
||||
if self.parent is not None and hasattr(self.parent, 'printInit'):
|
||||
self.parent.printInit()
|
||||
return
|
||||
print "%s %s %s" % ('='*22, self.name, '='*22)
|
||||
print "iter\tJc\t\tnorm(dJ)\tLS"
|
||||
print "%s" % '-'*57
|
||||
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)
|
||||
|
||||
def printIter(self):
|
||||
def printIter(self, inLS=False):
|
||||
"""
|
||||
**printIter** is called directly after function evaluations.
|
||||
|
||||
@@ -182,13 +271,11 @@ class Minimize(object):
|
||||
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'):
|
||||
self.parent.printIter()
|
||||
return
|
||||
print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS)
|
||||
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)
|
||||
|
||||
def printDone(self):
|
||||
def printDone(self, inLS=False):
|
||||
"""
|
||||
**printDone** is called at the end of the optimization routine.
|
||||
|
||||
@@ -196,31 +283,19 @@ class Minimize(object):
|
||||
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()
|
||||
return
|
||||
print "%s STOP! %s" % ('-'*25,'-'*25)
|
||||
# TODO: put controls on gradient value, min model update, and function value
|
||||
if self._iter > 0:
|
||||
print "%d : |fc-fOld| = %1.4e <= tolF*(1+|fStop|) = %1.4e" % (self._STOP[0], abs(self.f-self.fOld), self.tolF*(1+abs(self.fStop)))
|
||||
print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (self._STOP[1], norm(self.xc-self.xOld), self.tolX*(1+norm(self.x0)))
|
||||
print "%d : |g| = %1.4e <= tolG*(1+|fStop|) = %1.4e" % (self._STOP[2], norm(self.g), self.tolG*(1+abs(self.fStop)))
|
||||
print "%d : |g| = %1.4e <= 1e3*eps = %1.4e" % (self._STOP[3], norm(self.g), 1e3*self.eps)
|
||||
print "%d : iter = %3d\t <= maxIter\t = %3d" % (self._STOP[4], self._iter, self.maxIter)
|
||||
print "%s DONE! %s\n" % ('='*25,'='*25)
|
||||
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):
|
||||
|
||||
def stoppingCriteria(self, inLS=False):
|
||||
if self._iter == 0:
|
||||
self.fStop = self.f # Save this for stopping criteria
|
||||
self.f0 = self.f
|
||||
self.g0 = self.g
|
||||
return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
|
||||
|
||||
# check stopping rules
|
||||
self._STOP[0] = self._iter > 0 and (abs(self.f-self.fOld) <= self.tolF*(1+abs(self.fStop)))
|
||||
self._STOP[1] = self._iter > 0 and (norm(self.xc-self.xOld) <= self.tolX*(1+norm(self.x0)))
|
||||
self._STOP[2] = norm(self.g) <= self.tolG*(1+abs(self.fStop))
|
||||
self._STOP[3] = norm(self.g) <= 1e3*self.eps
|
||||
self._STOP[4] = self._iter >= self.maxIter
|
||||
return all(self._STOP[0:3]) | any(self._STOP[3:])
|
||||
|
||||
def projection(self, p):
|
||||
"""
|
||||
@@ -278,6 +353,8 @@ class Minimize(object):
|
||||
p = self.maxStep*p/np.abs(p.max())
|
||||
return p
|
||||
|
||||
nameLS = "Armijo linesearch"
|
||||
|
||||
def modifySearchDirection(self, p):
|
||||
"""
|
||||
**modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria.
|
||||
@@ -296,20 +373,23 @@ class Minimize(object):
|
||||
:rtype: numpy.ndarray,bool
|
||||
:return: (xt, passLS)
|
||||
"""
|
||||
# Armijo linesearch
|
||||
descent = np.inner(self.g, p)
|
||||
t = 1
|
||||
iterLS = 0
|
||||
while iterLS < self.maxIterLS:
|
||||
xt = self.projection(self.xc + t*p)
|
||||
ft = self.evalFunction(xt, return_g=False, return_H=False)
|
||||
if ft < self.f + t*self.LSreduction*descent:
|
||||
break
|
||||
iterLS += 1
|
||||
t = self.LSshorten*t
|
||||
# Projected Armijo linesearch
|
||||
self._LS_t = 1
|
||||
self._iterLS = 0
|
||||
while self._iterLS < self.maxIterLS:
|
||||
self._LS_xt = self.projection(self.xc + self._LS_t*p)
|
||||
self._LS_ft = self.evalFunction(self._LS_xt, return_g=False, return_H=False)
|
||||
self._LS_descent = np.inner(self.g, self._LS_xt - self.xc) # this takes into account multiplying by t, but is important for projection.
|
||||
if self.stoppingCriteria(inLS=True): break
|
||||
self._iterLS += 1
|
||||
self._LS_t = self.LSshorten*self._LS_t
|
||||
if self.debugLS:
|
||||
if self._iterLS == 1: self.printInit(inLS=True)
|
||||
self.printIter(inLS=True)
|
||||
|
||||
self._iterLS = iterLS
|
||||
return xt, iterLS < self.maxIterLS
|
||||
if self.debugLS and self._iterLS > 0: self.printDone(inLS=True)
|
||||
|
||||
return self._LS_xt, self._iterLS < self.maxIterLS
|
||||
|
||||
def modifySearchDirectionBreak(self, p):
|
||||
"""
|
||||
@@ -327,35 +407,227 @@ class Minimize(object):
|
||||
:rtype: numpy.ndarray,bool
|
||||
:return: (xt, breakCaught)
|
||||
"""
|
||||
self.printDone(inLS=True)
|
||||
print 'The linesearch got broken. Boo.'
|
||||
return p, False
|
||||
|
||||
def doEndIteration(self, xt):
|
||||
"""
|
||||
**doEndIteration** is called at the end of each minimize iteration.
|
||||
**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.
|
||||
By default, function values and x locations are shuffled to store 1 past iteration in memory.
|
||||
|
||||
self.xc must be updated in this code.
|
||||
self.xc must be updated in this code.
|
||||
|
||||
:param numpy.ndarray xt: tested new iterate that ensures a descent direction.
|
||||
:rtype: None
|
||||
:return: None
|
||||
|
||||
If you have things that also need to run at the end of every iteration, you can create a method::
|
||||
|
||||
def _doEndIteration*(self, xt):
|
||||
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)(xt)
|
||||
|
||||
# store old values
|
||||
self.fOld = self.f
|
||||
self.xOld, self.xc = self.xc, xt
|
||||
self.f_last = self.f
|
||||
self.x_last, self.xc = self.xc, xt
|
||||
self._iter += 1
|
||||
if self.debug: self.printDone()
|
||||
|
||||
|
||||
class GaussNewton(Minimize):
|
||||
name = 'GaussNewton'
|
||||
|
||||
class Remember(object):
|
||||
"""
|
||||
This mixin remembers all the things you tend to forget.
|
||||
|
||||
You can remember parameters directly, naming the str in Minimize,
|
||||
or pass a tuple with the name and the function that takes Minimize.
|
||||
|
||||
For Example::
|
||||
|
||||
opt.remember('f',('norm_g', lambda M: np.linalg.norm(M.g)))
|
||||
|
||||
opt.minimize(evalFunction, x0)
|
||||
|
||||
opt.recall('f')
|
||||
|
||||
The param name (str) can also be located in the parent (if no conflicts),
|
||||
and it will be looked up by default.
|
||||
"""
|
||||
|
||||
_rememberThese = []
|
||||
|
||||
def remember(self, *args):
|
||||
self._rememberThese = args
|
||||
|
||||
def recall(self, param):
|
||||
assert param in self._rememberList, "You didn't tell me to remember "+param+", you gotta tell me what to remember!"
|
||||
return self._rememberList[param]
|
||||
|
||||
def _startupRemember(self, x0):
|
||||
self._rememberList = {}
|
||||
for param in self._rememberThese:
|
||||
if type(param) is str:
|
||||
self._rememberList[param] = []
|
||||
elif type(param) is tuple:
|
||||
self._rememberList[param[0]] = []
|
||||
|
||||
def _doEndIterationRemember(self, *args):
|
||||
for param in self._rememberThese:
|
||||
if type(param) is str:
|
||||
if self.debug: print 'Remember is remembering: ' + param
|
||||
val = getattr(self, param, None)
|
||||
if val is None and getattr(self, 'parent', None) is not None:
|
||||
# Look to the parent for the param if not found here.
|
||||
val = getattr(self.parent, param, None)
|
||||
self._rememberList[param].append( val )
|
||||
elif type(param) is tuple:
|
||||
if self.debug: print 'Remember is remembering: ' + param[0]
|
||||
self._rememberList[param[0]].append( param[1](self) )
|
||||
|
||||
|
||||
|
||||
class ProjectedGradient(Minimize, Remember):
|
||||
name = 'Projected Gradient'
|
||||
|
||||
maxIterCG = 10
|
||||
tolCG = 1e-3
|
||||
|
||||
lower = -np.inf
|
||||
upper = np.inf
|
||||
|
||||
def __init__(self,**kwargs):
|
||||
super(ProjectedGradient, self).__init__(**kwargs)
|
||||
|
||||
self.stoppers.append(StoppingCriteria.bindingSet)
|
||||
self.stoppersLS.append(StoppingCriteria.bindingSet_LS)
|
||||
|
||||
self.printers.extend([ IterationPrinters.itType, IterationPrinters.aSet, IterationPrinters.bSet, IterationPrinters.comment ])
|
||||
|
||||
|
||||
def _startup(self, x0):
|
||||
# ensure bound vectors are the same size as the model
|
||||
if type(self.lower) is not np.ndarray:
|
||||
self.lower = np.ones_like(x0)*self.lower
|
||||
if type(self.upper) is not np.ndarray:
|
||||
self.upper = np.ones_like(x0)*self.upper
|
||||
|
||||
self.explorePG = True
|
||||
self.exploreCG = False
|
||||
self.stopDoingPG = False
|
||||
|
||||
self._itType = 'SD'
|
||||
self.projComment = ''
|
||||
|
||||
self.aSet_prev = self.activeSet(x0)
|
||||
|
||||
def projection(self, x):
|
||||
"""Make sure we are feasible."""
|
||||
return np.median(np.c_[self.lower,x,self.upper],axis=1)
|
||||
|
||||
def activeSet(self, x):
|
||||
"""If we are on a bound"""
|
||||
return np.logical_or(x == self.lower, x == self.upper)
|
||||
|
||||
def inactiveSet(self, x):
|
||||
"""The free variables."""
|
||||
return np.logical_not(self.activeSet(x))
|
||||
|
||||
def bindingSet(self, x):
|
||||
"""
|
||||
If we are on a bound and the negative gradient points away from the feasible set.
|
||||
|
||||
Optimality condition. (Satisfies Kuhn-Tucker) MoreToraldo91
|
||||
|
||||
"""
|
||||
bind_up = np.logical_and(x == self.lower, self.g >= 0)
|
||||
bind_low = np.logical_and(x == self.upper, self.g <= 0)
|
||||
return np.logical_or(bind_up, bind_low)
|
||||
|
||||
def findSearchDirection(self):
|
||||
self.aSet_prev = self.activeSet(self.xc)
|
||||
allBoundsAreActive = sum(self.aSet_prev) == self.xc.size
|
||||
|
||||
if self.debug: print 'findSearchDirection: stopDoingPG: ', self.stopDoingPG
|
||||
if self.debug: print 'findSearchDirection: explorePG: ', self.explorePG
|
||||
if self.debug: print 'findSearchDirection: exploreCG: ', self.exploreCG
|
||||
if self.debug: print 'findSearchDirection: aSet', np.sum(self.activeSet(self.xc))
|
||||
if self.debug: print 'findSearchDirection: bSet', np.sum(self.bindingSet(self.xc))
|
||||
if self.debug: print 'findSearchDirection: allBoundsAreActive: ', allBoundsAreActive
|
||||
|
||||
if self.explorePG or not self.exploreCG or allBoundsAreActive:
|
||||
if self.debug: print 'findSearchDirection.PG: doingPG'
|
||||
self._itType = 'SD'
|
||||
p = -self.g
|
||||
else:
|
||||
if self.debug: print 'findSearchDirection.CG: doingCG'
|
||||
# Reset the max decrease each time you do a CG iteration
|
||||
self.f_decrease_max = -np.inf
|
||||
|
||||
self._itType = '.CG.'
|
||||
|
||||
iSet = self.inactiveSet(self.xc) # The inactive set (free variables)
|
||||
bSet = self.bindingSet(self.xc)
|
||||
shape = (self.xc.size, np.sum(iSet))
|
||||
v = np.ones(shape[1])
|
||||
i = np.where(iSet)[0]
|
||||
j = np.arange(shape[1])
|
||||
if self.debug: print 'findSearchDirection.CG: Z.shape', shape
|
||||
Z = sp.csr_matrix((v, (i, j)), shape=shape)
|
||||
|
||||
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 )
|
||||
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
|
||||
|
||||
def _doEndIteration_ProjectedGradient(self, xt):
|
||||
aSet = self.activeSet(xt)
|
||||
bSet = self.bindingSet(xt)
|
||||
|
||||
self.explorePG = not np.all(aSet == self.aSet_prev) # explore proximal gradient
|
||||
self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient
|
||||
|
||||
f_current_decrease = self.f_last - self.f
|
||||
self.projComment = ''
|
||||
if self._iter < 1:
|
||||
# Note that this is reset on every CG iteration.
|
||||
self.f_decrease_max = -np.inf
|
||||
else:
|
||||
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.explorePG = False
|
||||
self.exploreCG = True
|
||||
# implement 3.8, MoreToraldo91
|
||||
#self.eta_2 * max_decrease where max decrease
|
||||
# if true go to CG
|
||||
# don't do too many steps of PG in a row.
|
||||
|
||||
if self.debug: print 'doEndIteration.ProjGrad, f_current_decrease: ', f_current_decrease
|
||||
if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max
|
||||
if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD
|
||||
|
||||
class GaussNewton(Minimize, Remember):
|
||||
name = 'Gauss Newton'
|
||||
def findSearchDirection(self):
|
||||
return Solver(self.H).solve(-self.g)
|
||||
|
||||
|
||||
class InexactGaussNewton(Minimize):
|
||||
name = 'InexactGaussNewton'
|
||||
class InexactGaussNewton(Minimize, Remember):
|
||||
name = 'Inexact Gauss Newton'
|
||||
|
||||
maxIterCG = 10
|
||||
tolCG = 1e-5
|
||||
@@ -366,8 +638,8 @@ class InexactGaussNewton(Minimize):
|
||||
return p
|
||||
|
||||
|
||||
class SteepestDescent(Minimize):
|
||||
name = 'SteepestDescent'
|
||||
class SteepestDescent(Minimize, Remember):
|
||||
name = 'Steepest Descent'
|
||||
def findSearchDirection(self):
|
||||
return -self.g
|
||||
|
||||
@@ -377,9 +649,9 @@ 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')
|
||||
# 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])
|
||||
|
||||
@@ -5,11 +5,17 @@ from SimPEG.utils import mkvc, sdiag
|
||||
from SimPEG import utils
|
||||
from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
import inspect
|
||||
|
||||
happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.']
|
||||
sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo :(', 'Testing is important. Do it again.']
|
||||
try:
|
||||
import getpass
|
||||
name = getpass.getuser()[0].upper() + getpass.getuser()[1:]
|
||||
except Exception, e:
|
||||
name = 'You'
|
||||
happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.','Not just a pretty face '+name,'You deserve a pat on the back!','Well done '+name+'!', 'Awesome, '+name+', just awesome.']
|
||||
sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo :(', 'Testing is important. Do it again.',"Did you put your clever trousers on today?",'Just think about a dancing dinosaur and life will get better!','You had so much promise '+name+', oh well...', name.upper()+' ERROR!','Get on it '+name+'!', 'You break it, you fix it.']
|
||||
|
||||
class OrderTest(unittest.TestCase):
|
||||
"""
|
||||
@@ -174,14 +180,14 @@ def Rosenbrock(x, return_g=True, return_H=True):
|
||||
|
||||
f = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
|
||||
g = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1), 200*(x[1]-x[0]**2)])
|
||||
H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]])
|
||||
H = sp.csr_matrix(np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]]))
|
||||
|
||||
out = (f,)
|
||||
if return_g:
|
||||
out += (g,)
|
||||
if return_H:
|
||||
out += (H,)
|
||||
return out
|
||||
return out if len(out) > 1 else out[0]
|
||||
|
||||
def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
|
||||
"""
|
||||
@@ -269,6 +275,28 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
|
||||
return passTest
|
||||
|
||||
|
||||
|
||||
def getQuadratic(A, b):
|
||||
"""
|
||||
Given A and b, this returns a quadratic, Q
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{Q( x ) = 0.5 x A x + b x}
|
||||
"""
|
||||
def Quadratic(x, return_g=True, return_H=True):
|
||||
f = 0.5 * x.dot( A.dot(x)) + b.dot( x )
|
||||
out = (f,)
|
||||
if return_g:
|
||||
g = A.dot(x) + b
|
||||
out += (g,)
|
||||
if return_H:
|
||||
H = A
|
||||
out += (H,)
|
||||
return out if len(out) > 1 else out[0]
|
||||
return Quadratic
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
def simplePass(x):
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
import TestUtils
|
||||
from TestUtils import checkDerivative, Rosenbrock, OrderTest
|
||||
from TestUtils import checkDerivative, Rosenbrock, OrderTest, getQuadratic
|
||||
|
||||
@@ -0,0 +1,54 @@
|
||||
import unittest
|
||||
from SimPEG import Solver
|
||||
from SimPEG.mesh import TensorMesh
|
||||
from SimPEG.utils import sdiag
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from SimPEG import inverse
|
||||
from SimPEG.tests import getQuadratic, Rosenbrock
|
||||
|
||||
TOL = 1e-2
|
||||
|
||||
class TestOptimizers(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.A = sp.identity(2).tocsr()
|
||||
self.b = np.array([-5,-5])
|
||||
|
||||
def test_GN_Rosenbrock(self):
|
||||
GN = inverse.GaussNewton()
|
||||
xopt = GN.minimize(Rosenbrock,np.array([0,0]))
|
||||
x_true = np.array([1.,1.])
|
||||
print 'xopt: ', xopt
|
||||
print 'x_true: ', x_true
|
||||
self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
|
||||
|
||||
def test_GN_quadratic(self):
|
||||
GN = inverse.GaussNewton()
|
||||
xopt = GN.minimize(getQuadratic(self.A,self.b),np.array([0,0]))
|
||||
x_true = np.array([5.,5.])
|
||||
print 'xopt: ', xopt
|
||||
print 'x_true: ', x_true
|
||||
self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
|
||||
|
||||
def test_ProjGradient_quadraticBounded(self):
|
||||
PG = inverse.ProjectedGradient()
|
||||
PG.lower, PG.upper = -2, 2
|
||||
xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0]))
|
||||
x_true = np.array([2.,2.])
|
||||
print 'xopt: ', xopt
|
||||
print 'x_true: ', x_true
|
||||
self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
|
||||
|
||||
def test_ProjGradient_quadratic1Bound(self):
|
||||
myB = np.array([-5,1])
|
||||
PG = inverse.ProjectedGradient()
|
||||
PG.lower, PG.upper = -2, 2
|
||||
xopt = PG.minimize(getQuadratic(self.A,myB),np.array([0,0]))
|
||||
x_true = np.array([2.,-1.])
|
||||
print 'xopt: ', xopt
|
||||
print 'x_true: ', x_true
|
||||
self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -93,7 +93,7 @@ class Solver(object):
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if backend is None: backend = DEFAULTS['scipy']
|
||||
if backend is None: backend = DEFAULTS['direct']
|
||||
|
||||
assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
|
||||
|
||||
|
||||
@@ -7,5 +7,55 @@ from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
|
||||
from sputils import spzeros, kron3, speye, sdiag
|
||||
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
|
||||
from interputils import interpmat
|
||||
from ipythonUtils import easyAnimate as animate
|
||||
import Solver
|
||||
from Solver import Solver
|
||||
|
||||
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:
|
||||
if hasattr(obj, attr):
|
||||
setattr(obj, attr, kwargs[attr])
|
||||
else:
|
||||
raise Exception('%s attr is not recognized' % attr)
|
||||
|
||||
def printTitles(obj, printers, name='Print Titles', pad=''):
|
||||
titles = ''
|
||||
widths = 0
|
||||
for printer in printers:
|
||||
titles += ('{:^%i}'%printer['width']).format(printer['title']) + ''
|
||||
widths += printer['width']
|
||||
print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name)
|
||||
print pad + titles
|
||||
print pad + "%s" % '-'*widths
|
||||
|
||||
def printLine(obj, printers, pad=''):
|
||||
values = ''
|
||||
for printer in printers:
|
||||
values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj))
|
||||
print pad + values
|
||||
|
||||
def checkStoppers(obj, stoppers):
|
||||
# check stopping rules
|
||||
optimal = []
|
||||
critical = []
|
||||
for stopper in stoppers:
|
||||
l = stopper['left'](obj)
|
||||
r = stopper['right'](obj)
|
||||
if stopper['stopType'] == 'optimal':
|
||||
optimal.append(l <= r)
|
||||
if stopper['stopType'] == 'critical':
|
||||
critical.append(l <= r)
|
||||
|
||||
if obj.debug: print 'checkStoppers.optimal: ', optimal
|
||||
if obj.debug: print 'checkStoppers.critical: ', critical
|
||||
|
||||
return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical))
|
||||
|
||||
def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
|
||||
print pad + "%s%s%s" % ('-'*25,stop,'-'*25)
|
||||
for stopper in stoppers:
|
||||
l = stopper['left'](obj)
|
||||
r = stopper['right'](obj)
|
||||
print pad + stopper['str'] % (l<=r,l,r)
|
||||
print pad + "%s%s%s" % ('-'*25,done,'-'*25)
|
||||
|
||||
@@ -0,0 +1,29 @@
|
||||
from tempfile import NamedTemporaryFile
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib import animation
|
||||
from IPython.display import HTML
|
||||
|
||||
# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
|
||||
# http://www.renevolution.com/how-to-install-ffmpeg-on-mac-os-x/
|
||||
|
||||
VIDEO_TAG = """<video controls loop>
|
||||
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
|
||||
Your browser does not support the video tag.
|
||||
</video>"""
|
||||
|
||||
def anim_to_html(anim):
|
||||
if not hasattr(anim, '_encoded_video'):
|
||||
with NamedTemporaryFile(suffix='.mp4') as f:
|
||||
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
|
||||
video = open(f.name, "rb").read()
|
||||
anim._encoded_video = video.encode("base64")
|
||||
|
||||
return VIDEO_TAG.format(anim._encoded_video)
|
||||
|
||||
def display_animation(anim):
|
||||
plt.close(anim._fig)
|
||||
return anim_to_html(anim)
|
||||
|
||||
animation.Animation._repr_html_ = display_animation
|
||||
|
||||
easyAnimate = animation.FuncAnimation
|
||||
+524
-331
File diff suppressed because it is too large
Load Diff
File diff suppressed because one or more lines are too long
Reference in New Issue
Block a user