Merge branch 'develop' of https://github.com/simpeg/simpeg into hdf5

Conflicts:
	SimPEG/utils/__init__.py

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)

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)

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

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()

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):

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()
 

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()

notebooks/3DRenderingWithvtkTools.ipynb

@@ -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": {}
-  }
- ]
-}

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": {}
+  }
+ ]
+}