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/__init__.py b/SimPEG/__init__.py index f4628b6e..41f7c431 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -5,3 +5,4 @@ import inverse import visulize import forward import regularization +import examples diff --git a/SimPEG/examples/__init__.py b/SimPEG/examples/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 132966a8..34046a52 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -5,7 +5,6 @@ from SimPEG.utils import ModelBuilder, sdiag, mkvc from SimPEG import Solver import numpy as np import scipy.sparse as sp -import scipy.sparse.linalg as linalg class DCProblem(ModelTransforms.LogModel, Problem): @@ -16,7 +15,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/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 6a9a8032..2647cf44 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -102,7 +102,7 @@ class BaseMesh(object): elif xType in ['F', 'E']: # This will only deal with components of fields, not full 'F' or 'E' xx = mkvc(xx) # unwrap it in case it is a matrix - nn = self.nF if xType == 'F' else self.nE + nn = self.nFv if xType == 'F' else self.nEv nn = np.r_[0, nn] nx = [0, 0, 0] @@ -331,7 +331,7 @@ class BaseMesh(object): return locals() nEz = property(**nEz()) - def nE(): + def nEv(): doc = """ Total number of edges in each direction @@ -346,6 +346,18 @@ class BaseMesh(object): """ fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None]) return locals() + nEv = property(**nEv()) + + def nE(): + doc = """ + Total number of edges. + + :rtype: int + :return: sum([prod(nEx), prod(nEy), prod(nEz)]) + + """ + fget = lambda self: np.sum(self.nEv) + return locals() nE = property(**nE()) def nFx(): @@ -391,7 +403,7 @@ class BaseMesh(object): return locals() nFz = property(**nFz()) - def nF(): + def nFv(): doc = """ Total number of faces in each direction @@ -406,6 +418,19 @@ class BaseMesh(object): """ fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None]) return locals() + nFv = property(**nFv()) + + + def nF(): + doc = """ + Total number of faces. + + :rtype: int + :return: sum([prod(nFx), prod(nFy), prod(nFz)]) + + """ + fget = lambda self: np.sum(self.nFv) + return locals() nF = property(**nF()) def normals(): @@ -418,13 +443,13 @@ class BaseMesh(object): def fget(self): if self.dim == 2: - nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0])] - nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1])] + nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])] + nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])] return np.r_[nX, nY] elif self.dim == 3: - nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0]), np.zeros(self.nF[0])] - nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1]), np.zeros(self.nF[1])] - nZ = np.c_[np.zeros(self.nF[2]), np.zeros(self.nF[2]), np.ones(self.nF[2])] + nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])] + nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])] + nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])] return np.r_[nX, nY, nZ] return locals() normals = property(**normals()) @@ -439,13 +464,13 @@ class BaseMesh(object): def fget(self): if self.dim == 2: - tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0])] - tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1])] + tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])] + tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])] return np.r_[tX, tY] elif self.dim == 3: - tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0]), np.zeros(self.nE[0])] - tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1]), np.zeros(self.nE[1])] - tZ = np.c_[np.zeros(self.nE[2]), np.zeros(self.nE[2]), np.ones(self.nE[2])] + tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])] + tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])] + tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])] return np.r_[tX, tY, tZ] return locals() tangents = property(**tangents()) diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 93b82b25..fa2bb729 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -110,6 +110,12 @@ class Cyl1DMesh(object): return locals() nFz = property(**nFz()) + def nFv(): + doc = "Total number of faces in each direction" + fget = lambda self: np.array([self.nFr, self.nFz]) + return locals() + nFv = property(**nFv()) + def nF(): doc = "Total number of faces" fget = lambda self: self.nFr + self.nFz diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 598b392a..fa3f42e0 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -1,11 +1,6 @@ import numpy as np from scipy import sparse as sp -from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros - - -def ddx(n): - """Define 1D derivatives, inner, this means we go from n+1 to n+1""" - return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") +from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av def checkBC(bc): @@ -39,11 +34,6 @@ def ddxCellGrad(n, bc): return D -def av(n): - """Define 1D averaging operator from cell-centres to nodes.""" - return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") - - class DiffOperators(object): """ Class creates the differential operators that you need! @@ -107,6 +97,38 @@ class DiffOperators(object): _nodalGrad = None nodalGrad = property(**nodalGrad()) + def nodalLaplacian(): + doc = "Construct laplacian operator (nodes to edges)." + + def fget(self): + if(self._nodalLaplacian is None): + print 'Warning: Laplacian has not been tested rigorously.' + # The number of cell centers in each direction + n = self.n + # Compute divergence operator on faces + if(self.dim == 1): + D1 = sdiag(1./self.hx) * ddx(mesh.nCx) + L = - D1.T*D1 + elif(self.dim == 2): + D1 = sdiag(1./self.hx) * ddx(n[0]) + D2 = sdiag(1./self.hy) * ddx(n[1]) + L1 = sp.kron(speye(n[1]+1), - D1.T * D1) + L2 = sp.kron(- D2.T * D2, speye(n[0]+1)) + L = L1 + L2 + elif(self.dim == 3): + D1 = sdiag(1./self.hx) * ddx(n[0]) + D2 = sdiag(1./self.hy) * ddx(n[1]) + D3 = sdiag(1./self.hz) * ddx(n[2]) + L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1) + L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1)) + L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1)) + L = L1 + L2 + L3 + self._nodalLaplacian = L + return self._nodalLaplacian + return locals() + _nodalLaplacian = None + nodalLaplacian = property(**nodalLaplacian()) + def setCellGradBC(self, BC): """ Function that sets the boundary conditions for cell-centred derivative operators. @@ -182,7 +204,6 @@ class DiffOperators(object): return locals() cellGradx = property(**cellGradx()) - def cellGrady(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): @@ -203,8 +224,6 @@ class DiffOperators(object): return locals() cellGrady = property(**cellGrady()) - - def cellGradz(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): @@ -222,7 +241,6 @@ class DiffOperators(object): return locals() cellGradz = property(**cellGradz()) - def edgeCurl(): doc = "Construct the 3D curl operator." @@ -265,6 +283,8 @@ class DiffOperators(object): _edgeCurl = None edgeCurl = property(**edgeCurl()) + # --------------- Averaging --------------------- + def aveF2CC(): doc = "Construct the averaging operator on cell faces to cell centers." @@ -274,10 +294,10 @@ class DiffOperators(object): if(self.dim == 1): self._aveF2CC = av(n[0]) elif(self.dim == 2): - self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])), + self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])), sp.kron(av(n[1]), speye(n[0]))), format="csr") elif(self.dim == 3): - self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), kron3(speye(n[2]), av(n[1]), speye(n[0])), kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") return self._aveF2CC @@ -295,10 +315,10 @@ class DiffOperators(object): if(self.dim == 1): raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') elif(self.dim == 2): - self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])), sp.kron(speye(n[1]), av(n[0]))), format="csr") elif(self.dim == 3): - self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), kron3(av(n[2]), speye(n[1]), av(n[0])), kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._aveE2CC @@ -316,37 +336,57 @@ class DiffOperators(object): if(self.dim == 1): self._aveN2CC = av(n[0]) elif(self.dim == 2): - self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") + self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr() elif(self.dim == 3): - self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") + self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr() return self._aveN2CC return locals() _aveN2CC = None aveN2CC = property(**aveN2CC()) - def aveN2CCv(): - doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." + def aveN2E(): + doc = "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate." def fget(self): - if(self._aveN2CCv is None): + if(self._aveN2E is None): # The number of cell centers in each direction n = self.n if(self.dim == 1): - self._aveN2CCv = av(n[0]) + self._aveN2E = av(n[0]) elif(self.dim == 2): - self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") + self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])), + sp.kron(av(n[1]), speye(n[0]+1))), format="csr") elif(self.dim == 3): - self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveN2CCv + self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])), + kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)), + kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr") + return self._aveN2E return locals() - _aveN2CCv = None - aveN2CCv = property(**aveN2CCv()) + _aveN2E = None + aveN2E = property(**aveN2E()) + + def aveN2F(): + doc = "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate." + + def fget(self): + if(self._aveN2F is None): + # The number of cell centers in each direction + n = self.n + if(self.dim == 1): + self._aveN2F = av(n[0]) + elif(self.dim == 2): + self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)), + sp.kron(speye(n[1]+1), av(n[0]))), format="csr") + elif(self.dim == 3): + self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)), + kron3(av(n[2]), speye(n[1]+1), av(n[0])), + kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") + return self._aveN2F + return locals() + _aveN2F = None + aveN2F = property(**aveN2F()) + + # --------------- Methods --------------------- def getMass(self, materialProp=None, loc='e'): """ Produces mass matricies. diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index f0ac4ab0..5c2e7b5f 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -174,8 +174,8 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): def Pxxx(pos): ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0] - ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1] + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nFv[0] + ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nFv[0] + mesh.nFv[1] IND = np.r_[ind1, ind2, ind3].flatten() @@ -286,7 +286,7 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False): def Pxx(pos): ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nF[0] + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nFv[0] IND = np.r_[ind1, ind2].flatten() @@ -387,8 +387,8 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False): def Pxxx(pos): ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0] - ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1] + ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nEv[0] + ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nEv[0] + mesh.nEv[1] IND = np.r_[ind1, ind2, ind3].flatten() @@ -499,7 +499,7 @@ def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): def Pxx(pos): ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nE[0] + ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nEv[0] IND = np.r_[ind1, ind2].flatten() diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index b510a754..5c4a73db 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -45,8 +45,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridCC is None: - ccV = (self.aveN2CCv*mkvc(self.gridN)) - self._gridCC = ccV.reshape((-1, self.dim), order='F') + self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F') return self._gridCC return locals() _gridCC = None # Store grid by default diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 6cc5d8a6..90f8fd10 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -396,7 +396,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: - nF_nE = self.nF if 'F' in locType else self.nE + nF_nE = self.nFv if 'F' in locType else self.nEv components = [spzeros(loc.shape[0], n) for n in nF_nE] components[ind] = interpmat(loc, *self.getTensor(locType)) Q = sp.hstack(components) 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/runTests.py b/SimPEG/tests/runTests.py index c44f9a1e..58d94303 100644 --- a/SimPEG/tests/runTests.py +++ b/SimPEG/tests/runTests.py @@ -4,47 +4,54 @@ import unittest import HTMLTestRunner # This code will run all tests in directory named test_*.py +def main(html=False): + TITLE = 'Test Results' + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) -TITLE = 'Test Results' -test_file_strings = glob.glob('test_*.py') -module_strings = [str[0:len(str)-3] for str in test_file_strings] -suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str - in module_strings] -testSuite = unittest.TestSuite(suites) -unittest.TextTestRunner(verbosity=2).run(testSuite) + if not html: + unittest.TextTestRunner(verbosity=2).run(testSuite) + return -outfile = open("report.html", "w") -runner = HTMLTestRunner.HTMLTestRunner( - stream=outfile, - title=TITLE, - description='SimPEG Test Report was automatically generated.' - ) + outfile = open("report.html", "w") + runner = HTMLTestRunner.HTMLTestRunner( + stream=outfile, + title=TITLE, + description='SimPEG Test Report was automatically generated.', + verbosity=2 + ) -runner.run(testSuite) -outfile.close() + runner.run(testSuite) + outfile.close() -reader = open("report.html", "r") -writer = open("../../docs/api_TestResults.rst", "w") + reader = open("report.html", "r") + writer = open("../../docs/api_TestResults.rst", "w") -writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n') + writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n') -go = False -for line in reader: - skip = False - if line == '