diff --git a/.gitignore b/.gitignore index bcfb3d73..5c39733b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ *.pyc -SimPEG.sublime-project -SimPEG.sublime-workspace -docs/_build/ \ No newline at end of file +*.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/setup.py b/SimPEG/setup.py new file mode 100644 index 00000000..c421cb4b --- /dev/null +++ b/SimPEG/setup.py @@ -0,0 +1,8 @@ +import os +print 'Compiling TriSolve.' +os.system('f2py -c utils/TriSolve.f -m TriSolve') +print 'TriSolve Compiled! yay.' +print 'Moving TriSolve into Utils.' +os.system('mv TriSolve.so utils/TriSolve.so') +print 'Thats it. Well Done Computer.' + 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/api_TestResults.rst b/SimPEG/tests/api_TestResults.rst deleted file mode 100644 index 0d63a328..00000000 --- a/SimPEG/tests/api_TestResults.rst +++ /dev/null @@ -1,355 +0,0 @@ -.. _api_TestResults: - -.. raw:: html - - - - - -
-

Test Report

-

Start Time: 2013-11-05 15:24:44

-

Duration: 0:00:00.007500

-

Status: Pass 22

- -

This demonstrates the report output by Prasanna.Yelsangikar.

-
- - - -

Show - Summary - Failed - All -

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Test Group/Test caseCountPassFailErrorView
test_basemesh.TestBaseMesh111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_basemesh.TestMeshNumbers2D111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
Total222200 
- diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py index 9b5cc0e7..f8c2964d 100644 --- a/SimPEG/tests/test_Solver.py +++ b/SimPEG/tests/test_Solver.py @@ -58,23 +58,39 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directLower_1(self): + def test_directLower_1_python(self): AL = sparse.tril(self.A) - solve = Solver(AL, doDirect=True, flag='L', options={}) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'}) e = np.ones(self.M.nC) rhs = AL.dot(e) x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directLower_M(self): + def test_directLower_M_python(self): AL = sparse.tril(self.A) - solve = Solver(AL, doDirect=True, flag='L', options={}) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'}) e = np.ones((self.M.nC,numRHS)) rhs = AL.dot(e) x = solve.solve(rhs) + + def test_directLower_1_fortran(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) + e = np.ones(self.M.nC) + rhs = AL.dot(e) + x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directUpper_1(self): + def test_directLower_M_fortran(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) + e = np.ones((self.M.nC,numRHS)) + rhs = AL.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_1_python(self): AU = sparse.triu(self.A) solve = Solver(AU, doDirect=True, flag='U', options={}) e = np.ones(self.M.nC) @@ -82,7 +98,7 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directUpper_M(self): + def test_directUpper_M_python(self): AU = sparse.triu(self.A) solve = Solver(AU, doDirect=True, flag='U', options={}) e = np.ones((self.M.nC,numRHS)) @@ -90,6 +106,23 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_1_fortran(self): + AU = sparse.triu(self.A) + solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'}) + e = np.ones(self.M.nC) + rhs = AU.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_M_fortran(self): + AU = sparse.triu(self.A) + solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'}) + e = np.ones((self.M.nC,numRHS)) + rhs = AU.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + def test_directDiagonal_1(self): AD = sdiag(self.A.diagonal()) solve = Solver(AD, doDirect=True, flag='D', options={}) 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 d9ac4a1d..a1194487 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -1,7 +1,22 @@ import numpy as np import scipy.sparse as sparse import scipy.sparse.linalg as linalg +from SimPEG.utils import mkvc +DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'} + +try: + import TriSolve +except Exception, e: + try: + import os + # Note: this may not work from SublimeText, if that is the case, just run the command in your shell. + os.system('f2py -c TriSolve.f -m TriSolve') + import TriSolve + except Exception, e: + print 'Warning: Python backend is being used for solver. Run setup.py from the command line.' + DEFAULTS['forward'] = 'python' + DEFAULTS['backward'] = 'python' class Solver(object): """ @@ -55,11 +70,11 @@ class Solver(object): elif self.flag is None and not self.doDirect: return self.solveIter(b, **self.options) elif self.flag == 'U': - return self.solveBackward(b) + return self.solveBackward(b, **self.options) elif self.flag == 'L': - return self.solveForward(b) + return self.solveForward(b, **self.options) elif self.flag == 'D': - return self.solveDiagonal(b) + return self.solveDiagonal(b, **self.options) else: raise Exception('Unknown flag.') pass @@ -69,7 +84,7 @@ class Solver(object): del self.dsolve self.dsolve = None - def solveDirect(self, b, factorize=False, backend='scipy'): + def solveDirect(self, b, factorize=False, backend=None): """ Use solve instead of this interface. @@ -78,6 +93,8 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['direct'] + assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' if factorize and self.dsolve is None: @@ -104,7 +121,7 @@ class Solver(object): def solveIter(self, b, M=None, iterSolver='CG'): pass - def solveBackward(self, b, backend='python'): + def solveBackward(self, b, backend=None): """ Use solve instead of this interface. @@ -114,21 +131,29 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['backward'] if type(self.A) is not sparse.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) vals = self.A.data rowptr = self.A.indptr colind = self.A.indices - x = np.empty_like(b) # empty() is faster than zeros(). - for i in reversed(xrange(self.A.shape[0])): - ith_row = vals[rowptr[i] : rowptr[i+1]] - cols = colind[rowptr[i] : rowptr[i+1]] - x_vals = x[cols] - x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0] + if backend == 'fortran': + if len(b.shape) == 1 or b.shape[1] == 1: + x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size, 1) + x = mkvc(x) + else: + x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1]) + elif backend == 'python': + x = np.empty_like(b) # empty() is faster than zeros(). + for i in reversed(xrange(self.A.shape[0])): + ith_row = vals[rowptr[i] : rowptr[i+1]] + cols = colind[rowptr[i] : rowptr[i+1]] + x_vals = x[cols] + x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0] return x - def solveForward(self, b, backend='python'): + def solveForward(self, b, backend=None): """ Use solve instead of this interface. @@ -138,21 +163,29 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['forward'] if type(self.A) is not sparse.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) vals = self.A.data rowptr = self.A.indptr colind = self.A.indices - x = np.empty_like(b) # empty() is faster than zeros(). - for i in xrange(self.A.shape[0]): - ith_row = vals[rowptr[i] : rowptr[i+1]] - cols = colind[rowptr[i] : rowptr[i+1]] - x_vals = x[cols] - x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1] + if backend == 'fortran': + if len(b.shape) == 1 or b.shape[1] == 1: + x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size, 1) + x = mkvc(x) + else: + x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1]) + elif backend == 'python': + x = np.empty_like(b) # empty() is faster than zeros(). + for i in xrange(self.A.shape[0]): + ith_row = vals[rowptr[i] : rowptr[i+1]] + cols = colind[rowptr[i] : rowptr[i+1]] + x_vals = x[cols] + x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1] return x - def solveDiagonal(self, b, backend='python'): + def solveDiagonal(self, b, backend=None): """ Use solve instead of this interface. @@ -162,6 +195,8 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['diagonal'] + diagA = self.A.diagonal() if len(b.shape) == 1 or b.shape[1] == 1: # Just one RHS @@ -205,3 +240,20 @@ if __name__ == '__main__': print np.linalg.norm(e-x,np.inf) + n = 6000 + A_dense = np.random.random((n,n)) + L = np.tril(np.dot(A_dense, A_dense)) # Positive definite is better conditioned. + e = np.ones(n) + b = np.dot(L, e) + + A = sparse.csr_matrix(L) + pSolve = Solver(A,flag='L',options={'backend':'python'}); + fSolve = Solver(A,flag='L',options={'backend':'fortran'}) + tic = time() + x = pSolve.solve(b) + toc = time() - tic + print 'Error Forward Python = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc + tic = time() + x = fSolve.solve(b) + toc = time() - tic + print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc diff --git a/SimPEG/utils/TriSolve.f b/SimPEG/utils/TriSolve.f new file mode 100644 index 00000000..ab9048d6 --- /dev/null +++ b/SimPEG/utils/TriSolve.f @@ -0,0 +1,64 @@ +c File TriSolve.f + subroutine forward(al, ial, jal, b, nv, n, nRHS, x) + double precision al(nv) + integer ial(n+1) + integer jal(nv) + double precision b(n,nRHS) + double precision x(n,nRHS) + integer nv + integer n + integer nRHS + integer rhs +cf2py intent(in) :: al +cf2py intent(in) :: ial +cf2py intent(in) :: jal +cf2py intent(in) :: b +cf2py intent(in) :: nv +cf2py intent(in) :: n +cf2py intent(in) :: nRHS +cf2py intent(out) :: x + real ( kind = 8 ) t + + do rhs = 1, nRHS + do k = 1, n + t = b(k,rhs) + do j = ial(k)+1, ial(k+1) + t = t - al(j) * x(jal(j)+1,rhs) + end do + x(k,rhs) = t/al(ial(k+1)) + end do + end do + end subroutine forward + + + subroutine backward(au,iau, jau, b, nv, n, nRHS, x) + double precision au(nv) + integer iau(n+1) + integer jau(nv) + double precision b(n,nRHS) + double precision x(n,nRHS) + integer nv + integer n + integer nRHS + integer rhs +cf2py intent(in) :: au +cf2py intent(in) :: iau +cf2py intent(in) :: jau +cf2py intent(in) :: b +cf2py intent(in) :: nv +cf2py intent(in) :: n +cf2py intent(in) :: nRHS +cf2py intent(out) :: x + real ( kind = 8 ) t + + do rhs = 1, nRHS + do k = n, 1, -1 + t = b(k,rhs) + do j = iau(k)+1, iau(k+1) + t = t - au(j) * x(jau(j)+1,rhs) + end do + x(k,rhs) = t/au(iau(k)+1) + end do + end do + + end subroutine backward diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 244c3103..d2170721 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,9 +3,59 @@ import sputils import lomutils import interputils import ModelBuilder -import Solver -from Solver import Solver from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag, ddx, av 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..7b66e131 --- /dev/null +++ b/SimPEG/utils/ipythonUtils.py @@ -0,0 +1,28 @@ +from tempfile import NamedTemporaryFile +import matplotlib.pyplot as plt +from matplotlib import animation + +# 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 8e432568..9e8b2df3 100644 --- a/docs/api_TestResults.rst +++ b/docs/api_TestResults.rst @@ -185,9 +185,9 @@ Test Results -->
-

