diff --git a/README.md b/README.md index e107fea3..c2eafa26 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png) -Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. +Simulation and Parameter Estimation in Geoscience - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features: diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 1ea0601c..f421bc82 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,10 +1,12 @@ +import numpy as np +import scipy.sparse as sp import utils from utils import Solver import mesh +import data +import forward import inverse import visualize -import forward -import regularization import examples import scipy.version as _v diff --git a/SimPEG/data/__init__.py b/SimPEG/data/__init__.py new file mode 100644 index 00000000..86123003 --- /dev/null +++ b/SimPEG/data/__init__.py @@ -0,0 +1,20 @@ +from SimPEG import utils + + +class SimPEGData(object): + """Data holds the observed data, and the standard deviations.""" + + __metaclass__ = utils.Save.Savable + + std = None #: Estimated Standard Deviations + dobs = None #: Observed data + dtrue = None #: True data, if data is synthetic + mtrue = None #: True model, if data is synthetic + + def __init__(self, prob, **kwargs): + utils.setKwargs(self, **kwargs) + self.prob = prob + + def isSynthetic(self): + "Check if the data is synthetic." + return self.mtrue is not None diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/examples/DC.py similarity index 98% rename from SimPEG/forward/DCProblem.py rename to SimPEG/examples/DC.py index 074cdb8b..ad98d6c2 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/examples/DC.py @@ -168,7 +168,6 @@ def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): if __name__ == '__main__': - from SimPEG.regularization import Regularization from SimPEG import inverse import matplotlib.pyplot as plt @@ -217,7 +216,7 @@ if __name__ == '__main__': m0 = mesh.gridCC[:,0]*0+sig2 opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - reg = Regularization(mesh) + reg = inverse.Regularization(mesh) inv = inverse.Inversion(problem, reg, opt, beta0=1e4) # Check Derivative diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/examples/Linear.py similarity index 51% rename from SimPEG/forward/LinearProblem.py rename to SimPEG/examples/Linear.py index ef92957b..dbee7dff 100644 --- a/SimPEG/forward/LinearProblem.py +++ b/SimPEG/examples/Linear.py @@ -1,14 +1,13 @@ -import numpy as np -from SimPEG.mesh import TensorMesh -from SimPEG.forward import Problem -from SimPEG.regularization import Regularization -from SimPEG.inverse import * +from SimPEG import mesh, forward, inverse, np import matplotlib.pyplot as plt -class LinearProblem(Problem): +class LinearProblem(forward.Problem): """docstring for LinearProblem""" + def __init__(self, *args, **kwargs): + forward.Problem.__init__(self, *args, **kwargs) + def dpred(self, m, u=None): return self.G.dot(m) @@ -21,7 +20,7 @@ class LinearProblem(Problem): def example(N): h = np.ones(N)/N - M = TensorMesh([h]) + M = mesh.TensorMesh([h]) nk = 20 jk = np.linspace(1.,20.,nk) @@ -35,35 +34,27 @@ def example(N): for i in range(nk): G[i,:] = g(i) - - m_true = np.zeros(M.nC) - m_true[M.vectorCCx > 0.3] = 1. - m_true[M.vectorCCx > 0.45] = -0.5 - m_true[M.vectorCCx > 0.6] = 0 - - - d_true = G.dot(m_true) - noise = 0.1 * np.random.rand(d_true.size) - - d_obs = d_true + noise + mtrue = np.zeros(M.nC) + mtrue[M.vectorCCx > 0.3] = 1. + mtrue[M.vectorCCx > 0.45] = -0.5 + mtrue[M.vectorCCx > 0.6] = 0 prob = LinearProblem(M) prob.G = G - prob.dobs = d_obs - prob.std = np.ones_like(d_obs)*0.1 + data = prob.createSyntheticData(mtrue, std=0.01) - return prob, m_true + return prob, data if __name__ == '__main__': - prob, m_true = example(100) + prob, data = 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) + reg = inverse.Regularization(M) + opt = inverse.InexactGaussNewton(maxIter=20) + inv = inverse.Inversion(prob,reg,opt,data) + m0 = np.zeros_like(data.mtrue) mrec = inv.run(m0) @@ -73,9 +64,7 @@ if __name__ == '__main__': plt.figure(2) - plt.plot(M.vectorCCx, m_true, 'b-') + plt.plot(M.vectorCCx, data.mtrue, 'b-') plt.plot(M.vectorCCx, mrec, 'r-') - - plt.show() diff --git a/SimPEG/examples/__init__.py b/SimPEG/examples/__init__.py index e69de29b..a5c37345 100644 --- a/SimPEG/examples/__init__.py +++ b/SimPEG/examples/__init__.py @@ -0,0 +1,2 @@ +import DC +import Linear diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 99747f7d..1412ff5b 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -1,6 +1,4 @@ -import numpy as np -from SimPEG.utils import mkvc, sdiag, count, timeIt -import scipy.sparse as sp +from SimPEG import utils, data, np, sp norm = np.linalg.norm @@ -37,9 +35,13 @@ class Problem(object): to (locally) find how model parameters change the data, and optimize! """ - counter = None + __metaclass__ = utils.Save.Savable - def __init__(self, mesh): + counter = None #: A SimPEG.utils.Counter object + + + def __init__(self, mesh, *args, **kwargs): + utils.setKwargs(self, **kwargs) self.mesh = mesh @property @@ -65,27 +67,7 @@ class Problem(object): def P(self, value): self._P = value - @property - def std(self): - """ - Estimated Standard Deviations. - """ - return self._std - @std.setter - def std(self, value): - self._std = value - - @property - def dobs(self): - """ - Observed data. - """ - return self._dobs - @dobs.setter - def dobs(self, value): - self._dobs = value - - @count + @utils.count def dpred(self, m, u=None): """ Predicted data. @@ -97,8 +79,8 @@ class Problem(object): u = self.field(m) return self.P*u - @count - def dataResidual(self, m, u=None): + @utils.count + def dataResidual(self, m, data, u=None): """ :param numpy.array m: geophysical model :param numpy.array u: fields @@ -115,9 +97,9 @@ class Problem(object): u is the field of interest; d_obs is the observed data. """ - return self.dpred(m, u=u) - self.dobs + return self.dpred(m, u=u) - data.dobs - @timeIt + @utils.timeIt def J(self, m, v, u=None): """ :param numpy.array m: model @@ -147,7 +129,7 @@ class Problem(object): """ raise NotImplementedError('J is not yet implemented.') - @timeIt + @utils.timeIt def Jt(self, m, v, u=None): """ :param numpy.array m: model @@ -161,7 +143,7 @@ class Problem(object): raise NotImplementedError('Jt is not yet implemented.') - @timeIt + @utils.timeIt def J_approx(self, m, v, u=None): """ @@ -176,7 +158,7 @@ class Problem(object): """ return self.J(m, v, u) - @timeIt + @utils.timeIt def Jt_approx(self, m, v, u=None): """ :param numpy.array m: model @@ -238,12 +220,11 @@ class Problem(object): Returns the observed data with random Gaussian noise and Wd which is the same size as data, and can be used to weight the inversion. """ - dobs = self.dpred(m,u=u) - noise = std*abs(dobs)*np.random.randn(*dobs.shape) - dobs = dobs+noise - eps = np.linalg.norm(mkvc(dobs),2)*1e-5 - Wd = 1/(abs(dobs)*std+eps) - return dobs, Wd + dtrue = self.dpred(m,u=u) + noise = std*abs(dtrue)*np.random.randn(*dtrue.shape) + dobs = dtrue+noise + stdev = dobs*0 + std + return data.SimPEGData(self, dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m) diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index 33c9a6b1..ce6f7e52 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1,4 +1,2 @@ from Problem import * -import DCProblem -from LinearProblem import LinearProblem import ModelTransforms diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 160678aa..926f719a 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,13 +1,14 @@ -import numpy as np -import scipy.sparse as sp import SimPEG -from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt, callHooks +from SimPEG import utils, sp, np from Optimize import Remember from BetaSchedule import Cooling from SimPEG.inverse import IterationPrinters, StoppingCriteria class BaseInversion(object): - """docstring for BaseInversion""" + """BaseInversion(prob, reg, opt, data, **kwargs) + """ + + __metaclass__ = utils.Save.Savable maxIter = 1 #: Maximum number of iterations name = 'BaseInversion' @@ -17,14 +18,15 @@ class BaseInversion(object): comment = '' #: Used by some functions to indicate what is going on in the algorithm counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things - beta0 = None #: The initial Beta (regularization parameter) + beta0 = None #: The initial Beta (regularization parameter) + beta0_ratio = 0.1 #: When beta0 is set to None, estimateBeta0 is used with this ratio - - def __init__(self, prob, reg, opt, **kwargs): - setKwargs(self, **kwargs) + def __init__(self, prob, reg, opt, data, **kwargs): + utils.setKwargs(self, **kwargs) self.prob = prob self.reg = reg self.opt = opt + self.data = data self.opt.parent = self self.stoppers = [StoppingCriteria.iteration] @@ -46,8 +48,8 @@ class BaseInversion(object): Standard deviation weighting matrix. """ if getattr(self,'_Wd',None) is None: - eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5 - self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) + eps = np.linalg.norm(utils.mkvc(self.data.dobs),2)*1e-5 + self._Wd = 1/(abs(self.data.dobs)*self.data.std+eps) return self._Wd @Wd.setter def Wd(self, value): @@ -63,14 +65,14 @@ class BaseInversion(object): Note that we do not set the target if it is None, but we return the default value. """ if getattr(self, '_phi_d_target', None) is None: - return self.prob.dobs.size # + return self.data.dobs.size # return self._phi_d_target @phi_d_target.setter def phi_d_target(self, value): self._phi_d_target = value - @timeIt + @utils.timeIt def run(self, m0): """run(m0) @@ -85,25 +87,19 @@ class BaseInversion(object): if self.stoppingCriteria(): break self.printDone() + self.finish() + return self.m + @utils.callHooks('startup') 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 """ - callHooks(self,'startup',m0) if not hasattr(self.reg, '_mref'): print 'Regularization has not set mref. SimPEG will set it to m0.' @@ -115,43 +111,25 @@ class BaseInversion(object): self.phi_d_last = np.nan self.phi_m_last = np.nan + @utils.callHooks('doStartIteration') def doStartIteration(self): """ **doStartIteration** 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 _doStartIteration*(self): - pass - - Where the * can be any string. If present, _doStartIteration* will be called at the start of the default doStartIteration call. - You may also completely overwrite this function. - :rtype: None :return: None """ - callHooks(self,'doStartIteration') - self._beta = self.getBeta() + @utils.callHooks('doEndIteration') 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): - 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. - :rtype: None :return: None """ - callHooks(self,'doEndIteration') - # store old values self.phi_d_last = self.phi_d self.phi_m_last = self.phi_m @@ -203,7 +181,7 @@ class BaseInversion(object): def stoppingCriteria(self): if self.debug: print 'checking stoppingCriteria' - return checkStoppers(self, self.stoppers) + return utils.checkStoppers(self, self.stoppers) def printDone(self): @@ -211,9 +189,17 @@ class BaseInversion(object): **printDone** is called at the end of the inversion routine. """ - printStoppers(self, self.stoppers) + utils.printStoppers(self, self.stoppers) - @timeIt + @utils.callHooks('finish') + def finish(self): + """finish() + + **finish** is called at the end of the optimization. + """ + pass + + @utils.timeIt def evalFunction(self, m, return_g=True, return_H=True): """evalFunction(m, return_g=True, return_H=True) @@ -223,7 +209,7 @@ class BaseInversion(object): u = self.prob.field(m) if self._iter is 0 and self._beta is None: - self._beta = self.beta0 = self.estimateBeta0(u=u) + self._beta = self.beta0 = self.estimateBeta0(u=u,ratio=self.beta0_ratio) phi_d = self.dataObj(m, u) phi_m = self.reg.modelObj(m) @@ -253,7 +239,7 @@ class BaseInversion(object): out += (operator,) return out if len(out) > 1 else out[0] - @timeIt + @utils.timeIt def dataObj(self, m, u=None): """dataObj(m, u=None) @@ -272,11 +258,11 @@ class BaseInversion(object): u is the field of interest; d_obs is the observed data; and W is the weighting matrix. """ # TODO: ensure that this is a data is vector and Wd is a matrix. - R = self.Wd*self.prob.dataResidual(m, u=u) - R = mkvc(R) + R = self.Wd*self.prob.dataResidual(m, self.data, u=u) + R = utils.mkvc(R) return 0.5*np.vdot(R, R) - @timeIt + @utils.timeIt def dataObjDeriv(self, m, u=None): """dataObjDeriv(m, u=None) @@ -312,13 +298,13 @@ class BaseInversion(object): if u is None: u = self.prob.field(m) - R = self.Wd*self.prob.dataResidual(m, u=u) + R = self.Wd*self.prob.dataResidual(m, self.data, u=u) dmisfit = self.prob.Jt(m, self.Wd * R, u=u) return dmisfit - @timeIt + @utils.timeIt def dataObj2Deriv(self, m, v, u=None): """dataObj2Deriv(m, v, u=None) @@ -357,7 +343,7 @@ class BaseInversion(object): if u is None: u = self.prob.field(m) - R = self.Wd*self.prob.dataResidual(m, u=u) + R = self.Wd*self.prob.dataResidual(m, self.data, u=u) # TODO: abstract to different norms a little cleaner. # \/ it goes here. in l2 it is the identity. @@ -365,13 +351,19 @@ class BaseInversion(object): return dmisfit + def save(self, group): + group.attrs['phi_d'] = self.phi_d + group.attrs['phi_m'] = self.phi_m + group.setArray('m', self.m) + group.setArray('dpred', self.dpred) + class Inversion(Cooling, Remember, BaseInversion): maxIter = 10 name = "SimPEG Inversion" - def __init__(self, prob, reg, opt, **kwargs): - BaseInversion.__init__(self, prob, reg, opt, **kwargs) + def __init__(self, prob, reg, opt, data, **kwargs): + BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) @@ -387,8 +379,8 @@ class TimeSteppingInversion(Remember, BaseInversion): maxIter = 1 name = "Time-Stepping SimPEG Inversion" - def __init__(self, prob, reg, opt, **kwargs): - BaseInversion.__init__(self, prob, reg, opt, **kwargs) + def __init__(self, prob, reg, opt, data, **kwargs): + BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 7e934e13..9d5edcc2 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -1,9 +1,6 @@ -import numpy as np +from SimPEG import Solver, utils, sp, np import matplotlib.pyplot as plt -from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt, callHooks norm = np.linalg.norm -import scipy.sparse as sp -from SimPEG import Solver __all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters'] @@ -85,6 +82,8 @@ class Minimize(object): Minimize is a general class for derivative based optimization. """ + __metaclass__ = utils.Save.Savable + name = "General Optimization Algorithm" #: The name of the optimization algorithm maxIter = 20 #: Maximum number of iterations @@ -110,9 +109,9 @@ class Minimize(object): 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) + utils.setKwargs(self, **kwargs) - @timeIt + @utils.timeIt def minimize(self, evalFunction, x0): """minimize(evalFunction, x0) @@ -155,6 +154,7 @@ class Minimize(object): doEndIteration(xt) printDone() + finish() return xc """ self.evalFunction = evalFunction @@ -175,6 +175,7 @@ class Minimize(object): self.doEndIteration(xt) self.printDone() + self.finish() return self.xc @@ -188,6 +189,7 @@ class Minimize(object): def parent(self, value): self._parent = value + @utils.callHooks('startup') def startup(self, x0): """ **startup** is called at the start of any new minimize call. @@ -198,19 +200,10 @@ 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 """ - callHooks(self,'startup',x0) self._iter = 0 self._iterLS = 0 @@ -221,7 +214,8 @@ class Minimize(object): self.f_last = np.nan self.x_last = x0 - @count + @utils.count + @utils.callHooks('doStartIteration') def doStartIteration(self): """doStartIteration() @@ -230,7 +224,8 @@ class Minimize(object): :rtype: None :return: None """ - callHooks(self,'doStartIteration') + pass + def printInit(self, inLS=False): """ @@ -242,8 +237,9 @@ class Minimize(object): """ 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) + utils.printTitles(self, self.printers if not inLS else self.printersLS, name, pad) + @utils.callHooks('printIter') def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -252,10 +248,8 @@ class Minimize(object): parent.printIter function and call that. """ - callHooks(self,'printIter',inLS) - pad = ' '*10 if inLS else '' - printLine(self, self.printers if not inLS else self.printersLS, pad=pad) + utils.printLine(self, self.printers if not inLS else self.printersLS, pad=pad) def printDone(self, inLS=False): """ @@ -268,15 +262,29 @@ class Minimize(object): 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) + utils.printStoppers(self, stoppers, pad='', stop=stop, done=done) + + + @utils.callHooks('finish') + def finish(self): + """finish() + + **finish** is called at the end of the optimization. + + :rtype: None + :return: None + + """ + pass def stoppingCriteria(self, inLS=False): if self._iter == 0: self.f0 = self.f self.g0 = self.g - return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) + return utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) - @timeIt + @utils.timeIt + @utils.callHooks('projection') def projection(self, p): """projection(p) @@ -288,10 +296,9 @@ class Minimize(object): :rtype: numpy.ndarray :return: p, projected search direction """ - callHooks(self,'projection',p) return p - @timeIt + @utils.timeIt def findSearchDirection(self): """findSearchDirection() @@ -322,7 +329,7 @@ class Minimize(object): """ return -self.g - @count + @utils.count def scaleSearchDirection(self, p): """scaleSearchDirection(p) @@ -341,7 +348,7 @@ class Minimize(object): nameLS = "Armijo linesearch" #: The line-search name - @timeIt + @utils.timeIt def modifySearchDirection(self, p): """modifySearchDirection(p) @@ -379,7 +386,7 @@ class Minimize(object): return self._LS_xt, self._iterLS < self.maxIterLS - @count + @utils.count def modifySearchDirectionBreak(self, p): """modifySearchDirectionBreak(p) @@ -401,7 +408,8 @@ class Minimize(object): print 'The linesearch got broken. Boo.' return p, False - @count + @utils.count + @utils.callHooks('doEndIteration') def doEndIteration(self, xt): """doEndIteration(xt) @@ -411,21 +419,10 @@ class Minimize(object): self.xc must be updated in this code. - - 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 """ - callHooks(self,'doEndIteration',xt) - # store old values self.f_last = self.f self.x_last, self.xc = self.xc, xt @@ -433,6 +430,19 @@ class Minimize(object): if self.debug: self.printDone() + def save(self, group): + group.setArray('searchDirection', self.searchDirection) + + if getattr(self,'parent',None) is None: + group.setArray('x', self.xc) + else: # Assume inversion is the parent + group.attrs['phi_d'] = self.parent.phi_d + group.attrs['phi_m'] = self.parent.phi_m + group.setArray('m', self.xc) + group.setArray('dpred', self.parent.dpred) + + + class Remember(object): """ This mixin remembers all the things you tend to forget. @@ -517,7 +527,7 @@ class ProjectedGradient(Minimize, Remember): self.aSet_prev = self.activeSet(x0) - @count + @utils.count def projection(self, x): """projection(x) @@ -526,7 +536,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.median(np.c_[self.lower,x,self.upper],axis=1) - @count + @utils.count def activeSet(self, x): """activeSet(x) @@ -535,7 +545,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.logical_or(x == self.lower, x == self.upper) - @count + @utils.count def inactiveSet(self, x): """inactiveSet(x) @@ -544,7 +554,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.logical_not(self.activeSet(x)) - @count + @utils.count def bindingSet(self, x): """bindingSet(x) @@ -557,7 +567,7 @@ class ProjectedGradient(Minimize, Remember): bind_low = np.logical_and(x == self.upper, self.g <= 0) return np.logical_or(bind_up, bind_low) - @timeIt + @utils.timeIt def findSearchDirection(self): """findSearchDirection() @@ -602,7 +612,7 @@ class ProjectedGradient(Minimize, Remember): # aSet_after = self.activeSet(self.xc+p) return p - @timeIt + @utils.timeIt def _doEndIteration_ProjectedGradient(self, xt): """_doEndIteration_ProjectedGradient(xt)""" aSet = self.activeSet(xt) @@ -630,13 +640,16 @@ class ProjectedGradient(Minimize, Remember): 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 + if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingPG class BFGS(Minimize, Remember): name = 'BFGS' nbfgs = 10 + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + @property def bfgsH0(self): """ @@ -702,7 +715,10 @@ class BFGS(Minimize, Remember): class GaussNewton(Minimize, Remember): name = 'Gauss Newton' - @timeIt + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + + @utils.timeIt def findSearchDirection(self): return Solver(self.H).solve(-self.g) @@ -749,7 +765,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value - @timeIt + @utils.timeIt def findSearchDirection(self): Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG}) p = Hinv.solve(-self.g) @@ -759,7 +775,10 @@ class InexactGaussNewton(BFGS, Minimize, Remember): class SteepestDescent(Minimize, Remember): name = 'Steepest Descent' - @timeIt + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + + @utils.timeIt def findSearchDirection(self): return -self.g @@ -792,15 +811,26 @@ class NewtonRoot(object): doLS = True def __init__(self, **kwargs): - setKwargs(self, **kwargs) + utils.setKwargs(self, **kwargs) def root(self, fun, x): + """root(fun, x) + + Function Should have the form:: + + def evalFunction(x, return_g=False): + out = (f,) + if return_g: + out += (g,) + return out if len(out) > 1 else out[0] + + """ if self.comments: print 'Newton Method:\n' self._iter = 0 while True: - [r,J] = fun(x); + r, J = fun(x, return_g=True) if self.solveTol == 0: Jinv = Solver(J) dh = - Jinv.solve(r) @@ -812,13 +842,12 @@ class NewtonRoot(object): muLS = 1. LScnt = 1 xt = x + dh - rt, Jt = fun(xt) # TODO: get rid of Jt + rt = fun(xt, return_g=False) - if self.comments: print '\tLinesearch:\n' + if self.comments and self.doLS: print '\tLinesearch:\n' # Enter Linesearch while True and self.doLS: - if self.comments: - print '\t\tResid: %e\n'%norm(rt) + if self.comments: print '\t\tResid: %e\n'%norm(rt) if norm(rt) <= norm(r) or norm(rt) < self.tol: break @@ -827,14 +856,16 @@ class NewtonRoot(object): print '.' if LScnt > self.maxLS: print 'Newton Method: Line search break.' - root = NaN - return + return None xt = x + muLS*dh - rt, Jt = fun(xt) # TODO: get rid of Jt + rt = fun(xt, return_g=False) x = xt self._iter += 1 - if norm(rt) < self.tol or self._iter > self.maxIter: + if norm(rt) < self.tol: + break + if self._iter > self.maxIter: + print 'NewtonRoot stopped by maxIters. norm: %4.4e' % norm(rt) break return x @@ -854,7 +885,7 @@ if __name__ == '__main__': print 'test the newtonRoot finding.' - fun = lambda x: (np.sin(x), sdiag(np.cos(x))) + fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), utils.sdiag( np.cos(x) ) ) x = np.array([np.pi-0.3, np.pi+0.1, 0]) - pnt = NewtonRoot(comments=False).root(fun,x) + pnt = NewtonRoot(comments=True).root(fun,x) print pnt diff --git a/SimPEG/inverse/Regularization.py b/SimPEG/inverse/Regularization.py new file mode 100644 index 00000000..8c6f19ef --- /dev/null +++ b/SimPEG/inverse/Regularization.py @@ -0,0 +1,219 @@ +from SimPEG import utils, np, sp + +class Regularization(object): + """**Regularization** + + Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same: + + .. math:: + + R(m) = \int_\Omega \\frac{\\alpha_x}{2}\left(\\frac{\partial m}{\partial x}\\right)^2 + \\frac{\\alpha_y}{2}\left(\\frac{\partial m}{\partial y}\\right)^2 \partial v + + Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this: + + .. math:: + + R(m) \\approx \sum_{ij} \left[\\frac{\\alpha_x}{2}\left[\left(\\frac{m_{i+1,j} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i-1,j}}{h}\\right)^2\\right] + + \\frac{\\alpha_y}{2}\left[\left(\\frac{m_{i,j+1} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i,j-1}}{h}\\right)^2\\right] + \\right]h^2 + + If we let D_1 be the derivative matrix in the x direction + + .. math:: + + \mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1 + + .. math:: + + \mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1 + + Where d_1 is the one dimensional derivative: + + .. math:: + + \mathbf{d}_1 = \\frac{1}{h} \left[ \\begin{array}{cccc} + -1 & 1 & & \\\\ + & \ddots & \ddots&\\\\ + & & -1 & 1\end{array} \\right] + + .. math:: + + R(m) \\approx \mathbf{v}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \\frac{\\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \\right] + + Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise) + + .. math:: + + \mathbf{a\odot b} = \\text{diag}(\mathbf{a})\mathbf{b} = \\text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a} + + and the transpose also is true (but the sizes have to make sense...): + + .. math:: + + \mathbf{a}^\\top\\text{diag}(\mathbf{b}) = \mathbf{b}^\\top\\text{diag}(\mathbf{a}) + + So R(m) can simplify to: + + .. math:: + + R(m) \\approx \mathbf{m}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{D}_1^\\top \\text{diag}(\mathbf{A}_1^\\top\mathbf{v}) \mathbf{D}_1 + \\frac{\\alpha_y}{2}\mathbf{D}_2^\\top \\text{diag}(\mathbf{A}_2^\\top \mathbf{v}) \mathbf{D}_2 \\right] \mathbf{m} + + We will define W_x as: + + .. math:: + + \mathbf{W}_x = \sqrt{\\alpha_x}\\text{diag}\left(\sqrt{\mathbf{A}_1^\\top\mathbf{v}}\\right) \mathbf{D}_1 + + + And then W as a tall matrix of all of the different regularization terms: + + .. math:: + + \mathbf{W} = \left[ \\begin{array}{c} + \mathbf{W}_s\\\\ + \mathbf{W}_x\\\\ + \mathbf{W}_y\end{array} \\right] + + Then we can write + + .. math:: + + R(m) \\approx \\frac{1}{2}\mathbf{m^\\top W^\\top W m} + + + """ + + __metaclass__ = utils.Save.Savable + + alpha_s = utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight") + alpha_x = utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + alpha_xx = utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") + alpha_yy = utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") + alpha_zz = utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") + + counter = None + + def __init__(self, mesh, **kwargs): + utils.setKwargs(self, **kwargs) + self.mesh = mesh + + @property + def mref(self): + if getattr(self, '_mref', None) is None: + return np.zeros(self.mesh.nC); + return self._mref + @mref.setter + def mref(self, value): + self._mref = value + + @property + def Ws(self): + """Regularization matrix Ws""" + if getattr(self,'_Ws', None) is None: + self._Ws = utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + return self._Ws + + @property + def Wx(self): + """Regularization matrix Wx""" + if getattr(self, '_Wx', None) is None: + Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFv[0]].T*self.mesh.vol + self._Wx = utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx + return self._Wx + + @property + def Wy(self): + """Regularization matrix Wy""" + if getattr(self, '_Wy', None) is None: + Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2])].T*self.mesh.vol + self._Wy = utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady + return self._Wy + + @property + def Wz(self): + """Regularization matrix Wz""" + if getattr(self, '_Wz', None) is None: + Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.nFv[:2]):].T*self.mesh.vol + self._Wz = utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz + return self._Wz + + @property + def Wxx(self): + """Regularization matrix Wxx""" + if getattr(self, '_Wxx', None) is None: + self._Wxx = utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + return self._Wxx + + @property + def Wyy(self): + """Regularization matrix Wyy""" + if getattr(self, '_Wyy', None) is None: + self._Wyy = utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + return self._Wyy + + @property + def Wzz(self): + """Regularization matrix Wzz""" + if getattr(self, '_Wzz', None) is None: + self._Wzz = utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + return self._Wzz + + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Ws, self.Wx, self.Wxx) + if self.mesh.dim > 1: + wlist += (self.Wy, self.Wyy) + if self.mesh.dim > 2: + wlist += (self.Wz, self.Wzz) + self._W = sp.vstack(wlist) + return self._W + + + @utils.timeIt + def modelObj(self, m): + r = self.W * (m - self.mref) + return 0.5*r.dot(r) + + @utils.timeIt + def modelObjDeriv(self, m): + """ + + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} + + """ + return self.W.T * ( self.W * (m - self.mref) ) + + @utils.timeIt + def modelObj2Deriv(self): + """ + + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the second derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W} + + """ + return self.W.T * self.W + diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py index 14bddce7..d83a8269 100644 --- a/SimPEG/inverse/__init__.py +++ b/SimPEG/inverse/__init__.py @@ -1,3 +1,4 @@ from Optimize import * from Inversion import * +from Regularization import Regularization import BetaSchedule diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 62ec4251..6151d20b 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -1,6 +1,5 @@ import numpy as np -from SimPEG.utils import mkvc - +from SimPEG import utils class BaseMesh(object): @@ -12,6 +11,7 @@ class BaseMesh(object): :param numpy.array,list x0: Origin of the mesh (dim, ) """ + def __init__(self, n, x0=None): # Check inputs @@ -78,7 +78,7 @@ class BaseMesh(object): x_array = np.ones((x.size, len(x))) # Unwrap it and put it in a np array for i, xi in enumerate(x): - x_array[:, i] = mkvc(xi) + x_array[:, i] = utils.mkvc(xi) x = x_array assert type(x) == np.ndarray, "x must be a numpy array" @@ -91,7 +91,7 @@ class BaseMesh(object): if format == 'M': return xx.reshape(nn, order='F') elif format == 'V': - return mkvc(xx) + return utils.mkvc(xx) def switchKernal(xx): """Switches over the different options.""" @@ -101,7 +101,7 @@ class BaseMesh(object): return outKernal(xx, nn) 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 + xx = utils.mkvc(xx) # unwrap it in case it is a matrix nn = self.nFv if xType == 'F' else self.nEv nn = np.r_[0, nn] @@ -308,7 +308,7 @@ class BaseMesh(object): """ fget = lambda self: np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None]) return locals() - nNv = property(**nNv()) + nNv = property(**nNv()) def nEx(): doc = """ diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index 5c2e7b5f..45853071 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -81,9 +81,9 @@ class InnerProducts(object): def getFaceInnerProduct(self, mu=None, returnP=False): """Wrapper function, - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct` - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D` """ if self.dim == 2: return getFaceInnerProduct2D(self, mu, returnP) @@ -93,9 +93,9 @@ class InnerProducts(object): def getEdgeInnerProduct(self, sigma=None, returnP=False): """Wrapper function, - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct` - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D` """ if self.dim == 2: return getEdgeInnerProduct2D(self, sigma, returnP) diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index 5c4a73db..b3dcd095 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -1,15 +1,14 @@ -import numpy as np +from SimPEG import utils, np from BaseMesh import BaseMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts from LomView import LomView -from SimPEG.utils import mkvc, ndgrid, volTetra, indexCube, faceInfo # Some helper functions. length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 -normalize2D = lambda x: x/np.kron(np.ones((1, 2)), mkvc(length2D(x), 2)) -normalize3D = lambda x: x/np.kron(np.ones((1, 3)), mkvc(length3D(x), 2)) +normalize2D = lambda x: x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2)) +normalize3D = lambda x: x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): @@ -21,6 +20,9 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): .. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py """ + + __metaclass__ = utils.Save.Savable + _meshType = 'LOM' def __init__(self, nodes): @@ -38,7 +40,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): # Save nodes to private variable _gridN as vectors self._gridN = np.ones((nodes[0].size, self.dim)) for i, node_i in enumerate(nodes): - self._gridN[:, i] = mkvc(node_i.astype(float)) + self._gridN[:, i] = utils.mkvc(node_i.astype(float)) def gridCC(): doc = "Cell-centered grid." @@ -69,10 +71,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridFx is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] + XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] self._gridFx = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] + XYZ = [utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridFx return locals() @@ -86,10 +88,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridFy is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] + XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] self._gridFy = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] + XYZ = [utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridFy return locals() @@ -102,7 +104,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridFz is None and self.dim == 3: N = self.r(self.gridN, 'N', 'N', 'M') - XYZ = [mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] + XYZ = [utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridFz return locals() @@ -116,10 +118,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridEx is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] + XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] self._gridEx = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] + XYZ = [utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridEx return locals() @@ -133,10 +135,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if self._gridEy is None: N = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - XY = [mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] + XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] self._gridEy = np.c_[XY[0], XY[1]] elif self.dim == 3: - XYZ = [mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] + XYZ = [utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridEy return locals() @@ -149,7 +151,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridEz is None and self.dim == 3: N = self.r(self.gridN, 'N', 'N', 'M') - XYZ = [mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] + XYZ = [utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]] return self._gridEz return locals() @@ -192,25 +194,25 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if(self._vol is None): if self.dim == 2: - A, B, C, D = indexCube('ABCD', self.n+1) - normal, area = faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D) + A, B, C, D = utils.indexCube('ABCD', self.n+1) + normal, area = utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D) self._vol = area elif self.dim == 3: # Each polyhedron can be decomposed into 5 tetrahedrons # However, this presents a choice so we may as well divide in two ways and average. - A, B, C, D, E, F, G, H = indexCube('ABCDEFGH', self.n+1) + A, B, C, D, E, F, G, H = utils.indexCube('ABCDEFGH', self.n+1) - vol1 = (volTetra(self.gridN, A, B, D, E) + # cutted edge top - volTetra(self.gridN, B, E, F, G) + # cutted edge top - volTetra(self.gridN, B, D, E, G) + # middle - volTetra(self.gridN, B, C, D, G) + # cutted edge bottom - volTetra(self.gridN, D, E, G, H)) # cutted edge bottom + vol1 = (utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top + utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top + utils.volTetra(self.gridN, B, D, E, G) + # middle + utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom + utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom - vol2 = (volTetra(self.gridN, A, F, B, C) + # cutted edge top - volTetra(self.gridN, A, E, F, H) + # cutted edge top - volTetra(self.gridN, A, H, F, C) + # middle - volTetra(self.gridN, C, H, D, A) + # cutted edge bottom - volTetra(self.gridN, C, G, H, F)) # cutted edge bottom + vol2 = (utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top + utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top + utils.volTetra(self.gridN, A, H, F, C) + # middle + utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom + utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom self._vol = (vol1 + vol2)/2 return self._vol @@ -226,30 +228,30 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): # Compute areas of cell faces if(self.dim == 2): xy = self.gridN - A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) + A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) edge1 = xy[B, :] - xy[A, :] normal1 = np.c_[edge1[:, 1], -edge1[:, 0]] area1 = length2D(edge1) - A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) + A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) # Note that we are doing A-D to make sure the normal points the right way. # Think about it. Look at the picture. Normal points towards C iff you do this. edge2 = xy[A, :] - xy[D, :] normal2 = np.c_[edge2[:, 1], -edge2[:, 0]] area2 = length2D(edge2) - self._area = np.r_[mkvc(area1), mkvc(area2)] + self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)] self._normals = [normalize2D(normal1), normalize2D(normal2)] elif(self.dim == 3): - A, E, F, B = indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz])) - normal1, area1 = faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False) + A, E, F, B = utils.indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz])) + normal1, area1 = utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False) - A, D, H, E = indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz])) - normal2, area2 = faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False) + A, D, H, E = utils.indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz])) + normal2, area2 = utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False) - A, B, C, D = indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz])) - normal3, area3 = faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False) + A, B, C, D = utils.indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz])) + normal3, area3 = utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False) - self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] + self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)] self._normals = [normal1, normal2, normal3] return self._area return locals() @@ -289,21 +291,21 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): if(self._edge is None or self._tangents is None): if(self.dim == 2): xy = self.gridN - A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) + A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) edge1 = xy[D, :] - xy[A, :] - A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) + A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) edge2 = xy[B, :] - xy[A, :] - self._edge = np.r_[mkvc(length2D(edge1)), mkvc(length2D(edge2))] + self._edge = np.r_[utils.mkvc(length2D(edge1)), utils.mkvc(length2D(edge2))] self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge] elif(self.dim == 3): xyz = self.gridN - A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz])) + A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz])) edge1 = xyz[D, :] - xyz[A, :] - A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz])) + A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz])) edge2 = xyz[B, :] - xyz[A, :] - A, E = indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz])) + A, E = utils.indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz])) edge3 = xyz[E, :] - xyz[A, :] - self._edge = np.r_[mkvc(length3D(edge1)), mkvc(length3D(edge2)), mkvc(length3D(edge3))] + self._edge = np.r_[utils.mkvc(length3D(edge1)), utils.mkvc(length3D(edge2)), utils.mkvc(length3D(edge3))] self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge] return self._edge return locals() @@ -329,10 +331,10 @@ if __name__ == '__main__': h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)]) dee3 = True if dee3: - X, Y, Z = ndgrid(h1, h2, h3, vector=False) + X, Y, Z = utils.ndgrid(h1, h2, h3, vector=False) M = LogicallyOrthogonalMesh([X, Y, Z]) else: - X, Y = ndgrid(h1, h2, vector=False) + X, Y = utils.ndgrid(h1, h2, vector=False) M = LogicallyOrthogonalMesh([X, Y]) print M.r(M.normals, 'F', 'Fx', 'V') diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 4c0841b8..d4ab86b0 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -1,11 +1,8 @@ -import numpy as np -import scipy.sparse as sp +from SimPEG import utils, np, sp from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -from SimPEG.utils import ndgrid, mkvc, spzeros, interpmat - class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): """ @@ -24,13 +21,20 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): Example of a padded tensor mesh: - .. plot:: examples/mesh/plot_TensorMesh.py + .. plot:: + + from SimPEG import mesh, utils + M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() For a quick tensor mesh on a (10x12x15) unit cube:: mesh = TensorMesh([10, 12, 15]) """ + + __metaclass__ = utils.Save.Savable + _meshType = 'TENSOR' def __init__(self, h_in, x0=None): @@ -48,7 +52,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" # Ensure h contains 1D vectors - self._h = [mkvc(x.astype(float)) for x in h] + self._h = [utils.mkvc(x.astype(float)) for x in h] def __str__(self): outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) @@ -166,7 +170,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridCC is None: - self._gridCC = ndgrid(self.getTensor('CC')) + self._gridCC = utils.ndgrid(self.getTensor('CC')) return self._gridCC return locals() _gridCC = None # Store grid by default @@ -177,7 +181,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridN is None: - self._gridN = ndgrid(self.getTensor('N')) + self._gridN = utils.ndgrid(self.getTensor('N')) return self._gridN return locals() _gridN = None # Store grid by default @@ -188,7 +192,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFx is None: - self._gridFx = ndgrid(self.getTensor('Fx')) + self._gridFx = utils.ndgrid(self.getTensor('Fx')) return self._gridFx return locals() _gridFx = None # Store grid by default @@ -199,7 +203,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFy is None and self.dim > 1: - self._gridFy = ndgrid(self.getTensor('Fy')) + self._gridFy = utils.ndgrid(self.getTensor('Fy')) return self._gridFy return locals() _gridFy = None # Store grid by default @@ -210,7 +214,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFz is None and self.dim > 2: - self._gridFz = ndgrid(self.getTensor('Fz')) + self._gridFz = utils.ndgrid(self.getTensor('Fz')) return self._gridFz return locals() _gridFz = None # Store grid by default @@ -221,7 +225,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEx is None: - self._gridEx = ndgrid(self.getTensor('Ex')) + self._gridEx = utils.ndgrid(self.getTensor('Ex')) return self._gridEx return locals() _gridEx = None # Store grid by default @@ -232,7 +236,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEy is None and self.dim > 1: - self._gridEy = ndgrid(self.getTensor('Ey')) + self._gridEy = utils.ndgrid(self.getTensor('Ey')) return self._gridEy return locals() _gridEy = None # Store grid by default @@ -243,7 +247,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEz is None and self.dim > 2: - self._gridEz = ndgrid(self.getTensor('Ez')) + self._gridEz = utils.ndgrid(self.getTensor('Ez')) return self._gridEz return locals() _gridEz = None # Store grid by default @@ -258,13 +262,13 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): vh = self.h # Compute cell volumes if(self.dim == 1): - self._vol = mkvc(vh[0]) + self._vol = utils.mkvc(vh[0]) elif(self.dim == 2): # Cell sizes in each direction - self._vol = mkvc(np.outer(vh[0], vh[1])) + self._vol = utils.mkvc(np.outer(vh[0], vh[1])) elif(self.dim == 3): # Cell sizes in each direction - self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2])) + self._vol = utils.mkvc(np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) return self._vol return locals() _vol = None @@ -285,12 +289,12 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): elif(self.dim == 2): area1 = np.outer(np.ones(n[0]+1), vh[1]) area2 = np.outer(vh[0], np.ones(n[1]+1)) - self._area = np.r_[mkvc(area1), mkvc(area2)] + self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)] elif(self.dim == 3): - area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2]))) - area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] + area1 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)] return self._area return locals() _area = None @@ -307,16 +311,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): n = self.n # Compute edge lengths if(self.dim == 1): - self._edge = mkvc(vh[0]) + self._edge = utils.mkvc(vh[0]) elif(self.dim == 2): l1 = np.outer(vh[0], np.ones(n[1]+1)) l2 = np.outer(np.ones(n[0]+1), vh[1]) - self._edge = np.r_[mkvc(l1), mkvc(l2)] + self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2)] elif(self.dim == 3): - l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) - l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)] + l1 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2), utils.mkvc(l3)] return self._edge return locals() _edge = None @@ -406,11 +410,11 @@ 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.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)) + components = [utils.spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = utils.interpmat(loc, *self.getTensor(locType)) Q = sp.hstack(components) elif locType in ['CC', 'N']: - Q = interpmat(loc, *self.getTensor(locType)) + Q = utils.interpmat(loc, *self.getTensor(locType)) else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) return Q diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 9956565d..a745c1cd 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -372,7 +372,7 @@ class TensorView(object): return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)]) - def video(mesh,var,function,figsize=(10,8),colorbar=True): + def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1): """ Call a function for a list of models to create a video. @@ -393,9 +393,11 @@ class TensorView(object): if colorbar: plt.colorbar(function(var[0],ax,clim,tlt,0)) - def animateFrame(i): + frames = np.arange(0,len(var),skip) + def animateFrame(j): + i = frames[j] function(var[i],ax,clim,tlt,i) - return animate(fig, animateFrame, frames=len(var)) + return animate(fig, animateFrame, frames=len(frames)) diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py deleted file mode 100644 index 41938152..00000000 --- a/SimPEG/regularization/Regularization.py +++ /dev/null @@ -1,119 +0,0 @@ -from SimPEG.utils import sdiag, count, timeIt, setKwargs -import numpy as np - -class Regularization(object): - """docstring for Regularization""" - - @property - def mref(self): - if getattr(self, '_mref', None) is None: - return np.zeros(self.mesh.nC); - return self._mref - @mref.setter - def mref(self, value): - self._mref = value - - @property - def Ws(self): - if getattr(self,'_Ws', None) is None: - self._Ws = sdiag(self.mesh.vol) - return self._Ws - - @property - def Wx(self): - if getattr(self, '_Wx', None) is None: - self._Wx = self.mesh.cellGradx*sdiag(self.mesh.vol) - return self._Wx - - @property - def Wy(self): - if getattr(self, '_Wy', None) is None: - self._Wy = self.mesh.cellGrady*sdiag(self.mesh.vol) - return self._Wy - - @property - def Wz(self): - if getattr(self, '_Wz', None) is None: - self._Wz = self.mesh.cellGradz*sdiag(self.mesh.vol) - return self._Wz - - alpha_s = 1e-6 - alpha_x = 1.0 - alpha_y = 1.0 - alpha_z = 1.0 - - counter = None - - def __init__(self, mesh, **kwargs): - setKwargs(self, **kwargs) - self.mesh = mesh - - - def pnorm(self, r): - return 0.5*r.dot(r) - - @timeIt - def modelObj(self, m): - mresid = m - self.mref - - mobj = self.alpha_s * self.pnorm( self.Ws * mresid ) - - mobj += self.alpha_x * self.pnorm( self.Wx * mresid ) - - if self.mesh.dim > 1: - mobj += self.alpha_y * self.pnorm( self.Wy * mresid ) - if self.mesh.dim > 2: - mobj += self.alpha_z * self.pnorm( self.Wz * mresid ) - - return mobj - - @timeIt - def modelObjDeriv(self, m): - """ - - In 1D: - - .. math:: - - m_{\\text{obj}} = {1 \over 2}\\alpha_s \left\| W_s (m- m_{\\text{ref}})\\right\|^2_2 - + {1 \over 2}\\alpha_x \left\| W_x (m- m_{\\text{ref}})\\right\|^2_2 - - \\frac{ \partial m_{\\text{obj}} }{\partial m} = - \\alpha_s W_s^{\\top} W_s (m - m_{\\text{ref}}) + - \\alpha_x W_x^{\\top} W_x (m - m_{\\text{ref}}) - - - \\frac{ \partial^2 m_{\\text{obj}} }{\partial m^2} = - \\alpha_s W_s^{\\top} W_s + - \\alpha_x W_x^{\\top} W_x - - """ - - mresid = m - self.mref - - mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid) - - mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid) - - if self.mesh.dim > 1: - mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid) - if self.mesh.dim > 2: - mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid) - - return mobjDeriv - - - @timeIt - def modelObj2Deriv(self): - - mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws - - mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx - - if self.mesh.dim > 1: - mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy - if self.mesh.dim > 2: - mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz - - return mobj2Deriv - diff --git a/SimPEG/regularization/__init__.py b/SimPEG/regularization/__init__.py deleted file mode 100644 index 0230f8c3..00000000 --- a/SimPEG/regularization/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from Regularization import Regularization diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index 847d1055..90312fdf 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -3,10 +3,9 @@ import unittest from SimPEG.mesh import TensorMesh from SimPEG.utils import ModelBuilder, sdiag from SimPEG.forward import Problem -from SimPEG.forward.DCProblem import * +from SimPEG.examples.DC import * from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve -from SimPEG.regularization import Regularization from SimPEG import inverse @@ -44,23 +43,19 @@ class DCProblemTests(unittest.TestCase): problem = DCProblem(mesh) problem.P = P problem.RHS = q - dobs, Wd = problem.createSyntheticData(mSynth, std=0.05) + data = problem.createSyntheticData(mSynth, std=0.05) # Now set up the problem to do some minimization - problem.W = Wd - problem.dobs = dobs - problem.std = dobs*0 + 0.05 - opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - reg = Regularization(mesh) - inv = inverse.Inversion(problem, reg, opt, beta0=1e4) + reg = inverse.Regularization(mesh) + inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4) self.inv = inv self.reg = reg self.p = problem self.mesh = mesh self.m0 = mSynth - self.dobs = dobs + self.data = data def test_misfit(self): derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)] @@ -71,7 +66,7 @@ class DCProblemTests(unittest.TestCase): # Adjoint Test u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1]) v = np.random.rand(self.mesh.nC) - w = np.random.rand(self.dobs.shape[0]) + w = np.random.rand(self.data.dobs.shape[0]) wtJv = w.dot(self.p.J(self.m0, v, u=u)) vtJtw = v.dot(self.p.Jt(self.m0, w, u=u)) passed = (wtJv - vtJtw) < 1e-10 diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py index ba49efd8..1d4de045 100644 --- a/SimPEG/tests/test_forward_problem.py +++ b/SimPEG/tests/test_forward_problem.py @@ -1,8 +1,6 @@ import numpy as np import unittest -from SimPEG.mesh import TensorMesh -from SimPEG.forward import Problem -from SimPEG.regularization import Regularization +from SimPEG import mesh, forward, inverse from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve @@ -14,9 +12,9 @@ class ProblemTests(unittest.TestCase): a = np.array([1, 1, 1]) b = np.array([1, 2]) c = np.array([1, 4]) - self.mesh2 = TensorMesh([a, b], np.array([3, 5])) - self.p2 = Problem(self.mesh2) - self.reg = Regularization(self.mesh2) + self.mesh2 = mesh.TensorMesh([a, b], np.array([3, 5])) + self.p2 = forward.Problem(self.mesh2) + self.reg = inverse.Regularization(self.mesh2) def test_modelTransform(self): print 'SimPEG.forward.Problem: Testing Model Transform' diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/tests/test_optimizers.py index 6fff75d3..0253852c 100644 --- a/SimPEG/tests/test_optimizers.py +++ b/SimPEG/tests/test_optimizers.py @@ -32,7 +32,7 @@ class TestOptimizers(unittest.TestCase): self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) def test_ProjGradient_quadraticBounded(self): - PG = inverse.ProjectedGradient() + PG = inverse.ProjectedGradient(debug=True) PG.lower, PG.upper = -2, 2 xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0])) x_true = np.array([2.,2.]) @@ -50,5 +50,15 @@ class TestOptimizers(unittest.TestCase): print 'x_true: ', x_true self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) + def test_NewtonRoot(self): + fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), sdiag( np.cos(x) ) ) + x = np.array([np.pi-0.3, np.pi+0.1, 0]) + xopt = inverse.NewtonRoot(comments=False).root(fun,x) + x_true = np.array([np.pi,np.pi,0]) + print 'Newton Root Finding' + 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/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py index 56df6250..969551ed 100644 --- a/SimPEG/utils/ModelBuilder.py +++ b/SimPEG/utils/ModelBuilder.py @@ -172,7 +172,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]): if len(shape) is 1: smth = np.array([1,10.,1],dtype=float) elif len(shape) is 2: - smth = np.array([[1,2,1],[7,10,7],[1,2,1]],dtype=float) + smth = np.array([[1,7,1],[2,10,2],[1,7,1]],dtype=float) elif len(shape) is 3: kernal = np.array([1,4,1], dtype=float).reshape((1,3)) smth = np.array(sp.kron(sp.kron(kernal,kernal.T).todense()[:],kernal).todense()).reshape((3,3,3)) diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py new file mode 100644 index 00000000..4ad550de --- /dev/null +++ b/SimPEG/utils/Save.py @@ -0,0 +1,348 @@ +import numpy as np +import time +import re + +try: + import h5py +except Exception, e: + print 'Warning: SimPEG.utils.Save needs h5py to be installed.' + + +SAVEABLES = {} + +def natural_keys(text): + ''' + alist.sort(key=natural_keys) sorts in human order + http://nedbatchelder.com/blog/200712/human_sorting.html + (See Toothy's implementation in the comments) + ''' + atoi = lambda text: int(text) if text.isdigit() else text + return [ atoi(c) for c in re.split('(\d+)', text) ] + + +def preIteration(group): + group.attrs['complete'] = False + group.attrs['time'] = time.time() + +def postIteration(group): + group.attrs['time'] = time.time() - group.attrs['time'] + group.attrs['date'] = time.ctime() + group.attrs['complete'] = True + +class SimPEGTable: + """ + This is a wrapper class on the HDF5 file. + """ + def __init__(self, name, mode='a'): + if '.hdf5' not in name: + name += '.hdf5' + self.f = h5py.File(name, mode) + self.root = hdf5Group(self,self.f) + + self.inversions = hdf5InversionGroup(self,self.root.addGroup('inversions',soft=True)) + + def show(self): self.root.show() + + def saveInversion(self, invObj): + + # Create a new inversion anytime this is run. + def _startup_hdf5_inv(invObj, m0): + node = self.inversions.addGroup('%d'%self.inversions.numChildren) + saveSavable(invObj,node.addGroup('rebuild')) + results = node.addGroup('results') + preIteration(results) + invObj._invNode = results + self.f.flush() + invObj.hook(_startup_hdf5_inv, overwrite=True) + + # At the start of every iteration we will create a inversion iteration node. + def _doStartIteration_hdf5_inv(invObj): + invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) + preIteration(invObj._invNodeIt) + invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) + + def _doEndIteration_hdf5_inv(invObj): + invObj.save(invObj._invNodeIt) + postIteration(invObj._invNodeIt) + self.f.flush() + invObj.hook(_doEndIteration_hdf5_inv, overwrite=True) + + # Delete all iterates that did not finish. + def _finish_hdf5_inv(invObj): + postIteration(invObj._invNode) + for it in invObj._invNode: + if not it.attrs['complete']: + del self.f[it.path] + del invObj._invNode + self.f.flush() + invObj.hook(_finish_hdf5_inv, overwrite=True) + + def _doStartIteration_hdf5_opt(optObj): + optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) + preIteration(optObj._optNodeIt) + invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True) + + def _doEndIteration_hdf5_opt(optObj, xt): + optObj.save(optObj._optNodeIt) + postIteration(optObj._optNodeIt) + self.f.flush() + invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True) + + + +class hdf5Group(object): + """Has some low level support for wrapping the native HDF5-Group class""" + + def __init__(self, T, groupNode): + self.T = T + # check if you are inputing a hdf5Group rather than a raw node, and act accordingly + if issubclass(groupNode.__class__, hdf5Group): + self.node = groupNode.node + else: + self.node = groupNode + + self.childClass = hdf5Group + self.parentClass = hdf5Group + + @property + def children(self): + """Children names in a list + + Use obj[name] to return the actual node. + """ + myChildren = [c for c in self.node] + myChildren.sort(key=natural_keys) + return myChildren + + @property + def numChildren(self): + """Returns the len(children)""" + return len(self.children) + + @property + def parent(self): + """Returns parent node""" + return self.parentClass(self.T, self.node.parent) + + @property + def name(self): + return self.path.split('/')[-1] + + @property + def path(self): + """Returns the root path of the group""" + return self.node.name + + @property + def attrs(self): + """Returns a list of attributes in the group""" + return self.node.attrs + + def addGroup(self, name, soft=False): + """Adds a child group to the current node.""" + if name in self.children and soft: + return self[name] + assert name not in self.children, 'Already a child called: '+self.path+'/'+name + return self.childClass(self.T, self.node.create_group(name)) + + def setArray(self, name, data): + a = self.node.create_dataset(name, data.shape) + a[...] = data + return a + + def __getitem__(self, val): + if type(val) is int: + val = self.children[val] + child = self.node[val] + if type(child) is h5py.Group: + child = self.childClass(self.T, child) + return child + + def __contains__(self, key): + return key in self.children + + def show(self, pad='', maxDepth=1, depth=0): + """ + Recursively show the structure of the database. + + For example:: + + + - + - + - + - + - + """ + s = self.__str__() + pad += ' '*4 + if maxDepth <= 0: print s + if depth >= maxDepth: return s + + for c in self.children: + if issubclass(self[c].__class__, hdf5Group): + s += '\n%s- %s' % (pad, self[c].show(pad=pad,depth=depth+1,maxDepth=maxDepth)) + else: + s += '\n%s- %s' % (pad, self[c].__str__()) + if depth is 0: + print s + else: + return s + + def __str__(self): + return '<%s "%s" (%d member%s, attrs=[%s])>' % (self.__class__.__name__, self.path, self.numChildren, '' if self.numChildren == 1 else 's',', '.join([a for a in self.attrs])) + + + +class hdf5InversionGroup(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.childClass = hdf5Inversion + +class hdf5Inversion(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.parentClass = hdf5InversionGroup + self.childClass = hdf5InversionResults + + def rebuild(self): + return loadSavable(self['rebuild']) + +class hdf5InversionResults(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.parentClass = hdf5Inversion + self.childClass = hdf5InversionIteration + +class hdf5InversionIteration(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.parentClass = hdf5InversionResults + + + +class Savable(type): + def __new__(cls, name, bases, attrs): + __init__ = attrs['__init__'] + def init_wrapper(self, *args, **kwargs): + self._args_init = args + self._kwargs_init = kwargs + return __init__(self, *args,**kwargs) + attrs['__init__'] = init_wrapper + + newClass = super(Savable, cls).__new__(cls, name, bases, attrs) + SAVEABLES[name] = newClass + return newClass + + +def saveSavable(obj, group, debug=False): + """ + This creates softlinks if _savable exists in children object. + + The first object is always created. + """ + assert type(obj.__class__) is Savable, 'Can only save objects that are Savable objects.' + + def doSave(grp, name, val): + if debug: print name, val + if type(val.__class__) is Savable: + link = getattr(val,'_savable',None) + if link is not None: + group.node[name] = h5py.SoftLink(link.path) + if debug: 'Created a softlink path to %s' % link.path + else: + subgrp = grp.addGroup(name) + saveSavable(val, subgrp, debug=debug) + elif type(val) in [list, tuple]: + # Split up, and save each element + for i, v in enumerate(val): + doSave(grp, name + '[%d]'%i, v) + elif type(val) is np.ndarray: + grp.setArray(name, val) + elif val is None: + grp.attrs[name] = 'None' + else: + # just try saving it as an attr + try: + grp.attrs[name] = val + except Exception, e: + print 'Warning: Could not save %s, problems may arise in loading.' % name + + group.attrs['__class__'] = obj.__class__.__name__ + for arg in obj._kwargs_init: + doSave(group, '_kwarg_'+arg, obj._kwargs_init[arg]) + for i, arg in enumerate(obj._args_init): + doSave(group, '_arg%d'%i, arg) + obj._savable = group + + +def loadSavable(node, pointers=None): + """ + pointers allow things that point to the same node in the h5py file to + be returned as the same object, if they have already been created. + """ + + if pointers is None: pointers = [] + for pointer in pointers: + if pointer._savable.node == node.node: return pointer + + args = ([a for a in node.attrs if '_arg' in a] + [a for a in node.children if '_arg' in a]) + kwargs = ([a for a in node.attrs if '_kwarg' in a] + [a for a in node.children if '_kwarg' in a]) + args.sort(key=natural_keys) + kwargs.sort(key=natural_keys) + + def get(node,key): + if key in node.children: return node[key] + elif key in node.attrs: return node.attrs[key] + + ARGS = [] + for name in args: + val = get(node, name) + if val.__class__ is h5py.Dataset: val = val[:] + if val is 'None': val = None + if '[' in name: # We are reloading a list + ind = int(name[4:name.index('[')]) + if len(ARGS) is ind: # Create the list + ARGS.append([val]) + else: + ARGS[ind].append(val) + elif issubclass(val.__class__,hdf5Group): + ARGS.append(loadSavable(val,pointers=pointers)) + else: + ind = int(name[4:]) + ARGS.append(val) + + KWARGS = {} + for name in kwargs: + val = get(node, name) + if val.__class__ is h5py.Dataset: val = val[:] + if val is 'None': val = None + if '[' in name: # We are reloading a list + key = name[7:name.index('[')] + if key not in KWARGS: # Create the list + KWARGS[key] = [val] + else: + KWARGS[key].append(val) + elif issubclass(val.__class__,hdf5Group): + key = name[7:] + KWARGS[key] = loadSavable(val,pointers=pointers) + else: + key = name[7:] + KWARGS[key] = val + + cls = get(node, '__class__') + if cls in SAVEABLES: + try: + out = SAVEABLES[cls](*ARGS, **KWARGS) + out._savable = node + pointers.append(out) # Because this is recursive. + return out + except Exception, e: + print 'Warning: %s Class could not be initiated.' % cls + print 'ARGS: ', ARGS + print 'KWARGS: ', KWARGS + return (cls, ARGS, KWARGS, node) + else: + print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls + return (cls, ARGS, KWARGS, node) + diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index c5151837..e0f1e017 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -73,11 +73,9 @@ class Solver(object): Jacobi = sdiag(1.0/M[1].diagonal()) options['M'] = Jacobi elif M[0] is 'GS': - LL = sp.tril(M[1]) - UU = sp.triu(M[1]) DD = sdiag(M[1].diagonal()) - Uinv = Solver(UU, flag='U') - Linv = Solver(LL, flag='L') + Uinv = Solver(M[1], flag='U') + Linv = Solver(M[1], flag='L') def GS(f): return Uinv.solve(DD*Linv.solve(f)) options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype ) diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 622557d6..8d9b0504 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,13 +3,16 @@ import sputils import lomutils import interputils import ModelBuilder +import meshutils from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap -from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird +from meshutils import exampleLomGird, meshTensors +from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube from interputils import interpmat from ipythonUtils import easyAnimate as animate import Solver from Solver import Solver +import Save import Geophysics import types @@ -23,7 +26,7 @@ def hook(obj, method, name=None, overwrite=False, silent=False): If name is None, the name of the method is used. """ - if name is None: + if name is None: name = method.__name__ if name == '': raise Exception('Must provide name to hook lambda functions.') @@ -42,7 +45,7 @@ def setKwargs(obj, **kwargs): setattr(obj, attr, kwargs[attr]) else: raise Exception('%s attr is not recognized' % attr) - hook(obj,callHooks, silent=True) + hook(obj,hook, silent=True) hook(obj,setKwargs, silent=True) @@ -87,11 +90,48 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'): print pad + stopper['str'] % (l<=r,l,r) print pad + "%s%s%s" % ('-'*25,done,'-'*25) -def callHooks(obj, match, *args, **kwargs): - for method in [posible for posible in dir(obj) if ('_'+match) in posible]: - if getattr(obj,'debug',False): print (match+' is calling self.'+method) - getattr(obj,method)(*args, **kwargs) +def callHooks(match): + """ + Use this to wrap a funciton:: + @callHooks('doEndIteration') + def doEndIteration(self): + pass + + This will call everything named _doEndIteration* at the beginning of the function call. + """ + def callHooksWrap(f): + @wraps(f) + def wrapper(self,*args,**kwargs): + + for method in [posible for posible in dir(self) if ('_'+match) in posible]: + if getattr(self,'debug',False): print (match+' is calling self.'+method) + getattr(self,method)(*args, **kwargs) + + return f(self,*args,**kwargs) + + extra = """ + If you have things that also need to run in the method %s, you can create a method:: + + def _%s*(self, ... ): + pass + + Where the * can be any string. If present, _%s* will be called at the start of the default %s call. + You may also completely overwrite this function. + """ % (match, match, match, match) + + wrapper.__doc__ += extra + return wrapper + return callHooksWrap + +def dependentProperty(name, value, children, doc): + def fget(self): return getattr(self,name,value) + def fset(self, val): + for child in children: + if hasattr(self, child): + delattr(self, child) + setattr(self, name, val) + return property(fget=fget, fset=fset, doc=doc) class Counter(object): diff --git a/SimPEG/utils/lomutils.py b/SimPEG/utils/lomutils.py index 9020ce12..653cfdde 100644 --- a/SimPEG/utils/lomutils.py +++ b/SimPEG/utils/lomutils.py @@ -4,29 +4,6 @@ from matutils import mkvc, ndgrid, sub2ind from sputils import sdiag -def exampleLomGird(nC, exType): - assert type(nC) == list, "nC must be a list containing the number of nodes" - assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" - exType = exType.lower() - - possibleTypes = ['rect', 'rotate'] - assert exType in possibleTypes, "Not a possible example type." - - if exType == 'rect': - return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) - elif exType == 'rotate': - if len(nC) == 2: - X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) - amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2) - amt[amt < 0] = 0 - return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt - elif len(nC) == 3: - X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) - amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2) - amt[amt < 0] = 0 - return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt - - def volTetra(xyz, A, B, C, D): """ Returns the volume for tetrahedras volume specified by the indexes A to D. diff --git a/SimPEG/utils/meshutils.py b/SimPEG/utils/meshutils.py new file mode 100644 index 00000000..5c0a753c --- /dev/null +++ b/SimPEG/utils/meshutils.py @@ -0,0 +1,58 @@ +import numpy as np +from scipy import sparse as sp +from matutils import mkvc, ndgrid, sub2ind +from sputils import sdiag + +def exampleLomGird(nC, exType): + assert type(nC) == list, "nC must be a list containing the number of nodes" + assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" + exType = exType.lower() + + possibleTypes = ['rect', 'rotate'] + assert exType in possibleTypes, "Not a possible example type." + + if exType == 'rect': + return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + elif exType == 'rotate': + if len(nC) == 2: + X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2) + amt[amt < 0] = 0 + return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt + elif len(nC) == 3: + X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2) + amt[amt < 0] = 0 + return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt + + +def meshTensors(*args): + """ + **meshTensors** takes any number of tuples that have the form:: + + h1 = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCode), (numPad, sizeStart [, increaseFactor]) ) + + .. plot:: + + from SimPEG import mesh, utils + M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() + + """ + def padding(num, start, factor=1.3, reverse=False): + pad = ((np.ones(num)*factor)**np.arange(num))*start + if reverse: pad = pad[::-1] + return pad + tensors = tuple() + for i, arg in enumerate(args): + tensors += (np.r_[padding(*arg[0],reverse=True),np.ones(arg[1][0])*arg[1][1],padding(*arg[2])],) + + return list(tensors) if len(tensors) > 1 else tensors[0] + +if __name__ == '__main__': + from SimPEG import mesh + import matplotlib.pyplot as plt + M = mesh.TensorMesh(meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() + plt.gca().axis('tight') + plt.show() diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index 9853d089..c3d9f4a7 100644 --- a/SimPEG/visualize/vtk/vtkTools.py +++ b/SimPEG/visualize/vtk/vtkTools.py @@ -46,7 +46,7 @@ class vtkTools(object): vZ = mesh.vectorNz zD = mesh.nNz # Use rectilinear VTK grid. - # Asaign the spatial information. + # Assign the spatial information. vtkObj = vtk.vtkRectilinearGrid() vtkObj.SetDimensions(xD,yD,zD) vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) @@ -61,6 +61,7 @@ class vtkTools(object): vtkObj.GetCellData().AddArray(vtkDoubleArr) vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() return vtkObj @staticmethod @@ -104,7 +105,7 @@ class vtkTools(object): # Cells -cell array FCellArr = vtk.vtkCellArray() FCellArr.SetNumberOfCells(mesh.nF) - FCellArr.SetCells(mesh.nF*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + FCellArr.SetCells(mesh.nF,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) # Cell type FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1) # Cell location @@ -166,7 +167,7 @@ class vtkTools(object): # Cells -cell array ECellArr = vtk.vtkCellArray() ECellArr.SetNumberOfCells(mesh.nE) - ECellArr.SetCells(mesh.nE*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) + ECellArr.SetCells(mesh.nE,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) # Cell type ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1) # Cell location @@ -186,6 +187,7 @@ class vtkTools(object): vtkObj.GetCellData().AddArray(vtkDoubleArr) vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() return vtkObj @staticmethod @@ -255,8 +257,7 @@ class vtkTools(object): cellThres = vtk.vtkThreshold() cellThres.AllScalarsOn() cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(limits[0]) - cellThres.ThresholdByLower(limits[1]) + cellThres.ThresholdBetween(limits[0],limits[1]) cellThres.Update() return cellThres.GetOutput(), cellCore.GetOutput() @@ -271,8 +272,7 @@ class vtkTools(object): cellThres = vtk.vtkThreshold() cellThres.AllScalarsOn() cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(limits[0]) - cellThres.ThresholdByLower(limits[1]) + cellThres.ThresholdBetween(limits[0],limits[1]) cellThres.Update() return cellThres.GetOutput(), cellCore.GetOutput() diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index abedd0bf..3b52b992 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -1,6 +1,6 @@ -import numpy as np +import numpy as np, matplotlib as mpl try: - import vtk + import vtk, vtk.util.numpy_support as npsup #import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import except Exception, e: print 'VTK import error. Please ensure you have VTK installed to use this visualization package.' @@ -8,16 +8,16 @@ import SimPEG as simpeg class vtkView(object): """ - Class for storing and view of SimPEG models in VTK (visulization toolkit). + Class for storing and view of SimPEG models in VTK (visualization toolkit). Inputs: :param mesh, SimPEG mesh. :param propdict, dictionary of property models. Can have these dictionary names: - 'cell' - cell model; 'face' - face model; 'edge' - edge model - The dictionary properties are given as dictionaries with: + 'C' - cell model; 'F' - face model; 'E' - edge model; ('V' - vector field : NOT SUPPORTED) + The dictionary values are given as dictionaries with: {'NameOfThePropertyModel': np.array of the properties}. - The property array has to be ordered in compliance with SimPEG standards. + The property np.array has to be ordered in compliance with SimPEG standards. :: Example of usages. @@ -30,22 +30,7 @@ class vtkView(object): """ """ - # ToDo: Set the properties up so that there are set/get methods - self.name = 'VTK figure of SimPEG model' - self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] - self.limits = [0, 1e12] - self.viewprop = {'cell':0} # Name of the tyep and Int order of the array or name of the vector. - self._mesh = mesh - - - # Set vtk object containers - self._cell = None - self._faces = None - self._edges = None - - self._readPropertyDictionary(propdict) - - # Setup hidden properties + # Setup hidden properties, used for the visualization self._ren = None self._iren = None self._renwin = None @@ -56,6 +41,164 @@ class vtkView(object): self._widget = None self._actor = None self._lut = None + # Set vtk object containers + self._cells = None + self._faces = None + self._edges = None + self._vectors = None # Not implemented + # Set default values + self.name = 'VTK figure of SimPEG model' + + + + # Error check the input mesh + if type(mesh).__name__ != 'TensorMesh': + raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh)) + # Set the mesh + self._mesh = mesh + + # Read the property dictionary + self._readPropertyDictionary(propdict) + + + + + # Set/Get properties + @property + def cmap(self): + ''' Colormap to use in vtkView. Colormap is a matplotlib cmap(cm) array, has to be uint8(use flag bytes=True during cmap generation).''' + if getattr(self,'_cmap',None) is None: + # Set default + self._cmap = mpl.cm.hsv(np.arange(0.,1.,0.05),bytes=True) + return self._cmap + @cmap.setter + def cmap(self,value): + if value.min() > 0 or value.max() < 255 or value.shape[1] != 4 or value.dtype != np.uint8: + raise Exception('Input not an allowed array.\n Use matplotlib.cm to generate an array of size [nrColors,4] and dtype = uint8(flag bytes=True).') + self._cmap = value + + @property + def range(self): + ''' Range of the colors in vtkView.''' + if getattr(self,'_range',None) is None: + self._range = np.array(self._getActiveVTKobj().GetArray(self.viewprop.values()[0]).GetRange()) + return self._range + @range.setter + def range(self,value): + if type(value) not in [tuple, list, np.ndarray] or len(value) != 2 or np.array(value).dtype is not np.dtype('float'): + raise Exception('Input not in correct format. \n Has to be a list, tuple or np.arry of 2 floats.') + self._range = np.array(value) + + @property + def extent(self): + ''' Extent of the sub-domain of the model to view''' + if getattr(self,'_extent',None) is None: + self._extent = [0,self._mesh.nCx-1,0,self._mesh.nCy-1,0,self._mesh.nCz-1] + return self._extent + @extent.setter + def extent(self,value): + + import warnings + # Error check + valnp = np.array(value,dtype=int) + if valnp.dtype != int or len(valnp) != 6: + raise Exception('.extent has to be list or nparray of 6 integers.') + # Test the range of the values + loB = np.zeros(3,dtype=int) + upB = np.array(self._mesh.nCv - np.ones(3),dtype=int) + # Test the bounds + change = 0 + # Test for lower bounds, can't be smaller the 0 + tlb = valnp[::2] < loB + if tlb.any(): + valnp[::2][tlb] = loB[tlb] + change = 1 + warnings.warn('Lower bounds smaller then 0') + # Test for lower bounds, can't be larger then upB + tlub = valnp[::2] > upB + if tlub.any(): + valnp[::2][tlub] = upB[tlub] - 1 + change = 1 + warnings.warn('Lower bounds larger then uppermost bounds') + # Test for upper bounds, can't be larger the extent of the mesh + tub = valnp[1::2] > upB + if tub.any(): + valnp[1::2][tub] = upB[tub] + change = 1 + warnings.warn('Upper bounds greater then number of cells') + # Test if lower is smaller the upper + tgt = valnp[::2] > valnp[1::2] + if tgt.any(): + valnp[1::2][tgt] = valnp[::2][tgt] + 1 + change = 1 + warnings.warn('Lower bounds greater the Upper bounds') + # Print a warning + if change: + warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist())) + + # Set extent + self._extent = valnp + + @property + def limits(self): + ''' Lower and upper limits (cutoffs) of the values to view. ''' + return getattr(self,'_limits',None) + @limits.setter + def limits(self,value): + if value is None: + self._limits = None + else: + valnp = np.array(value) + if valnp.dtype != float or len(valnp) != 2: + raise Exception('.limits has to be list or numpy array of 2 floats.') + self._limits = valnp + + + @property + def viewprop(self): + ''' Controls the property that will be viewed.''' + + if getattr(self,'_viewprop',None) is None: + self._viewprop = {'C':0} # Name of the type and Int order of the array or name of the vector. + return self._viewprop + @viewprop.setter + def viewprop(self,value): + if type(value) != dict: + raise Exception('{:s} has to be a python dictionary containing property type and name index. ') + if len(value) > 1: + raise Exception('Too many input items in the viewprop dictionary') + if value.keys()[0] not in ['C','F','E']: + raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) + if not(type(self.viewprop.values()[0]) is int or type(self.viewprop.values()[0]) is str): + raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be integer or a string with the index.') + + + self._viewprop = value + + def _getActiveVTKobj(self): + """ + Finds the active VTK object. + """ + + if self.viewprop.keys()[0] is 'C': + vtkCellData = self._cells.GetCellData() + elif self.viewprop.keys()[0] is 'F': + vtkCellData = self._faces.GetCellData() + elif self.viewprop.keys()[0] is 'E': + vtkCellData = self._edges.GetCellData() + + return vtkCellData + + def _getActiveArrayName(self): + """ + Finds the name of the active array. + """ + actArr = self.viewprop.values()[0] + if type(actArr) is str: + activeName = actArr + elif type(actArr) is int: + activeName = self._getActiveVTKobj().GetArrayName(actArr) + return activeName def _readPropertyDictionary(self,propdict): """ @@ -64,18 +207,20 @@ class vtkView(object): import SimPEG.visualize.vtk.vtkTools as vtkSP # Test the property dictionary - if len(propdict) > 3: - raise(Exception,'Too many input items in the property dictionary') + if type(propdict) != dict: + raise Exception('{:s} has to be a python dictionary containing property models. ') + if len(propdict) > 4: + raise Exception('Too many input items in the property dictionary') for propitem in propdict.iteritems(): - if propitem[0] in ['cell','face','edge']: - if propitem[0] == 'cell': - self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'face': - self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'edge': - self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) + if propitem[0] in ['C','F','E']: + if propitem[0] == 'C': + self._cells = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) + if propitem[0] == 'F': + self._faces = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) + if propitem[0] == 'E': + self._edges = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) else: - raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) + raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) def Show(self): """ @@ -89,27 +234,34 @@ class vtkView(object): # Make renderwindow. Returns the interactor. self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) - imageType = self.viewprop.keys()[0] - # Sort out the actor - if imageType == 'cell': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) - elif imageType == 'face': - extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] - self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._face,extent,self.limits) - elif imageType == 'edge': - extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] - self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edge,extent,self.limits) - else: - raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) - + # Set the active scalar. if type(self.viewprop.values()[0]) == int: - actScalar = self._vtkobj.GetCellData().GetArrayName(self.viewprop.values()[0]) + actScalar = self._getActiveVTKobj().GetArrayName(self.viewprop.values()[0]) elif type(self.viewprop.values()[0]) == str: actScalar = self.viewprop.values()[0] else : raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be interger or a string.') - self._vtkobj.GetCellData().SetActiveScalars(actScalar) + self._getActiveVTKobj().SetActiveScalars(actScalar) + # Sort out the actor + imageType = self.viewprop.keys()[0] + if imageType == 'C': + if self.limits is None: + self.limits = self._cells.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cells,self.extent,self.limits) + elif imageType == 'F': + if self.limits is None: + self.limits = self._faces.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() + extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] + self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._faces,extent,self.limits) + elif imageType == 'E': + if self.limits is None: + self.limits = self._edges.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() + extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] + self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edges,extent,self.limits) + else: + raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType)) + #self._vtkobj.GetCellData().SetActiveScalars(actScalar) # Set up the plane, clipper and the user interaction. global intPlane, intActor self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) @@ -125,15 +277,27 @@ class vtkView(object): self._widget.AddObserver("InteractionEvent",movePlane) lut = vtk.vtkLookupTable() - lut.SetNumberOfColors(256) - lut.SetHueRange(0,0.66667) + lut.SetNumberOfColors(len(self.cmap)) + lut.SetTable(npsup.numpy_to_vtk(self.cmap)) lut.Build() self._lut = lut + scalarBar = vtk.vtkScalarBarActor() + scalarBar.SetLookupTable(lut) + scalarBar.SetTitle(self._getActiveArrayName()) + scalarBar.GetPositionCoordinate().SetCoordinateSystemToNormalizedViewport() + scalarBar.GetPositionCoordinate().SetValue(0.1,0.01) + scalarBar.SetOrientationToHorizontal() + scalarBar.SetWidth(0.8) + scalarBar.SetHeight(0.17) + + self._actor.GetMapper().SetScalarRange(self.range) self._actor.GetMapper().SetLookupTable(lut) # Set renderer options self._ren.SetBackground(.5,.5,.5) self._ren.AddActor(self._actor) + self._ren.AddActor2D(scalarBar) + self._renwin.SetSize(450,450) # Start the render Window vtkSP.startRenderWindow(self._iren) @@ -144,17 +308,43 @@ class vtkView(object): if __name__ == '__main__': + + #Make a mesh and model x0 = np.zeros(3) - h1 = np.ones(20)*50 - h2 = np.ones(10)*100 - h3 = np.ones(5)*200 + h1 = np.ones(60)*50 + h2 = np.ones(60)*100 + h3 = np.ones(50)*200 mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0) # Make a models that correspond to the cells, faces and edges. - models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} + t = np.ones(mesh.nC) + t[10000:50000] = 100 + t[100000:120000] = 100 + t[100000:120000] = 50 + models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t, 'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,mesh.nF),'AllOnce':np.ones(mesh.nF)},'E':{'Test':np.arange(0,mesh.nE),'AllOnce':np.ones(mesh.nE)}} # Make the vtk viewer object. vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) + # Set the .viewprop for which model to view + vtkViewer.viewprop = {'F':'Test'} # Show the image vtkViewer.Show() + + # Set subset of the mesh to view (remove padding) + vtkViewer.extent = [4,14,0,7,0,3] + vtkViewer.Show() + + # Change viewing property + vtkViewer.viewprop = {'C':'Model'} + # Set the color range + # Reset extent. + vtkViewer.extent = [-1,1000,-1,1000,-1,1000] + vtkViewer.range = [0.,100.] + vtkViewer.Show() + # Change color scale, has to be set to bytes=True. + vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True) + vtkViewer.Show() + # Set limits of values to view + vtkViewer.limits = [5.0,100.0] + vtkViewer.Show() \ No newline at end of file diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst index 8ca74e27..a30e55f7 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Optimize.rst @@ -24,3 +24,12 @@ Beta Schedule .. automodule:: SimPEG.inverse.BetaSchedule :members: :undoc-members: + + +Regularization +************** + +.. automodule:: SimPEG.inverse.Regularization + :show-inheritance: + :members: + :undoc-members: diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst index b4c96c4a..daa75454 100644 --- a/docs/api_Problem.rst +++ b/docs/api_Problem.rst @@ -14,7 +14,7 @@ Problem DCProblem ********* -.. automodule:: SimPEG.forward.DCProblem +.. automodule:: SimPEG.examples.DC :show-inheritance: :members: :undoc-members: @@ -25,7 +25,7 @@ DCProblem Linear Problem ************** -.. automodule:: SimPEG.forward.LinearProblem +.. automodule:: SimPEG.examples.Linear :show-inheritance: :members: :undoc-members: diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 07903e3c..7960b3e1 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -37,6 +37,13 @@ LOM Utilities :members: :undoc-members: +Mesh Utilities +************** + +.. automodule:: SimPEG.utils.meshutils + :members: + :undoc-members: + Model Builder Utilities *********************** diff --git a/docs/examples/mesh/plot_TensorMesh.py b/docs/examples/mesh/plot_TensorMesh.py deleted file mode 100644 index e773ed13..00000000 --- a/docs/examples/mesh/plot_TensorMesh.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from SimPEG.mesh import TensorMesh - -pad = 7 -padfactor = 1.4 -xpad = (np.ones(pad)*padfactor)**np.arange(pad) -ypad = (np.ones(pad)*padfactor)**np.arange(pad) - -core = 15 -xcore = np.ones(core) -ycore = np.ones(core) - -h1 = np.r_[xpad[::-1],xcore,xpad] -h2 = np.r_[ypad[::-1],ycore,ypad] - -mesh = TensorMesh([h1, h2]) -mesh.plotGrid() -plt.axis('tight') - -plt.show() diff --git a/docs/index.rst b/docs/index.rst index b008e963..b2b9806f 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,14 +3,23 @@ :alt: SimPEG :align: center -SimPEG (Simulation and Parameter Estimation in Geophysics) is a python package for simulation and gradient based parameter estimation in the context of geophysical applications. +SimPEG (Simulation and Parameter Estimation in Geoscience) is a python +package for simulation and gradient based parameter estimation in the +context of geoscience applications. -The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features: +The vision is to create a package for finite volume simulation and parameter estimation with +applications to geophysical imaging and subsurface flow. To enable +these goals, this package has the following features: -* modular with respect to the spacial discretization -* is built with the inverse problem in mind -* supports different hyperbolic solvers (Euler, Semi-Lagrangian, Lagrangian) -* supports 2D and 3D problems +* is modular with respect to discretization, physics, optimization, and regularization +* is built with the (large-scale) inverse problem in mind +* provides a framework for geophysical and hydrogeologic problems +* supports 1D, 2D and 3D problems + + +.. raw:: html + + Meshing & Operators diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb deleted file mode 100644 index 1ad88317..00000000 --- a/notebooks/3DRenderingWithvtkTools.ipynb +++ /dev/null @@ -1,70 +0,0 @@ -{ - "metadata": { - "name": "" - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ - { - "cells": [ - { - "cell_type": "code", - "collapsed": false, - "input": [ - "import numpy as np, vtk\n", - "import SimPEG as simpeg" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#Make a mesh and model\n", - "x0 = np.zeros(3)\n", - "h1 = np.ones(20)*50\n", - "h2 = np.ones(10)*100\n", - "h3 = np.ones(5)*200\n", - "\n", - "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 6 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "# Make a models that correspond to the cells, faces and edges.\n", - "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n", - "# Make the vtk viewer object.\n", - "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) \n", - " " - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 7 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "# Show the image \n", - "vtkViewer.Show()\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": "*" - } - ], - "metadata": {} - } - ] -} diff --git a/notebooks/VisualizeWithvtkView-updated.ipynb b/notebooks/VisualizeWithvtkView-updated.ipynb new file mode 100644 index 00000000..f4e8a3a7 --- /dev/null +++ b/notebooks/VisualizeWithvtkView-updated.ipynb @@ -0,0 +1,142 @@ +{ + "metadata": { + "name": "VisualizeWithvtkView-updated" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import SimPEG as simpeg, matplotlib as mpl" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The history saving thread hit an unexpected error (OperationalError('disk I/O error',)).History will not be written to the database.\n", + "Warning: mumps solver not available." + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Simple notebook of how to use vtkView to visualize SimPEG models. It will pop-up external vtk windows." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make a mesh and model\n", + "x0 = np.zeros(3)\n", + "h1 = np.ones(60)*50\n", + "h2 = np.ones(60)*100\n", + "h3 = np.ones(50)*200\n", + "\n", + "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "\n", + "# Make a models that correspond to the cells, faces and edges.\n", + "t = np.ones(mesh.nC)\n", + "t[10000:50000] = 100\n", + "t[100000:120000] = 100\n", + "t[100000:120000] = 50\n", + "# Make models called 'Test' for all with a range. \n", + "models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t},'F':{'Test':np.arange(0,mesh.nF)},'E':{'Test':np.arange(0,mesh.nE)}}\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make the vtk viewer object.\n", + "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)\n", + "# Set the .viewprop for which model to view\n", + "vtkViewer.viewprop = {'F':'Test'}\n", + "# Show the image\n", + "vtkViewer.Show()\n", + "\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Set subset of the mesh to view (remove padding)\n", + "vtkViewer.extent = [4,14,0,7,0,3]\n", + "vtkViewer.Show()\n", + "\n", + "# Change viewing property \n", + "vtkViewer.viewprop = {'C':'Model'}\n", + "# Set the color range\n", + "# Reset extent. Error check will reset the limits correctly.\n", + "vtkViewer.extent = [-1,1000,-1,1000,-1,1000]\n", + "# Set the range\n", + "vtkViewer.range = [0.,100.]\n", + "# Show\n", + "vtkViewer.Show()\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:116: UserWarning: Lower bounds smaller then 0\n", + " warnings.warn('Lower bounds smaller then 0')\n", + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:128: UserWarning: Upper bounds greater then number of cells\n", + " warnings.warn('Upper bounds greater then number of cells')\n", + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:137: UserWarning: Changed given extent from [-1, 1000, -1, 1000, -1, 1000] to [0, 59, 0, 59, 0, 49]\n", + " warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))\n" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Change color scale, has to be set to bytes=True.\n", + "vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)\n", + "vtkViewer.Show()\n", + "# Set limits of values to view \n", + "vtkViewer.limits = [5.0,100.0]\n", + "vtkViewer.Show()" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file