diff --git a/.gitignore b/.gitignore
index abbec73c..5c39733b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
*.pyc
-SimPEG.sublime-project
-SimPEG.sublime-workspace
+*.so
+*.sublime-project
+*.sublime-workspace
docs/_build/
myNotebooks/*
diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py
index 132966a8..9ddb5332 100644
--- a/SimPEG/forward/DCProblem.py
+++ b/SimPEG/forward/DCProblem.py
@@ -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):
diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py
index d30a5b4d..ef92957b 100644
--- a/SimPEG/forward/LinearProblem.py
+++ b/SimPEG/forward/LinearProblem.py
@@ -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)
diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py
index 2e6831f7..cf22baae 100644
--- a/SimPEG/forward/Problem.py
+++ b/SimPEG/forward/Problem.py
@@ -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):
diff --git a/SimPEG/inverse/BetaSchedule.py b/SimPEG/inverse/BetaSchedule.py
index fe197340..af4d883d 100644
--- a/SimPEG/inverse/BetaSchedule.py
+++ b/SimPEG/inverse/BetaSchedule.py
@@ -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
diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py
index 3aa8ce58..e3d500b1 100644
--- a/SimPEG/inverse/Inversion.py
+++ b/SimPEG/inverse/Inversion.py
@@ -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)
diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py
index 6221246f..37c8b296 100644
--- a/SimPEG/inverse/Optimize.py
+++ b/SimPEG/inverse/Optimize.py
@@ -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])
diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py
index 140b56fa..1cc2bbae 100644
--- a/SimPEG/tests/TestUtils.py
+++ b/SimPEG/tests/TestUtils.py
@@ -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):
diff --git a/SimPEG/tests/__init__.py b/SimPEG/tests/__init__.py
index f896cfd0..d082562d 100644
--- a/SimPEG/tests/__init__.py
+++ b/SimPEG/tests/__init__.py
@@ -1,2 +1,2 @@
import TestUtils
-from TestUtils import checkDerivative, Rosenbrock, OrderTest
+from TestUtils import checkDerivative, Rosenbrock, OrderTest, getQuadratic
diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/tests/test_optimizers.py
new file mode 100644
index 00000000..6fff75d3
--- /dev/null
+++ b/SimPEG/tests/test_optimizers.py
@@ -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()
diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py
index 3efaf10b..a1194487 100644
--- a/SimPEG/utils/Solver.py
+++ b/SimPEG/utils/Solver.py
@@ -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'
diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py
index d3955170..3f82ba6d 100644
--- a/SimPEG/utils/__init__.py
+++ b/SimPEG/utils/__init__.py
@@ -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)
diff --git a/SimPEG/utils/ipythonUtils.py b/SimPEG/utils/ipythonUtils.py
new file mode 100644
index 00000000..aa1eefdd
--- /dev/null
+++ b/SimPEG/utils/ipythonUtils.py
@@ -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 = """"""
+
+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
diff --git a/docs/api_TestResults.rst b/docs/api_TestResults.rst
index 7b276006..8c16f36e 100644
--- a/docs/api_TestResults.rst
+++ b/docs/api_TestResults.rst
@@ -185,9 +185,9 @@ Test Results
-->
-
Start Time: 2013-11-05 16:10:07
-
Duration: 0:00:31.280082
-
Status: Pass 99
+
Start Time: 2013-11-12 11:29:23
+
Duration: 0:00:32.452108
+
Status: Pass 107
SimPEG Test Report was automatically generated.
@@ -377,15 +377,15 @@ Test Results
pt3.2: ==================== checkDerivative ====================
iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order
---------------------------------------------------------
- 0 1.00e-01 3.555e+01 3.750e+01 nan
- 1 1.00e-02 1.711e-01 3.662e-01 2.010
- 2 1.00e-03 1.586e-02 3.653e-03 2.001
- 3 1.00e-04 1.915e-03 3.652e-05 2.000
- 4 1.00e-05 1.948e-04 3.652e-07 2.000
- 5 1.00e-06 1.951e-05 3.662e-09 1.999
- 6 1.00e-07 1.951e-06 3.158e-11 2.064
+ 0 1.00e-01 2.025e+01 2.686e+01 nan
+ 1 1.00e-02 3.953e-01 2.660e-01 2.004
+ 2 1.00e-03 6.347e-02 2.657e-03 2.000
+ 3 1.00e-04 6.587e-03 2.657e-05 2.000
+ 4 1.00e-05 6.610e-04 2.657e-07 2.000
+ 5 1.00e-06 6.613e-05 2.657e-09 2.000
+ 6 1.00e-07 6.613e-06 2.700e-11 1.993
========================= PASS! =========================
- Testing is important.
+ You deserve a pat on the back!
@@ -414,15 +414,15 @@ Test Results
pt3.3: ==================== checkDerivative ====================
iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order
---------------------------------------------------------
- 0 1.00e-01 5.409e-02 4.593e-03 nan
- 1 1.00e-02 5.429e-03 4.574e-05 2.002
- 2 1.00e-03 5.431e-04 4.570e-07 2.000
- 3 1.00e-04 5.431e-05 4.570e-09 2.000
- 4 1.00e-05 5.431e-06 4.570e-11 2.000
- 5 1.00e-06 5.431e-07 4.553e-13 2.002
- 6 1.00e-07 5.431e-08 4.751e-15 1.982
+ 0 1.00e-01 5.087e-02 5.479e-03 nan
+ 1 1.00e-02 5.240e-03 5.471e-05 2.001
+ 2 1.00e-03 5.256e-04 5.468e-07 2.000
+ 3 1.00e-04 5.258e-05 5.467e-09 2.000
+ 4 1.00e-05 5.258e-06 5.467e-11 2.000
+ 5 1.00e-06 5.258e-07 5.467e-13 2.000
+ 6 1.00e-07 5.258e-08 5.672e-15 1.984
========================= PASS! =========================
- That was easy!
+ Testing is important.
@@ -451,15 +451,15 @@ Test Results
pt3.4: ==================== checkDerivative ====================
iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order
---------------------------------------------------------
- 0 1.00e-01 8.756e+00 7.946e+00 nan
- 1 1.00e-02 1.605e-01 7.946e-02 2.000
- 2 1.00e-03 8.894e-03 7.946e-04 2.000
- 3 1.00e-04 8.179e-04 7.946e-06 2.000
- 4 1.00e-05 8.108e-05 7.946e-08 2.000
- 5 1.00e-06 8.100e-06 7.946e-10 2.000
- 6 1.00e-07 8.100e-07 7.947e-12 2.000
+ 0 1.00e-01 6.461e+00 7.059e+00 nan
+ 1 1.00e-02 1.073e-02 7.059e-02 2.000
+ 2 1.00e-03 5.280e-03 7.059e-04 2.000
+ 3 1.00e-04 5.915e-04 7.059e-06 2.000
+ 4 1.00e-05 5.979e-05 7.059e-08 2.000
+ 5 1.00e-06 5.985e-06 7.059e-10 2.000
+ 6 1.00e-07 5.986e-07 7.059e-12 2.000
========================= PASS! =========================
- The test be workin!
+ That was easy!
@@ -498,15 +498,15 @@ Test Results
==================== checkDerivative ====================
iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order
---------------------------------------------------------
- 0 1.00e-01 1.309e-01 4.745e-17 nan
- 1 1.00e-02 1.309e-02 4.825e-17 -0.007
- 2 1.00e-03 1.309e-03 5.993e-17 -0.094
- 3 1.00e-04 1.309e-04 6.648e-17 -0.045
- 4 1.00e-05 1.309e-05 6.931e-17 -0.018
- 5 1.00e-06 1.309e-06 6.247e-17 0.045
- 6 1.00e-07 1.309e-07 5.708e-17 0.039
+ 0 1.00e-01 3.421e-01 3.925e-17 nan
+ 1 1.00e-02 3.421e-02 2.357e-17 0.222
+ 2 1.00e-03 3.421e-03 4.849e-17 -0.313
+ 3 1.00e-04 3.421e-04 5.925e-17 -0.087
+ 4 1.00e-05 3.421e-05 6.316e-17 -0.028
+ 5 1.00e-06 3.421e-06 4.061e-17 0.192
+ 6 1.00e-07 3.421e-07 6.370e-17 -0.196
========================= PASS! =========================
- You are awesome.
+ You get a gold star!
@@ -535,15 +535,15 @@ Test Results
pt4.2: ==================== checkDerivative ====================
iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order
---------------------------------------------------------
- 0 1.00e-01 3.133e-01 3.168e-02 nan
- 1 1.00e-02 3.418e-02 3.168e-04 2.000
- 2 1.00e-03 3.447e-03 3.168e-06 2.000
- 3 1.00e-04 3.449e-04 3.168e-08 2.000
- 4 1.00e-05 3.450e-05 3.168e-10 2.000
- 5 1.00e-06 3.450e-06 3.169e-12 2.000
- 6 1.00e-07 3.450e-07 3.481e-14 1.959
+ 0 1.00e-01 2.565e+00 1.802e-01 nan
+ 1 1.00e-02 2.727e-01 1.802e-03 2.000
+ 2 1.00e-03 2.743e-02 1.802e-05 2.000
+ 3 1.00e-04 2.745e-03 1.802e-07 2.000
+ 4 1.00e-05 2.745e-04 1.802e-09 2.000
+ 5 1.00e-06 2.745e-05 1.802e-11 2.000
+ 6 1.00e-07 2.745e-06 1.800e-13 2.000
========================= PASS! =========================
- That was easy!
+ Yay passed!
@@ -583,22 +583,22 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 2.76e-01 |
- 16 | 6.10e-02 | 4.5221 | 2.1770
- 32 | 1.44e-02 | 4.2290 | 2.0803
+ 8 | 2.58e-01 |
+ 16 | 6.41e-02 | 4.0278 | 2.0100
+ 32 | 1.77e-02 | 3.6279 | 1.8591
---------------------------------------------
- Happy little convergence test!
+ That was easy!
randomTensorMesh: Interpolation 1D: CC
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 2.46e-01 |
- 16 | 1.11e-01 | 2.2066 | 0.9651
- 32 | 4.07e-02 | 2.7360 | 3.1871
+ 8 | 4.97e-01 |
+ 16 | 1.59e-01 | 3.1240 | 1.6281
+ 32 | 4.04e-02 | 3.9430 | 2.3561
---------------------------------------------
- Testing is important.
+ Happy little convergence test!
@@ -629,22 +629,22 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 2.39e-01 |
- 16 | 6.89e-02 | 3.4636 | 1.7923
- 32 | 1.78e-02 | 3.8797 | 1.9559
+ 8 | 3.09e-01 |
+ 16 | 7.89e-02 | 3.9199 | 1.9708
+ 32 | 1.83e-02 | 4.3119 | 2.1083
---------------------------------------------
- You get a gold star!
+ Well done Rowan!
randomTensorMesh: Interpolation 1D: N
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 5.97e-01 |
- 16 | 1.38e-01 | 4.3354 | 1.7681
- 32 | 3.33e-02 | 4.1339 | 2.1579
+ 8 | 4.76e-01 |
+ 16 | 1.11e-01 | 4.2662 | 1.9952
+ 32 | 4.59e-02 | 2.4283 | 1.2771
---------------------------------------------
- Yay passed!
+ Once upon a time, a happy little test passed.
@@ -684,24 +684,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.06e-02 |
- 16 | 1.78e-02 | 3.9781 | 1.9921
- 32 | 4.72e-03 | 3.7595 | 1.9106
- 64 | 1.18e-03 | 3.9981 | 1.9993
+ 8 | 7.60e-02 |
+ 16 | 1.91e-02 | 3.9815 | 1.9933
+ 32 | 4.72e-03 | 4.0446 | 2.0160
+ 64 | 1.18e-03 | 4.0053 | 2.0019
---------------------------------------------
- The test be workin!
+ Happy little convergence test!
randomTensorMesh: Interpolation 2D: CC
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 9.16e-02 |
- 16 | 3.39e-02 | 2.7044 | 1.6128
- 32 | 9.31e-03 | 3.6391 | 1.9612
- 64 | 4.10e-03 | 2.2687 | 1.5130
+ 8 | 1.75e-01 |
+ 16 | 2.69e-02 | 6.5233 | 2.3005
+ 32 | 1.12e-02 | 2.4031 | 1.2481
+ 64 | 3.53e-03 | 3.1726 | 1.8392
---------------------------------------------
- Testing is important.
+ You are awesome.
@@ -732,24 +732,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.04e-02 |
- 16 | 1.88e-02 | 3.7329 | 1.9003
- 32 | 4.77e-03 | 3.9517 | 1.9825
- 64 | 1.19e-03 | 4.0081 | 2.0029
+ 8 | 7.01e-02 |
+ 16 | 1.88e-02 | 3.7205 | 1.8955
+ 32 | 4.59e-03 | 4.1006 | 2.0358
+ 64 | 1.16e-03 | 3.9440 | 1.9797
---------------------------------------------
- Testing is important.
+ Well done Rowan!
randomTensorMesh: Interpolation 2D: Ex
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 1.96e-01 |
- 16 | 4.30e-02 | 4.5625 | 2.7996
- 32 | 2.01e-02 | 2.1406 | 1.0860
- 64 | 3.32e-03 | 6.0522 | 2.4217
+ 8 | 1.12e-01 |
+ 16 | 8.66e-02 | 1.2909 | 0.8416
+ 32 | 2.00e-02 | 4.3264 | 1.8545
+ 64 | 2.47e-03 | 8.1124 | 2.6759
---------------------------------------------
- You get a gold star!
+ You deserve a pat on the back!
@@ -780,10 +780,10 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.01e-02 |
- 16 | 1.88e-02 | 3.7175 | 1.8943
- 32 | 4.59e-03 | 4.1014 | 2.0361
- 64 | 1.19e-03 | 3.8653 | 1.9506
+ 8 | 6.84e-02 |
+ 16 | 1.60e-02 | 4.2800 | 2.0976
+ 32 | 4.79e-03 | 3.3359 | 1.7381
+ 64 | 1.20e-03 | 3.9828 | 1.9938
---------------------------------------------
The test be workin!
@@ -792,12 +792,12 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 1.92e-01 |
- 16 | 5.13e-02 | 3.7457 | 2.6593
- 32 | 1.09e-02 | 4.7164 | 2.0227
- 64 | 4.07e-03 | 2.6707 | 1.4891
+ 8 | 6.00e-02 |
+ 16 | 4.16e-02 | 1.4412 | 0.4902
+ 32 | 7.60e-03 | 5.4787 | 2.8421
+ 64 | 2.60e-03 | 2.9195 | 1.2969
---------------------------------------------
- Go Test Go!
+ Not just a pretty face Rowan
@@ -828,24 +828,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.06e-02 |
- 16 | 1.78e-02 | 3.9781 | 1.9921
- 32 | 4.72e-03 | 3.7595 | 1.9106
- 64 | 1.18e-03 | 3.9981 | 1.9993
+ 8 | 7.60e-02 |
+ 16 | 1.91e-02 | 3.9815 | 1.9933
+ 32 | 4.72e-03 | 4.0446 | 2.0160
+ 64 | 1.18e-03 | 4.0053 | 2.0019
---------------------------------------------
- Yay passed!
+ Happy little convergence test!
randomTensorMesh: Interpolation 2D: Fx
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 9.82e-02 |
- 16 | 5.58e-02 | 1.7592 | 6.6696
- 32 | 6.57e-03 | 8.4887 | 2.2061
- 64 | 4.28e-03 | 1.5342 | 0.6081
+ 8 | 7.56e-02 |
+ 16 | 3.64e-02 | 2.0777 | 0.8083
+ 32 | 1.07e-02 | 3.3892 | 1.8786
+ 64 | 3.58e-03 | 3.0004 | 1.6103
---------------------------------------------
- The test be workin!
+ Go Test Go!
@@ -876,10 +876,10 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.57e-02 |
- 16 | 1.88e-02 | 4.0222 | 2.0080
- 32 | 4.72e-03 | 3.9905 | 1.9966
- 64 | 1.18e-03 | 4.0110 | 2.0040
+ 8 | 7.49e-02 |
+ 16 | 1.80e-02 | 4.1557 | 2.0551
+ 32 | 4.00e-03 | 4.5049 | 2.1715
+ 64 | 1.20e-03 | 3.3419 | 1.7407
---------------------------------------------
That was easy!
@@ -888,12 +888,12 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 6.72e-02 |
- 16 | 3.28e-02 | 2.0482 | 1.6433
- 32 | 1.13e-02 | 2.9135 | 2.1786
- 64 | 3.35e-03 | 3.3628 | 1.6438
+ 8 | 9.96e-02 |
+ 16 | 5.79e-02 | 1.7215 | 0.5101
+ 32 | 1.21e-02 | 4.7877 | 2.1054
+ 64 | 2.60e-03 | 4.6488 | 2.4950
---------------------------------------------
- Yay passed!
+ Happy little convergence test!
@@ -924,10 +924,10 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.04e-02 |
- 16 | 1.88e-02 | 3.7329 | 1.9003
- 32 | 4.77e-03 | 3.9517 | 1.9825
- 64 | 1.19e-03 | 4.0081 | 2.0029
+ 8 | 7.01e-02 |
+ 16 | 1.88e-02 | 3.7205 | 1.8955
+ 32 | 4.59e-03 | 4.1006 | 2.0358
+ 64 | 1.16e-03 | 3.9440 | 1.9797
---------------------------------------------
Once upon a time, a happy little test passed.
@@ -936,12 +936,12 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 1.24e-01 |
- 16 | 6.30e-02 | 1.9636 | 1.2918
- 32 | 2.15e-02 | 2.9310 | 1.6864
- 64 | 4.96e-03 | 4.3359 | 1.9158
+ 8 | 1.13e-01 |
+ 16 | 4.02e-02 | 2.8137 | 1.7733
+ 32 | 1.12e-02 | 3.5726 | 1.4513
+ 64 | 3.25e-03 | 3.4586 | 1.4744
---------------------------------------------
- Happy little convergence test!
+ You get a gold star!
@@ -981,24 +981,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.36e-02 |
- 16 | 1.74e-02 | 4.2291 | 2.0803
- 32 | 4.16e-03 | 4.1900 | 2.0669
- 64 | 9.48e-04 | 4.3858 | 2.1328
+ 8 | 7.59e-02 |
+ 16 | 1.90e-02 | 3.9963 | 1.9987
+ 32 | 4.69e-03 | 4.0474 | 2.0170
+ 64 | 1.18e-03 | 3.9935 | 1.9977
---------------------------------------------
- That was easy!
+ Once upon a time, a happy little test passed.
randomTensorMesh: Interpolation CC
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 2.30e-01 |
- 16 | 1.70e-02 | 13.5276 | 3.5152
- 32 | 6.70e-03 | 2.5405 | 1.2651
- 64 | 1.88e-03 | 3.5689 | 1.8184
+ 8 | 7.88e-02 |
+ 16 | 2.49e-02 | 3.1639 | 2.5664
+ 32 | 1.40e-02 | 1.7739 | 0.8582
+ 64 | 2.03e-03 | 6.9317 | 2.8790
---------------------------------------------
- That was easy!
+ And then everyone was happy.
@@ -1029,24 +1029,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 6.98e-02 |
- 16 | 1.79e-02 | 3.8902 | 1.9598
- 32 | 4.41e-03 | 4.0681 | 2.0244
- 64 | 1.18e-03 | 3.7418 | 1.9037
+ 8 | 7.03e-02 |
+ 16 | 1.88e-02 | 3.7326 | 1.9002
+ 32 | 4.77e-03 | 3.9468 | 1.9807
+ 64 | 1.12e-03 | 4.2670 | 2.0932
---------------------------------------------
- Happy little convergence test!
+ You are awesome.
randomTensorMesh: Interpolation Ex
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 2.92e-01 |
- 16 | 8.20e-02 | 3.5597 | 4.3217
- 32 | 9.49e-03 | 8.6460 | 1.9353
- 64 | 2.89e-03 | 3.2854 | 1.7236
+ 8 | 9.15e-02 |
+ 16 | 7.21e-02 | 1.2700 | 0.5458
+ 32 | 1.22e-02 | 5.9140 | 2.6646
+ 64 | 2.75e-03 | 4.4260 | 1.8486
---------------------------------------------
- Testing is important.
+ That was easy!
@@ -1077,24 +1077,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.04e-02 |
- 16 | 1.87e-02 | 3.7726 | 1.9155
- 32 | 4.61e-03 | 4.0477 | 2.0171
- 64 | 1.09e-03 | 4.2392 | 2.0838
+ 8 | 7.03e-02 |
+ 16 | 1.88e-02 | 3.7324 | 1.9001
+ 32 | 4.32e-03 | 4.3577 | 2.1236
+ 64 | 1.19e-03 | 3.6424 | 1.8649
---------------------------------------------
- You are awesome.
+ You get a gold star!
randomTensorMesh: Interpolation Ey
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 1.30e-01 |
- 16 | 4.52e-02 | 2.8770 | 1.4488
- 32 | 1.36e-02 | 3.3196 | 2.2159
- 64 | 2.81e-03 | 4.8481 | 1.8199
+ 8 | 1.53e-01 |
+ 16 | 6.69e-02 | 2.2859 | 1.4178
+ 32 | 2.05e-02 | 3.2684 | 2.0771
+ 64 | 2.74e-03 | 7.4881 | 2.3296
---------------------------------------------
- You are awesome.
+ Testing is important.
@@ -1125,24 +1125,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.04e-02 |
- 16 | 1.87e-02 | 3.7711 | 1.9150
- 32 | 4.61e-03 | 4.0489 | 2.0175
- 64 | 1.15e-03 | 4.0151 | 2.0054
+ 8 | 7.03e-02 |
+ 16 | 1.88e-02 | 3.7376 | 1.9021
+ 32 | 4.57e-03 | 4.1226 | 2.0435
+ 64 | 1.16e-03 | 3.9523 | 1.9827
---------------------------------------------
- Yay passed!
+ You are awesome.
randomTensorMesh: Interpolation Ez
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 3.81e-01 |
- 16 | 3.89e-02 | 9.7974 | 2.2747
- 32 | 2.09e-02 | 1.8614 | 2.2906
- 64 | 3.97e-03 | 5.2691 | 1.7820
+ 8 | 1.48e-01 |
+ 16 | 9.32e-02 | 1.5867 | 0.4554
+ 32 | 2.27e-02 | 4.1106 | 1.5265
+ 64 | 2.68e-03 | 8.4640 | 3.0885
---------------------------------------------
- Once upon a time, a happy little test passed.
+ Testing is important.
@@ -1173,24 +1173,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.36e-02 |
- 16 | 1.74e-02 | 4.2291 | 2.0803
- 32 | 4.16e-03 | 4.1900 | 2.0669
- 64 | 9.48e-04 | 4.3858 | 2.1328
+ 8 | 7.59e-02 |
+ 16 | 1.90e-02 | 3.9963 | 1.9987
+ 32 | 4.69e-03 | 4.0474 | 2.0170
+ 64 | 1.18e-03 | 3.9935 | 1.9977
---------------------------------------------
- Once upon a time, a happy little test passed.
+ You deserve a pat on the back!
randomTensorMesh: Interpolation Fx
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 1.59e-01 |
- 16 | 7.66e-02 | 2.0691 | 1.2662
- 32 | 1.01e-02 | 7.5592 | 3.0198
- 64 | 1.82e-03 | 5.5842 | 2.5778
+ 8 | 8.49e-02 |
+ 16 | 2.72e-02 | 3.1225 | 1.2154
+ 32 | 1.12e-02 | 2.4354 | 1.9171
+ 64 | 3.67e-03 | 3.0455 | 1.3837
---------------------------------------------
- That was easy!
+ You deserve a pat on the back!
@@ -1221,24 +1221,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.61e-02 |
- 16 | 1.92e-02 | 3.9682 | 1.9885
- 32 | 4.77e-03 | 4.0172 | 2.0062
- 64 | 1.16e-03 | 4.1081 | 2.0385
+ 8 | 7.60e-02 |
+ 16 | 1.91e-02 | 3.9765 | 1.9915
+ 32 | 4.72e-03 | 4.0514 | 2.0184
+ 64 | 1.18e-03 | 4.0055 | 2.0020
---------------------------------------------
- Once upon a time, a happy little test passed.
+ Happy little convergence test!
randomTensorMesh: Interpolation Fy
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 1.45e-01 |
- 16 | 5.01e-02 | 2.8858 | 1.2456
- 32 | 8.71e-03 | 5.7525 | 3.4043
- 64 | 1.85e-03 | 4.7102 | 1.9579
+ 8 | 7.84e-02 |
+ 16 | 2.69e-02 | 2.9183 | 1.3569
+ 32 | 1.20e-02 | 2.2396 | 1.1828
+ 64 | 2.31e-03 | 5.1989 | 2.2465
---------------------------------------------
- The test be workin!
+ You get a gold star!
@@ -1269,24 +1269,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.41e-02 |
- 16 | 1.78e-02 | 4.1730 | 2.0611
- 32 | 4.55e-03 | 3.9004 | 1.9636
- 64 | 1.11e-03 | 4.0930 | 2.0332
+ 8 | 7.61e-02 |
+ 16 | 1.92e-02 | 3.9650 | 1.9873
+ 32 | 4.79e-03 | 4.0044 | 2.0016
+ 64 | 1.18e-03 | 4.0544 | 2.0195
---------------------------------------------
- And then everyone was happy.
+ Testing is important.
randomTensorMesh: Interpolation Fz
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 7.63e-02 |
- 16 | 4.19e-02 | 1.8225 | 1.1195
- 32 | 9.07e-03 | 4.6176 | 1.4900
- 64 | 2.70e-03 | 3.3657 | 2.0946
+ 8 | 2.11e-01 |
+ 16 | 3.01e-02 | 7.0210 | 2.8286
+ 32 | 1.37e-02 | 2.1966 | 1.3516
+ 64 | 2.88e-03 | 4.7606 | 1.9822
---------------------------------------------
- Go Test Go!
+ Not just a pretty face Rowan
@@ -1317,24 +1317,24 @@ Test Results
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 6.98e-02 |
- 16 | 1.79e-02 | 3.8902 | 1.9598
- 32 | 4.41e-03 | 4.0681 | 2.0244
- 64 | 1.18e-03 | 3.7418 | 1.9037
+ 8 | 7.03e-02 |
+ 16 | 1.88e-02 | 3.7326 | 1.9002
+ 32 | 4.77e-03 | 3.9468 | 1.9807
+ 64 | 1.12e-03 | 4.2670 | 2.0932
---------------------------------------------
- You are awesome.
+ Yay passed!
randomTensorMesh: Interpolation N
_____________________________________________
h | error | e(i-1)/e(i) | order
~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
- 8 | 2.76e-01 |
- 16 | 4.47e-02 | 6.1758 | 2.9429
- 32 | 1.13e-02 | 3.9576 | 1.7499
- 64 | 4.59e-03 | 2.4595 | 2.1993
+ 8 | 4.98e-02 |
+ 16 | 3.64e-02 | 1.3674 | 0.3550
+ 32 | 1.19e-02 | 3.0611 | 2.4683
+ 64 | 3.98e-03 | 2.9868 | 1.4933
---------------------------------------------
- And then everyone was happy.
+ Yay passed!
@@ -1426,7 +1426,7 @@ Test Results
16 | 2.42e-03 |
32 | 6.06e-04 | 4.0001 | 2.0000
---------------------------------------------
- Once upon a time, a happy little test passed.
+ That was easy!
uniformLOM: Edge Inner Product - Isotropic
@@ -1446,7 +1446,7 @@ Test Results
16 | 2.81e-03 |
32 | 7.11e-04 | 3.9432 | 1.9794
---------------------------------------------
- The test be workin!
+ Testing is important.
@@ -1480,7 +1480,7 @@ Test Results
16 | 6.29e-04 |
32 | 1.57e-04 | 3.9978 | 1.9992
---------------------------------------------
- You are awesome.
+ Yay passed!
uniformLOM: Face Inner Product - Isotropic
@@ -1490,7 +1490,7 @@ Test Results
16 | 6.29e-04 |
32 | 1.57e-04 | 3.9978 | 1.9992
---------------------------------------------
- Yay passed!
+ Once upon a time, a happy little test passed.
rotateLOM: Face Inner Product - Isotropic
@@ -1500,7 +1500,7 @@ Test Results
16 | 3.08e-04 |
32 | 7.07e-05 | 4.3564 | 2.1231
---------------------------------------------
- And then everyone was happy.
+ The test be workin!
@@ -1534,7 +1534,7 @@ Test Results
16 | 6.99e-03 |
32 | 1.75e-03 | 3.9996 | 1.9998
---------------------------------------------
- Yay passed!
+ Awesome, Rowan, just awesome.
uniformLOM: Edge Inner Product - Anisotropic
@@ -1544,7 +1544,7 @@ Test Results
16 | 6.99e-03 |
32 | 1.75e-03 | 3.9996 | 1.9998
---------------------------------------------
- Happy little convergence test!
+ You deserve a pat on the back!
rotateLOM: Edge Inner Product - Anisotropic
@@ -1588,7 +1588,7 @@ Test Results
16 | 2.68e-03 |
32 | 6.69e-04 | 3.9982 | 1.9993
---------------------------------------------
- You are awesome.
+ Go Test Go!
uniformLOM: Face Inner Product - Anisotropic
@@ -1598,7 +1598,7 @@ Test Results
16 | 2.68e-03 |
32 | 6.69e-04 | 3.9982 | 1.9993
---------------------------------------------
- Testing is important.
+ Yay passed!
rotateLOM: Face Inner Product - Anisotropic
@@ -1608,7 +1608,7 @@ Test Results
16 | 2.15e-03 |
32 | 5.25e-04 | 4.0845 | 2.0302
---------------------------------------------
- The test be workin!
+ That was easy!
@@ -1642,7 +1642,7 @@ Test Results
16 | 6.79e-03 |
32 | 1.70e-03 | 3.9996 | 1.9998
---------------------------------------------
- And then everyone was happy.
+ Not just a pretty face Rowan
uniformLOM: Edge Inner Product - Full Tensor
@@ -1652,7 +1652,7 @@ Test Results
16 | 6.79e-03 |
32 | 1.70e-03 | 3.9996 | 1.9998
---------------------------------------------
- You get a gold star!
+ Yay passed!
rotateLOM: Edge Inner Product - Full Tensor
@@ -1662,7 +1662,7 @@ Test Results
16 | 7.49e-03 |
32 | 1.89e-03 | 3.9617 | 1.9861
---------------------------------------------
- The test be workin!
+ That was easy!
@@ -1696,7 +1696,7 @@ Test Results
16 | 3.10e-03 |
32 | 7.74e-04 | 3.9981 | 1.9993
---------------------------------------------
- Testing is important.
+ That was easy!
uniformLOM: Face Inner Product - Full Tensor
@@ -1706,7 +1706,7 @@ Test Results
16 | 3.10e-03 |
32 | 7.74e-04 | 3.9981 | 1.9993
---------------------------------------------
- Go Test Go!
+ Once upon a time, a happy little test passed.
rotateLOM: Face Inner Product - Full Tensor
@@ -1716,7 +1716,7 @@ Test Results
16 | 2.54e-03 |
32 | 6.23e-04 | 4.0741 | 2.0265
---------------------------------------------
- You are awesome.
+ Go Test Go!
@@ -1763,7 +1763,7 @@ Test Results
64 | 2.18e-02 | 4.0001 | 2.0000
128 | 5.46e-03 | 4.0000 | 2.0000
---------------------------------------------
- Go Test Go!
+ Once upon a time, a happy little test passed.
uniformLOM: 2D Edge Inner Product - Isotropic
@@ -1777,7 +1777,7 @@ Test Results
64 | 2.18e-02 | 4.0001 | 2.0000
128 | 5.46e-03 | 4.0000 | 2.0000
---------------------------------------------
- And then everyone was happy.
+ Yay passed!
rotateLOM: 2D Edge Inner Product - Isotropic
@@ -1791,7 +1791,7 @@ Test Results
64 | 2.00e-02 | 4.0155 | 2.0056
128 | 5.00e-03 | 4.0038 | 2.0014
---------------------------------------------
- You get a gold star!
+ You are awesome.
@@ -1829,7 +1829,7 @@ Test Results
64 | 2.53e-02 | 4.0000 | 2.0000
128 | 6.32e-03 | 4.0000 | 2.0000
---------------------------------------------
- And then everyone was happy.
+ Not just a pretty face Rowan
uniformLOM: 2D Face Inner Product - Isotropic
@@ -1843,7 +1843,7 @@ Test Results
64 | 2.53e-02 | 4.0000 | 2.0000
128 | 6.32e-03 | 4.0000 | 2.0000
---------------------------------------------
- Yay passed!
+ Happy little convergence test!
rotateLOM: 2D Face Inner Product - Isotropic
@@ -1857,7 +1857,7 @@ Test Results
64 | 2.30e-02 | 4.0132 | 2.0048
128 | 5.74e-03 | 4.0009 | 2.0003
---------------------------------------------
- And then everyone was happy.
+ Happy little convergence test!
@@ -1895,7 +1895,7 @@ Test Results
64 | 1.30e-01 | 4.0000 | 2.0000
128 | 3.24e-02 | 4.0000 | 2.0000
---------------------------------------------
- You are awesome.
+ And then everyone was happy.
uniformLOM: 2D Face Inner Product - Anisotropic
@@ -1909,7 +1909,7 @@ Test Results
64 | 1.30e-01 | 4.0000 | 2.0000
128 | 3.24e-02 | 4.0000 | 2.0000
---------------------------------------------
- And then everyone was happy.
+ Yay passed!
rotateLOM: 2D Face Inner Product - Anisotropic
@@ -1923,7 +1923,7 @@ Test Results
64 | 1.28e-01 | 4.0007 | 2.0003
128 | 3.19e-02 | 3.9975 | 1.9991
---------------------------------------------
- And then everyone was happy.
+ You get a gold star!
@@ -1961,7 +1961,7 @@ Test Results
64 | 3.66e-02 | 4.0000 | 2.0000
128 | 9.14e-03 | 4.0000 | 2.0000
---------------------------------------------
- And then everyone was happy.
+ Not just a pretty face Rowan
uniformLOM: 2D Edge Inner Product - Anisotropic
@@ -2027,7 +2027,7 @@ Test Results
64 | 1.52e-01 | 4.0000 | 2.0000
128 | 3.80e-02 | 4.0000 | 2.0000
---------------------------------------------
- Yay passed!
+ You get a gold star!
uniformLOM: 2D Face Inner Product - Full Tensor
@@ -2041,7 +2041,7 @@ Test Results
64 | 1.52e-01 | 4.0000 | 2.0000
128 | 3.80e-02 | 4.0000 | 2.0000
---------------------------------------------
- Once upon a time, a happy little test passed.
+ Go Test Go!
rotateLOM: 2D Face Inner Product - Full Tensor
@@ -2055,7 +2055,7 @@ Test Results
64 | 1.45e-01 | 4.0074 | 2.0027
128 | 3.62e-02 | 3.9984 | 1.9994
---------------------------------------------
- And then everyone was happy.
+ You get a gold star!
@@ -2093,7 +2093,7 @@ Test Results
64 | 3.78e-03 | 4.0001 | 2.0000
128 | 9.46e-04 | 4.0000 | 2.0000
---------------------------------------------
- You get a gold star!
+ The test be workin!
uniformLOM: 2D Edge Inner Product - Full Tensor
@@ -2107,7 +2107,7 @@ Test Results
64 | 3.78e-03 | 4.0001 | 2.0000
128 | 9.46e-04 | 4.0000 | 2.0000
---------------------------------------------
- Testing is important.
+ That was easy!
rotateLOM: 2D Edge Inner Product - Full Tensor
@@ -2121,7 +2121,7 @@ Test Results
64 | 1.98e-02 | 3.8708 | 1.9526
128 | 5.03e-03 | 3.9418 | 1.9789
---------------------------------------------
- Testing is important.
+ You are awesome.
@@ -2166,7 +2166,7 @@ Test Results
16 | 3.95e-02 | 3.7462 | 1.9054
32 | 1.00e-02 | 3.9364 | 1.9769
---------------------------------------------
- You are awesome.
+ And then everyone was happy.
uniformLOM: Curl
@@ -2178,7 +2178,7 @@ Test Results
16 | 3.95e-02 | 3.7462 | 1.9054
32 | 1.00e-02 | 3.9364 | 1.9769
---------------------------------------------
- The test be workin!
+ Happy little convergence test!
rotateLOM: Curl
@@ -2190,7 +2190,7 @@ Test Results
16 | 1.70e-02 | 5.2040 | 2.3796
32 | 3.77e-03 | 4.5126 | 2.1740
---------------------------------------------
- Happy little convergence test!
+ Go Test Go!
@@ -2234,7 +2234,7 @@ Test Results
16 | 1.19e-01 | 3.7462 | 1.9054
32 | 3.01e-02 | 3.9364 | 1.9769
---------------------------------------------
- That was easy!
+ Go Test Go!
uniformLOM: Face Divergence
@@ -2245,7 +2245,7 @@ Test Results
16 | 1.19e-01 | 3.7462 | 1.9054
32 | 3.01e-02 | 3.9364 | 1.9769
---------------------------------------------
- Once upon a time, a happy little test passed.
+ Well done Rowan!
rotateLOM: Face Divergence
@@ -2256,7 +2256,7 @@ Test Results
16 | 9.53e-04 | 9.5374 | 3.2536
32 | 2.75e-04 | 3.4594 | 1.7905
---------------------------------------------
- That was easy!
+ Testing is important.
@@ -2313,7 +2313,7 @@ Test Results
32 | 2.01e-02 | 3.9364 | 1.9769
64 | 5.04e-03 | 3.9841 | 1.9943
---------------------------------------------
- The test be workin!
+ Go Test Go!
rotateLOM: Face Divergence 2D
@@ -2325,7 +2325,7 @@ Test Results
32 | 2.01e-02 | 3.9364 | 1.9769
64 | 5.57e-03 | 3.6062 | 1.8505
---------------------------------------------
- Go Test Go!
+ And then everyone was happy.
@@ -2370,7 +2370,7 @@ Test Results
16 | 1.34e-04 | 3.9116 | 1.9678
32 | 3.39e-05 | 3.9578 | 1.9847
---------------------------------------------
- Yay passed!
+ The test be workin!
uniformLOM: Nodal Gradient
@@ -2382,7 +2382,7 @@ Test Results
16 | 1.34e-04 | 3.9116 | 1.9678
32 | 3.39e-05 | 3.9578 | 1.9847
---------------------------------------------
- You get a gold star!
+ Happy little convergence test!
rotateLOM: Nodal Gradient
@@ -2394,7 +2394,7 @@ Test Results
16 | 1.80e-04 | 3.6283 | 1.8593
32 | 4.66e-05 | 3.8703 | 1.9525
---------------------------------------------
- You get a gold star!
+ The test be workin!
@@ -2439,7 +2439,7 @@ Test Results
16 | 1.34e-04 | 3.9116 | 1.9678
32 | 3.39e-05 | 3.9578 | 1.9847
---------------------------------------------
- Go Test Go!
+ Awesome, Rowan, just awesome.
uniformLOM: Nodal Gradient 2D
@@ -2451,7 +2451,7 @@ Test Results
16 | 1.34e-04 | 3.9116 | 1.9678
32 | 3.39e-05 | 3.9578 | 1.9847
---------------------------------------------
- You are awesome.
+ You get a gold star!
rotateLOM: Nodal Gradient 2D
@@ -2463,10 +2463,183 @@ Test Results
16 | 1.80e-04 | 3.6343 | 1.8617
32 | 4.64e-05 | 3.8804 | 1.9562
---------------------------------------------
- Go Test Go!
+ You deserve a pat on the back!
+
+
+
+
+
+
+
+
+ | test_optimizers.TestOptimizers |
+ 4 |
+ 4 |
+ 0 |
+ 0 |
+ Detail |
+
+
+
+ test_GN_Rosenbrock |
+
+
+
+
+
+
+
+
+ |
+
+
+
+ test_GN_quadratic |
+
+
+
+
+
+
+
+
+ |
+
+
+
+ test_ProjGradient_quadratic1Bound |
+
+
+
+
+
+
+
+
+ |
+
+
+
+ test_ProjGradient_quadraticBounded |
+
+
+
+
+
+
@@ -2476,60 +2649,80 @@ Test Results
|
| test_Solver.TestSolver |
- 10 |
- 10 |
+ 14 |
+ 14 |
0 |
0 |
- Detail |
+ Detail |
-
+
test_directDiagonal_1 |
pass |
-
+
test_directDiagonal_M |
pass |
-
+
test_directFactored_1 |
pass |
-
+
test_directFactored_M |
pass |
-
- test_directLower_1 |
+
+ test_directLower_1_fortran |
pass |
-
- test_directLower_M |
+
+ test_directLower_1_python |
pass |
-
+
+ test_directLower_M_fortran |
+ pass |
+
+
+
+ test_directLower_M_python |
+ pass |
+
+
+
test_directSpsolve_1 |
pass |
-
+
test_directSpsolve_M |
pass |
-
- test_directUpper_1 |
+
+ test_directUpper_1_fortran |
pass |
-
- test_directUpper_M |
+
+ test_directUpper_1_python |
+ pass |
+
+
+
+ test_directUpper_M_fortran |
+ pass |
+
+
+
+ test_directUpper_M_python |
pass |
@@ -2539,40 +2732,40 @@ Test Results
7 |
0 |
0 |
- Detail |
+ Detail |
-
+
test_area_3D |
pass |
-
+
test_edge_2D |
pass |
-
+
test_edge_3D |
pass |
-
+
test_vectorCC_2D |
pass |
-
+
test_vectorN_2D |
pass |
-
+
test_vol_2D |
pass |
-
+
test_vol_3D |
pass |
@@ -2583,25 +2776,25 @@ Test Results
2 |
0 |
0 |
- Detail |
+ Detail |
-
+
test_orderBackward |
-
- |
-
+
test_orderForward |
-
- |
-
+
test_simpleFail |
-
- |
-
+
test_simpleFunction |
-
- |
-
+
test_simplePass |
-
- |
-
+
test_indexCube_2D |
pass |
-
+
test_indexCube_3D |
pass |
-
+
test_invXXXBlockDiagonal |
pass |
-
+
test_mkvc1 |
pass |
-
+
test_mkvc2 |
pass |
-
+
test_mkvc3 |
pass |
-
+
test_ndgrid_2D |
pass |
-
+
test_ndgrid_3D |
pass |
| Total |
- 99 |
- 99 |
+ 107 |
+ 107 |
0 |
0 |
|
diff --git a/notebooks/Bound Constraint Problem.ipynb b/notebooks/Bound Constraint Problem.ipynb
new file mode 100644
index 00000000..efaa487c
--- /dev/null
+++ b/notebooks/Bound Constraint Problem.ipynb
@@ -0,0 +1,620 @@
+{
+ "metadata": {
+ "name": ""
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
+ {
+ "cells": [
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "import SimPEG\n",
+ "from SimPEG.mesh import TensorMesh\n",
+ "from SimPEG.regularization import Regularization\n",
+ "import SimPEG.inverse as inverse\n",
+ "from SimPEG.inverse import Minimize\n",
+ "import numpy as np\n",
+ "import scipy.sparse as sp"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [],
+ "prompt_number": 1
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "from SimPEG.forward.LinearProblem import example as LinExample\n",
+ "\n",
+ "prob, m_true = LinExample(1000)\n",
+ "M = prob.mesh\n",
+ "\n",
+ "reg = Regularization(M)\n",
+ "opt = inverse.InexactGaussNewton(maxIter=100,maxStep=0.2,debug=False,LSreduction=0.3,maxIterLS=50)\n",
+ "opt.remember('f', ('norm_g', lambda M:np.linalg.norm(M.g)))\n",
+ "inv = inverse.Inversion(prob,reg,opt,beta0=1e-4,maxIter=2,debug=False)\n",
+ "m0 = np.zeros_like(m_true)\n",
+ "\n",
+ "mrec = inv.run(m0)\n",
+ "\n",
+ "plt.plot(M.vectorCCx, m_true, 'b-')\n",
+ "plt.plot(M.vectorCCx, mrec, 'r-')\n",
+ "\n",
+ "figure()\n",
+ "plot(opt.recall('f'))\n"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "====================== Inexact Gauss Newton ======================\n",
+ " # beta phi_d phi_m f |proj(x-g)-x| LS \n",
+ "-----------------------------------------------------------------\n",
+ " 0 1.00e-04 9.62e+02 0.00e+00 9.62e+02 1.02e+04 0 \n",
+ " 1 1.00e-04 6.63e+02 7.73e+02 6.63e+02 8.47e+03 0 \n",
+ " 2 1.00e-04 4.15e+02 4.39e+03 4.15e+02 6.63e+03 0 \n",
+ " 3 1.00e-04 2.26e+02 6.85e+03 2.27e+02 4.80e+03 0 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 4 1.00e-04 9.74e+01 1.07e+04 9.84e+01 3.02e+03 0 \n",
+ " 5 1.00e-04 2.42e+01 1.64e+04 2.58e+01 1.35e+03 0 \n",
+ " 6 1.00e-04 1.06e+00 2.34e+04 3.39e+00 9.21e+00 0 \n",
+ "------------------------- STOP! -------------------------\n",
+ "1 : |fc-fOld| = 2.2421e+01 <= tolF*(1+|f0|) = 9.6325e+01\n",
+ "0 : |xc-x_last| = 2.1351e+00 <= tolX*(1+|x0|) = 1.0000e-01\n",
+ "0 : |proj(x-g)-x| = 9.2073e+00 <= tolG = 1.0000e-01\n",
+ "0 : |proj(x-g)-x| = 9.2073e+00 <= 1e3*eps = 1.0000e-02\n",
+ "0 : maxIter = 100 <= iter = 6\n",
+ "1 : phi_d = 1.0595e+00 <= phi_d_target = 2.0000e+01 \n",
+ "------------------------- DONE! -------------------------\n",
+ "-------------------------STOP!-------------------------\n",
+ "0 : maxIter = 2 <= iter = 1\n",
+ "1 : phi_d = 1.0595e+00 <= phi_d_target = 2.0000e+01 \n",
+ "-------------------------DONE!-------------------------\n"
+ ]
+ },
+ {
+ "metadata": {},
+ "output_type": "pyout",
+ "prompt_number": 3,
+ "text": [
+ "[]"
+ ]
+ },
+ {
+ "metadata": {},
+ "output_type": "display_data",
+ "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX+P/DXYcncwgVzSUXBBVAQBEGvgriAGmndopLW\nq9VXLZUyy6/lfRhere79ZWXWVatrt+916RZa7imYgEkC7hum4oJaLkjiigLz+f3xkVGcGWZfz+v5\neMxDmDlzzpvjzGs+8zmf8zmKEEKAiIhUw8vZBRARkWMx+ImIVIbBT0SkMgx+IiKVYfATEakMg5+I\nSGWsDv7Ro0ejZcuWCAsL0/t4dnY2/Pz8EBkZicjISMycOdPaTRIRkRV8rF3BqFGjMGHCBDz33HMG\nl+nfvz9Wrlxp7aaIiMgGrG7xx8XFoWnTpnUuw3PEiIhch937+BVFQV5eHiIiIjBp0iQUFxfbe5NE\nRFQHuwd/z549cfLkSRQWFiI0NBRpaWn23iQREdVF2MCxY8dE9+7djS6n0WjE/fffLyoqKnQeCwoK\nEgB444033ngz4xYUFGR2Ztu9xX/27FltH/+qVasQHh6OevXq6SxXXFwMIQRvQmD69OlOr8FVbtwX\n3BfcF3XfLOk+t3pUT2pqKnJyclBaWop27dohPT0dlZWVAIAxY8YgIyMD8+bNg4+PD8LDwzF79mxr\nN0lERFawOviXLl1a5+OvvPIKXnnlFWs3Q0RENsIzd11QQkKCs0twGdwXt3Ff3MZ9YR1FCCGcXQQg\nh326SClERG7Dkuxki5+ISGUY/EREKsPgJyJSGQY/EZHKMPiJiFSGwU9EpDIMfiIilWHwExGpDIOf\niEhlGPxERCrD4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqQyDn4hIZRj8REQqw+AnIlIZ\nBj8Rkcow+ImIVMbH2QWQyly5AmRmAhcuAD17ApGRgKI4uyoiVWGLnxxn0SIgKAiYPx/YsgV4/HEg\nKQk4dcrZlRGpClv85BjvvQf861/A+vVARIS8r6pK3t+vH7BpE9Cxo3NrJFIJBj/Z37ffAp9/Dvzy\nC9Cq1e37fXyAv/4VaNwYGDFCPt6okfPqJFIJRQghnF0EACiKAhcphWzpyBGgTx9gwwbZn6+PEMDz\nzwP33Qd8+qlj6yNyc5ZkJ4Of7EcIYMgQ2Y8/eXLdy/7xB9C9O5CRIT8oiMgklmQnD+6S/Xz3HXD2\nLJCWZnzZpk2BWbOAKVPkBwYR2Q2Dn+yjqgqYNg2YPRvw9TXtOc8+C5SWAj/+aN/aiFSOwU/2sXgx\n0KYNMGiQ6c/x9gZmzJA3IrIb9vGT7Wk0QHAw8MUXQP/+5j23uhro1An45hsgNtY+9RF5EPbxk2tY\nt06O0ImPN/+53t7AxInAnDm2r4uIALDFT/YwZAjw9NPAc89Z9vzyciAgAPj1V6BlS9vWRuRh2OIn\n5ysqAnbvBp580vJ1+PkBDz8MLFliu7qISMvq4B89ejRatmyJsLAwg8tMnToVgYGBiIqKwsGDB63d\nJLmyBQuAF18E6tWzbj2jRgFffcWhnUR2YHXwjxo1Cj/WMfyuoKAAmzdvxrZt2zB58mRMNnYiD7mv\nmzdlK33UKOvXFR8vZ/LcudP6dRFRLVYHf1xcHJo2bWrw8fz8fKSkpKBZs2ZITU1FUVGRtZskV7Vu\nnRzNExRk/bq8vIDUVDnPDxHZlN0naSsoKMCzzz6r/b1FixYoLi5GkC3CgexCowF+/hmorDTveWEf\nfI3S2Ofx+0bb1NGoXQrC/vY4fhn8nsVz9gcFAR062KYeIk9h9+AXQugccVYMvInfeecd7c8JCQlI\nSEiwY2VkSFER8OCD5g2jv6/yAv6z9SdM8fo3rtmqd0ZE4D9lAv99azeKG0eY/fTSUjk4aOVKG9VD\n5AKys7ORnZ1t1TpsMpzz+PHjGD58OPbu3avz2Ny5c1FVVYXXXnsNABAUFITi4mLdQjic02Xs3i1H\nYu7ebcaTvvxSzsBp666ZKVPklA8zZ5r91LVr5WSfa9fatiQiV+KSwzljY2OxbNkyXLhwAUuWLEFI\nSIi9N0lWEsKCnpWMDHlFLVtLSZHrJiKbsbqrJzU1FTk5OSgtLUW7du2Qnp6Oyludw2PGjEFMTAz6\n9euH6OhoNGvWDIsWLbK6aLIvs4O/rExeRGXZMtsXEx0NXL4MHD4MdO5s1lMVhaNBifSxOviXLl1q\ndJn3338f77//vrWbIgcxO/h/+AFITAQaNrR9MYoiDzisXWva9M53PZXBT6SLZ+6SDrODPyNDdsnY\ny7BhcqgoEdkEg590mBX8Fy8CW7YAycn2K2jwYLmNq1fNehpb/ET6MfhJh1lhuW6dPMu2cWO71YP7\n7gN69QI2bTL7qQx+Il0MftJhVot/9Wpg+HC71gNAdveYOS7TwnO+iDweg590mBz8VVXyMon27Oap\nkZQEZGWZ9RR29RDpx+AnHSYH/y+/yFNjH3jA7jUhLEwOGz150v7bIvJwDH7SYXLwO6qbB5CTtg0Y\nYFY/P1v8RPox+EmHycG/ahXw0EN2r0dr0CDgp59MXpzBT6Qfg590mBT8xcWy6yUqyiE1AQAGDpTB\nzzQnsgqDn3SYFPxr1sgzar0c+BLq3FnOGa1nkj992OIn0o/BTzpMCv7MTGDoUIfUo6Uot1v9Ji7O\n4CfSxeAnHUaDv7ISyM2VIexoZgQ/EenH4CcdRoM/Px/o1Anw93dYTVpxccDmzSY15dniJ9KPwU86\njAZ/VpacjdMZAgNlP/+JEyYtzuAn0sXgJx1Ggz8zU06c5gyKAvTrJy8KbMKiRKSLwU866gz+S5eA\nPXuAvn0dWlMtZgQ/W/xEuhj8pKPO4M/OBnr3BurXd2RJtZkY/ACDn0gfBj/pqDP4s7Kc181To0cP\noKQE+OOPOhdjVw+Rfgx+0mE0+J11YLeGjw8QEwPk5dW5GLt6iPRj8JNeeoP/1Cng3DkgIsLh9egw\no7uHiGpj8JMOgy3+jRvlCVSOnKbBEBOCny1+Iv1c4B1MrsZg8GdmOr+bp0bv3sCOHcCNGwYXYfAT\n6cfgJx16g18I1ziwW6NRI3n28O7dzq6EyO0w+EmH3uDfv1+GbceOTqlJr9hYOX2EAWzxE+nH4Ccd\neoPfmWfrGsLgJ7IIg5906A1+V+rmqWEk+IlIPwY/6dBpJd+8KWfEdMY0zHUJCQHOngUuXND7MFv8\nRPox+EmHTot/61aga1egWTOn1aSXtzcQHQ0UFBhchMFPpIvBTzp0gt+VhnHerY7uHk7ZQKQfg590\n6AS/K/bv1zAS/GzxE+li8JOOWsFfXg7s2wf86U9Orcmg2FjZ1WMg4Rn8RLoY/KSjVvBv2iRD/957\nnVqTQa1bAw0bAkeO6DzErh4i/Rj8pKNW8LtyN08NA9097Ooh0o/BTzpqBb8rH9itwfH8RGZh8JMO\nbfCXlABlZUB4uLNLqltMDFv8RGawOvhzc3MREhKCzp07Y+7cuTqPZ2dnw8/PD5GRkYiMjMTMmTOt\n3STZmTb4N24EBg1yjWmY6xIVJQ9A3zVTJ4OfSD8fa1eQlpaGBQsWICAgAEOGDEFqair8/f1rLdO/\nf3+sXLnS2k2Rg2iD3x26eQB5cLdTJ2DvXnlCFxHVyaqmXHl5OQAgPj4eAQEBSEpKQr6er9yCzS63\nIgTgBY1s8bv6gd0a0dHAtm217mKLn0g/q4K/sLAQwcHB2t9DQ0OxdevWWssoioK8vDxERERg0qRJ\nKC4utmaT5ABCAAGX9gL33QcEBDi7HNNERwPbt+vczeAn0mX3ztuePXvi5MmTKCwsRGhoKNLS0uy9\nSbKSEED4WTfp5qkRFaW3xU9Euqzq4+/VqxfeeOMN7e/79+/H0KFDay3TuHFj7c8vvPAC3n77bdy4\ncQP16tXTWd8777yj/TkhIQEJCQnWlEcWEgIIO5sJJI1zdimmCw8Hfv0VqKjQnmzGrh7yRNnZ2cjO\nzrZqHVYFv5+fHwA5sqd9+/bIzMzE9OnTay1z9uxZ3H///VAUBatWrUJ4eLje0AdqBz85j3KjAl1K\n84AB3zq7FNPVry9nEN2zRw7vvIXBT57m7kZxenq62euwelTPxx9/jDFjxqCyshITJ06Ev78/FixY\nAAAYM2YMMjIyMG/ePPj4+CA8PByzZ8+2dpNkZy1+/Rknm4Shy60PdrdR091zK/jZ1UOknyJcZMiN\noigc/eMi9j00BYdK7sWje8xvSTjV/PlywraFCwHIof0jR8p/iTyVJdnp4mfmkDO03peJ/a3d6MBu\nDT0je9iWINLF4Kfazp1D43NHcbRFrLMrMV9YGHD4MHDtGgB29RAZwuCn2jZuxJmu/aHx9nV2Jear\nV09eh3f3bgAc1UNkCIOfasvMxG/dEt23tXzXGbwMfiJdDH66TQhgwwacCnXz4L/Vz++2fwORnTH4\n6baDBwEvL5Tf38XZlVjurjN42eIn0sXgp9tuzcYpoLhva7l7d+DoUeDqVff9G4jsjMFPt61bB9ya\ncsNtQ/Oee2T479rFg7tEBjD4Sbp2Dfj5ZyApqfalF93RHd09DH4iXQx+kn76SQamn5/7B/+tkT1u\n/TcQ2RGDn6S1a4HkZADwjODfvp1dPUQGMPhJpuOaNcCDD2p/devgDw0FTpyAcuUyg59IDwY/AQcO\nyKQPDQXgAcHv6wuEheHeop3OroTIJTH46XY3z620d/vgB4DoaNy7fztb/ER6MPipVjcP4DnBX29v\nIYOfSA8Gv9qVlwM7dgADBmjv8pjg36d78XUiYvDT2rVAfDzQoIH2Lo8I/uBg+Jw9jUbV5c6uhMjl\nMPjVbvly4NFHa93lEcHv44ObIRHoVsFWP9HdGPxqdv06sGEDMGJErbs9IvgB3AyPRvjNbcYXJFIZ\nBr+abdggT3by9691t6cE/42waIQx+Il0MPjVTE83D+A5wX+zRy+E32DwE92Nwa9WlZXA6tXAI4/o\nPOQpwV8d2BlNNBeACxecXQqRS2Hwq9WmTUCXLsADD+g85CnBDy8v7LunZ60LsxARg1+9vvtObzcP\n4DnBryjAHt9oBj/RXXycXQA5QUWF7N/ftUvvw54U/Lt8ewHbljq7FCKXwha/Gq1dC/ToAbRrp/dh\nTwl+ANjNFj+RDga/Gi1aBDzzjMGHPSX4FQUo8e4IXL0KnDnj7HKIXAaDX23KyoCNG4HHHjO4iKcE\nPwAIKNorchGRxOBXm4wMYMgQwM/P4CKeEvzaK3D16sXgJ7oDg19tvvwSeP75OhfxlKmMtR9ebPET\n1cLgV5MdO4CzZ4GhQ+tczFNa/MCtD7Ga4PeUTzQiKzH41WTBAuCllwBv7zoX85Tg13b1tG0rfzh1\nytklEbkEjuNXi8uXgW+/BfbvN7qoJwW/9oeaVr+BIaxEasIWv1osXAgMHgy0aWN0UU8JfuCO3h0e\n4CXSYovfWjduyDNgT5wASksBX1+geXOga1cgONhot4pDVFYCH34oR/SYwFOCv9bfEB0NzJ3rtFqI\nXIlrBf+mTUCnTq7/dfzcOTnXzbffylZk587y1rw5UFUlHy8qAv74Qw6dfPFFICHBeWn67bdAYKBs\n9ZrAk4Jf2+KPirp9gNcT/jgiK1jd1ZObm4uQkBB07twZcw20qKZOnYrAwEBERUXh4MGDhlc2cCDQ\nvr18Y374IXD0qLXl2daePcBf/iJb81u3Am+8IUfJ7NolPwjmz5fDJVeuBA4flvf36QNMmAB06yZb\n3I4eWaLRAH//O/DmmyY/xZOyUbu7W7cG6tcHjh1zaj1ErsDqFn9aWhoWLFiAgIAADBkyBKmpqfC/\n44pOBQUF2Lx5M7Zt24b169dj8uTJWL16tf6VFRUB27cD33wDfPIJ8PrrQEQE8K9/AT17Wluq5Y4c\nAd56C9i8WYZ4cTHQrJnx57VtC4wfD7zyirza1f/+LzB7NvD550BYmP3rBuS+rF/f6BDOOwkBeHnA\n0R+dD6+afv7AQKfU43Rnz8prMGzeDBw6dHsai0aNgA4dZIOmXz95a97cqaXaxbVr8m++elVOVHjz\nJnDPPfL90aAB0KqV/FcFrAr+8vJyAEB8fDwAICkpCfn5+UhOTtYuk5+fj5SUFDRr1gypqamYNm2a\n4RUGB8vb00/L9Hn5ZdmKjoqSb9Yff5RdKo5y7hwwY4YMz0mTgK++Aho2NH89iiK7fBIT5YfYwIHy\nw+Ctt+QLz15u3gSmTZMHds1owntkix+4PbLniSecVo9T7NkDpKfLqTqGDQP695fdj61ayU/4S5eA\n48eBffuAf/4TeO45OYnf448DKSny25I7OX9efiPfu1f+7UVFwMmTMvhbtZLv4fr1gXr15DG669dv\nz+dUv74cANGpk/yW3q2bbKSFhgI+rtUzbg2r2nWFhYUIDg7W/h4aGoqtW7fWWqagoAChoaHa31u0\naIHi4mLjK1cUYN48ObfM11/L1kqXLrKr5fp1a8o27to1YObM2//ZBw/KkLYk9O/k5SXH0e/cCRQW\nAnFx9u16mD9ftuISEsx+qicEf60+fkAGf0GB0+pxuCtXgDFjgKQkID4eKCkBli4Fxo6VrfpOnWSD\nKiJCXolt2jRg/XrZ4HnzTfkh2a2bfGztWqC62tl/kX7V1TLop08HYmLk3/Xpp8DFi0BysmywHTok\nc+P4cTmkeds2YMsW+e/+/fL+69flt/ulS4Fnn5UfDCtWAE8+Kb/hDxwo99GaNcCtRq+7svtHmBAC\n4q5+bcVAqvj66ru3KYDn4CWexkten2P216+i3tdfY5LXx5jrlWbTWr1ENZ4R/0G6ZhrylL6Y5lWA\nY58FAp/ZdDMA2gJiNdLEx3gzMBYvey3ACq8/23QLbcRpbKv+GwZ7Z+OA3v1qWHU1MGeOTctxCp2X\nWUyM7EqsqnK71ltqqsmDsgAAPcQuLK1+HD8r8Zjs9SsuTfYDJpv67HoAHgLwEBqKK3hi1X/x4op0\ntMJYLPR6EV8ro3BKce4AjBbiHJLEegwR65AoNuAMWuNHZRh+VP6OPKUvKn+6B/jJ3LUqAJrfuoUD\nSNE+0lSUISYnH32y8/AnzEa0GIl9CEOWkoiNXonIRyyqFDPfaDbQpYtlz1PE3alshvLyciQkJGDn\nzp0AgAkTJmDo0KG1unrmzp2LqqoqvPbaawCAoKAgvS1+RVEwbdp07e/9+yegf/8E3Y1WVMDr4w/h\n9fFsiPAe0Lw5FWLQYOuaqBoNlBU/wHtmOkTjxtD8/QOI2N6Wr88MSkE+vJ8ZCU3ycGje/3+ylWEt\nIeD9+KMQ3cOgeWeGRavw8XH/Vn9ZGRAUJAdXaYWGAkuWyFauGxkwQH7pNOXLm7L+R3iPfg7VH34C\n8eRI2xWxaxe8/vU5vL77L0Rsb2hGvQiR/JChFpttVVdDKciH8uM6KOvXQSk+ApEwEJqhwyCShjp+\nJOD161DytkDJyoTXxkzgaDFEv3iIxCRoEofILmk7vYFycrKRk5Ot/X3mzHSdxrVRwkoREREiJydH\nHDt2THTt2lWcP3++1uP5+fmib9++orS0VCxevFgkJyfrXY/Zpdy8KcTXXwsRGipEz55CLFkixLVr\n5q3jyhUhFi4UIiREiF69hPjhByE0GvPWYQtlZUI88oj8Ow4ftn59//ynEOHhQly/bv263FhZmRB+\nfnfdOWqU3D9uJj5eiOxsExb84Qch7r9fiJ9/tl8xV6/K915cnBAtWwrxxhtCFBba/r1z8qQQX30l\nxJNPCtGsmXxNT5kid8SNG7bdlrXOnRPim2+EGD1aiAceECIgQIiXXhIiI0O+EO3Ikhi3Ovizs7NF\ncHCwCAoKEnPmzBFCCDF//nwxf/587TJTpkwRHTp0ED179hQHDhzQX4iln0HV1UKsWCHE4MFCNG0q\nd/x//yvE77/rLqvRyBfTokVCPPusEE2aCJGcLERWlnMC/+7aPvlECH9/Wb+l8vLkOn791Xa1uak/\n/tAT/J9/LsRzzzmlHmv06ydEbq6RhTIzhWjRQoawoxQVyTDu3FmIdu2EmDBBiOXLhTh71rz1aDRC\nFBfL1/748UIEB8uwT0kR4ssvhTh1yj7124NGI8T+/UJ89JEQQ4cK0aiREL17CzF9uhBbtghRWWnT\nzVmSnVZ19diSoijmf1252+nT8mSln36SB26qq+WQSl9fOcKlpEQO1+rbV46wefRReZTflWzfLg8m\nDRoEfPAB0Lix6c/du1dOy/DVV8CDD9qvRjdx8SIQEHDXcbi9e+VFaA4dclpdlujbF/jHP+S/eh08\nKA/gLlsmBw04mhDAgQPyHJbNm4G8PHlAtEsX2e3RurV8LTdsKN+L167JPrgTJ+QAh/375XszOlqe\n+5KYKLvjPGFccUWFzKP16+Ww7hMn5IHipCR569jRqtVbkp2eFfx3EkK+80+dkh8APj7yQ6BJE9tt\nw14uXQLS0oCsLHmU9c9/Nt5fmJkpL6f4ySfyg4NQXi67fi9duuPO6mqgaVMZNm40Vr1PH3lOY58+\neh68eBGIjQWmTAFGj3Z4bXpVV8sRMocPy9u5c3KiwKtX5RDmhg3lxYACAuQ5BMHBrtcIs5czZ+T7\ndcMGefPzkw22gQPlQZw7zoMyBYPf0+TkyBPA6teXb+rkZODee2svc/IkMGuWbGktXSrHaBMAGfht\n294V/IB8k732mtyfbiI2Vn6mx8be9YAQcrhl+/aci8gdaTTA7t2yl2LTJvltqWlTOQItLk5+A+rW\nDbjvPoOrsCQ73WtMm9r07y9fFMuXy3HJL74o3/kdOsgW1YED8iv+88/Ln93h24yD6X0/9O4tx327\nUfBrNAZ6Pb74Qn6r/e47h9dENuDlBURGytvrr8uhxjt2yHMFMjPl+RQVFfKD/dFH5Yl10dHyJDM/\nP9ltZgG2+N3J2bPyxK+SEjnrZ5cu8oNAJaeZm+vyZdm1fOXKXQ+sWQN89JHsSnMTUVFypo+oqDvu\nPHwY+NOfgNxcICTEabWRnWVny5M+T52SXUP79t1+rEcPKLt3s6uHqIbB4C8tlQP8y8pcY9psE0RG\nypk3IiNv3aHRyIO5TzwBTJzo1NrIwS5flo0/Hx+gUSMobduyq4eohs6UDTX8/YGWLWX3mKMmy7OS\nRnPX8f2FC2V33/jxTquJnKRxY9nvbwUPGCtFpF+dA6Fq+vndRK0ZU0tLgbfflnNZecJwR3I4vmrI\noxn8Bty7N/DLLw6txRq1Du6++Sbw1FNuN+0EuQ529ZDHMtjVA8gzodxoJjpt8BcUyBOBioqcXRK5\nMbb4yWPV2dXTvbs8qajmYiQuTqMBvBQhW/vp6XWO6yYyhsFPHs1gi9/bW85Jv3mzQ+uxlEYDNM5d\nIy8y8pe/OLsccnMMfvJYdXb1AHI4ZG6uw+qxhlJdBf8PpsjrJ7vZtQTI9TD4yWMZnQ7djYL/z5f+\nDU3zFm51tjG5Lp7ARR7rxg055NngWe2VlXIGyRMn5L+u6sYN/NawE/BdBtr8+e7JekjtLMlOtvjJ\nYxnt6vH1lcM6t2xxWE0W+b//wyHfbqiKYuiTbTD4Sd1cvbunqgp47z3MaTSN52qRzfClRB7LaIsf\ncP3gX7oUaN8eBff0c/trIJPrYPCTxzIp+GNi5NWfLl92SE1mqa4G3n0XmDbN8LTMRBbgS4nUrX59\noFcv12z1L18u51wfNIjBTzbFlxJ5LJNa/IC8vmtmpt3rMYsQwMyZwLRpgKIw+Mmm+FIij2Vyn3hi\noutdlGX1avkH3Bq3z+AnW+JLiTyaSS3+nj2B336TN1dwV2sfYPCTbfGlRB7L5Ba/tzcwcKDrtPqz\nsuTB5kcf1d7F4Cdb4kuJPJZZwx9dqbtn5kx5oZU7kp7BT7bElxJ5PJMP8GZlmbiwHeXmAqdPA08+\nWetuBj/ZEl9K5PFMyvLAQKBBA2D3brvXU6dZs4CpU3Vm4GTwky3xpURUY/hwYNUq522/oEBeWevZ\nZ3Ue0rnYOpEVGPzk0Uweyw8AI0Y4N/hnzQKmTAHuuUfnoVoXWyeyEqdlJo/m7S2nZzbp2iWVlUDL\nlsC+fUCbNnavrZbdu4Fhw4DiYnk28V18fICKCl6DhXRxWmYia/j6AkOHAmvWOH7b774LvP663tAH\n2MdPtsWXEnk0s7p6ANnPv3Kl3erRq6gI2LQJGDNG78NCyBv7+MlW2NVDHs3XF7h6VW+3uX5//AEE\nBAC//w40bGjX2rSefhro1g146y29D2s0sotHo3FMOeRe2NVDZK2mTYE+feRcOY7w66/Ahg3A+PEG\nF2E3D9kaX07k0czu6gGA1FTgm2/sUo+OWbOAiROB++4zuAiDn2yNXT3k0erVAy5dkv+arLwcaN9e\nXoS9SRO71YYjR+Q1f4uL5bz7BlRUyDIqKuxXCrkvdvUQ6WF2e8LPT07a9v33dqlHa9Ys2cVTR+gD\nbPGT7Vn8crp8+TIefvhhtG/fHo888giuXLmid7kOHTogPDwckZGRiImJsbhQIktY1NUDyO6exYtt\nXo/WkSNy9FBamtFFedYu2ZrFwT9v3jy0b98ehw8fRtu2bTF//ny9yymKguzsbOzcuRMFBQUWF0rk\nUCNGyJOqjh61z/rffhuYNEkeTDaCZ+2SrVn8ciooKMALL7yAevXqYfTo0cjPzze4LPvuyVksbvHf\ne6+cM+eLL2xeEwoLgZ9/Bl591aTF2dVDtmbxCeCFhYUIDg4GAAQHBxtszSuKgoEDB6Jjx44YPXo0\nRowYYekmicymKEBGhpkHd29p3PZ/MGBGAlZ3T4fGx9QTAYwQAgkzp6Dkwek4utq08wSuXmXwk23V\nGfyJiYk4c+aMzv2zZs0yuRW/ZcsWtG7dGkVFRRg+fDhiYmLQqlUrvcu+88472p8TEhKQkJBg0jaI\nDPmf/7FmBoZgBNwTjFOf/YCt7Z6wST09fl+PuJLT+DBkNDTLTX/eiy/aZPPkAbKzs5GdnW3VOiwe\nzvnYY49h2rRpiIyMxPbt2/Hee+8hIyOjzudMmjQJISEheOmll3QL4XBOckXffy/n0SkosP4I640b\nQHg48MEHcmoIIhtw6HDO2NhYLFy4ENevX8fChQvRu3dvnWWuXbuGy5cvAwDOnz+P9evXY+jQoZZu\nksjxHn5wKB3JAAAHuElEQVRY9rX89JP16/rgAyA4mKFPTmdx8I8bNw4lJSXo2rUrTp8+jbFjxwIA\nfvvtNyQnJwMAzpw5g7i4OERERGDkyJF4/fXX0a5dO9tUTuQIXl5yjvy//c26yzIeOwZ89BEwZ47t\naiOyEM/cJTKmqkp20fzjH8BDD5n/fI0GGDwYGDJEfogQ2RDP3CWyBx8f2U0zebK8WIu5Zs+Wz5s8\n2fa1EVmALX4iUwgh++ajooD0dNOft3WrfF5hIdChg93KI/WyJDsZ/ESm+v13ICICWLFCTq5mTEmJ\nnOJ5wQLLuoiITMCuHiJ7at0a+PJL4LHHjE/lcPIkMGgQ8MYbDH1yOQx+InMMHw789a8y1Pft079M\nXh7Qrx8wbpzJ0zIQOZLFUzYQqdbYsfKyjP37y+vkPvWU/DZw6BDw1VeyK+iLL+REb0QuiH38RJYq\nKZGjfdatA0pL5bV6H3kEmDABaN7c2dWRSvDgLhGRyvDgLhERGcXgJyJSGQY/EZHKMPiJiFSGwU9E\npDIMfiIilWHwExGpDIOfiEhlGPxERCrD4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqQyD\nn4hIZRj8REQqw+AnIlIZBj8Rkcow+ImIVIbBT0SkMgx+IiKVYfATEakMg5+ISGUY/EREKsPgJyJS\nGQY/EZHKWBz83333Hbp16wZvb2/s2LHD4HK5ubkICQlB586dMXfuXEs3R0RENmJx8IeFheH7779H\nfHx8nculpaVhwYIFyMrKwmeffYbS0lJLN6ka2dnZzi7BZXBf3MZ9cRv3hXUsDv7g4GB06dKlzmXK\ny8sBAPHx8QgICEBSUhLy8/Mt3aRq8EV9G/fFbdwXt3FfWMeuffyFhYUIDg7W/h4aGoqtW7fac5NE\nRGSET10PJiYm4syZMzr3v/vuuxg+fLjdiiIiIjsSVkpISBDbt2/X+9jFixdFRESE9vfx48eL1atX\n6102KChIAOCNN954482MW1BQkNm5XWeL31RCCL33+/n5AZAje9q3b4/MzExMnz5d77JHjhyxRSlE\nRGSExX3833//Pdq1a4etW7ciOTkZw4YNAwD89ttvSE5O1i738ccfY8yYMRg8eDBefvll+Pv7W181\nERFZTBGGmutEROSRHHrmriknc02dOhWBgYGIiorCwYMHHVmeQxnbF4sXL0aPHj3Qo0cPPPXUUzh0\n6JATqnQMU0/yKywshI+PD5YvX+7A6hzLlH1RWFiIXr16ISQkBAkJCY4t0IGM7Yvr16/j+eefR2Rk\nJPr3748VK1Y4oUr7Gz16NFq2bImwsDCDy5idm2YfFbBCRESEyMnJEcePHxddu3YV58+fr/V4fn6+\n6Nu3r7hw4YJYsmSJSE5OdmR5DmVsX+Tl5YmLFy8KIYT497//LZ555hlnlOkQxvaFEEJUVVWJAQMG\niOTkZJGRkeGEKh3D2L7QaDSie/fuIjMzUwgh9O4rT2FsX8ybN0+MGzdOCCHE8ePHRWBgoNBoNM4o\n1a5yc3PFjh07RPfu3fU+bkluOqzFb8rJXPn5+UhJSUGzZs2QmpqKoqIiR5XnUKbsiz59+mgPjicn\nJyMnJ8fhdTqCqSf5zZ07FykpKWjRooWjS3QYU/bFtm3bEB4ejsGDBwOAxx4zM2Vf+Pn54fLly6is\nrERZWRkaNGgARVGcUa5dxcXFoWnTpgYftyQ3HRb8ppzMVVBQgNDQUO3vLVq0QHFxsaNKdBhzT2z7\n/PPPPfa8CVP2xenTp7FixQqMGzcOADzyzQ2Yti/Wr18PRVEQFxeH4cOHY/369Y4u0yFM2Repqamo\nrq6Gv78/+vXrh8WLFzu6TJdgSW7aZDinrQghdIaGeuqb3FRZWVlYtGgR8vLynF2K07z66qt4//33\noSiK3teImlRUVGDXrl3IysrCtWvXkJiYiH379qF+/frOLs3hPv30U/j4+OD333/H3r17kZycjBMn\nTsDLS12TDluSmw7bQ7169ap10GH//v3o3bt3rWViY2Nx4MAB7e/nz59HYGCgo0p0GFP2BQDs2bMH\nY8eOxcqVK9GkSRNHlugwpuyL7du3Y+TIkejYsSOWLVuGl19+GStXrnR0qXZnyr7o06cPhg0bhlat\nWiEwMBDR0dHIzc11dKl2Z8q+yM3NxdNPP40GDRogNjYWbdq08ehBEIZYkpsOC/47T+Y6fvw4MjMz\nERsbW2uZ2NhYLFu2DBcuXMCSJUsQEhLiqPIcypR9UVJSgsceewyLFy9Gp06dnFGmQ5iyL44ePYpj\nx47h2LFjSElJwbx58zBixAhnlGtXpuyL3r17IycnB9euXUNZWRl27tyJvn37OqNcuzJlXwwaNAir\nVq2CRqPB0aNHUVZWVqt7SC0syU2HdvXUnMxVWVmJiRMnwt/fHwsWLAAAjBkzBjExMejXrx+io6PR\nrFkzLFq0yJHlOZSxfTFjxgyUlZVh7NixAABfX18UFBQ4s2S7MbYv1MTYvmjevDlGjRqF6OhotGjR\nAjNmzECjRo2cXLV9GNsXI0eOxIEDB7T7Ys6cOU6u2D5SU1ORk5OD0tJStGvXDunp6aisrARgeW7y\nBC4iIpVR11EQIiJi8BMRqQ2Dn4hIZRj8REQqw+AnIlIZBj8Rkcow+ImIVIbBT0SkMv8fsrmNQhBd\nxXsAAAAASUVORK5CYII=\n",
+ "text": [
+ ""
+ ]
+ },
+ {
+ "metadata": {},
+ "output_type": "display_data",
+ "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHo1JREFUeJzt3Xt01OWdx/H3BCQUuRRFgrswgpDmAkoC5oIYElAhoiEB\nRMoKWBOtRm1A4NSzqy3Y3SMFbSGgCdjdWC3HWo/dKiIKJDZAUJJwKSiEkggt4HIRRJ0IKMLsH49M\nQLkkk8k8c/m8zpkzYTK/zHfmcD558v09v+dxuN1uNyIiEtIibBcgIiItT2EvIhIGFPYiImFAYS8i\nEgYU9iIiYUBhLyISBi4a9rm5uURFRXHdddd5HnO5XGRnZ+N0OsnJyaG+vt7zvQULFhAdHU18fDwV\nFRWex2tqahgwYADXXnstjz/+eAu8DRERuZiLhv29997LO++8c85jxcXFOJ1Oamtr6d69O4sWLQLg\n0KFDFBUVUVZWRnFxMQUFBZ5jpk+fzmOPPUZ1dTWrV69mw4YNLfBWRETkQi4a9mlpaXTu3Pmcx6qq\nqsjLyyMyMpLc3FwqKysBqKysJDMzE6fTSXp6Om632zPq//vf/8748eO58sorGTNmjOcYERHxjyb3\n7Kurq4mNjQUgNjaWqqoqwIR9XFyc53kxMTFUVlZSV1dH165dPY/Hx8ezfv365tYtIiJN0OSwb8rq\nCg6Ho1nHi4iIb7Ru6gFJSUnU1NSQmJhITU0NSUlJAKSkpFBaWup53o4dO0hKSqJDhw4cPHjQ8/j2\n7dtJTU0978/u06cPH330UVNLEhEJa71796auru6iz2nyyD4lJYWSkhKOHz9OSUmJJ7iTk5NZsWIF\ne/bsoby8nIiICDp06ACYds8rr7zC4cOH+ctf/kJKSsp5f/ZHH32E2+3Wze1m5syZ1msIlJs+C30W\n+iwufmvMIPmiYT9hwgRuvPFGdu7cSY8ePXjhhRfIz89nz549xMTE8PHHH/Pggw8CEBUVRX5+PsOG\nDeOhhx6isLDQ83OeeeYZ5s6dS1JSEmlpadxwww2N+80iIiI+cdE2zh//+MfzPv7GG2+c9/EpU6Yw\nZcqU7z0eHx/Ppk2bvChPRER8QVfQBqiMjAzbJQQMfRYN9Fk00GfRNA632x0w02McDgcBVI6ISFBo\nTHZqZC8iEgYU9iIiYSDgwn7XLtsViIiEnoAL+2HDFPgiIr4WcGH/2GMwdKgCX0TEl5q8XEJLy88H\nh8ME/rvvQu/etisSEQl+ARf2AN9elMuwYQp8ERFfCMiwBxP4DocCX0TEFwI27AEeeODclk6fPrYr\nEhEJTgEd9gA//am5PzPCV+CLiDRdwIc9KPBFRJorKMIeTOCf6eGXlUF0tO2KRESCR9CEPcD99597\n0laBLyLSOEEV9gD33WfuFfgiIo0XdGEPJvDPbun86Ee2KxIRCWxBGfYAeXkm8G++WYEvInIpQRv2\nALm55l6BLyJycUEd9mAC/+yWTkyM7YpERAJP0Ic9wL33mvszI3wFvojIuUIi7MEE/tk9fAW+iEiD\nkAl7gJ/8pCHwS0shNtZ2RSIigSGkwh7gnnvM/ZkRvgJfRCQEwx4U+CIi3xWSYQ8m8M9u6cTF2a5I\nRMSekA17gMmTTeDfcgusWgXx8bYrEhGxI6TDHmDSJBP4t96qwBeR8BXyYQ8wcaK5v+UW09JR4ItI\nuAmLsAcT+Ge3dPr2tV2RiIj/hE3YA9x9t7k/09JR4ItIuAirsAcT+Gf38BX4IhIOwi7sAf7t38z9\nrbfCypXQr5/dekREWlpYhj2YwD97hK/AF5FQFrZhDzBhgrlX4ItIqAvrsAcT+GdG+CtXwnXX2a5I\nRMT3Irw98He/+x033ngjAwcOZOrUqQC4XC6ys7NxOp3k5ORQX1/vef6CBQuIjo4mPj6eioqK5lfu\nQz/+McyfD8OHw9attqsREfE9r8L+008/5amnnmLVqlVUV1ezc+dOVqxYQXFxMU6nk9raWrp3786i\nRYsAOHToEEVFRZSVlVFcXExBQYFP34QvjB9vAn/ECAW+iIQer8L+Bz/4AW63m88//5zjx49z7Ngx\nfvjDH1JVVUVeXh6RkZHk5uZSWVkJQGVlJZmZmTidTtLT03G73bhcLp++EV8YPx4KCxX4IhJ6vA77\n4uJievbsSbdu3Rg8eDApKSlUV1cT++16wrGxsVRVVQEm7OPOWnYyJibG871Ac9ddJvCHD4ctW2xX\nIyLiG16doP3kk0/Iz89n+/btdO7cmXHjxrFs2TLcbnejf4bD4Tjv47NmzfJ8nZGRQUZGhjclNstd\nd5mTtiNGwIoV0L+/30sQEbmg8vJyysvLm3SMV2FfVVVFamoqffr0AWDcuHGsXbuWpKQkampqSExM\npKamhqSkJABSUlIoLS31HL9jxw7P977r7LC3adw4c6/AF5FA892B8JNPPnnJY7xq46SlpbFhwwY+\n/fRTvvrqK95++22GDx9OSkoKJSUlHD9+nJKSElJTUwFITk5mxYoV7Nmzh/LyciIiIujQoYM3L+1X\n48bBs8+awP/b32xXIyLiPa9G9h07duSJJ55g9OjRHDt2jMzMTIYOHUpycjITJ04kJiaGAQMGMGfO\nHACioqLIz89n2LBhtGnThsWLF/v0TbSkO+8095mZ8M47kJBgtx4REW843E1ptLcwh8PRpL6/P/35\nz/Dww/D225CYaLsaEZEGjcnOsL+CtrHGjjX3t92mwBeR4KOwb4KxY80sHQW+iAQbhX0TjRnTEPjL\nl8OAAbYrEhG5NIW9F0aPNvdnRvgKfBEJdAp7LynwRSSYKOybYfToc1s6AwfarkhE5PwU9s2Uk2MC\nf+RIeOstuOEG2xWJiHyfwt4HsrPN/e23K/BFJDAp7H0kO9uM8BX4IhKIFPY+NGpUQ+AvWwYXWOtN\nRMTvvN6WUM4vKwv++79N4Afokv0iEoYU9i0gKwv+53/gjjsU+CISGBT2LSQrC0pKzL0CX0RsU9i3\noDvuaBjhf7sdr4iIFQr7FnbHHfDCC2aEr8AXEVsU9n5w++3w+9+bwF+/3nY1IhKOFPZ+MnKkCfxR\noxT4IuJ/Cns/UuCLiC0Kez8bORJefNEE/vvv265GRMKFwt6C226Dl14ySyy8957takQkHCjsLcnM\nNIGfk6PAF5GWp7C3KDMT/vAHBb6ItDyFvWUjRjQE/rp1tqsRkVClsA8AZwJ/9GioqLBdjYiEIoV9\ngBgxApYsgTFjFPgi4nsK+wAyfHhD4K9aZbsaEQklCvsAM3w4vPYaTJoECxeC2227IhEJBQ63O3Di\nxOFwEEDlWLVrl7nwavBgE/pt2tiuSEQCVWOyUyP7AHXttWY65v/9nxntHz5suyIRCWYK+wDWsSO8\n/jqkpkJyMnz4oe2KRCRYqY0TJJYsgUcfNZuhjBpluxoRCSSNyc7WfqpFmmniRIiONjN1amrg5z8H\nh8N2VSISLDSyDzL79pkF1OLj4Xe/g7ZtbVckIrbpBG0I6t4d1q6Fr7+GjAzYv992RSISDBT2Qahd\nO3jlFbPdYXIybNxouyIRCXReh/2XX37JPffcw49+9CPi4+OprKzE5XKRnZ2N0+kkJyeH+vp6z/MX\nLFhAdHQ08fHxVGg9gGZzOOAXv4D5883qmX/6k+2KRCSQeR32M2fOxOl0snXrVrZu3UpsbCzFxcU4\nnU5qa2vp3r07ixYtAuDQoUMUFRVRVlZGcXExBQUFPnsD4W7sWLO0ws9/Dr/8JZw+bbsiEQlEXod9\naWkp//Ef/0Hbtm1p3bo1nTp1oqqqiry8PCIjI8nNzaWyshKAyspKMjMzcTqdpKen43a7cblcPnsT\n4S4hAaqqoKwMxo2DL7+0XZGIBBqvwn7fvn2cOHGC/Px8UlJSmDNnDsePH6e6uprY2FgAYmNjqaqq\nAkzYx8XFeY6PiYnxfE98IyoK3n3XXIg1eDD885+2KxKRQOJV2J84cYKdO3cyduxYysvL2bZtG6++\n+mqTpk06NEnc5yIjoaTELKI2aJA2QxGRBl5dVNWnTx9iYmLIysoCYMKECbz00kskJSVRU1NDYmIi\nNTU1JCUlAZCSkkJpaann+B07dni+912zZs3yfJ2RkUFGRoY3JYYthwOmTzfz8EePhrlz4Sc/sV2V\niPhSeXk55eXlTTrG64uqRo0axeOPP05SUhIFBQUkJiZy5MgR9u7dy9y5c5kxYwa9evVixowZHDx4\nkPT0dFauXMmuXbuYNm0amzZt+n4xuqjKp2pqICvLbHk4Zw60amW7IhFpCY3JTq/DfufOnUyePJkT\nJ05wyy238OSTT3L69GkmTpzI5s2bGTBgAEuWLKF9+/YAFBYWsnDhQtq0acPixYtJS0vzqmBpmiNH\n4K67TIvnj3+ETp1sVyQivtaiYd8SFPYt4+RJs4haWRksXWrW2BGR0KHlEgSAyy6DZ5+FKVPgpptM\n6ItIeFHYh5EHHzTLLNx9Nzz3nO1qRMSf1MYJQx99ZNbEHzIEFiwwI38RCV5q48h59e4N778Pe/ea\nLQ+PHLFdkYi0NIV9mOrYEd54A5KSzMqZ27bZrkhEWpLCPoy1amUuupo506yNv2yZ7YpEpKWoZy+A\naevceSdMnQozZmjLQ5Fgonn20iR795otD/v1g+ef15aHIsFCJ2ilSXr0MFseHj8OQ4fCgQO2KxIR\nX1HYyzkuv9zsepWZaU7cnmcJIxEJQmrjyAW99hrk50NRkdkURUQCk3r20mybN5s+fm6u2fYwQn8L\nigQchb34xIEDZm38f/1XePFF0+oRkcChE7TiE926wV//akL+pptgzx7bFYlIUynspVHatoXf/94s\nopaaaubli0jwUBtHmuytt+Dee+Hpp+Gee2xXIyLq2UuL2b7dbHk4Zgz8+tfa8lDEJoW9tKgjR8wS\nC+3amS0PO3a0XZFIeNIJWmlRV14JK1fCNdeYPn5dne2KRORCFPbSLJddZi66euQRGDzYzNoRkcCj\nsBefeOgh08r58Y+huNh2NSLyXerZi0/V1ZkTt8OGwfz52vJQxB/Usxe/69MH1q+H3bthxAhteSgS\nKBT24nOdOsGbb8LAgZCSAjU1tisSEYW9tIhWrcxFV7/4BaSnw/LltisSCW/q2UuLe+89Mx9/2jSY\nPl1bHor4mi6qkoCxZw+MGgUJCbB4MURG2q5IJHToBK0EDKcT1q2D+npteShig8Je/Obyy+HVV+HW\nW82J282bbVckEj7UxhErXn0VHn4YFi2CsWNtVyMS3NSzl4C2cSPk5MD995tZOzpxK+Idhb0EvP37\nzZaHTqfZHKVdO9sViQQfnaCVgHf11VBebnbCSkuDfftsVyQSmhT2Yl3btmYj8/HjzYnb9ettVyQS\netTGkYCybJnZ8vC3v4VJk2xXIxIc1LOXoPThh+YCrHHj4KmntOWhyKW0eM/+1KlTJCYmkpWVBYDL\n5SI7Oxun00lOTg719fWe5y5YsIDo6Gji4+OpqKhozstKiOvXD6qqoLISsrPh6FHbFYkEv2aFfWFh\nIfHx8Ti+nTNXXFyM0+mktraW7t27s2jRIgAOHTpEUVERZWVlFBcXU1BQ0PzKJaR16WK2POzTB/r2\nNRuj6I8+Ee95Hfb79u1j+fLl3HfffZ4/H6qqqsjLyyMyMpLc3FwqKysBqKysJDMzE6fTSXp6Om63\nG5fL5Zt3ICGrTRuzAcqf/wyzZ8Ntt8GuXbarEglOXof9o48+ytNPP01ERMOPqK6uJjY2FoDY2Fiq\nqqoAE/ZxcXGe58XExHi+J3IpgwaZC7CGDoXkZJgzB06etF2VSHDxKuyXLVtG165dSUxMPOekQFNO\nrjp0uaQ0wWWXwWOPmV7+X/9qNkbRFE2RxmvtzUHvvfceS5cuZfny5Zw4cYIvvviCSZMmkZSURE1N\nDYmJidTU1JCUlARASkoKpaWlnuN37Njh+d53zZo1y/N1RkYGGRkZ3pQoIeraa+Htt+FPf4IxY8xy\nC7Nnm92xRMJFeXk55eXlTTqm2VMvV69ezTPPPMObb77J3Llz2bt3L3PnzmXGjBn06tWLGTNmcPDg\nQdLT01m5ciW7du1i2rRpbNq06fvFaOqlNMHRo2a0/9Zbprd/551aX0fCk9+WSzjTksnPz2fPnj3E\nxMTw8ccf8+CDDwIQFRVFfn4+w4YN46GHHqKwsNAXLythrnNneP55M8qfNQuysuCf/7RdlUhg0kVV\nEhK+/trseTtvHvz7v8OUKdDaqyalSPDRFbQSdurq4MEH4cgRM+q/wKkhkZCiVS8l7PTpA6tWmc3N\ns7KgoAC++MJ2VSL2Kewl5DgcZhG1bdvgyy/NFbivv267KhG71MaRkLd6NTzwAMTGwsKF0KOH7YpE\nfEttHBEgPR22bIHERHMrLIRTp2xXJeJfGtlLWNmxw5zAra83J3AHDLBdkUjzaWQv8h2xsWa5hYcf\nNgurTZ9ugl8k1CnsJew4HGY3rA8/hE8+MSdwly2zXZVIy1IbR8JeWZlp7SQkmH7+v/yL7YpEmkZt\nHJFGuPlm2LrVtHj694fnntMJXAk9GtmLnGX7dvjpT+Gbb8wJ3Ouvt12RyKVpZC/SRPHxsGYN5OXB\nLbeYVTWPHbNdlUjzKexFviMiAu6/Hz74APbuNRugv/OO7apEmkdtHJFLWLECHnrILKo2fz5062a7\nIpFzqY0j4gMjRphRfq9epof//PNw+rTtqkSaRiN7kSb44ANzArdVK1i82MzRF7FNI3sRH7vuOli3\nDu6+GzIy4PHH4fhx21WJXJrCXqSJIiIgP98srlZba34BlJbarkrk4tTGEWmmt94ya+2kpcFvfgNd\nu9quSMKN2jgifnD77WajlKgoM8ovKQGNWSTQaGQv4kObN5uNUn7wA3MCNzbWdkUSDjSyF/GzxER4\n/32480646SaYORNOnLBdlYjCXsTnWrWCn/0M/vY3M1Wzf3+zhr6ITWrjiLSwN94w4X/zzfD009Cl\ni+2KJNSojSMSALKzzQncTp3MOjsvvaQTuOJ/GtmL+NGGDeYK3M6dYdEiiI62XZGEAo3sRQLMDTdA\nVRXccQcMGgT/9V/w9de2q5JwoLAX8bPWreHRR2HjRqisNNshrl1ruyoJdWrjiFjkdsP//i9MmQK3\n3QZz5sAVV9iuSoKN2jgiAc7hgLFjzQncyEiziubLL+sErvieRvYiAWT9enMFbrduUFQEvXvbrkiC\ngUb2IkEmNdXM2Ln5ZkhJgV//Gk6etF2VhAKN7EUC1O7dZjvEffvM7liDBtmuSAJVY7JTYS8SwNxu\nePVVM3snOxt++Uu4+mrbVUmgURtHJMg5HDB+PGzfDm3aQHy82SWrstJ2ZRJsNLIXCSKffWbWy3/u\nObPGzs9+BuPGmZk8Er5abGS/d+9ehg4dSt++fcnIyODll18GwOVykZ2djdPpJCcnh/r6es8xCxYs\nIDo6mvj4eCoqKrx5WZGw98MfwrRpsHMnPPEEvPgi9OwJs2bB/v22q5NA5lXYX3bZZcybN49t27bx\n2muv8cQTT+ByuSguLsbpdFJbW0v37t1ZtGgRAIcOHaKoqIiysjKKi4spKCjw6ZsQCTetWkFWFqxa\nBWVlcOiQWjxycV6Ffbdu3UhISACgS5cu9O3bl+rqaqqqqsjLyyMyMpLc3Fwqv/1fV1lZSWZmJk6n\nk/T0dNxuNy6Xy3fvQiSMxcebOfm7d8PAgTBhgpm2uWQJfPWV7eokUDT7BG1dXR3btm0jOTmZ6upq\nYr/dhy02NpaqqirAhH1cXJznmJiYGM/3RMQ3zrR4amvV4pHva92cg10uF+PHj2fevHm0b9++SSdX\nHQ7HeR+fNWuW5+uMjAwyMjKaU6JI2DnT4snKMrN4nn3WjP5HjoSCAjPql+BWXl5OeXl50w5ye+nr\nr79233rrre558+Z5HhszZox706ZNbrfb7d6wYYN77Nixbrfb7V66dKm7oKDA87z+/fu7v/jii+/9\nzGaUIyIXcfSo2/2b37jdvXq53cnJbvcf/uB2nzhhuyrxlcZkp1dtHLfbTV5eHv369WPq1Kmex1NS\nUigpKeH48eOUlJSQmpoKQHJyMitWrGDPnj2Ul5cTERFBhw4dvHlpEfGCWjzi1Tz7iooKhgwZwvXX\nX+9px8yePZvBgwczceJENm/ezIABA1iyZAnt27cHoLCwkIULF9KmTRsWL15MWlra94vRPHsRv9m+\nHRYuhFdeUYsn2Gm5BBG5pKNH4YUXTG//qqt0oVYwUtiLSKOdOgXLl8OCBfDhh2ap5Qce0Fo8wUBr\n44hIo333Qq2DB3WhVijRyF5ELkgtnuCgNo6I+MSpU/DWW+aErlo8gUdtHBHxiVatYNQo0+IpLVWL\nJxhpZC8iXlGLJ3CojSMiLU4tHvvUxhGRFqcWT3DQyF5EfE4tHv9SG0dErFKLxz/UxhERqy7U4pk4\nUS0ef9PIXkT86ujRhk3Tr7rKLMA2bhy0aWO7suClNo6IBKwzLZ4FC2DbNrV4mkNtHBEJWGdaPKWl\navH4g0b2IhIw1OLxjto4IhKU1OJpGrVxRCQoqcXjexrZi0hQOLvF06WL+WUwZAgkJ0Pbtrars0tt\nHBEJOadOwcqVZsS/Zg3U1MCAASb4hwyBQYOgQwfbVfqXwl5EQp7LBe+/b4J/zRrYuBH69m0I/5tu\ngiuusF1ly1LYi0jYOXECqqoawv/996FXr4bwT0sLvRO9CnsRCXsnT8LmzQ3hX1Fhev5nwn/IELjm\nGnA4bFfqPYW9iMh3nD5tFmU7E/5r1ph5/GeHf0xMcIW/wl5E5BLcbqitPTf8jx0z7Z4z4X/99WY6\naKBS2IuIeGHPHli7tiH89++HwYMbwn/gwMC6qldhLyLiAwcPml7/mfCvqzPz+8+Ef0oKtGtnrz6F\nvYhIC/jsM1i3riH8t26FhISG8L/xRujUyX/1KOxFRPzgyy9h/fqG8K+uNid5z57rf9VVLff6CnsR\nEQu++go2bDDBv3at+Suge/dzT/p27+6711PYi4gEgG++Ma2es2f8dOx47nTP3r29n+6psBcRCUCn\nT8OOHQ3Bv3q1mQJ6dvjHx0NEI9clVtiLiAQBtxt27z535H/06Lltn4QEaN36/Mcr7EVEgtTHHzfM\n9V+71sz9HzSoIfyTkiAy0jxXYS8iEiKOHDl3rn9NDdxwgwn+//xPhb2ISEj64guzoufq1TB7doBt\nS7hmzRri4uKIjo5m4cKF/nxpEZGQ0rEjjBgBTz3VuOf7NeynTJnC4sWLKS0t5bnnnuPw4cP+fPmg\nUl5ebruEgKHPooE+iwb6LJrGb2H/+eefAzBkyBCuueYahg8fTqV2Dr4g/UduoM+igT6LBvosmsZv\nYV9dXU1sbKzn3/Hx8axfv95fLy8iEtb82sYRERFL3H7y2WefuRMSEjz/fuSRR9zLli075zm9e/d2\nA7rppptuujXh1rt370tm8AWux/K9Tt+u97lmzRqcTierVq1i5syZ5zynrq7OX+WIiIQVv4U9wPz5\n83nggQc4efIkBQUFdOnSxZ8vLyIStgLqoioREWkZAXGCVhdbNcjNzSUqKorrrrvOdilW7d27l6FD\nh9K3b18yMjJ4+eWXbZdkzYkTJ0hJSSEhIYHU1FTmzZtnuyTrTp06RWJiIllZWbZLsapnz55cf/31\nJCYmkpycfNHnBsTIPjExkcLCQq655hpGjBhBRUVF2LZ41q5dS/v27Zk8eTIffPCB7XKsOXDgAAcO\nHCAhIYHDhw+TnJzMli1b6NChg+3SrDh27Bjt2rXjq6++YuDAgbz++uv06dPHdlnW/Pa3v2Xjxo24\nXC6WLl1quxxrevXqxcaNG7niiisu+VzrI3tdbHWutLQ0OnfubLsM67p160ZCQgIAXbp0oW/fvmzY\nsMFyVfa0+3Y36/r6er755hsizyx3GIb27dvH8uXLue+++7SWFjT6M7Ae9rrYSi6lrq6Obdu2XfLP\n1FB2+vRp+vfvT1RUFI888gg9evSwXZI1jz76KE8//TQRjd3ZI4Q5HA6GDRtGTk7OJf/C0aclAc3l\ncjF+/HjmzZvH5ZdfbrscayIiItiyZQt1dXUUFRWxefNm2yVZsWzZMrp27UpiYqJG9cC6devYsmUL\ns2fPZtq0aRw4cOCCz7Ue9klJSezYscPz723btpGammqxIgkUJ0+eZOzYsUyaNIns7Gzb5QSEnj17\nMnLkyLBtdb733nssXbqUXr16MWHCBN59910mT55suyxrrr76agDi4uIYNWoUb7755gWfaz3sz77Y\n6h//+AerVq0iJSXFclVim9vtJi8vj379+jF16lTb5Vh1+PBhPvvsMwCOHDnCypUrw/aX31NPPcXe\nvXvZvXs3r7zyCsOGDeOll16yXZYVx44dw+VyAfDJJ5+wYsUKMjMzL/h8v15UdSG62KrBhAkTWL16\nNUeOHKFHjx786le/4t5777Vdlt+tW7eOJUuWeKaVAcyePfui/5lD1f79+7nnnns4deoU3bp1Y8aM\nGZ4RXbhzOBy2S7Dm4MGDjB49GoArr7yS6dOnX/RcTkBMvRQRkZZlvY0jIiItT2EvIhIGFPYiImFA\nYS8iEgYU9iIiYUBhLyISBhT2IiJhQGEvIhIG/h+tMGlORAtokQAAAABJRU5ErkJggg==\n",
+ "text": [
+ ""
+ ]
+ }
+ ],
+ "prompt_number": 3
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "FUN = SimPEG.tests.getQuadratic(sp.identity(2),np.array([-5,-5]))\n",
+ "# FUN = SimPEG.tests.Rosenbrock\n",
+ "\n",
+ "f = FUN(np.array([1,2]))\n",
+ "print f\n",
+ "n,l,u = 50,-10,10\n",
+ "I = np.zeros((n,n))\n",
+ "X = np.linspace(l,u,n)\n",
+ "for i, x in enumerate(X):\n",
+ " for j, y in enumerate(X):\n",
+ " f, g, H = FUN(np.array([x,y]))\n",
+ " I[i,j] = f\n",
+ "\n",
+ "colorbar(contourf(X,X, I.T))\n",
+ "\n",
+ "GN = inverse.GaussNewton()\n",
+ "xopt = GN.minimize(FUN,np.array([0,0]))\n",
+ "print xopt"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "(-12.5, array([-4., -3.]), <2x2 sparse matrix of type ''\n",
+ "\twith 2 stored elements (1 diagonals) in DIAgonal format>)\n",
+ "=========== Gauss Newton ==========="
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " # f |proj(x-g)-x| LS \n",
+ "-----------------------------------\n",
+ " 0 0.00e+00 7.07e+00 0 \n",
+ " 1 -2.50e+01 0.00e+00 0 \n",
+ "------------------------- STOP! -------------------------\n",
+ "0 : |fc-fOld| = 2.5000e+01 <= tolF*(1+|f0|) = 1.0000e-01\n",
+ "0 : |xc-x_last| = 7.0711e+00 <= tolX*(1+|x0|) = 1.0000e-01\n",
+ "1 : |proj(x-g)-x| = 0.0000e+00 <= tolG = 1.0000e-01\n",
+ "1 : |proj(x-g)-x| = 0.0000e+00 <= 1e3*eps = 1.0000e-02\n",
+ "0 : maxIter = 20 <= iter = 1\n",
+ "------------------------- DONE! -------------------------\n",
+ "[ 5. 5.]\n"
+ ]
+ },
+ {
+ "metadata": {},
+ "output_type": "display_data",
+ "png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAEACAYAAACEfgxsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH/lJREFUeJzt3X1wVOX5N/DvRppQJQ8FbZfMJAszGLIbiLDxSTbtr4QQ\neQnQkFTaqbQVh+RxYLEEBdJfLWgPTJVB7aAwEtLWtNNamF+lIjgtBEJdMsGyG95E01ChQGMZNCZ9\nHljeqinn+QN3zXvOvpy973PO9zOzM+7Zk3Nf7sKXK9ees2tTVVUFEREJlyS6ACIiuo2BTEQkCQYy\nEZEkGMhERJJgIBMRSYKBTEQkiZgDuaKiAna7HTk5OeFtwWAQZWVlcDgcKC8vx9WrV2NdhohImA8+\n+ADTp0/HxIkTUVRUhO3btwMAXnvtNUycOBF33HEHjh8/3uNnNm/ejMzMTGRnZ6OpqUnTOjEH8uLF\ni7Fv374e22pqauBwOHDmzBmkp6dj27ZtsS5DRCTMF77wBWzatAktLS3YuXMn1q5di2AwiJycHOza\ntQuFhYU99m9vb8fWrVtx8OBB1NTUoKqqStM6MQfy1KlTMWrUqB7bAoEAKisrkZKSgoqKCvj9/liX\nISISZsyYMZgyZQoA4J577sHEiRNx9OhROJ1OTJgwoc/+fr8fJSUlcDgcmDZtGlRVRTAYHHIdXWbI\nzc3NcDqdAACn04lAIKDHMkRECXf27Fm0tLQgPz9/wH0CgQBcLlf4flZWlqYc1CWQeTU2EZlRMBjE\nd77zHWzatAl33XXXgPv1l4E2m23I4w+LqboB5OXlobW1FW63G62trcjLy+uzz5h7v4iP/n5Tj+WJ\nyGTGjx+Ps2fPxnSM/2WzYeihwW0jRozoM2L49NNPsWDBAjz88MMoKysb9Oc9Hg8aGhrC90+fPt1v\nDvamSyB7PB7U1dXhueeeQ11dHQoKCvrs89Hfb2K3OmvQ42zDEj3K02Rv44OJW+xXCrBYGXo/DbsY\n2lt+AL8A8KjoSiTB5yLk73/vmyGRCgLQdq4D8PVeZ4apqorKykpMmjQJjz/+eL8/070rzs/PR3V1\nNdra2nDu3DkkJSUhNTV1yHVjHlksXLgQX/va1/D+++8jIyMDv/rVr+D1etHW1oasrCxcvHgRS5cu\njXWZhEpoGNNtb/GNX5LX4cOH8eqrr+LPf/4z3G433G439u7dizfeeAMZGRk4cuQI5s2bhzlz5gAA\n7HY7vF4viouLsWzZMrz00kua1om5Q96xY0e/23fv3h3roak7RXQBOmIYk+S+/vWv49atW/0+Vl5e\n3u/2FStWYMWKFRGtI+2VeqLGFUK64ylFiV9TFn3COFdIGXLic2E10gaypbiLBn9cSUQRAvTbGd+f\n8DLkxefCahjIJAbHFER9SBnIlhpXDEURXYAOGMZE/ZIykOkziugCdMAwJhqQdIHM7tjEGMZEg5Iu\nkOkziugC4oxhTDQkBjIRkSQYyJBwXKGILiDO2B0TaSJVIIv87ArSCcOYSDOpAlkEdsc6YhgTRcTy\ngUw6YRgTRUyaQOa4goisTppAFoHjCp2wOyaKiqUDmXTAMCaKmhSBzHEFzNMdE1HUpAhkEaQbV5gB\nu2OimFg2kImIZMNAloEiuoA4YHdMFDPhgcz5sQkwjMnkKioqYLfbkZOTE97217/+Fd/4xjcwZcoU\nlJaWorW1NfzY5s2bkZmZiezsbDQ1af2uawkCWQTOj4koEosXL8a+fft6bFu/fj0WLVqEkydP4rvf\n/S7Wr18PAGhvb8fWrVtx8OBB1NTUoKqqSvM6MX/rNMVIEV1AjNgdkwVMnToVFy5c6LFt5MiR6Ozs\nxK1bt9DZ2YlRo0YBAPx+P0pKSuBwOOBwOKCqKoLBIFJTU4dch4FMRBSF559/Hvn5+fjRj36EtLQ0\nHD16FAAQCATgcrnC+2VlZSEQCOCBBx4Y8piWHFlQnLA7JgurqKjA8uXL0dnZCa/Xi4qKCgCAqqp9\n9rXZbJqOKbRDFvGGnlTzY0V0AUTW8l8L+9/u+wjwtXfb8N7Qx2pqasJvf/tbDBs2DJWVldiwYQMA\nwOPxoKGhIbzf6dOnkZeXp6k+dshEZHlFdkDJ+fymxfTp07Fnzx4AwO7duzFz5kwAQH5+Purr69HW\n1gafz4ekpCRN82OAM2SKFscVZCELFy7EoUOH0NHRgYyMDKxfvx5r167FT3/6Uzz77LOYNGkSnnrq\nKQCA3W6H1+tFcXExkpOTUVtbq3kdm9rfwCMBbDYb5qh/SPi60owsFNEFxIiBTAlV0O9sNhI2mw3q\nACOLPvvu6H8WrDdLjSykCWOjYxgT6cJSgUxEJDMGMhGRJBjIFBmOK4h0w0AWQRFdABHJiIFMRCQJ\nywQyz7CIA44riHRlmUAmIpIdA5mISBIM5ERTRBdARLJiIBMRSYKBTNrwDT0i3TGQiYgkwUAmIpIE\nA5mISBIMZCIiSVgikHmVHhEZgSUCmWLEMyyIEoKBTEQkCV0Dedy4cbjvvvvgdruRn5+v51LGoIgu\ngIiiUVFRAbvdjpycz7+SWlEUpKenw+12w+12Y+/eveHHNm/ejMzMTGRnZ6OpqUnzOrp+67TNZoPP\n58Po0aP1XIaISFeLFy/G8uXLsWjRovA2m82GlStXYuXKlT32bW9vx9atW3Hw4EGcP38eVVVVOH78\nuKZ1dA1kQMw3txIRxdPUqVNx4cKFPtv7yze/34+SkhI4HA44HA6oqopgMIjU1NQh19F1ZGGz2VBc\nXIzy8nLs2bNHz6WIiBJuy5YtKCgowMaNGxEMBgEAgUAALpcrvE9WVhYCgYCm4+naIR8+fBhpaWlo\nbW1FaWkp8vPzMWbMmPDjZ5T/Cf/36KKJuLtokp7lEJFhHAOg7df8ePB9BPjaI/sZr9eLp59+Gleu\nXEF1dTVqa2uxevXqfrtmm82m6Zi6BnJaWhoAwOVyYf78+XjzzTfx6KOPhh/PVL6j5/JEZFj3f3YL\neSU+h/3v/jcXfXYLWTdl6EN95StfAQCMHDkSjz32GJYtW4bVq1fD4/GgoaEhvN/p06eRl5enqTzd\nRhbXr18Pt/Aff/wx6uvrUVJSotdyREQJdenSJQBAV1cXtm/fjrlz5wIA8vPzUV9fj7a2Nvh8PiQl\nJWmaHwM6dsgfffQRvvnNbwIA7r77bqxatQoZGRl6LUdEpJuFCxfi0KFD6OjoQEZGBtatWwefz4eT\nJ08iOTkZhYWF8Hq9AAC73Q6v14vi4mIkJyejtrZW8zo2VdBpEDabDXPUPyRkLakunVZEFxAFXqlH\nwhXEfMaWzWaDelLjvlPEnCHGK/WIiCTBQKahTfeIroDIEhjIRESSsEQgzyl8XXQJRERDskQgExEZ\nAQOZiEgSDGQiIkkwkImIJMFAJm146huR7hjIRESSYCAnmiK6ACKSFQOZiEgSDGTSjnNkIl1ZJpB5\ntR4Ryc4ygUxEJDsGsgiK6AKISEYMZIoM58hEumEgExFJwlKBzDf2iEhmlgpkqSiiC4gBxxZkMRUV\nFbDb7cjJyQlvq66uhsvlQm5uLh5//HHcuHEj/NjmzZuRmZmJ7OxsNDU1aV6HgUxENITFixdj3759\nPbbNmjULLS0tOHr0KK5du4bt27cDANrb27F161YcPHgQNTU1qKqq0rwOA5miwy6ZLGTq1KkYNWpU\nj20zZ85EUlISkpKSMHv2bBw6dAgA4Pf7UVJSAofDgWnTpkFVVQSDQU3rMJCJiGL0i1/8AqWlpQCA\nQCAAl8sVfiwrKwuBQEDTcYQG8lLUJnxNqd7YU0QXECN2yWQSvmZAqfn8Fon169cjNTUV3/72twEA\nqqr22cdms2k61rDIliYiMq49k2f1/8BkIPf/dLtfu1/T8X7961+jvr4eBw8eDG/zeDxoaGgI3z99\n+jTy8vI0HY8jC4oNu2SyqH379uH555/Hnj17MHz48PD2/Px81NfXo62tDT6fD0lJSUhNTdV0THbI\noikw/uiCyOQWLlyIQ4cOoaOjAxkZGVi3bh02bNiATz75BDNmzAAAfPWrX8XWrVtht9vh9XpRXFyM\n5ORk1NZqH83a1P4GHglgs9mwW73968M2LEn4+nsbH0z4mgNSRBcQB2/5RVdAplbQ72w2Et0zZyhl\ntv0xrxcNjiyIiCTBQJaBIrqAOOAsmShmlg1kqU5/IyKCJIEs4nxk0gG7ZKKYSBHIBHOMLQCGMlEM\nGMhERJKwdCBLN0dWRBcQJ+ySiaIiTSBzjmwyDGWiiEkTyPQZRXQBRCSK5QNZurGFmbBLJoqIVIHM\nsYUJMZSJNJMqkEWRrktWRBcQZwxlIk0YyJQYDGWiIUkXyBxbfEYRXQARJZp0gSyKdGMLwHyhzC6Z\naFAMZEoshjLRgKQMZFFjC3bJCcJQJuqXlIFMFsBQJuqDgdwLu+QEYigT9aBbIDc2NsLlciEzMxNb\ntmyJ+Od5toVFMJTJILZv345p06Zh4sSJ+OUvfwkACAaDKCsrg8PhQHl5Oa5evRrTGroF8ooVK1Bb\nW4uGhga8/PLL6Ojo0Gspa1BEF6AjhjJJ7vLly1i3bh3eeOMN+P1+/PznP8fly5dRU1MDh8OBM2fO\nID09Hdu2bYtpHV0C+fLlywCAwsJCjB07FrNmzYLfH/m3EvPNvV4U0QXoiKFMEnv77beRm5uLUaNG\nYcSIEZg+fTr+8pe/IBAIoLKyEikpKaioqIgq57rTJZCbm5vhdDrD97Ozs3HkyBE9liIzme5hMJOU\nCgsLEQgEcP78eVy6dAl/+tOf8Pbbb/fIOqfTiUAgENM6w+JRbLR2KGfD/z2paDRyikYLrKanOYWv\nY2/jg6LL6EuBuTtl4HYovxVbp0FGdwzA8dv/OS4duKDvau/6/oX3fP8a8PG77roLL774Ih577DFc\nvnwZOTk5SElJgaqqca1Dl0DOy8tDdXV1+H5LSwtKSkr67LdQuXfIYy1FLbZhSVzrMzwFDGUyuftv\n30K/MV1YF5ejDpglRZ/dQtYt6LNLaWkpSktLAQAPPfQQSkpKcPz4cbS2tsLtdqO1tRV5eXkx1afL\nyGLkyJEAbp9pceHCBRw4cAAej/F+FZV2lmwVHF9Yl4Tjq/b2dgBAQ0MD3nvvPeTm5sLj8aCurg43\nbtxAXV0dCgoKYlpDt7MsXnzxRSxZsgQzZszAsmXLcM8990R9LJGnwEkbyoroAhJEwr+YpDNJX+9v\nfetbcDqdePLJJ1FXVwebzQav14u2tjZkZWXh4sWLWLp0aUxr2NR4D0G0LmyzYbc6S/P+IscWUs6S\nQxTRBSQQRxjmN1AYv2WLeV5rs9kwR/2Dpn332hbEfT6shWGu1GOXTOyWTYyvLQADBTINQBFdgAD8\ni2sufD3DDBXI7JIHoIguQAB2VMbH17APQwUyDUIRXYAg/AttPAziARkukNklD0IRXYAg/AtuDHyd\nhmS4QBZN+lC2Mv6FlxdfF00MGcj8aM5BKKILkACDWR58LSJiyEAWTfouWRFdgCQYBuLwuY+KYQOZ\nXfIQFNEFSIThkBih55nPddQMG8iiSd8lAwzl3hgW+uDzGjeGDmTRXTJD2aDYycUHn8O4E/p5yJQg\nChjMA+keKPysjKExgHVl6A4ZYJesmSK6AANg59w/Pi8JY/hABhjKmimiCzAQq4eQ1f//BeHIIk6k\n/cqn3hQwmCPVO5TMONpg8ErBNIHMr3qKgAKGcizMENAMYCmZJpBlYJguGWAox9OAH6ouQVAzeA3F\nMN8YopUMXbJhQhlgKIsWS2hbKWwt8o0hpgtkgKEcMUV0AURDsEggm+IsC4qRIroAIvldu3YNjzzy\nCCZMmIDs7Gz4/X4Eg0GUlZXB4XCgvLwcV69ejWkNUway6NPgAAOdCheiiC6ASG4/+clP4HA4cOrU\nKZw6dQpOpxM1NTVwOBw4c+YM0tPTsW3btpjWMGUgAwzlqCiiCyCSV0NDA3784x9j+PDhGDZsGEaO\nHIlAIIDKykqkpKSgoqICfn9sb+SaNpBlYchQVgTXQCSZf/7zn7h58ya8Xi88Hg82btyIGzduoLm5\nGU6nEwDgdDoRCARiWsfUgSxDlwwYMJQBhjJZSqfvPZxR/id86+3mzZt4//33sWDBAvh8PrS0tOD3\nv/993N/4M+VZFr3JcNYFYLAzL0IU0QUQIW5nWeCQxmNM67uey+VCa2srAGDv3r34zW9+g08++QRr\n166F2+3GsWPHsGHDBuzcuTPqGk3dIVMcKKILIJJDZmYm/H4/bt26hT/+8Y+YMWMGPB4P6urqcOPG\nDdTV1aGgoCCmNSwRyBxdxEgBg5ks74UXXsCKFSuQm5uL4cOH46GHHoLX60VbWxuysrJw8eJFLF26\nNKY1LDGyCOHoIg4U0QWQJUkwskgES3TIsjFspwwwkIl0ZKlAlmV0AZgglBXBNRCZkNBAnv/O/oSv\nyVCOI0V0AUTmIrxDZiibIJQVwTUQmYTwQBaFoRxniugCiIxPikAW0SXLxjShrAiugcjApAhkgKML\nwCShDDCYiaIkTSCLwlDWkSK6ACJjkSqQRY0uGMo6UsBgJtJIqkAGGMohpgplgMFMpIF0gQwwlENM\nF8oAg5loEFIGMn1uTuHrDGYii5A2kNkl92TKUAYYzETdSBvIAEO5N9OGMsBgJoLkgQwwlHszdSgD\nDGbqSYGl/jxIH8giMZQFUmC5v4zUjQJLvvZCP6BePal9/z2TE/th9t3J8sH2/TH0h91HShFdAOlK\nGeSxOHxgvBE+oN4wgQwwlAdiqVAOUUQXQHGjaNiHgazzwlEEMsBQHoglQzlEEV0ARUyJcH+LBLIu\nM2RFUZCeng632w232419+/bF7dgiPxlO1pkyYOLzlbVQwFA2AgWGfa1u3rwJj8eDKVOmoKCgAJs2\nbQIABINBlJWVweFwoLy8HFevXo1pHV065HXr1iE1NRUrV64ceOEoO2RAbJcMyN0pAxbvlrtTRBdA\ncXsNJOiQr1+/jjvvvBP//ve/cf/992PXrl3YtWsXPvjgA7zwwgtYtWoVxo0bh9WrV0ddo25nWejZ\n7ov+/GSZO2XAImdhaKHAsB2ZoSkw5fN+5513AgCuXr2Krq4upKSkIBAIoLKyEikpKaioqIDf749p\nDd0CecuWLSgoKMDGjRsRDAbjfnyG8uAsPcLojwJThoQUFFji+b116xYmT54Mu92OH/zgB3A4HGhu\nbobT6QQAOJ1OBAKBmNaIemQxc+ZMfPjhh322P/PMMygoKMCXv/xlXLlyBdXV1ZgwYUKfNt5ms+En\n3X7zL/rfQFFe5HVwfDE0jjCGoIguwGCUBKxxwgec9H1+/9fr4jOymD7AMf6vD/h/3da7MPB6Fy5c\nwNy5c/G73/0OZWVleP/99zF8+HBcv34dLpcL//jHP6KvUe+zLN555x0sW7YMhw8f7rlwDDPk3hjK\n2jCYI6CILkAiiugCEL8Z8kCB3Ntbg6+3evVq3HvvvThw4ADWrl0Lt9uNY8eOYcOGDdi5c2fUNQ6L\n+icHcenSJaSlpaGrqwvbt2/H3Llz9VgmbP47+4WG8lLUGiKU5xS+zlDWStG4zUwU0QXIq6OjA8OG\nDcOXvvQldHZ2Yv/+/Vi1ahWuXLmCuro6PPfcc6irq0NBQUFM6+jSIS9atAgnT55EcnIyCgsLsXbt\nWowePbrnwnHskEPYKWvHYI4zRXQBEVBEFxAFwR3yu+++i0ceeQT/+c9/MGbMGHzve9/DokWLEAwG\n8f3vfx8nTpxAbm4uXn31VYwYMSL6Go12YchQRIcyYJxgZigLoBj02KJJNrLQi+kCGWAoR4rBTNKz\nSCCb8tPeRJ8SB8h/Wlx3PEWOSA6mDGSAoRwNhjKRWKYNZIChHA12y0TimDqQAYZytBjMRIlnyjf1\n+iPDG32Asd7s645v/JFQfFPPXGTolAFjdssAO2aiRLBMIAMM5XhgMBPpR+zIYiGA/0782rKMLwDj\njjBCOMqghODIIkE2Jn5JWTplwNjdMsCOmSiexAeyIAzl+GIwE8VO/MgiRMDoApBrfAEYf4QRwlEG\nxRVHFgkmYHQByNUpA+bolgF2zETRkKdDDmGnHGaWbjmEXTNFjR2yIOyUw8zSLYewayYanHwdcgg7\n5R7M1i2HsGsmTSzSIcsbyICwUAbkDGazhnIIw5kGxEDWeWEtgQwwlPvBYCYrmVP4OvbaFjCQdV1Y\nayADDOV+mD2UQxjO1tT7vQYZArmxsRFLlixBV1cXqqqqsHz58pjq6bdGQwQywFAegFWCOYQBbV6D\nveErQyC73W689NJLGDt2LGbPno2mpibcc889MdXUp0bDBDLAUB6E1YIZYDibgdazbkQH8uXLl1FU\nVIQTJ04AAKqqqjB79mzMmzcvppp6GxbXo+ltI4SFcui0OFmDeSlqLRfKfX6tZUBLz6inPTY3N8Pp\ndIbvZ2dn48iRIxYPZEBoKAO3g1nmUAas2S0DDGgZGTWARTHWyKI3jjAGZdVgHggDWn96BXDcRhY4\nMsCjxwAc73b/lUFHFsuXL0dJSUncO2RjBzIgNJQBBrPRMaSjl8juV/9A7q1gwDf1HA4HSkpK+Kbe\ngBjKQ2Ioa8eQ7kv06EGGQD506BCWLl2KTz/9FFVVVaiqqoqpnn5rNEUgAwxljRjMsTFzWIsO3cHI\nEMiJYJ5ABoSHMsBgtjpZA1vmsNWCgaz3wnoEMsBQjhCDmYyAgaz3wnoFMiBFKAMMZqJ4sUogy/d5\nyPGwEcI+V7k7GT9jeSBLUWu6z18mMhpzdsjdsVuOGrtmkgU7ZLOQoFMGjNUth7BrJkos83fIIZJ0\nyoAxu2WAHTOJww7ZbCSZKwPG7JYBdsxEehPaITcB+K9Edskh7Jbjhl0zJQI75AQ5vEPAopJ0yoBx\nu+UQds1E8SO8Qw6xeqcMGL9bDmHXTPFmlQ5ZmkAGBIUywGDWEcOZ4oGBrPfC/QQywFDuzkzBDDCc\nKXoMZL0XHiCQAYZyd2YL5RCGM0WCgaz3woMEcgiD+XNmDWaA4UxDYyDrvbCGQAYYyr2ZOZgBhjP1\nj4Gs98IaAxkQGMoAg1kwBrR1dT+dssy2n4Gs68IRBDLAUB6IVYIZYDib3WDnszOQ9V44wkAGBIcy\nwGCWDAPa2CK5oMgqgRz1lXqvvfYaJk6ciDvuuAPHjx/v8djmzZuRmZmJ7OxsNDVFGrsDO7xD0JV9\nIRJd4dfb/Hf2G/6qv0iFrhLsfiN5mfW1euqppzB58mRMmTIFDz/8MDo7O8OPRZqFUXfIp0+fRlJS\nEpYsWYKf/exnyM3NBQC0t7ejsLAQ+/fvx/nz5/HEE0/0CWwgug65OzN1y75moCgvfscDjNsxv+v7\nF3KKRsfteEbuojt97+Huokmiy4hKvANX5g45GAwiNTUVALB+/Xp0dXVh/fr1mrOwu2Eaq+vD6XT2\nu93v96OkpAQOhwMOhwOqqvYoOF4O7xAcyhsRt1D2HY1/IIe6ZaMF83txDuSBgsEIQf0vX4shAtlM\n3W40QtnW1dWFa9euYeTIkQCiy8KoA3kggUAALpcrfD8rKwuBQAAPPPBAvJcKjy+EBXNohCHxbNmo\nwaw3Iwe1KFYP3sGsWbMGtbW1yMrKgs/nAxBdFg4ayDNnzsSHH37YZ/uzzz6L0tLSfn+mvzb/9q8K\n+jFTt6wXBrM2DGoGb3+GysJnnnkGa9aswZo1a/DDH/4QmzZtii4L1RgVFRWpx44dC9/fs2ePWlVV\nFb4/efJk9cqVK31+bvz48SoA3njjjbchb+PHj481qiJab8SIEVGtcerUKdXj8aiqqj0Lu4vLyELt\n9i9Bfn4+qqur0dbWhnPnziEpKanfmcnZs2fjsTQRkSaqTqexnTlzBpmZmejq6sKOHTvw4IMPAtCe\nhd1FHci7du1CVVUVOjo6MG/ePLjdbuzduxd2ux1erxfFxcVITk5GbS1//SEi83ryySfxt7/9DV/8\n4hdRVFSERx99FACiykJhF4YQEVFPCf8KJxEXlBiBoihIT0+H2+2G2+3Gvn37RJeUcI2NjXC5XMjM\nzMSWLVtElyPUuHHjcN9998HtdiM/P190OQlVUVEBu92OnJyc8LZgMIiysjI4HA6Ul5fj6tWrAivU\nT8IDOScnB7t27UJhYWGP7e3t7di6dSsOHjyImpoaVFVVJbo0oWw2G1auXIkTJ07gxIkTKCkpEV1S\nwq1YsQK1tbVoaGjAyy+/jI6ODtElCWOz2eDz+XDixAkEAgHR5STU4sWL+zQkNTU1cDgcOHPmDNLT\n07Ft2zZB1ekr4YHsdDoxYcKEPtu7n0Q9bdq08EnUVmLl6dHly5cBAIWFhRg7dixmzZoFv98vuCqx\nrPrnYerUqRg1alSPbYFAAJWVlUhJSUFFRYVp/2wI/9bpkIFOoraSLVu2oKCgABs3brTcP0bNzc09\nrv7Mzs7GkSNaL3M1H5vNhuLiYpSXl2PPnj2iyxGu+58Pp9Np2myI+5V6gHEuKEm0gZ6XZ555Bl6v\nF08//TSuXLmC6upq1NbWYvXq1QKqJBkcPnwYaWlpaG1tRWlpKfLz8zFmzBjRZQljld8WdAnkAwcO\nRPwzHo8HDQ0N4funT59GXl6cP+BBMC3Py8iRI/HYY49h2bJllgrkvLw8VFdXh++3tLRYco4ekpaW\nBgBwuVyYP38+3nzzzfDpVFaUl5eH1tZWuN1utLa2mi4bQoSOLHpfUFJfX4+2tjb4fD5NJ1GbyaVL\nlwDc/oCS7du3Y+7cuYIrSqzQB7I0NjbiwoULOHDgADwej+CqxLh+/Xp4ZPXxxx+jvr7e0v84Abcb\ntrq6Oty4cQN1dXUoKCgQXZI+oro+MAavv/66mp6erg4fPly12+1qSUlJ+LEXX3xRHT9+vOpyudTG\nxsZElybUww8/rObk5Kj333+/+sQTT6idnZ2iS0o4n8+nOp1Odfz48epLL70kuhxhzp07p06ePFmd\nPHmyWlxcrL7yyiuiS0qohx56SE1LS1OTk5PV9PR0ta6uTr1y5Yo6f/58NSMjQy0rK1ODwaDoMnXB\nC0OIiCQhzVkWRERWx0AmIpIEA5mISBIMZCIiSTCQiYgkwUAmIpIEA5mISBIMZCIiSfx/MpMmbrDt\nBFYAAAAASUVORK5CYII=\n",
+ "text": [
+ ""
+ ]
+ }
+ ],
+ "prompt_number": 5
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "opt = inverse.ProjectedGradient(maxIter=50,maxStep=np.inf, maxIterLS=20, debug=False)\n",
+ "opt.remember('f')\n",
+ "opt.lower = -2\n",
+ "opt.upper = 2\n",
+ "opt.minimize(FUN,np.array([0,0]))"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "======================= Projected Gradient =======================\n",
+ " # f |proj(x-g)-x| LS itType aSet bSet Comment\n",
+ "------------------------------------------------------------------\n",
+ " 0 0.00e+00 2.83e+00 0 SD 0 0 \n",
+ " 1 -1.60e+01 0.00e+00 0 SD 2 2 \n",
+ "------------------------- STOP! -------------------------\n",
+ "0 : |fc-fOld| = 1.6000e+01 <= tolF*(1+|f0|) = 1.0000e-01\n",
+ "0 : |xc-x_last| = 2.8284e+00 <= tolX*(1+|x0|) = 1.0000e-01\n",
+ "1 : |proj(x-g)-x| = 0.0000e+00 <= tolG = 1.0000e-01\n",
+ "1 : |proj(x-g)-x| = 0.0000e+00 <= 1e3*eps = 1.0000e-02\n",
+ "0 : maxIter = 50 <= iter = 1\n",
+ "1 : probSize = 2 <= bindingSet = 2\n",
+ "------------------------- DONE! -------------------------\n"
+ ]
+ },
+ {
+ "metadata": {},
+ "output_type": "pyout",
+ "prompt_number": 7,
+ "text": [
+ "array([ 2., 2.])"
+ ]
+ }
+ ],
+ "prompt_number": 7
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "opt = inverse.ProjectedGradient(maxIter=50, maxStep=0.3,debug=False, maxIterLS=20, tolCG=1e-5, maxIterCG=1000)\n",
+ "opt.lower, opt.upper = -0.4, 0.9\n",
+ "opt.remember('f', 'xc', ('norm_g', lambda M: np.linalg.norm(M.g)), 'phi_d', 'phi_m')\n",
+ "inv = inverse.Inversion(prob,reg,opt,beta0=1e-3,debug=False)\n",
+ "inv.remember(('phi_d',lambda I:I.opt.recall('phi_d')),('phi_m',lambda I:I.opt.recall('phi_m')))\n",
+ "m0 = np.zeros_like(m_true)\n",
+ "mrecB = inv.run(m0)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "====================================== Projected Gradient ======================================\n",
+ " # beta phi_d phi_m f |proj(x-g)-x| LS itType aSet bSet Comment\n",
+ "------------------------------------------------------------------------------------------------\n",
+ " 0 1.00e-03 9.62e+02 0.00e+00 9.62e+02 2.14e+01 0 SD 0 0 \n",
+ " 1 1.00e-03 9.53e+02 8.49e-03 9.53e+02 2.18e+01 9 SD 0 0 \n",
+ " 2 1.00e-03 5.32e+02 1.42e+03 5.33e+02 2.20e+01 0 .CG. 0 0 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 3 1.00e-03 2.32e+02 5.66e+03 2.38e+02 2.26e+01 0 .CG. 0 0 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 4 1.00e-03 1.29e+02 8.85e+03 1.38e+02 2.30e+01 1 .CG. 0 0 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 5 1.00e-03 1.10e+02 9.74e+03 1.20e+02 2.15e+01 3 .CG. 22 22 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 6 1.00e-03 1.09e+02 9.74e+03 1.19e+02 2.36e+01 10 SD 22 0 Stop SD\n",
+ " 7 1.00e-03 7.65e+01 1.33e+04 8.98e+01 2.21e+01 1 .CG. 85 75 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 8 1.00e-03 5.11e+01 1.33e+04 6.44e+01 2.38e+01 9 SD 75 3 \n",
+ " 9 1.00e-03 5.07e+01 1.33e+04 6.41e+01 2.11e+01 12 SD 4 4 \n",
+ " 10 1.00e-03 5.06e+01 1.33e+04 6.40e+01 2.20e+01 12 SD 17 9 Stop SD\n",
+ " 11 1.00e-03 5.04e+01 1.34e+04 6.37e+01 2.05e+01 8 .CG. 85 81 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 12 1.00e-03 5.03e+01 1.34e+04 6.37e+01 1.65e+01 13 SD 81 26 \n",
+ " 13 1.00e-03 5.03e+01 1.34e+04 6.37e+01 1.67e+01 15 SD 30 30 \n",
+ " 14 1.00e-03 5.03e+01 1.34e+04 6.36e+01 1.42e+01 15 SD 84 84 Stop SD\n",
+ " 15 1.00e-03 2.65e+01 2.26e+04 4.91e+01 2.52e+01 0 .CG. 190 50 "
+ ]
+ },
+ {
+ "output_type": "stream",
+ "stream": "stdout",
+ "text": [
+ "\n",
+ " 16 1.00e-03 1.51e+01 2.26e+04 3.77e+01 2.24e+01 10 SD 50 1 \n",
+ "------------------------- STOP! -------------------------\n",
+ "1 : |fc-fOld| = 1.1392e+01 <= tolF*(1+|f0|) = 9.6325e+01\n",
+ "1 : |xc-x_last| = 3.4809e-03 <= tolX*(1+|x0|) = 1.0000e-01\n",
+ "0 : |proj(x-g)-x| = 2.2369e+01 <= tolG = 1.0000e-01\n",
+ "0 : |proj(x-g)-x| = 2.2369e+01 <= 1e3*eps = 1.0000e-02\n",
+ "0 : maxIter = 50 <= iter = 16\n",
+ "0 : probSize = 1000 <= bindingSet = 1\n",
+ "1 : phi_d = 1.5077e+01 <= phi_d_target = 2.0000e+01 \n",
+ "------------------------- DONE! -------------------------\n",
+ "-------------------------STOP!-------------------------\n",
+ "0 : maxIter = 10 <= iter = 1\n",
+ "1 : phi_d = 1.5077e+01 <= phi_d_target = 2.0000e+01 \n",
+ "-------------------------DONE!-------------------------\n"
+ ]
+ }
+ ],
+ "prompt_number": 8
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "# First set up the figure, the axis, and the plot element we want to animate\n",
+ "fig = plt.figure()\n",
+ "ax = plt.axes()\n",
+ "ax.plot(M.vectorCCx, m_true, 'b-')\n",
+ "txt = plt.text(0.8,0.9,'')\n",
+ "line, = ax.plot([], [], 'r-', lw=1)\n",
+ "\n",
+ "def animate(i):\n",
+ " txt.set_text('iteration %d'%i)\n",
+ " line.set_data(M.vectorCCx, np.array(opt.recall('xc')[i]))\n",
+ "\n",
+ "SimPEG.utils.animate(fig, animate, frames=len(opt.recall('xc')))"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [
+ {
+ "html": [
+ ""
+ ],
+ "metadata": {},
+ "output_type": "pyout",
+ "prompt_number": 9,
+ "text": [
+ ""
+ ]
+ }
+ ],
+ "prompt_number": 9
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "plot(opt.recall('phi_d'))"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": [
+ {
+ "metadata": {},
+ "output_type": "pyout",
+ "prompt_number": 10,
+ "text": [
+ "[]"
+ ]
+ },
+ {
+ "metadata": {},
+ "output_type": "display_data",
+ "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHg5JREFUeJzt3XtwVGWCxuFfIlcFlAIJOqEhQkgnoCRAp4MsJDBymXFj\ngiyFzHqpSdjSOEwEYcpdpUqc2nKKi8WANQSWmTjlWmKxO+OKyDXMNBgv3cEwKBDkbkCFAEHoQIIY\nzv5xTHMLkDSdnO4+71PVReju0/1Gi/d0f+f7zokxDMNARERsI9bqACIi0rpU/CIiNqPiFxGxGRW/\niIjNqPhFRGxGxS8iYjM3LP68vDzi4uK4//77A/f5/X5ycnJwOBzk5uZSU1MTeGzx4sUkJiaSkpJC\naWlp4P6KigoGDx7Mfffdx0svvdQCv4aIiDTVDYv/l7/8JevWrbvivqKiIhwOB3v37iU+Pp6lS5cC\nUFVVxZIlS9i0aRNFRUUUFhYGtpk5cyYvvPACZWVlbN68ma1bt7bAryIiIk1xw+IfMWIEXbt2veI+\nn89Hfn4+7du3Jy8vD6/XC4DX62X8+PE4HA4yMzMxDCPwbeDLL79k8uTJdOvWjUcffTSwjYiItL5m\nj/GXlZXhdDoBcDqd+Hw+wCz+5OTkwPOSkpLwer3s27ePHj16BO5PSUnh008/vdXcIiISpGYXf3PO\n8BATE3NL24uISOi1ae4GLpeLiooK0tLSqKiowOVyAeB2uykpKQk8b/fu3bhcLjp37syxY8cC9+/a\ntYuMjIxGX7tfv37s37+/uZFERGytb9++7Nu3r8nPb/YnfrfbTXFxMbW1tRQXFwdKPD09nfXr11NZ\nWYnH4yE2NpbOnTsD5pDQO++8w4kTJ3j33Xdxu92Nvvb+/fsxDCPsby+//LLlGaIlZyRkVE7lDPdb\ncz8w37D4p0yZwoMPPsiePXvo1asXb7zxBgUFBVRWVpKUlMTXX3/NM888A0BcXBwFBQWMHj2aZ599\nlkWLFgVeZ8GCBcybNw+Xy8WIESMYOnRos0KKiEjo3HCoZ8WKFY3e/9577zV6/3PPPcdzzz13zf0p\nKSmUl5cHEU9EREJNK3eDkJWVZXWEJomEnJGQEZQz1JTTWjGGYYTNNJuYmBjCKI6ISERobnfqE7+I\niM2o+EVEbEbFLyJiMyp+ERGbUfGLiNiMil9ExGZU/CIiNhN2xX/ihNUJRESiW9gVv07VLyLSssKu\n+D/5xOoEIiLRTcUvImIzYXeuns6dDaqroU2zLxEjImJPEX+unp/8BHbssDqFiEj0CrviHzZMwz0i\nIi1JxS8iYjMqfhERmwm74k9JgePHzZuIiIRe2BV/bCykp2shl4hISwm74gcN94iItCQVv4iIzYTd\nAi7DMDh1ChwOOHVKC7lERG4m4hdwAXTtCr16wRdfWJ1ERCT6hGXxg4Z7RERaiopfRMRmVPwiIjYT\ntsWfnGxejauqyuokIiLRJWyLPzYW3G4t5BIRCbWwLX7QcI+ISEtQ8YuI2ExYLuBq8N135nx+LeQS\nEbm+qFjA1eCuu8wVvJ9/bnUSEZHoEdbFDxruEREJNRW/iIjNqPhFRGwm6OJfvnw5Dz74IEOGDGH6\n9OkA+P1+cnJycDgc5ObmUlNTE3j+4sWLSUxMJCUlhdLS0ia/j9MJ1dVayCUiEipBFX91dTWvvvoq\nGzdupKysjD179rB+/XqKiopwOBzs3buX+Ph4li5dCkBVVRVLlixh06ZNFBUVUVhY2PSAPy7k0qd+\nEZHQCKr4O3bsiGEYnD59mtraWs6dO8ddd92Fz+cjPz+f9u3bk5eXh9frBcDr9TJ+/HgcDgeZmZkY\nhoHf72/y+2m4R0QkdIIu/qKiIvr06UPPnj0ZPnw4brebsrIynE4nAE6nE5/PB5jFn5ycHNg+KSkp\n8FhTqPhFREInqGVRx48fp6CggF27dtG1a1cmTZrE6tWrm7WAICYmptH758yZE/g5KyuLrKws3G74\n7DO4cAHatg0msYhI9PB4PHg8nqC3D6r4fT4fGRkZ9OvXD4BJkybx4Ycf4nK5qKioIC0tjYqKClwu\nFwBut5uSkpLA9rt37w48drXLi7/BnXdCnz7mQq4hQ4JJLCISPRo+FDd45ZVXmrV9UEM9I0aMYOvW\nrVRXV3P+/HnWrl3L2LFjcbvdFBcXU1tbS3FxMRkZGQCkp6ezfv16Kisr8Xg8xMbG0rlz52a9p4Z7\nRERCI6hP/F26dGH27NlMmDCBc+fOMX78eEaNGkV6ejqPP/44SUlJDB48mLlz5wIQFxdHQUEBo0eP\npl27dixbtqzZ7zlsGGzaBNOmBZNYREQahPVJ2i5XUQEPPwwHDrRyKBGRMBdVJ2m7XFKSebbOY8es\nTiIiEtkipvi1kEtEJDQipvhBB3hFREJBxS8iYjMRc3AX4MwZuPde84pcWsglImKK2oO7AF26QEIC\nbN9udRIRkcgVUcUPGu4REblVKn4REZtR8YuI2EzEFX///nD6NBw9anUSEZHIFHHFHxsLGRn61C8i\nEqyIK37QcI+IyK1Q8YuI2ExELeBq0LCQq7oa2rVrhWAiImEsqhdwNejSBe67Twu5RESCEZHFDxru\nEREJlopfRMRmVPwiIjYTscXfvz/4/fDtt1YnERGJLBFb/DExWsglIhKMiC1+0HCPiEgwVPwiIjYT\nkQu4Gvj9cM89WsglIvZmiwVcDTp3hr594R//sDqJiEjkiOjiBw33iIg0l4pfRMRmVPwiIjYT8cWf\nmAhnz8I331idREQkMkR88Wshl4hI80R88YOGe0REmkPFLyJiMxG9gKtBTQ3ExcGpU1rIJSL2Y6sF\nXA06dTIP8m7bZnUSEZHwFxXFDxruERFpKhW/iIjNBF38Z8+e5amnnqJ///6kpKTg9Xrx+/3k5OTg\ncDjIzc2lpqYm8PzFixeTmJhISkoKpaWlIQl/ORW/iEjTBF38L7/8Mg6Hg88//5zPP/8cp9NJUVER\nDoeDvXv3Eh8fz9KlSwGoqqpiyZIlbNq0iaKiIgoLC0P2CzTo1w9qa+Hrr0P+0iIiUSXo4i8pKeHF\nF1+kQ4cOtGnThjvvvBOfz0d+fj7t27cnLy8Pr9cLgNfrZfz48TgcDjIzMzEMA7/fH7JfArSQS0Sk\nqYIq/iNHjlBXV0dBQQFut5u5c+dSW1tLWVkZTqcTAKfTic/nA8ziT05ODmyflJQUeCyUNNwjInJz\nQRV/XV0de/bsYeLEiXg8Hnbu3MnKlSubNY80JiYmmLe+IRW/iMjNtQlmo379+pGUlER2djYAU6ZM\n4c0338TlclFRUUFaWhoVFRW4XC4A3G43JSUlge13794deOxqc+bMCfyclZVFVlZWk3O5XLB9O5w/\nD+3bN//3EhGJBB6PB4/HE/T2Qa/cfeSRR3jppZdwuVwUFhaSlpbGyZMnOXz4MPPmzWPWrFkkJCQw\na9Ysjh07RmZmJhs2bODAgQM8//zzlJeXXxsmyJW7l0tLg6Iic7xfRMQOmtudQX3iB1iwYAFPPvkk\ndXV1PPTQQzz22GNcvHiRxx9/nKSkJAYPHszcuXMBiIuLo6CggNGjR9OuXTuWLVsW7NveVMNwj4pf\nRKRxUXGunsv993/D++/DypUhCiUiEuZsea6ey+kAr4jIjUVd8fftC3V1cOSI1UlERMJT1BV/TIw+\n9YuI3EjUFT+o+EVEbkTFLyJiM1E3qwfg7Fno0QOqq7WQS0Sin+1n9QDccQckJUEja8RERGwvKosf\nNNwjInI9Kn4REZtR8YuI2EzUFv9998H338Phw1YnEREJL1Fb/DExkJkJGzdanUREJLxEbfED5OTA\ne+9ZnUJEJLxE5Tz+BqdOQe/e8O235hRPEZFopHn8l+naFdLTYf16q5OIiISPqC5+gAkT4P/+z+oU\nIiLhI6qHesCc1ZOaCkePQtu2IX1pEZGwoKGeq/TqZU7t3LLF6iQiIuEh6osfNNwjInK5qB/qAdi1\nC8aNg8pKc36/iEg00VBPI5KToWNHna1TRARsUvwxMZCbC+++a3USERHr2aL4QeP8IiINbFP8bjec\nOAF791qdRETEWrYp/thYnbtHRARsVPygcX4REbDJdM4G589DXBzs3g09e7bY24iItCpN57yB9u1h\n/Hh4/32rk4iIWMdWxQ/mcI9m94iIndlqqAfgzBmIj4cjR6BLlxZ9KxGRVqGhnpvo0gWGD4d166xO\nIiJiDdsVP2i4R0TszXZDPWBeijElBY4dg3btWvztRERalIZ6muCee8wTt/3971YnERFpfbYsftBw\nj4jYly2HegD27IGsLHN2T6xtd38iEg1adainvr6etLQ0srOzAfD7/eTk5OBwOMjNzaWmpibw3MWL\nF5OYmEhKSgqlpaW38rYh0b8/3HUX+HxWJxERaV23VPyLFi0iJSWFmB8va1VUVITD4WDv3r3Ex8ez\ndOlSAKqqqliyZAmbNm2iqKiIwsLCW08eAjpVs4jYUdDFf+TIEdasWcPUqVMDXzF8Ph/5+fm0b9+e\nvLw8vF4vAF6vl/Hjx+NwOMjMzMQwDPx+f2h+g1ugcX4RsaOgi3/GjBnMnz+f2MsGyMvKynA6nQA4\nnU58P46jeL1ekpOTA89LSkoKPGalIUOgpsY8aZuIiF0EVfyrV6+mR48epKWlXXFAoTkHF2LC4Krn\nsbE6VbOI2E+bYDb6+OOPWbVqFWvWrKGuro4zZ87wxBNP4HK5qKioIC0tjYqKClwuFwBut5uSkpLA\n9rt37w48drU5c+YEfs7KyiIrKyuYiE2WmwsvvQT/8R8t+jYiIiHj8XjweDxBb3/L0zk3b97MggUL\neP/995k3bx6HDx9m3rx5zJo1i4SEBGbNmsWxY8fIzMxkw4YNHDhwgOeff57y8vJrw7TidM4GFy6Y\n5+j/4gv4yU9a9a1FRELCkpW7DcM2BQUFVFZWkpSUxNdff80zzzwDQFxcHAUFBYwePZpnn32WRYsW\nheJtQ6JtW/j5z2HVKquTiIi0Dtsu4LrcX/4Cy5bBhg2t/tYiIresud2p4sec2XPvvVBZaS7qEhGJ\nJDpJWxA6dYLMTFizxuokIiItT8X/Iy3mEhG70FDPj44fh8REOHoUOnSwJIKISFA01BOku++GBx6A\nTZusTiIi0rJU/JfRcI+I2IGGei5z4AAMGwbffAO33WZZDBGRZtFQzy247z7o2RM++cTqJCIiLUfF\nfxUN94hItFPxX6Wh+MNnAExEJLRU/FdJTYUffoAdO6xOIiLSMlT8V4mJ0XCPiEQ3FX8jVPwiEs00\nnbMRP/xgzu4pLweHw+o0IiI3pumcIdCmDWRn61O/iEQnFf91aLhHRKKVhnqu49w5c7jn4EHo1s3q\nNCIi16ehnhC5/Xb46U/hgw+sTiIiEloq/huYMAHefdfqFCIioaWhnhs4eRISEsxz9N9+u9VpREQa\np6GeEOrWDYYOhY0brU4iIhI6Kv6bmDBBs3tEJLpoqOcmKith8GBzuKdNG6vTiIhcS0M9IeZwQO/e\nUFpqdRIRkdBQ8TeBFnOJSDTRUE8T7NgBDz8Mhw6ZZ+8UEQknGuppAQMGQNu28I9/WJ1EROTWqfib\nQOfoF5FoouJvIhW/iEQLFX8TDRtmTuk8cMDqJCIit0bF30S33QaPPKJP/SIS+VT8zaDhHhGJBprO\n2Qx1deY5+rdvNxd1iYiEA03nbEEdOsDs2fDP/wynTlmdRkQkOPrE30yGATNmmBdiX78eOna0OpGI\n2F1zu1PFH4SLF+Ff/9Uc+vmf/9HJ20TEWq0y1HP48GFGjRrFgAEDyMrK4u233wbA7/eTk5ODw+Eg\nNzeXmpqawDaLFy8mMTGRlJQUSiP8jGexsfDnP0NNDfzqV+a3ABGRSBHUJ/6jR49y9OhRUlNTOXHi\nBOnp6Wzfvp2ioiIOHz7MggULmDlzJn369GHWrFlUVVUxcuRINmzYwMGDB5kxYwbl5eXXhomQT/wN\n/H7IyoLsbJgzx+o0ImJXrfKJv2fPnqSmpgLQvXt3BgwYQFlZGT6fj/z8fNq3b09eXh5erxcAr9fL\n+PHjcTgcZGZmYhgGfr8/mLcOK507w5o18NZbsHSp1WlERJrmlmf17Nu3j507d5Kenk5ZWRlOpxMA\np9OJz+cDzOJPTk4ObJOUlBR4LNLFxZkHeX/7W/jrX61OIyJyc7d0WNLv9zN58mQWLlxIp06dmndU\n+TrnN55z2ZhJVlYWWVlZtxKxVfTtC6tXw/jx5nV6MzOtTiQi0czj8eDxeIJ/ASNI33//vTFmzBhj\n4cKFgfseffRRo7y83DAMw9i6dasxceJEwzAMY9WqVUZhYWHgeYMGDTLOnDlzzWveQpywUFJiGD16\nGMb27VYnERE7aW53BjXUYxgG+fn5DBw4kOnTpwfud7vdFBcXU1tbS3FxMRkZGQCkp6ezfv16Kisr\n8Xg8xMbG0rlz5+D3VmHqpz+FxYvh5z83L9oiIhKOgprVU1paysiRI3nggQcCQza/+93vGD58OI8/\n/jjbtm1j8ODBvPXWW3Tq1AmARYsW8frrr9OuXTuWLVvGiBEjrg0TYbN6rmfxYvjDH+Cjj6B7d6vT\niEi00wKuMPHii7BpE/ztb3DHHVanEZFopuIPE4YB+fnmOfzfe8+8dKOISEvQSdrCREwM/Nd/mat8\np07V6l4RCR8q/hbUpg2sXAl79sC//7vVaURETCr+Fnb77eYc//ffh4ULrU4jInKLC7ikabp1g3Xr\nYPhwc6XvL35hdSIRsTMVfytxOGDtWnOu/913w5gxVicSEbvSUE8rGjgQ/vIX81z+W7danUZE7ErF\n38r+6Z9g+XJ45BHYu9fqNCJiRxrqsUBODlRVmSd1++gj8wLuIiKtRcVvkX/7N3Nx189+Bps3Q5cu\nVicSEbvQUI+FZs+GBx80h302boTqaqsTiYgd6JQNFquvh//8T/OcPtu2mVM/hwy58tatm9UpRSSc\n6Vw9EeziRfOA72efXbqVl2tnICI3puKPMo3tDLZtg65dL+0Ehg7VzkDEzlT8NnDxIuzbd+03g8t3\nBoMGmaeDjo2F22678s/G7mvKnx06QMeOVv/2InI1Fb9NXbwI+/ebC8M++wy++ALOnzePIVy8GJo/\nv/8eRoyAKVNgwgS4806rf2sRARW/tKCzZ80Tzr39Nng88NBD5k7g4Yf1TUDESip+aRWnTsFf/wor\nVpjfMLKzzZ3AQw/pojMirU3FL63u6FHzugMrVpjDTf/yL+ZOYPhw8/iAiLQsFb9Y6sABeOcdcydw\n+jRMnmzuBNLSzKuSiUjoqfglbOzYYe4AVqyAdu3gscfMnUBSktXJRKKLil/CjmGA12vuAFauhHvv\nNXcAkydDr15WpxOJfCp+CWv19eaMoBUr4N13ITnZPEvpmDHmQrTbbrM6oUjkUfFLxDh/3jxH0caN\nsGEDfPMNjB5t7gTGjoWEBKsTikQGFb9ErG++gZIScydQUgKdOl3aCYwaBXfdZXVCkfCk4peoYBjm\n6uOGbwMff2xeunLsWHNn4HZrvYBIAxW/RKW6OigtNXcEGzea6wWyssydwJgx0L+/pouKfan4xRaq\nqmDTpkvfCG677dJOICPD/DYQExPcDW78WMPjofhZJBRU/GI7hgG7d1/aCWzbZp5YzjCCuzW8ZmP3\nNzwW7M+NuXpnE4qfr/ceobovHCUmmseCRo0yV43fcYfViVqPil8kQjS2ownFz429TyjvC0cXL5oL\nBv/+d/O2bZu5WrxhRzBsmHla8Wil4hcR2zt7Fj766NKOYMcOcLku7QjcbnM1ebRQ8YuIXOXMGXNy\nQMOO4MsvzWNBDTuCoUMje5aYil9E5CZOnYItWy7tCA4eNI8LNOwIBg+OrFXkKn4RkWY6cQI2b760\nI/j6a/NqcwMHQu/eV97C8aCxil9E5BYdO2buCPbsga++gkOHzD8rK80V5X36XLtDaLh17dr6M6HC\nuvi3bNnC008/zQ8//EBhYSG//vWvrwyj4heRMHbxormG5KuvGr8dOmQ+73o7hd69IS4u9BcoCuvi\nT0tLY9GiRfTu3Ztx48ZRWlpK9+7dL4WJkOL3eDxkZWVZHeOmIiFnJGQE5Qy1aM1pGPDdd9ffMXz1\nFaxbZx5DCKXmdmerXRjv9OnTAIwcOZLevXszduxYvF5va719SHk8HqsjNEkk5IyEjKCcoRatOWNi\nzKGe1FTIyYHCQnjtNfjf/4WyMvPbQqhLPxitVvxlZWU4nc7A31NSUvj0009b6+1FRORHuhS2iIjd\nGK3ku+++M1JTUwN/nzZtmrF69eorntO3b18D0E033XTTrRm3vn37NquP29BK7rzzTsCc2eNwONi4\ncSMvv/zyFc/Zt29fa8UREbGtVit+gN///vc8/fTTXLhwgcLCwitm9IiISOsIqwVcIiLS8sLi4O6W\nLVtITk4mMTGR119/3eo4jTp8+DCjRo1iwIABZGVl8fbbb1sd6Ybq6+tJS0sjOzvb6ijXdfbsWZ56\n6in69+8f1rO8li9fzoMPPsiQIUOYPn261XEC8vLyiIuL4/777w/c5/f7ycnJweFwkJubS01NjYUJ\nTY3l/M1vfkNycjKDBw9m+vTp1NbWWpiw8YwNXnvtNWJjY6murrYg2ZWul/ONN94gOTmZAQMG8MIL\nL9z8hUJw3PaWpaamGps3bzYOHTpkJCUlGcePH7c60jW+/fZbY9u2bYZhGMbx48eNhIQE48yZMxan\nur7XXnvN+MUvfmFkZ2dbHeW6Zs6cacyePduora01Lly4YHz33XdWR7rGyZMnjT59+hg1NTVGfX29\n8bOf/cxYt26d1bEMwzCMLVu2GOXl5cbAgQMD982dO9eYNm2aUVdXZ/zqV78y5s+fb2FCU2M5N2zY\nYNTX1xv19fXG1KlTjT/+8Y8WJmw8o2EYRmVlpTFu3DijT58+xsmTJy1Kd0ljOb/44gsjIyPD2LNn\nj2EYhlFVVXXT17H8E3+kLOzq2bMnqampAHTv3p0BAwawdetWi1M17siRI6xZs4apU6eG9UrokpIS\nXnzxRTp06ECbNm0CEwDCSceOHTEMg9OnT1NbW8u5c+fo2rWr1bEAGDFixDVZfD4f+fn5tG/fnry8\nvLD4t9RYzjFjxhAbG0tsbCzjxo1j8+bNFqUzNZYR4Pnnn2fevHkWJGpcYznXrl1Lfn4+iYmJANx9\n9903fR3Liz8SF3bt27ePnTt3kp6ebnWURs2YMYP58+cTG+oTgoTQkSNHqKuro6CgALfbzdy5c6mr\nq7M61jU6duxIUVERffr0oWfPngwfPjxs/7/Dlf+enE4nPp/P4kQ3t3z58rAcknzvvfeIj4/ngQce\nsDrKDW3YsIEdO3YwdOhQpk6dyq5du266Tfg2Q5jy+/1MnjyZhQsXckcYnp919erV9OjRg7S0tLD+\ntF9XV8eePXuYOHEiHo+HnTt3snLlSqtjXeP48eMUFBSwa9cuDh06xCeffMIHH3xgdazrCuf/5435\n7W9/S+fOnZk0aZLVUa5w7tw5Xn31VV555ZXAfeH637auro7q6mo+/PBDcnJymDZt2k23sbz4XS4X\nu3fvDvx9586dZGRkWJjo+i5cuMDEiRN54oknyMnJsTpOoz7++GNWrVpFQkICU6ZM4W9/+xtPPvmk\n1bGu0a9fP5KSksjOzqZjx45MmTKFtWvXWh3rGj6fj4yMDPr160e3bt2YNGkSW7ZssTrWdblcLioq\nKgCoqKjA5XJZnOj6/vznP7N+/Xreeustq6NcY//+/Rw6dIhBgwaRkJDAkSNHGDJkCFVVVVZHu0ZG\nRgaTJ0+mY8eOZGdns3v37pt+e7a8+C9f2HXo0CE2btyI2+22ONW1DMMgPz+fgQMHhtXMjqu9+uqr\nHD58mIMHD/LOO+8wevRo3nzzTatjNSoxMRGv18vFixf54IMPeOihh6yOdI0RI0awdetWqqurOX/+\nPGvXrmXs2LFWx7out9tNcXExtbW1FBcXh+2HqHXr1jF//nxWrVpFhzC8Cvr999/PsWPHOHjwIAcP\nHiQ+Pp7y8nJ69OhhdbRrDBs2jLVr12IYBl6vl759+978v2nojzs3n8fjMZxOp9G3b19j0aJFVsdp\n1IcffmjExMQYgwYNMlJTU43U1FRj7dq1Vse6IY/HE9azer788kvD7XYbgwYNMmbOnGnU1NRYHalR\nb7zxhjFy5Ehj6NChxuzZs436+nqrIxmGYRiPPfaYcc899xjt2rUz4uPjjeLiYuPMmTPGI488YvTq\n1cvIyckx/H6/1TEDOdu2bWvEx8cbf/rTn4x+/foZDocj8G+poKAgLDJe/t/ycgkJCWExq6exnD/8\n8IPx9NNPG06n08jNzTV8Pt9NX0cLuEREbMbyoR4REWldKn4REZtR8YuI2IyKX0TEZlT8IiI2o+IX\nEbEZFb+IiM2o+EVEbOb/AQsFDtCkTc8sAAAAAElFTkSuQmCC\n",
+ "text": [
+ ""
+ ]
+ }
+ ],
+ "prompt_number": 10
+ }
+ ],
+ "metadata": {}
+ }
+ ]
+}
\ No newline at end of file