Start Time: 2013-11-06 12:40:59

-

Duration: 0:00:32.198864

-

Status: Pass 108 Failure 1

+

Start Time: 2013-11-12 13:34:53

+

Duration: 0:00:31.732498

+

Status: Pass 117

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 4.588e+01 2.799e+01 nan - 1 1.00e-02 2.072e+00 2.831e-01 1.995 - 2 1.00e-03 1.817e-01 2.835e-03 1.999 - 3 1.00e-04 1.792e-02 2.835e-05 2.000 - 4 1.00e-05 1.789e-03 2.835e-07 2.000 - 5 1.00e-06 1.789e-04 2.831e-09 2.001 - 6 1.00e-07 1.789e-05 2.124e-11 2.125 + 0 1.00e-01 5.517e+01 4.463e+01 nan + 1 1.00e-02 1.516e+00 4.609e-01 1.986 + 2 1.00e-03 1.101e-01 4.624e-03 1.999 + 3 1.00e-04 1.059e-02 4.626e-05 2.000 + 4 1.00e-05 1.055e-03 4.626e-07 2.000 + 5 1.00e-06 1.055e-04 4.623e-09 2.000 + 6 1.00e-07 1.055e-05 4.504e-11 2.011 ========================= PASS! ========================= - That was easy! + Awesome, Rowan, just awesome. @@ -414,15 +414,15 @@ Test Results pt3.3: ==================== checkDerivative ==================== iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order --------------------------------------------------------- - 0 1.00e-01 9.439e-02 8.315e-03 nan - 1 1.00e-02 9.344e-03 8.324e-05 2.000 - 2 1.00e-03 2.787e+00 2.787e+00 -4.525 - 3 1.00e-04 9.336e-05 8.323e-09 8.525 - 4 1.00e-05 9.336e-06 8.323e-11 2.000 - 5 1.00e-06 9.336e-07 8.335e-13 1.999 - 6 1.00e-07 9.336e-08 9.031e-15 1.965 + 0 1.00e-01 1.180e-01 6.488e-03 nan + 1 1.00e-02 1.175e-02 6.401e-05 2.006 + 2 1.00e-03 1.175e-03 6.390e-07 2.001 + 3 1.00e-04 1.175e-04 6.389e-09 2.000 + 4 1.00e-05 1.175e-05 6.389e-11 2.000 + 5 1.00e-06 1.175e-06 6.392e-13 2.000 + 6 1.00e-07 1.175e-07 7.710e-15 1.919 ========================= PASS! ========================= - You are awesome. + Go Test Go! @@ -451,13 +451,13 @@ Test Results pt3.4: ==================== checkDerivative ==================== iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order --------------------------------------------------------- - 0 1.00e-01 8.866e+00 8.270e+00 nan - 1 1.00e-02 1.423e-01 8.270e-02 2.000 - 2 1.00e-03 6.787e-03 8.270e-04 2.000 - 3 1.00e-04 6.042e-04 8.270e-06 2.000 - 4 1.00e-05 5.968e-05 8.270e-08 2.000 - 5 1.00e-06 5.960e-06 8.270e-10 2.000 - 6 1.00e-07 5.960e-07 8.271e-12 2.000 + 0 1.00e-01 6.898e+00 7.074e+00 nan + 1 1.00e-02 5.316e-02 7.074e-02 2.000 + 2 1.00e-03 1.051e-03 7.074e-04 2.000 + 3 1.00e-04 1.688e-04 7.074e-06 2.000 + 4 1.00e-05 1.752e-05 7.074e-08 2.000 + 5 1.00e-06 1.758e-06 7.074e-10 2.000 + 6 1.00e-07 1.759e-07 7.076e-12 2.000 ========================= PASS! ========================= Testing is important. @@ -498,15 +498,15 @@ Test Results ==================== checkDerivative ==================== iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order --------------------------------------------------------- - 0 1.00e-01 2.817e-01 5.296e-17 nan - 1 1.00e-02 2.817e-02 6.505e-17 -0.089 - 2 1.00e-03 2.817e-03 4.768e-17 0.135 - 3 1.00e-04 2.817e-04 5.112e-17 -0.030 - 4 1.00e-05 2.817e-05 7.422e-17 -0.162 - 5 1.00e-06 2.817e-06 8.012e-17 -0.033 - 6 1.00e-07 2.817e-07 4.930e-17 0.211 + 0 1.00e-01 2.762e-01 4.105e-17 nan + 1 1.00e-02 2.762e-02 2.532e-17 0.210 + 2 1.00e-03 2.762e-03 4.066e-17 -0.206 + 3 1.00e-04 2.762e-04 4.297e-17 -0.024 + 4 1.00e-05 2.762e-05 5.494e-17 -0.107 + 5 1.00e-06 2.762e-06 4.696e-17 0.068 + 6 1.00e-07 2.762e-07 2.834e-17 0.219 ========================= PASS! ========================= - Once upon a time, a happy little test passed. + You deserve a pat on the back! @@ -535,15 +535,15 @@ Test Results pt4.2: ==================== checkDerivative ==================== iter h |J0-Jt| |J0+h*dJ'*dx-Jt| Order --------------------------------------------------------- - 0 1.00e-01 4.194e-01 8.200e-02 nan - 1 1.00e-02 3.456e-02 8.200e-04 2.000 - 2 1.00e-03 3.382e-03 8.200e-06 2.000 - 3 1.00e-04 3.375e-04 8.200e-08 2.000 - 4 1.00e-05 3.374e-05 8.200e-10 2.000 - 5 1.00e-06 3.374e-06 8.204e-12 2.000 - 6 1.00e-07 3.374e-07 8.252e-14 1.997 + 0 1.00e-01 5.835e-01 2.875e-01 nan + 1 1.00e-02 3.247e-02 2.875e-03 2.000 + 2 1.00e-03 2.989e-03 2.875e-05 2.000 + 3 1.00e-04 2.963e-04 2.875e-07 2.000 + 4 1.00e-05 2.960e-05 2.875e-09 2.000 + 5 1.00e-06 2.960e-06 2.875e-11 2.000 + 6 1.00e-07 2.960e-07 2.879e-13 1.999 ========================= PASS! ========================= - Testing is important. + Awesome, Rowan, just awesome. @@ -583,22 +583,22 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 2.73e-01 | - 16 | 7.11e-02 | 3.8392 | 1.9408 - 32 | 1.56e-02 | 4.5463 | 2.1847 + 8 | 2.14e-01 | + 16 | 6.31e-02 | 3.3843 | 1.7589 + 32 | 1.72e-02 | 3.6712 | 1.8763 --------------------------------------------- - You get a gold star! + Well done Rowan! randomTensorMesh: Interpolation 1D: CC _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 2.38e-01 | - 16 | 9.97e-02 | 2.3888 | 1.3742 - 32 | 3.33e-02 | 2.9985 | 1.6242 + 8 | 3.35e-01 | + 16 | 1.18e-01 | 2.8469 | 1.8016 + 32 | 2.27e-02 | 5.1806 | 3.9741 --------------------------------------------- - Yay passed! + Happy little convergence test! @@ -629,22 +629,22 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 2.22e-01 | - 16 | 5.39e-02 | 4.1229 | 2.0437 - 32 | 1.59e-02 | 3.3980 | 1.7647 + 8 | 2.98e-01 | + 16 | 6.69e-02 | 4.4487 | 2.1534 + 32 | 1.55e-02 | 4.3044 | 2.1058 --------------------------------------------- - That was easy! + Well done Rowan! randomTensorMesh: Interpolation 1D: N _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 4.51e-01 | - 16 | 9.69e-02 | 4.6544 | 2.0364 - 32 | 4.23e-02 | 2.2891 | 1.9120 + 8 | 7.16e-01 | + 16 | 1.60e-01 | 4.4731 | 1.8851 + 32 | 3.72e-02 | 4.3056 | 2.5981 --------------------------------------------- - And then everyone was happy. + Not just a pretty face Rowan @@ -655,11 +655,11 @@ Test Results - + test_interpolation.TestInterpolation2d 6 - 5 - 1 + 6 + 0 0 Detail @@ -684,24 +684,24 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 7.61e-02 | - 16 | 1.92e-02 | 3.9713 | 1.9896 - 32 | 4.75e-03 | 4.0298 | 2.0107 - 64 | 1.14e-03 | 4.1621 | 2.0573 + 8 | 7.48e-02 | + 16 | 1.79e-02 | 4.1840 | 2.0649 + 32 | 4.38e-03 | 4.0830 | 2.0296 + 64 | 1.16e-03 | 3.7604 | 1.9109 --------------------------------------------- - You get a gold star! + Testing is important. randomTensorMesh: Interpolation 2D: CC _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 1.54e-01 | - 16 | 3.87e-02 | 3.9740 | 2.8437 - 32 | 7.65e-03 | 5.0512 | 2.1271 - 64 | 3.67e-03 | 2.0838 | 0.9869 + 8 | 4.76e-02 | + 16 | 6.08e-02 | 0.7830 | -1.0896 + 32 | 1.14e-02 | 5.3162 | 1.7097 + 64 | 1.19e-03 | 9.6340 | 3.5674 --------------------------------------------- - Testing is important. + Well done Rowan! @@ -732,24 +732,24 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 6.94e-02 | - 16 | 1.84e-02 | 3.7623 | 1.9116 - 32 | 4.71e-03 | 3.9144 | 1.9688 - 64 | 1.17e-03 | 4.0345 | 2.0124 + 8 | 6.98e-02 | + 16 | 1.85e-02 | 3.7771 | 1.9173 + 32 | 4.76e-03 | 3.8827 | 1.9571 + 64 | 1.20e-03 | 3.9676 | 1.9883 --------------------------------------------- - And then everyone was happy. + That was easy! randomTensorMesh: Interpolation 2D: Ex _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 1.24e-01 | - 16 | 3.44e-02 | 3.6027 | 2.3303 - 32 | 1.47e-02 | 2.3355 | 1.0793 - 64 | 4.05e-03 | 3.6347 | 2.3831 + 8 | 2.61e-01 | + 16 | 3.89e-02 | 6.7177 | 2.4144 + 32 | 1.41e-02 | 2.7587 | 1.3921 + 64 | 4.59e-03 | 3.0696 | 1.9255 --------------------------------------------- - And then everyone was happy. + That was easy! @@ -780,24 +780,24 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 7.03e-02 | - 16 | 1.88e-02 | 3.7352 | 1.9012 - 32 | 4.77e-03 | 3.9519 | 1.9826 - 64 | 1.20e-03 | 3.9618 | 1.9862 + 8 | 6.97e-02 | + 16 | 1.87e-02 | 3.7164 | 1.8939 + 32 | 4.74e-03 | 3.9516 | 1.9825 + 64 | 1.14e-03 | 4.1658 | 2.0586 --------------------------------------------- - The test be workin! + That was easy! randomTensorMesh: Interpolation 2D: Ey _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 2.09e-01 | - 16 | 4.41e-02 | 4.7441 | 1.6758 - 32 | 1.04e-02 | 4.2565 | 2.2490 - 64 | 2.73e-03 | 3.8032 | 2.2280 + 8 | 1.82e-01 | + 16 | 3.85e-02 | 4.7305 | 2.2574 + 32 | 9.33e-03 | 4.1268 | 2.3056 + 64 | 3.29e-03 | 2.8348 | 1.1080 --------------------------------------------- - Yay passed! + Happy little convergence test! @@ -828,24 +828,24 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 7.61e-02 | - 16 | 1.92e-02 | 3.9713 | 1.9896 - 32 | 4.75e-03 | 4.0298 | 2.0107 - 64 | 1.14e-03 | 4.1621 | 2.0573 + 8 | 7.48e-02 | + 16 | 1.79e-02 | 4.1840 | 2.0649 + 32 | 4.38e-03 | 4.0830 | 2.0296 + 64 | 1.16e-03 | 3.7604 | 1.9109 --------------------------------------------- - The test be workin! + You are awesome. randomTensorMesh: Interpolation 2D: Fx _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 1.34e-01 | - 16 | 2.08e-02 | 6.4246 | 3.1408 - 32 | 9.86e-03 | 2.1089 | 1.1456 - 64 | 2.32e-03 | 4.2458 | 1.9464 + 8 | 9.78e-02 | + 16 | 4.87e-02 | 2.0098 | 1.1486 + 32 | 6.39e-03 | 7.6229 | 2.5269 + 64 | 3.06e-03 | 2.0873 | 1.0490 --------------------------------------------- - Yay passed! + You are awesome. @@ -876,10 +876,10 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 7.50e-02 | - 16 | 1.81e-02 | 4.1461 | 2.0517 - 32 | 4.70e-03 | 3.8447 | 1.9429 - 64 | 1.17e-03 | 4.0103 | 2.0037 + 8 | 7.59e-02 | + 16 | 1.90e-02 | 3.9914 | 1.9969 + 32 | 4.62e-03 | 4.1142 | 2.0406 + 64 | 1.15e-03 | 4.0248 | 2.0089 --------------------------------------------- Happy little convergence test! @@ -888,12 +888,12 @@ Test Results _____________________________________________ h | error | e(i-1)/e(i) | order ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~ - 8 | 7.60e-02 | - 16 | 3.01e-02 | 2.5304 | 1.4734 - 32 | 1.20e-02 | 2.5084 | 1.0636 - 64 | 2.94e-03 | 4.0744 | 1.9937 + 8 | 1.34e-01 | + 16 | 2.73e-02 | 4.9171 | 2.4681 + 32 | 1.10e-02 | 2.4900 | 1.2825 + 64 | 3.06e-03 | 3.5774 | 1.9074 --------------------------------------------- - You get a gold star! + Testing is important. @@ -904,52 +904,45 @@ Test Results - -
test_orderN
+ +
test_orderN
- - fail + + pass -