diff --git a/README.md b/README.md index 172cf62b..e107fea3 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,13 @@ -simpeg -====== +![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. +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: + +* modular with respect to the spacial discretization, optimization routine, and geophysical problem +* built with the inverse problem in mind +* provides a framework for geophysical and hydrogeologic problems +* supports 1D, 2D and 3D problems +* designed for large-scale inversions [![Build Status](https://travis-ci.org/simpeg/simpeg.png)](https://travis-ci.org/simpeg/simpeg) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 160678aa..a2ef2513 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -85,25 +85,19 @@ class BaseInversion(object): if self.stoppingCriteria(): break self.printDone() + self.finish() + return self.m + @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 +109,25 @@ class BaseInversion(object): self.phi_d_last = np.nan self.phi_m_last = np.nan + @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() + @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 @@ -213,6 +189,14 @@ class BaseInversion(object): """ printStoppers(self, self.stoppers) + @callHooks('finish') + def finish(self): + """finish() + + **finish** is called at the end of the optimization. + """ + pass + @timeIt def evalFunction(self, m, return_g=True, return_H=True): """evalFunction(m, return_g=True, return_H=True) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 7e934e13..80a03632 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -155,6 +155,7 @@ class Minimize(object): doEndIteration(xt) printDone() + finish() return xc """ self.evalFunction = evalFunction @@ -175,6 +176,7 @@ class Minimize(object): self.doEndIteration(xt) self.printDone() + self.finish() return self.xc @@ -188,6 +190,7 @@ class Minimize(object): def parent(self, value): self._parent = value + @callHooks('startup') def startup(self, x0): """ **startup** is called at the start of any new minimize call. @@ -198,19 +201,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 @@ -222,6 +216,7 @@ class Minimize(object): self.x_last = x0 @count + @callHooks('doStartIteration') def doStartIteration(self): """doStartIteration() @@ -230,7 +225,8 @@ class Minimize(object): :rtype: None :return: None """ - callHooks(self,'doStartIteration') + pass + def printInit(self, inLS=False): """ @@ -244,6 +240,7 @@ class Minimize(object): name = self.name if not inLS else self.nameLS printTitles(self, self.printers if not inLS else self.printersLS, name, pad) + @callHooks('printIter') def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -252,8 +249,6 @@ 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) @@ -270,6 +265,11 @@ class Minimize(object): stoppers = self.stoppers if not inLS else self.stoppersLS printStoppers(self, stoppers, pad='', stop=stop, done=done) + + def finish(self): + pass + + def stoppingCriteria(self, inLS=False): if self._iter == 0: self.f0 = self.f @@ -277,6 +277,7 @@ class Minimize(object): return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) @timeIt + @callHooks('projection') def projection(self, p): """projection(p) @@ -288,7 +289,6 @@ class Minimize(object): :rtype: numpy.ndarray :return: p, projected search direction """ - callHooks(self,'projection',p) return p @timeIt @@ -402,6 +402,7 @@ class Minimize(object): return p, False @count + @callHooks('doEndIteration') def doEndIteration(self, xt): """doEndIteration(xt) @@ -411,21 +412,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 @@ -630,7 +620,7 @@ 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): @@ -795,12 +785,23 @@ class NewtonRoot(object): 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 +813,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' # 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,10 +827,9 @@ 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 @@ -854,7 +853,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), 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/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/__init__.py b/SimPEG/utils/__init__.py index 05cf68a5..7762824c 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -88,20 +88,39 @@ 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): +def callHooks(match): """ - If you have things that also need to run at the end of every iteration, you can create a method:: + Use this to wrap a funciton:: - def _doEndIteration*(self, xt): - pass + @callHooks('doEndIteration') + 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. + This will call everything named _doEndIteration* at the beginning of the function call. """ - 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 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 class Counter(object): 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/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