From 6a101483a914e03c6c6ef34019ecb0be220ce1c9 Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Mon, 25 Nov 2013 17:22:46 -0800 Subject: [PATCH 01/35] Drasl --- SimPEG/InteractiveView.py | 335 ++++++++++++++++++++++ SimPEG/visulize/__init__.py | 2 + SimPEG/visulize/vtk/__init__.py | 2 + SimPEG/visulize/vtk/vtkTools.py | 361 ++++++++++++++++++++++++ SimPEG/visulize/vtk/vtkView.py | 137 +++++++++ notebooks/3DRenderingWithvtkTools.ipynb | 59 ++++ 6 files changed, 896 insertions(+) create mode 100644 SimPEG/InteractiveView.py create mode 100644 SimPEG/visulize/__init__.py create mode 100644 SimPEG/visulize/vtk/__init__.py create mode 100644 SimPEG/visulize/vtk/vtkTools.py create mode 100644 SimPEG/visulize/vtk/vtkView.py create mode 100644 notebooks/3DRenderingWithvtkTools.ipynb diff --git a/SimPEG/InteractiveView.py b/SimPEG/InteractiveView.py new file mode 100644 index 00000000..ba4f1b63 --- /dev/null +++ b/SimPEG/InteractiveView.py @@ -0,0 +1,335 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +from mpl_toolkits.mplot3d import Axes3D +from utils import mkvc + + +class TensorView(object): + """ + Provides viewing functions for TensorMesh + + This class is inherited by TensorMesh + """ + def __init__(self): + pass + + def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): + """ + Mesh.plotImage(I) + + Plots scalar fields on the given mesh. + + Input: + + :param numpy.array I: scalar field + + Optional Input: + + :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz + :param int figNum: number of figure to plot to + :param matplotlib.axes.Axes ax: axis to plot to + :param str direction: slice dimensions, 3D only ('x', 'y', 'z') + :param bool numbering: show numbering of slices, 3D only + :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' + :param bool showIt: call plt.show() + + .. plot:: examples/mesh/plot_image_2D.py + :include-source: + + .. plot:: examples/mesh/plot_image_3D.py + :include-source: + """ + assert type(I) == np.ndarray, "I must be a numpy array" + assert type(numbering) == bool, "numbering must be a bool" + assert direction in ["x", "y","z"], "direction must be either x,y, or z" + + + if imageType == 'CC': + assert I.size == self.nC, "Incorrect dimensions for CC." + elif imageType == 'N': + assert I.size == self.nN, "Incorrect dimensions for N." + elif imageType == 'Fx': + if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') + elif imageType == 'Fy': + if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') + elif imageType == 'Fz': + if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') + elif imageType == 'Ex': + if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') + elif imageType == 'Ey': + if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') + elif imageType == 'Ez': + if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') + elif imageType[0] == 'E': + plotAll = len(imageType) == 1 + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + fig = plt.figure(figNum) + # Determine the subplot number: 131, 121 + numPlots = 130 if plotAll else len(imageType)/2*10+100 + pltNum = 1 + ex, ey, ez = self.r(I,'E','E','M') + if plotAll or 'Ex' in imageType: + ax_x = plt.subplot(numPlots+pltNum) + self.plotImage(ex, imageType='Ex', ax=ax_x, **options) + pltNum +=1 + if plotAll or 'Ey' in imageType: + ax_y = plt.subplot(numPlots+pltNum) + self.plotImage(ey, imageType='Ey', ax=ax_y, **options) + pltNum +=1 + if plotAll or 'Ez' in imageType: + ax_z = plt.subplot(numPlots+pltNum) + self.plotImage(ez, imageType='Ez', ax=ax_z, **options) + pltNum +=1 + return + elif imageType[0] == 'F': + plotAll = len(imageType) == 1 + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + fig = plt.figure(figNum) + # Determine the subplot number: 131, 121 + numPlots = 130 if plotAll else len(imageType)/2*10+100 + pltNum = 1 + fx, fy, fz = self.r(I,'F','F','M') + if plotAll or 'Fx' in imageType: + ax_x = plt.subplot(numPlots+pltNum) + self.plotImage(fx, imageType='Fx', ax=ax_x, **options) + pltNum +=1 + if plotAll or 'Fy' in imageType: + ax_y = plt.subplot(numPlots+pltNum) + self.plotImage(fy, imageType='Fy', ax=ax_y, **options) + pltNum +=1 + if plotAll or 'Fz' in imageType: + ax_z = plt.subplot(numPlots+pltNum) + self.plotImage(fz, imageType='Fz', ax=ax_z, **options) + pltNum +=1 + return + else: + raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") + + + if ax is None: + fig = plt.figure(figNum) + fig.clf() + ax = plt.subplot(111) + else: + assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" + fig = ax.figure + + if self.dim == 1: + if imageType == 'CC': + ph = ax.plot(self.vectorCCx, I, '-ro') + elif imageType == 'N': + ph = ax.plot(self.vectorNx, I, '-bs') + ax.set_xlabel("x") + ax.axis('tight') + elif self.dim == 2: + if imageType == 'CC': + C = I[:].reshape(self.n, order='F') + elif imageType == 'N': + C = I[:].reshape(self.n+1, order='F') + C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) + elif imageType == 'Fx': + C = I[:].reshape(self.nFx, order='F') + C = 0.5*(C[:-1, :] + C[1:, :] ) + elif imageType == 'Fy': + C = I[:].reshape(self.nFy, order='F') + C = 0.5*(C[:, :-1] + C[:, 1:] ) + elif imageType == 'Ex': + C = I[:].reshape(self.nEx, order='F') + C = 0.5*(C[:,:-1] + C[:,1:] ) + elif imageType == 'Ey': + C = I[:].reshape(self.nEy, order='F') + C = 0.5*(C[:-1,:] + C[1:,:] ) + + ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) + ax.axis('tight') + ax.set_xlabel("x") + ax.set_ylabel("y") + + elif self.dim == 3: + if direction == 'z': + + # get copy of image and average to cell-centres is necessary + if imageType == 'CC': + Ic = I[:].reshape(self.n, order='F') + elif imageType == 'N': + Ic = I[:].reshape(self.n+1, order='F') + Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) + elif imageType == 'Fx': + Ic = I[:].reshape(self.nFx, order='F') + Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) + elif imageType == 'Fy': + Ic = I[:].reshape(self.nFy, order='F') + Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) + elif imageType == 'Fz': + Ic = I[:].reshape(self.nFz, order='F') + Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) + elif imageType == 'Ex': + Ic = I[:].reshape(self.nEx, order='F') + Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) + elif imageType == 'Ey': + Ic = I[:].reshape(self.nEy, order='F') + Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) + elif imageType == 'Ez': + Ic = I[:].reshape(self.nEz, order='F') + Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) + + # determine number oE slices in x and y dimension + nX = np.ceil(np.sqrt(self.nCz)) + nY = np.ceil(self.nCz/nX) + + # allocate space for montage + nCx = self.nCx + nCy = self.nCy + + C = np.zeros((nX*nCx,nY*nCy)) + + for iy in range(int(nY)): + for ix in range(int(nX)): + iz = ix + iy*nX + if iz < self.nCz: + C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] + else: + C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan + + C = np.ma.masked_where(np.isnan(C), C) + xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] + yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] + # Plot the mesh + ph = ax.pcolormesh(xx, yy, C.T) + # Plot the lines + gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) + gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) + # Repeat and seperate with NaN + gxX = np.c_[gx, gx, gx+np.nan].ravel() + gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() + gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() + gyY = np.c_[gy, gy, gy+np.nan].ravel() + ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) + ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) + ax.axis('tight') + + if numbering: + pad = np.sum(self.hx)*0.04 + for iy in range(int(nY)): + for ix in range(int(nX)): + iz = ix + iy*nX + if iz < self.nCz: + ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, + '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') + + ax.set_title(imageType) + if showIt: plt.show() + return ph + + def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. + + :param bool nodes: plot nodes + :param bool faces: plot faces + :param bool centers: plot centers + :param bool edges: plot edges + :param bool lines: plot lines connecting nodes + :param bool showIt: call plt.show() + + .. plot:: examples/mesh/plot_grid_2D.py + :include-source: + + .. plot:: examples/mesh/plot_grid_3D.py + :include-source: + """ + if self.dim == 1: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + ax.hold(True) + ax.plot(xn, np.ones(np.shape(xn)), 'bs') + ax.plot(xc, np.ones(np.shape(xc)), 'ro') + ax.plot(xn, np.ones(np.shape(xn)), 'k--') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + if showIt: plt.show() + elif self.dim == 2: + fig = plt.figure(2) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + xs1 = self.gridFx + xs2 = self.gridFy + + ax.hold(True) + if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') + if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') + if faces: + ax.plot(xs1[:, 0], xs1[:, 1], 'g>') + ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + + # Plot the grid lines + if lines: + NN = self.r(self.gridN, 'N', 'N', 'M') + X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() + X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] + plt.plot(X, Y) + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + if showIt: plt.show() + elif self.dim == 3: + fig = plt.figure(3) + fig.clf() + ax = fig.add_subplot(111, projection='3d') + xn = self.gridN + xc = self.gridCC + xfs1 = self.gridFx + xfs2 = self.gridFy + xfs3 = self.gridFz + + xes1 = self.gridEx + xes2 = self.gridEy + xes3 = self.gridEz + + ax.hold(True) + if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) + if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) + if faces: + ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) + ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) + ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) + if edges: + ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) + ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) + ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) + + # Plot the grid lines + if lines: + NN = self.r(self.gridN, 'N', 'N', 'M') + X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() + Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() + X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() + Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() + X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() + Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() + Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() + X = np.r_[X1, X2, X3] + Y = np.r_[Y1, Y2, Y3] + Z = np.r_[Z1, Z2, Z3] + plt.plot(X, Y, 'b-', zs=Z) + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + ax.set_zlabel('x3') + if showIt: plt.show() diff --git a/SimPEG/visulize/__init__.py b/SimPEG/visulize/__init__.py new file mode 100644 index 00000000..485f9782 --- /dev/null +++ b/SimPEG/visulize/__init__.py @@ -0,0 +1,2 @@ +import vtk +#import mpl diff --git a/SimPEG/visulize/vtk/__init__.py b/SimPEG/visulize/vtk/__init__.py new file mode 100644 index 00000000..6d60ed5c --- /dev/null +++ b/SimPEG/visulize/vtk/__init__.py @@ -0,0 +1,2 @@ +from vtkTools import vtkTools +from vtkView import vtkView \ No newline at end of file diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py new file mode 100644 index 00000000..50387a14 --- /dev/null +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -0,0 +1,361 @@ +import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb +from SimPEG.utils import mkvc + + +class vtkTools(object): + """ + Class that interacts with VTK visulization toolkit. + + """ + + def __init__(self): + """ Initializes the VTK vtkTools. + + """ + + pass + + @staticmethod + def makeCellVTKObject(mesh,model): + """ + Make and return a cell based VTK object for a simpeg mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells + + Output: + :rtype: vtkRecilinearGrid object + :return: vtkObj + """ + + # Deal with dimensionalities + if mesh.dim >= 1: + vX = mesh.vectorNx + xD = mesh.nNx + yD,zD = 1,1 + vY, vZ = np.array([0,0]) + if mesh.dim >= 2: + vY = mesh.vectorNy + yD = mesh.nNy + if mesh.dim == 3: + vZ = mesh.vectorNz + zD = mesh.nNz + # Use rectilinear VTK grid. + # Asaign the spatial information. + vtkObj = vtk.vtkRectilinearGrid() + vtkObj.SetDimensions(xD,yD,zD) + vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) + vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) + vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1)) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + return vtkObj + + @staticmethod + def makeFaceVTKObject(mesh,model): + """ + Make and return a face based VTK object for a simpeg mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). + Property array must be order hstack(Fx,Fy,Fz) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def faceR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTFx = np.prod(mesh.nFx) + FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) + FyCellBlock = np.array([],dtype='int64') + FzCellBlock = np.array([],dtype='int64') + # Second direction + if mesh.dim >= 2: + nTFy = np.prod(mesh.nFy) + FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) + # Third direction + if mesh.dim == 3: + nTFz = np.prod(mesh.nFz) + FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) + # Cells -cell array + FCellArr = vtk.vtkCellArray() + FCellArr.SetNumberOfCells(np.sum(mesh.nF)) + FCellArr.SetCells(np.sum(mesh.nF)*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + # Cell type + FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) + # Cell location + FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF)*5,5,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(FCellType,FCellLoc,FCellArr) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() + return vtkObj + + @staticmethod + def makeEdgeVTKObject(mesh,model): + """ + Make and return a edge based VTK object for a simpeg mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). + Property array must be order hstack(Ex,Ey,Ez) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def edgeR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTEx = np.prod(mesh.nEx) + ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)]) + # Second direction + if mesh.dim >= 2: + nTEy = np.prod(mesh.nEy) + EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)]) + # Third direction + if mesh.dim == 3: + nTEz = np.prod(mesh.nEz) + EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) + # Cells -cell array + ECellArr = vtk.vtkCellArray() + ECellArr.SetNumberOfCells(np.sum(mesh.nE)) + ECellArr.SetCells(np.sum(mesh.nE)*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) + # Cell type + ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(np.sum(mesh.nE),dtype='uint8'),deep=1) + # Cell location + ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nE)*3,3,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(ECellType,ECellLoc,ECellArr) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + return vtkObj + + @staticmethod + def makeRenderWindow(ren): + renwin = vtk.vtkRenderWindow() + renwin.AddRenderer(ren) + iren = vtk.vtkRenderWindowInteractor() + iren.SetRenderWindow(renwin) + + return iren, renwin + + + @staticmethod + def closeRenderWindow(iren): + renwin = iren.GetRenderWindow() + renwin.Finalize() + iren.TerminateApp() + + del iren, renwin + + @staticmethod + def makeVTKActor(vtkObj): + """ Makes a vtk mapper and Actor""" + mapper = vtk.vtkDataSetMapper() + mapper.SetInput(vtkObj) + actor = vtk.vtkActor() + actor.SetMapper(mapper) + actor.GetProperty().SetColor(0,0,0) + actor.GetProperty().SetRepresentationToWireframe() + return actor + + @staticmethod + def makeVTKLODActor(vtkObj,clipper): + """Make LOD vtk Actor""" + selectMapper = vtk.vtkDataSetMapper() + selectMapper.SetInputConnection(clipper.GetOutputPort()) + selectMapper.SetScalarVisibility(1) + selectMapper.SetColorModeToMapScalars() + selectMapper.SetScalarModeToUseCellData() + selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange()) + + selectActor = vtk.vtkLODActor() + selectActor.SetMapper(selectMapper) + selectActor.GetProperty().SetEdgeColor(1,0.5,0) + selectActor.GetProperty().SetEdgeVisibility(0) + selectActor.VisibilityOn() + selectActor.SetScale(1.01, 1.01, 1.01) + return selectActor + + @staticmethod + def setScalar2View(vtkObj,scalarName): + """ Sets the sclar to view """ + useArr = vtkObj.GetCellData().GetArray(scalarName) + if useArr == None: + raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName)) + vtkObj.GetCellData().SetActiveScalars(scalarName) + + @staticmethod + def makeRectiVTKVOIThres(vtkObj): + """Make volume of interest and threshold for rectilinear grid.""" + cellCore = vtk.vtkExtractRectilinearGrid() + cellCore.SetInput(vtkObj) + cellCore.SetVOI(vtkObj.GetExtent()) + cellThres = vtk.vtkThreshold() + cellThres.AllScalarsOn() + cellThres.SetInputConnection(cellCore.GetOutputPort()) + cellThres.ThresholdByUpper(-1) + cellThres.Update() + return cellThres.GetOutput(), cellCore.GetOutput() + + @staticmethod + def makePlaneClipper(vtkObj): + """Makes a plane and clipper """ + plane = vtk.vtkPlane() + clipper = vtk.vtkClipDataSet() + clipper.SetInputConnection(vtkObj.GetProducerPort()) + clipper.SetClipFunction(plane) + clipper.InsideOutOff() + return clipper, plane + + @staticmethod + def makePlaneWidget(vtkObj,iren,plane,actor): + """Make an interactive planeWidget""" + + # Callback function + def movePlane(obj, events): + obj.GetPlane(intPlane) + intActor.VisibilityOn() + + # Associate the line widget with the interactor + planeWidget = vtk.vtkImplicitPlaneWidget() + planeWidget.SetInteractor(iren) + planeWidget.SetPlaceFactor(1.25) + planeWidget.SetInput(vtkObj) + planeWidget.PlaceWidget() + #planeWidget.AddObserver("InteractionEvent", movePlane) + planeWidget.SetScaleEnabled(0) + planeWidget.SetEnabled(1) + planeWidget.SetOutlineTranslation(0) + planeWidget.GetPlaneProperty().SetOpacity(0.1) + return planeWidget + + + @staticmethod + def startRenderWindow(iren): + """ Start a vtk rendering window""" + iren.Initialize() + renwin = iren.GetRenderWindow() + renwin.Render() + iren.Start() + + + # Simple write/read VTK xml model functions. + @staticmethod + def writeVTPFile(fileName,vtkPolyObject): + '''Function to write vtk polydata file (vtp).''' + polyWriter = vtk.vtkXMLPolyDataWriter() + polyWriter.SetInput(vtkPolyObject) + polyWriter.SetFileName(fileName) + polyWriter.Update() + + @staticmethod + def writeVTUFile(fileName,vtkUnstructuredGrid): + '''Function to write vtk unstructured grid (vtu).''' + Writer = vtk.vtkXMLUnstructuredGridWriter() + Writer.SetInput(vtkUnstructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTRFile(fileName,vtkRectilinearGrid): + '''Function to write vtk rectilinear grid (vtr).''' + Writer = vtk.vtkXMLRectilinearGridWriter() + Writer.SetInput(vtkRectilinearGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTSFile(fileName,vtkStructuredGrid): + '''Function to write vtk structured grid (vts).''' + Writer = vtk.vtkXMLStructuredGridWriter() + Writer.SetInput(vtkStructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def readVTSFile(fileName): + '''Function to read vtk structured grid (vts) and return a grid object.''' + Reader = vtk.vtkXMLStructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTUFile(fileName): + '''Function to read vtk structured grid (vtu) and return a grid object.''' + Reader = vtk.vtkXMLUnstructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTRFile(fileName): + '''Function to read vtk structured grid (vtr) and return a grid object.''' + Reader = vtk.vtkXMLRectilinearGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTPFile(fileName): + '''Function to read vtk structured grid (vtp) and return a grid object.''' + Reader = vtk.vtkXMLPolyDataReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py new file mode 100644 index 00000000..334561ba --- /dev/null +++ b/SimPEG/visulize/vtk/vtkView.py @@ -0,0 +1,137 @@ +import numpy as np, vtk +import SimPEG as simpeg +#import SimPEG.visulize.vtk.vtkTools as vtkSP # Always get an error for this import + +class vtkView(object): + """ + Class for storing and view of SimPEG models in VTK (visulization 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: + {'NameOfThePropertyModel': np.array of the properties}. + The property array has to be ordered in compliance with SimPEG standards. + + :: + Example of usages. + + ToDo + + """ + + def __init__(self,mesh,propdict): + """ + """ + + self.name = 'VTK figure of SimPEG model' + self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] + self.limits = [0, 10000] + self._mesh = mesh + + # Set vtk object containers + self._cell = None + self._faces = None + self._edges = None + + self._readPropertyDictionary(propdict) + + # Setup hidden properties + self._ren = None + self._iren = None + self._renwin = None + self._core = None + self._viewobj = None + self._plane = None + self._clipper = None + self._widget = None + self._actor = None + self._lut = None + + def _readPropertyDictionary(self,propdict): + """ + Reads the property and assigns to the object + """ + import SimPEG.visulize.vtk.vtkTools as vtkSP + + # Test the property dictionary + if len(propdict) > 3: + 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]) + else: + raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) + + def Show(self,imageType='cell'): + """ + Open the VTK figure window and show the mesh. + + Inputs: + param: str imageType: type of image {'cell','face','edge'} + + """ + #vtkSP = simpeg.visulize.vtk.vtkTools + import SimPEG.visulize.vtk.vtkTools as vtkSP + + # Make a renderer + self._ren = vtk.vtkRenderer() + # Make renderwindow. Returns the interactor. + self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) + + + # Sort out the actor + if imageType == 'cell': + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell) + elif imageType == 'face': + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._face) + elif imageType == 'edge': + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._edge) + else: + raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) + + + global intPlane, intActor + self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) + intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper) + self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor) + # Callback function + self._plane = intPlane + self._actor = intActor + def movePlane(obj, events): + global intPlane, intActor + obj.GetPlane(intPlane) + intActor.VisibilityOn() + + self._widget.AddObserver("InteractionEvent",movePlane) + lut = vtk.vtkLookupTable() + lut.SetNumberOfColors(256) + lut.SetHueRange(0,0.66667) + lut.Build() + self._lut = lut + self._actor.GetMapper().SetLookupTable(lut) + + # Set renderer options + self._ren.SetBackground(.5,.5,.5) + self._ren.AddActor(self._actor) + + # Start the render Window + vtkSP.startRenderWindow(self._iren) + # Close the window when exited + vtkSP.closeRenderWindow(self._iren) + del self._iren, self._renwin + + + + + + + + diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb new file mode 100644 index 00000000..7209b10e --- /dev/null +++ b/notebooks/3DRenderingWithvtkTools.ipynb @@ -0,0 +1,59 @@ +{ + "metadata": { + "name": "3D rendering with vtkTools" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": "import numpy as np, vtk\nimport SimPEG as simpeg", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'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.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + } + ], + "metadata": {} + } + ] +} From b29cc6e08d24aa3aea0bac455371e8436cfcc590 Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Mon, 25 Nov 2013 17:44:26 -0800 Subject: [PATCH 02/35] Deleted local 'visulize' folder --- SimPEG/visulize/__init__.py | 2 - SimPEG/visulize/vtk/__init__.py | 2 - SimPEG/visulize/vtk/vtkTools.py | 361 ------------------- SimPEG/visulize/vtk/vtkView.py | 137 ------- notebooks/3DRenderingWithvtkTools.ipynb.orig | 115 ++++++ 5 files changed, 115 insertions(+), 502 deletions(-) delete mode 100644 SimPEG/visulize/__init__.py delete mode 100644 SimPEG/visulize/vtk/__init__.py delete mode 100644 SimPEG/visulize/vtk/vtkTools.py delete mode 100644 SimPEG/visulize/vtk/vtkView.py create mode 100644 notebooks/3DRenderingWithvtkTools.ipynb.orig diff --git a/SimPEG/visulize/__init__.py b/SimPEG/visulize/__init__.py deleted file mode 100644 index 485f9782..00000000 --- a/SimPEG/visulize/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -import vtk -#import mpl diff --git a/SimPEG/visulize/vtk/__init__.py b/SimPEG/visulize/vtk/__init__.py deleted file mode 100644 index 6d60ed5c..00000000 --- a/SimPEG/visulize/vtk/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -from vtkTools import vtkTools -from vtkView import vtkView \ No newline at end of file diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py deleted file mode 100644 index 50387a14..00000000 --- a/SimPEG/visulize/vtk/vtkTools.py +++ /dev/null @@ -1,361 +0,0 @@ -import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb -from SimPEG.utils import mkvc - - -class vtkTools(object): - """ - Class that interacts with VTK visulization toolkit. - - """ - - def __init__(self): - """ Initializes the VTK vtkTools. - - """ - - pass - - @staticmethod - def makeCellVTKObject(mesh,model): - """ - Make and return a cell based VTK object for a simpeg mesh and model. - - Input: - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells - - Output: - :rtype: vtkRecilinearGrid object - :return: vtkObj - """ - - # Deal with dimensionalities - if mesh.dim >= 1: - vX = mesh.vectorNx - xD = mesh.nNx - yD,zD = 1,1 - vY, vZ = np.array([0,0]) - if mesh.dim >= 2: - vY = mesh.vectorNy - yD = mesh.nNy - if mesh.dim == 3: - vZ = mesh.vectorNz - zD = mesh.nNz - # Use rectilinear VTK grid. - # Asaign the spatial information. - vtkObj = vtk.vtkRectilinearGrid() - vtkObj.SetDimensions(xD,yD,zD) - vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) - vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) - vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1)) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - return vtkObj - - @staticmethod - def makeFaceVTKObject(mesh,model): - """ - Make and return a face based VTK object for a simpeg mesh and model. - - Input: - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). - Property array must be order hstack(Fx,Fy,Fz) - - Output: - :rtype: vtkUnstructuredGrid object - :return: vtkObj - """ - - ## Convert simpeg mesh to VTK properties - # Convert mesh nodes to vtkPoints - vtkPts = vtk.vtkPoints() - vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - - # Define the face "cells" - # Using VTK_QUAD cell for faces (see VTK file format) - nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') - def faceR(mat,length): - return mat.T.reshape((length,1)) - # First direction - nTFx = np.prod(mesh.nFx) - FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) - FyCellBlock = np.array([],dtype='int64') - FzCellBlock = np.array([],dtype='int64') - # Second direction - if mesh.dim >= 2: - nTFy = np.prod(mesh.nFy) - FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) - # Third direction - if mesh.dim == 3: - nTFz = np.prod(mesh.nFz) - FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) - # Cells -cell array - FCellArr = vtk.vtkCellArray() - FCellArr.SetNumberOfCells(np.sum(mesh.nF)) - FCellArr.SetCells(np.sum(mesh.nF)*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) - # Cell type - FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) - # Cell location - FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF)*5,5,dtype='int64'),deep=1) - - ## Make the object - vtkObj = vtk.vtkUnstructuredGrid() - # Set the objects properties - vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(FCellType,FCellLoc,FCellArr) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - vtkObj.Update() - return vtkObj - - @staticmethod - def makeEdgeVTKObject(mesh,model): - """ - Make and return a edge based VTK object for a simpeg mesh and model. - - Input: - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). - Property array must be order hstack(Ex,Ey,Ez) - - Output: - :rtype: vtkUnstructuredGrid object - :return: vtkObj - """ - - ## Convert simpeg mesh to VTK properties - # Convert mesh nodes to vtkPoints - vtkPts = vtk.vtkPoints() - vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - - # Define the face "cells" - # Using VTK_QUAD cell for faces (see VTK file format) - nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') - def edgeR(mat,length): - return mat.T.reshape((length,1)) - # First direction - nTEx = np.prod(mesh.nEx) - ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)]) - # Second direction - if mesh.dim >= 2: - nTEy = np.prod(mesh.nEy) - EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)]) - # Third direction - if mesh.dim == 3: - nTEz = np.prod(mesh.nEz) - EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) - # Cells -cell array - ECellArr = vtk.vtkCellArray() - ECellArr.SetNumberOfCells(np.sum(mesh.nE)) - ECellArr.SetCells(np.sum(mesh.nE)*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) - # Cell type - ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(np.sum(mesh.nE),dtype='uint8'),deep=1) - # Cell location - ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nE)*3,3,dtype='int64'),deep=1) - - ## Make the object - vtkObj = vtk.vtkUnstructuredGrid() - # Set the objects properties - vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(ECellType,ECellLoc,ECellArr) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - return vtkObj - - @staticmethod - def makeRenderWindow(ren): - renwin = vtk.vtkRenderWindow() - renwin.AddRenderer(ren) - iren = vtk.vtkRenderWindowInteractor() - iren.SetRenderWindow(renwin) - - return iren, renwin - - - @staticmethod - def closeRenderWindow(iren): - renwin = iren.GetRenderWindow() - renwin.Finalize() - iren.TerminateApp() - - del iren, renwin - - @staticmethod - def makeVTKActor(vtkObj): - """ Makes a vtk mapper and Actor""" - mapper = vtk.vtkDataSetMapper() - mapper.SetInput(vtkObj) - actor = vtk.vtkActor() - actor.SetMapper(mapper) - actor.GetProperty().SetColor(0,0,0) - actor.GetProperty().SetRepresentationToWireframe() - return actor - - @staticmethod - def makeVTKLODActor(vtkObj,clipper): - """Make LOD vtk Actor""" - selectMapper = vtk.vtkDataSetMapper() - selectMapper.SetInputConnection(clipper.GetOutputPort()) - selectMapper.SetScalarVisibility(1) - selectMapper.SetColorModeToMapScalars() - selectMapper.SetScalarModeToUseCellData() - selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange()) - - selectActor = vtk.vtkLODActor() - selectActor.SetMapper(selectMapper) - selectActor.GetProperty().SetEdgeColor(1,0.5,0) - selectActor.GetProperty().SetEdgeVisibility(0) - selectActor.VisibilityOn() - selectActor.SetScale(1.01, 1.01, 1.01) - return selectActor - - @staticmethod - def setScalar2View(vtkObj,scalarName): - """ Sets the sclar to view """ - useArr = vtkObj.GetCellData().GetArray(scalarName) - if useArr == None: - raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName)) - vtkObj.GetCellData().SetActiveScalars(scalarName) - - @staticmethod - def makeRectiVTKVOIThres(vtkObj): - """Make volume of interest and threshold for rectilinear grid.""" - cellCore = vtk.vtkExtractRectilinearGrid() - cellCore.SetInput(vtkObj) - cellCore.SetVOI(vtkObj.GetExtent()) - cellThres = vtk.vtkThreshold() - cellThres.AllScalarsOn() - cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(-1) - cellThres.Update() - return cellThres.GetOutput(), cellCore.GetOutput() - - @staticmethod - def makePlaneClipper(vtkObj): - """Makes a plane and clipper """ - plane = vtk.vtkPlane() - clipper = vtk.vtkClipDataSet() - clipper.SetInputConnection(vtkObj.GetProducerPort()) - clipper.SetClipFunction(plane) - clipper.InsideOutOff() - return clipper, plane - - @staticmethod - def makePlaneWidget(vtkObj,iren,plane,actor): - """Make an interactive planeWidget""" - - # Callback function - def movePlane(obj, events): - obj.GetPlane(intPlane) - intActor.VisibilityOn() - - # Associate the line widget with the interactor - planeWidget = vtk.vtkImplicitPlaneWidget() - planeWidget.SetInteractor(iren) - planeWidget.SetPlaceFactor(1.25) - planeWidget.SetInput(vtkObj) - planeWidget.PlaceWidget() - #planeWidget.AddObserver("InteractionEvent", movePlane) - planeWidget.SetScaleEnabled(0) - planeWidget.SetEnabled(1) - planeWidget.SetOutlineTranslation(0) - planeWidget.GetPlaneProperty().SetOpacity(0.1) - return planeWidget - - - @staticmethod - def startRenderWindow(iren): - """ Start a vtk rendering window""" - iren.Initialize() - renwin = iren.GetRenderWindow() - renwin.Render() - iren.Start() - - - # Simple write/read VTK xml model functions. - @staticmethod - def writeVTPFile(fileName,vtkPolyObject): - '''Function to write vtk polydata file (vtp).''' - polyWriter = vtk.vtkXMLPolyDataWriter() - polyWriter.SetInput(vtkPolyObject) - polyWriter.SetFileName(fileName) - polyWriter.Update() - - @staticmethod - def writeVTUFile(fileName,vtkUnstructuredGrid): - '''Function to write vtk unstructured grid (vtu).''' - Writer = vtk.vtkXMLUnstructuredGridWriter() - Writer.SetInput(vtkUnstructuredGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def writeVTRFile(fileName,vtkRectilinearGrid): - '''Function to write vtk rectilinear grid (vtr).''' - Writer = vtk.vtkXMLRectilinearGridWriter() - Writer.SetInput(vtkRectilinearGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def writeVTSFile(fileName,vtkStructuredGrid): - '''Function to write vtk structured grid (vts).''' - Writer = vtk.vtkXMLStructuredGridWriter() - Writer.SetInput(vtkStructuredGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def readVTSFile(fileName): - '''Function to read vtk structured grid (vts) and return a grid object.''' - Reader = vtk.vtkXMLStructuredGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTUFile(fileName): - '''Function to read vtk structured grid (vtu) and return a grid object.''' - Reader = vtk.vtkXMLUnstructuredGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTRFile(fileName): - '''Function to read vtk structured grid (vtr) and return a grid object.''' - Reader = vtk.vtkXMLRectilinearGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTPFile(fileName): - '''Function to read vtk structured grid (vtp) and return a grid object.''' - Reader = vtk.vtkXMLPolyDataReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py deleted file mode 100644 index 334561ba..00000000 --- a/SimPEG/visulize/vtk/vtkView.py +++ /dev/null @@ -1,137 +0,0 @@ -import numpy as np, vtk -import SimPEG as simpeg -#import SimPEG.visulize.vtk.vtkTools as vtkSP # Always get an error for this import - -class vtkView(object): - """ - Class for storing and view of SimPEG models in VTK (visulization 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: - {'NameOfThePropertyModel': np.array of the properties}. - The property array has to be ordered in compliance with SimPEG standards. - - :: - Example of usages. - - ToDo - - """ - - def __init__(self,mesh,propdict): - """ - """ - - self.name = 'VTK figure of SimPEG model' - self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] - self.limits = [0, 10000] - self._mesh = mesh - - # Set vtk object containers - self._cell = None - self._faces = None - self._edges = None - - self._readPropertyDictionary(propdict) - - # Setup hidden properties - self._ren = None - self._iren = None - self._renwin = None - self._core = None - self._viewobj = None - self._plane = None - self._clipper = None - self._widget = None - self._actor = None - self._lut = None - - def _readPropertyDictionary(self,propdict): - """ - Reads the property and assigns to the object - """ - import SimPEG.visulize.vtk.vtkTools as vtkSP - - # Test the property dictionary - if len(propdict) > 3: - 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]) - else: - raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) - - def Show(self,imageType='cell'): - """ - Open the VTK figure window and show the mesh. - - Inputs: - param: str imageType: type of image {'cell','face','edge'} - - """ - #vtkSP = simpeg.visulize.vtk.vtkTools - import SimPEG.visulize.vtk.vtkTools as vtkSP - - # Make a renderer - self._ren = vtk.vtkRenderer() - # Make renderwindow. Returns the interactor. - self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) - - - # Sort out the actor - if imageType == 'cell': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell) - elif imageType == 'face': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._face) - elif imageType == 'edge': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._edge) - else: - raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) - - - global intPlane, intActor - self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) - intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper) - self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor) - # Callback function - self._plane = intPlane - self._actor = intActor - def movePlane(obj, events): - global intPlane, intActor - obj.GetPlane(intPlane) - intActor.VisibilityOn() - - self._widget.AddObserver("InteractionEvent",movePlane) - lut = vtk.vtkLookupTable() - lut.SetNumberOfColors(256) - lut.SetHueRange(0,0.66667) - lut.Build() - self._lut = lut - self._actor.GetMapper().SetLookupTable(lut) - - # Set renderer options - self._ren.SetBackground(.5,.5,.5) - self._ren.AddActor(self._actor) - - # Start the render Window - vtkSP.startRenderWindow(self._iren) - # Close the window when exited - vtkSP.closeRenderWindow(self._iren) - del self._iren, self._renwin - - - - - - - - diff --git a/notebooks/3DRenderingWithvtkTools.ipynb.orig b/notebooks/3DRenderingWithvtkTools.ipynb.orig new file mode 100644 index 00000000..8d2c833d --- /dev/null +++ b/notebooks/3DRenderingWithvtkTools.ipynb.orig @@ -0,0 +1,115 @@ +{ + "metadata": { +<<<<<<< HEAD + "name": "3D rendering with vtkTools" +======= + "name": "" +>>>>>>> develop + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "import numpy as np, vtk\nimport SimPEG as simpeg", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 +======= + "input": [ + "import numpy as np, vtk\n", + "import SimPEG as simpeg" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 +>>>>>>> develop + }, + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 +======= + "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 +>>>>>>> develop + }, + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'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.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 +======= + "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 +>>>>>>> develop + }, + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 +======= + "input": [ + "# Show the image \n", + "vtkViewer.Show()\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": "*" +>>>>>>> develop + } + ], + "metadata": {} + } + ] +} From 4b292b574d68b2d4a63dd477b500743de95e9feb Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Tue, 26 Nov 2013 07:45:18 -0800 Subject: [PATCH 03/35] Updating vtkView --- SimPEG/visualize/vtk/vtkView.py | 60 ++++++++++++++++++++++++--------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index abedd0bf..f65c8e7b 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -14,7 +14,7 @@ class vtkView(object): :param mesh, SimPEG mesh. :param propdict, dictionary of property models. Can have these dictionary names: - 'cell' - cell model; 'face' - face model; 'edge' - edge model + 'C' - cell model; 'F' - face model; 'E' - edge model The dictionary properties are given as dictionaries with: {'NameOfThePropertyModel': np.array of the properties}. The property array has to be ordered in compliance with SimPEG standards. @@ -30,11 +30,16 @@ class vtkView(object): """ """ - # ToDo: Set the properties up so that there are set/get methods + # Error check inputs + if type(mesh).__name__ != 'TensorMesh': + raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh)) + + + # Set default values 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._extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] + self._limits = [0, 1e12] + self._viewprop = {'C':0} # Name of the tyep and Int order of the array or name of the vector. self._mesh = mesh @@ -57,6 +62,31 @@ class vtkView(object): self._actor = None self._lut = None + # Set/Get properties + @property + def extent(self): + return getattr(self,_extent) + + @extent.setter + def extent(self,value): + + # Error check + valnp = np.array(value) + if valnp.dtype != int or len(valnp) != 6: + raise Exception('.extent has to be list or nparray of 6 integers.') + + self._extent = valnp + + @property + def limits(self): + return getattr(self._limits) + + @limits.setter + def limits(self,value): + valnp = np.array(value) + if valnp.dtype != float or len(valnp) != 2: + raise Exception('.limits has to be list or nparray of 2 floats.') + def _readPropertyDictionary(self,propdict): """ Reads the property and assigns to the object @@ -67,15 +97,15 @@ class vtkView(object): if len(propdict) > 3: 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': + if propitem[0] in ['C','F','E']: + if propitem[0] == 'C': self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'face': + if propitem[0] == 'F': self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'edge': + if propitem[0] == 'E': self._edge = 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 dictonary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) def Show(self): """ @@ -91,16 +121,16 @@ class vtkView(object): imageType = self.viewprop.keys()[0] # Sort out the actor - if imageType == 'cell': + if imageType == 'C': self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) - elif imageType == 'face': + elif imageType == 'F': 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': + elif imageType == 'E': 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)) + raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType)) # Set the active scalar. if type(self.viewprop.values()[0]) == int: @@ -153,7 +183,7 @@ if __name__ == '__main__': 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))}} + models = {'C':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'E':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} # Make the vtk viewer object. vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) # Show the image From 1ec3c5cd0936d0b014571fdee1aca78c291dc6c0 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Tue, 26 Nov 2013 21:07:30 -0800 Subject: [PATCH 04/35] Fixed vtkTools and vtkView to work for all cell,faces, edges. Changed index for cells, faces and edges to be C,F,E respectively. Made viewprop, limits, extent as @properties, with error checks in the set methods. Need to finish the limits error checks. --- SimPEG/visualize/vtk/vtkTools.py | 6 +- SimPEG/visualize/vtk/vtkView.py | 172 +++++++++++++++++++++++-------- 2 files changed, 134 insertions(+), 44 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index 9853d089..7b405d4c 100644 --- a/SimPEG/visualize/vtk/vtkTools.py +++ b/SimPEG/visualize/vtk/vtkTools.py @@ -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 diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index f65c8e7b..dde60fd3 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -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: - 'C' - cell model; 'F' - face model; 'E' - 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,27 +30,7 @@ class vtkView(object): """ """ - # Error check inputs - if type(mesh).__name__ != 'TensorMesh': - raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh)) - - - # Set default values - 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 = {'C':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 @@ -61,31 +41,122 @@ 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 extent(self): - return getattr(self,_extent) + ''' 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) + 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 = np.prod(valnp[::2] < loB,axis=0,dtype=bool) + 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 + #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) + 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 = np.prod(valnp[1::2] > upB,axis=0,dtype=bool) + 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 = np.prod(valnp[::2] < valnp[1::2],axis=0,dtype=bool) + 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): - return getattr(self._limits) + ''' Lower and upper limits (cutoffs) of the values to view. ''' + return getattr(self,'_limits',None) + @limits.setter def limits(self,value): - valnp = np.array(value) - if valnp.dtype != float or len(valnp) != 2: - raise Exception('.limits has to be list or nparray of 2 floats.') + 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 _readPropertyDictionary(self,propdict): """ @@ -94,18 +165,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 ['C','F','E']: if propitem[0] == 'C': - self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) + self._cells = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) if propitem[0] == 'F': - self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) + self._faces = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) if propitem[0] == 'E': - self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) + self._edges = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) else: - raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) + raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) def Show(self): """ @@ -120,15 +193,22 @@ class vtkView(object): self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) imageType = self.viewprop.keys()[0] + # Sort out the actor if imageType == 'C': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) + 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._face,extent,self.limits) + 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._edge,extent,self.limits) + 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)) @@ -174,6 +254,8 @@ class vtkView(object): if __name__ == '__main__': + + #Make a mesh and model x0 = np.zeros(3) h1 = np.ones(20)*50 @@ -183,8 +265,14 @@ if __name__ == '__main__': mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0) # Make a models that correspond to the cells, faces and edges. - models = {'C':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'E':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} + models = {'C':{'Test':np.arange(0,mesh.nC),'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() From 6603fe77e132a2d84200531c164a8215635c19e8 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 26 Nov 2013 21:35:03 -0800 Subject: [PATCH 05/35] First attempt on saving to a database, have written a very light wrapper on h5py. --- SimPEG/utils/Save.py | 155 +++++++++++++++++++++++++++++++++++++++ SimPEG/utils/__init__.py | 12 ++- 2 files changed, 166 insertions(+), 1 deletion(-) create mode 100644 SimPEG/utils/Save.py diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py new file mode 100644 index 00000000..9f661b0d --- /dev/null +++ b/SimPEG/utils/Save.py @@ -0,0 +1,155 @@ +import numpy as np +import time +import h5py +import re + +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) ] + + +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)) + +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 + + +class hdf5InversionGroup(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.childClass = hdf5Inversion + + def saveInversion(self, invObj, dataPath): + invObj._invNode = self.addGroup('%d'%self.numChildren) + + # At the start of every iteration we will create a inversion iteration node. + def _doStartIteration_hdf5_inv(invObj): + invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) + invNodeIt.attrs['complete'] = False + invNodeIt.attrs['time'] = time.time() + invObj._invNodeIt = invNodeIt + invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) + + def _doEndIteration_hdf5_inv(invObj): + invNodeIt = invObj._invNodeIt + invNodeIt.attrs['time'] = time.time() - invNodeIt.attrs['time'] + invNodeIt.attrs['ctime'] = time.ctime() + invNodeIt.attrs['phi_d'] = invObj.phi_d + invNodeIt.attrs['phi_m'] = invObj.phi_m + invNodeIt.setArray('m', invObj.m) + invNodeIt.setArray('dpred', invObj.dpred) + invNodeIt.attrs['complete'] = True + invObj.hook(_doEndIteration_hdf5_inv, overwrite=True) + + def _doStartIteration_hdf5_opt(optObj): + optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) + optNodeIt.attrs['complete'] = False + optNodeIt.attrs['time'] = time.time() + optObj._optNodeIt = optNodeIt + + invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True) + + def _doEndIteration_hdf5_opt(optObj, xt): + optNodeIt = optObj._optNodeIt + optNodeIt.attrs['time'] = time.time() - optNodeIt.attrs['time'] + optNodeIt.attrs['ctime'] = time.ctime() + optNodeIt.setArray('m', xt) + optNodeIt.setArray('dpred', optObj.parent.dpred) + optNodeIt.setArray('searchDirection', optObj.searchDirection) + optNodeIt.attrs['complete'] = True + invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True) + + return invObj._invNode + + +class hdf5Inversion(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.parentClass = hdf5InversionGroup + self.childClass = hdf5InversionIteration + + +class hdf5InversionIteration(hdf5Group): + def __init__(self, T, groupNode): + hdf5Group.__init__(self, T, groupNode) + self.parentClass = hdf5Inversion diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 622557d6..05cf68a5 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -10,6 +10,7 @@ from interputils import interpmat from ipythonUtils import easyAnimate as animate import Solver from Solver import Solver +import Save import Geophysics import types @@ -23,7 +24,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.') @@ -88,6 +89,15 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'): print pad + "%s%s%s" % ('-'*25,done,'-'*25) def callHooks(obj, match, *args, **kwargs): + """ + 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. + """ 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) From 5fd4290927702e191b4900d09cf3c0e4dacc1759 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 09:11:02 -0800 Subject: [PATCH 06/35] Fixed errors. C,F,E works with limits, extent. Added .cmap for color mapping and .range to set the color range. Need to add colorbar, and more settings as it comes along. --- SimPEG/visualize/vtk/vtkTools.py | 8 ++-- SimPEG/visualize/vtk/vtkView.py | 78 ++++++++++++++++++++++---------- 2 files changed, 58 insertions(+), 28 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index 7b405d4c..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)) @@ -257,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() @@ -273,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 dde60fd3..666cde57 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.' @@ -64,13 +64,37 @@ class vtkView(object): # 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): @@ -85,28 +109,24 @@ class vtkView(object): # Test the bounds change = 0 # Test for lower bounds, can't be smaller the 0 - #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) 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 - #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) 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 = np.prod(valnp[1::2] > upB,axis=0,dtype=bool) 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 = np.prod(valnp[::2] < valnp[1::2],axis=0,dtype=bool) tgt = valnp[::2] > valnp[1::2] if tgt.any(): valnp[1::2][tgt] = valnp[::2][tgt] + 1 @@ -123,8 +143,6 @@ class vtkView(object): 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: @@ -143,7 +161,6 @@ class vtkView(object): 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: @@ -158,6 +175,20 @@ class vtkView(object): 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 _readPropertyDictionary(self,propdict): """ Reads the property and assigns to the object @@ -192,9 +223,17 @@ class vtkView(object): # Make renderwindow. Returns the interactor. self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) - imageType = self.viewprop.keys()[0] - + + # Set the active scalar. + if type(self.viewprop.values()[0]) == int: + 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._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() @@ -211,15 +250,7 @@ class vtkView(object): 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)) - - # Set the active scalar. - if type(self.viewprop.values()[0]) == int: - actScalar = self._vtkobj.GetCellData().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._vtkobj.GetCellData().SetActiveScalars(actScalar) # Set up the plane, clipper and the user interaction. global intPlane, intActor self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) @@ -235,10 +266,11 @@ 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 + self._actor.GetMapper().SetScalarRange(self.range) self._actor.GetMapper().SetLookupTable(lut) # Set renderer options From 4cdafc12c1b5a29a568b3bfc433fb9b9c21db090 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 09:29:26 -0800 Subject: [PATCH 07/35] Added a colorbar to the figure. --- SimPEG/visualize/vtk/vtkView.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index 666cde57..8143f6e2 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -189,6 +189,17 @@ class vtkView(object): 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): """ Reads the property and assigns to the object @@ -270,12 +281,23 @@ class vtkView(object): 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) From 9cca076651fe0487f93271ef66772cdd3598128e Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 09:52:25 -0800 Subject: [PATCH 08/35] Updated the example at the bottom of the vtkView. __name__ == __main__ --- SimPEG/visualize/vtk/vtkView.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index 8143f6e2..3b52b992 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -312,14 +312,18 @@ 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 = {'C':{'Test':np.arange(0,mesh.nC),'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)}} + 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 @@ -330,3 +334,17 @@ if __name__ == '__main__': # 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 From 674d89ed41663b53d0eb4f95eb93bb5dbe7e401e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 27 Nov 2013 12:36:28 -0700 Subject: [PATCH 09/35] update the requirements.txt file. --- requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements.txt b/requirements.txt index 6b7df820..77ac05b9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,3 +2,4 @@ numpy scipy ipython matplotlib +h5py From c9195f0651962fb2c332c55f14b38696c7862c46 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 12:51:16 -0800 Subject: [PATCH 10/35] Made a notebook that works with updated version of vtkView --- notebooks/3DRenderingWithvtkTools.ipynb | 115 ------------------------ 1 file changed, 115 deletions(-) delete mode 100644 notebooks/3DRenderingWithvtkTools.ipynb diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb deleted file mode 100644 index 8d2c833d..00000000 --- a/notebooks/3DRenderingWithvtkTools.ipynb +++ /dev/null @@ -1,115 +0,0 @@ -{ - "metadata": { -<<<<<<< HEAD - "name": "3D rendering with vtkTools" -======= - "name": "" ->>>>>>> develop - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ - { - "cells": [ - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "import numpy as np, vtk\nimport SimPEG as simpeg", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 1 -======= - "input": [ - "import numpy as np, vtk\n", - "import SimPEG as simpeg" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 -======= - "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 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'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.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 -======= - "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 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 -======= - "input": [ - "# Show the image \n", - "vtkViewer.Show()\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": "*" ->>>>>>> develop - } - ], - "metadata": {} - } - ] -} From e0e64a6f9e720c0d5caf39384f7204c850eedbf1 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 12:52:55 -0800 Subject: [PATCH 11/35] Commit the new notebook for vtkView --- notebooks/VisualizeWithvtkView-updated..ipynb | 142 ++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 notebooks/VisualizeWithvtkView-updated..ipynb diff --git a/notebooks/VisualizeWithvtkView-updated..ipynb b/notebooks/VisualizeWithvtkView-updated..ipynb new file mode 100644 index 00000000..6dabae72 --- /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 From 60c1fc7c426ab128f1b2e3dc9e6f09d03c9c1440 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 13:04:48 -0800 Subject: [PATCH 12/35] New notebook for vtkView example --- notebooks/VisualizeWithvtkView-updated.ipynb | 142 +++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 notebooks/VisualizeWithvtkView-updated.ipynb 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 From 1b58820b5e040071c93f6c5319de83b2779007e6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 27 Nov 2013 14:30:27 -0700 Subject: [PATCH 13/35] Delete 3DRenderingWithvtkTools.ipynb.orig --- notebooks/3DRenderingWithvtkTools.ipynb.orig | 115 ------------------- 1 file changed, 115 deletions(-) delete mode 100644 notebooks/3DRenderingWithvtkTools.ipynb.orig diff --git a/notebooks/3DRenderingWithvtkTools.ipynb.orig b/notebooks/3DRenderingWithvtkTools.ipynb.orig deleted file mode 100644 index 8d2c833d..00000000 --- a/notebooks/3DRenderingWithvtkTools.ipynb.orig +++ /dev/null @@ -1,115 +0,0 @@ -{ - "metadata": { -<<<<<<< HEAD - "name": "3D rendering with vtkTools" -======= - "name": "" ->>>>>>> develop - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ - { - "cells": [ - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "import numpy as np, vtk\nimport SimPEG as simpeg", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 1 -======= - "input": [ - "import numpy as np, vtk\n", - "import SimPEG as simpeg" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 -======= - "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 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'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.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 -======= - "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 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 -======= - "input": [ - "# Show the image \n", - "vtkViewer.Show()\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": "*" ->>>>>>> develop - } - ], - "metadata": {} - } - ] -} From f9dd00ca8851a4c55472bb64f72e627198891dac Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 27 Nov 2013 15:32:50 -0700 Subject: [PATCH 14/35] Delete InteractiveView.py --- SimPEG/InteractiveView.py | 335 -------------------------------------- 1 file changed, 335 deletions(-) delete mode 100644 SimPEG/InteractiveView.py diff --git a/SimPEG/InteractiveView.py b/SimPEG/InteractiveView.py deleted file mode 100644 index ba4f1b63..00000000 --- a/SimPEG/InteractiveView.py +++ /dev/null @@ -1,335 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -from utils import mkvc - - -class TensorView(object): - """ - Provides viewing functions for TensorMesh - - This class is inherited by TensorMesh - """ - def __init__(self): - pass - - def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): - """ - Mesh.plotImage(I) - - Plots scalar fields on the given mesh. - - Input: - - :param numpy.array I: scalar field - - Optional Input: - - :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz - :param int figNum: number of figure to plot to - :param matplotlib.axes.Axes ax: axis to plot to - :param str direction: slice dimensions, 3D only ('x', 'y', 'z') - :param bool numbering: show numbering of slices, 3D only - :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_image_2D.py - :include-source: - - .. plot:: examples/mesh/plot_image_3D.py - :include-source: - """ - assert type(I) == np.ndarray, "I must be a numpy array" - assert type(numbering) == bool, "numbering must be a bool" - assert direction in ["x", "y","z"], "direction must be either x,y, or z" - - - if imageType == 'CC': - assert I.size == self.nC, "Incorrect dimensions for CC." - elif imageType == 'N': - assert I.size == self.nN, "Incorrect dimensions for N." - elif imageType == 'Fx': - if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') - elif imageType == 'Fy': - if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') - elif imageType == 'Fz': - if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') - elif imageType == 'Ex': - if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') - elif imageType == 'Ey': - if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') - elif imageType == 'Ez': - if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') - elif imageType[0] == 'E': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - ex, ey, ez = self.r(I,'E','E','M') - if plotAll or 'Ex' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(ex, imageType='Ex', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Ey' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(ey, imageType='Ey', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Ez' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(ez, imageType='Ez', ax=ax_z, **options) - pltNum +=1 - return - elif imageType[0] == 'F': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - fx, fy, fz = self.r(I,'F','F','M') - if plotAll or 'Fx' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(fx, imageType='Fx', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Fy' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(fy, imageType='Fy', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Fz' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(fz, imageType='Fz', ax=ax_z, **options) - pltNum +=1 - return - else: - raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") - - - if ax is None: - fig = plt.figure(figNum) - fig.clf() - ax = plt.subplot(111) - else: - assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" - fig = ax.figure - - if self.dim == 1: - if imageType == 'CC': - ph = ax.plot(self.vectorCCx, I, '-ro') - elif imageType == 'N': - ph = ax.plot(self.vectorNx, I, '-bs') - ax.set_xlabel("x") - ax.axis('tight') - elif self.dim == 2: - if imageType == 'CC': - C = I[:].reshape(self.n, order='F') - elif imageType == 'N': - C = I[:].reshape(self.n+1, order='F') - C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) - elif imageType == 'Fx': - C = I[:].reshape(self.nFx, order='F') - C = 0.5*(C[:-1, :] + C[1:, :] ) - elif imageType == 'Fy': - C = I[:].reshape(self.nFy, order='F') - C = 0.5*(C[:, :-1] + C[:, 1:] ) - elif imageType == 'Ex': - C = I[:].reshape(self.nEx, order='F') - C = 0.5*(C[:,:-1] + C[:,1:] ) - elif imageType == 'Ey': - C = I[:].reshape(self.nEy, order='F') - C = 0.5*(C[:-1,:] + C[1:,:] ) - - ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) - ax.axis('tight') - ax.set_xlabel("x") - ax.set_ylabel("y") - - elif self.dim == 3: - if direction == 'z': - - # get copy of image and average to cell-centres is necessary - if imageType == 'CC': - Ic = I[:].reshape(self.n, order='F') - elif imageType == 'N': - Ic = I[:].reshape(self.n+1, order='F') - Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) - elif imageType == 'Fx': - Ic = I[:].reshape(self.nFx, order='F') - Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) - elif imageType == 'Fy': - Ic = I[:].reshape(self.nFy, order='F') - Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) - elif imageType == 'Fz': - Ic = I[:].reshape(self.nFz, order='F') - Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) - elif imageType == 'Ex': - Ic = I[:].reshape(self.nEx, order='F') - Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) - elif imageType == 'Ey': - Ic = I[:].reshape(self.nEy, order='F') - Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) - elif imageType == 'Ez': - Ic = I[:].reshape(self.nEz, order='F') - Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) - - # determine number oE slices in x and y dimension - nX = np.ceil(np.sqrt(self.nCz)) - nY = np.ceil(self.nCz/nX) - - # allocate space for montage - nCx = self.nCx - nCy = self.nCy - - C = np.zeros((nX*nCx,nY*nCy)) - - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] - else: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan - - C = np.ma.masked_where(np.isnan(C), C) - xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] - yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] - # Plot the mesh - ph = ax.pcolormesh(xx, yy, C.T) - # Plot the lines - gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) - gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) - # Repeat and seperate with NaN - gxX = np.c_[gx, gx, gx+np.nan].ravel() - gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() - gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() - gyY = np.c_[gy, gy, gy+np.nan].ravel() - ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) - ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) - ax.axis('tight') - - if numbering: - pad = np.sum(self.hx)*0.04 - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, - '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') - - ax.set_title(imageType) - if showIt: plt.show() - return ph - - def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. - - :param bool nodes: plot nodes - :param bool faces: plot faces - :param bool centers: plot centers - :param bool edges: plot edges - :param bool lines: plot lines connecting nodes - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_grid_2D.py - :include-source: - - .. plot:: examples/mesh/plot_grid_3D.py - :include-source: - """ - if self.dim == 1: - fig = plt.figure(1) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - ax.hold(True) - ax.plot(xn, np.ones(np.shape(xn)), 'bs') - ax.plot(xc, np.ones(np.shape(xc)), 'ro') - ax.plot(xn, np.ones(np.shape(xn)), 'k--') - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - if showIt: plt.show() - elif self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - xs1 = self.gridFx - xs2 = self.gridFy - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') - if faces: - ax.plot(xs1[:, 0], xs1[:, 1], 'g>') - ax.plot(xs2[:, 0], xs2[:, 1], 'g^') - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] - plt.plot(X, Y) - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - if showIt: plt.show() - elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') - xn = self.gridN - xc = self.gridCC - xfs1 = self.gridFx - xfs2 = self.gridFy - xfs3 = self.gridFz - - xes1 = self.gridEx - xes2 = self.gridEy - xes3 = self.gridEz - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) - if faces: - ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) - ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) - ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) - if edges: - ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) - ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) - ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() - Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() - Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() - X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() - X = np.r_[X1, X2, X3] - Y = np.r_[Y1, Y2, Y3] - Z = np.r_[Z1, Z2, Z3] - plt.plot(X, Y, 'b-', zs=Z) - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - ax.set_zlabel('x3') - if showIt: plt.show() From 2a79e31d1f011707e6867d4024a51df98d892753 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 27 Nov 2013 15:33:18 -0700 Subject: [PATCH 15/35] Delete VisualizeWithvtkView-updated..ipynb --- notebooks/VisualizeWithvtkView-updated..ipynb | 142 ------------------ 1 file changed, 142 deletions(-) delete mode 100644 notebooks/VisualizeWithvtkView-updated..ipynb diff --git a/notebooks/VisualizeWithvtkView-updated..ipynb b/notebooks/VisualizeWithvtkView-updated..ipynb deleted file mode 100644 index 6dabae72..00000000 --- a/notebooks/VisualizeWithvtkView-updated..ipynb +++ /dev/null @@ -1,142 +0,0 @@ -{ - "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 From a2259b78112c72c330cef20350aaebb6b0407284 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 11:48:17 -0800 Subject: [PATCH 16/35] updates to Newton Root Finding. --- SimPEG/inverse/Optimize.py | 27 ++++++++++++++++++--------- SimPEG/tests/test_optimizers.py | 10 ++++++++++ 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 7e934e13..6edacabd 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -795,12 +795,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 +823,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 +837,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 +863,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..654d1804 100644 --- a/SimPEG/tests/test_optimizers.py +++ b/SimPEG/tests/test_optimizers.py @@ -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() From 9f76d2c3d9ba249135f489f39f761f2e2bf907fc Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 13:23:58 -0800 Subject: [PATCH 17/35] try to get travis working with hdf5 --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 1918d175..8ef7869a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,7 +4,7 @@ python: virtualenv: system_site_packages: true before_install: - - sudo apt-get install -qq python-numpy python-scipy python-matplotlib + - sudo apt-get install -qq python-numpy python-scipy python-matplotlib libhdf5-dev - python SimPEG/setup.py # command to install dependencies install: "pip install -r requirements.txt --use-mirrors" From c9b98cc7c01642436f37d92ce4d9306ff19d01c7 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 14:12:37 -0800 Subject: [PATCH 18/35] Recursive printing to view the inversion. --- .travis.yml | 2 +- SimPEG/utils/Save.py | 126 ++++++++++++++++++++++++++++--------------- 2 files changed, 84 insertions(+), 44 deletions(-) diff --git a/.travis.yml b/.travis.yml index 8ef7869a..1918d175 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,7 +4,7 @@ python: virtualenv: system_site_packages: true before_install: - - sudo apt-get install -qq python-numpy python-scipy python-matplotlib libhdf5-dev + - sudo apt-get install -qq python-numpy python-scipy python-matplotlib - python SimPEG/setup.py # command to install dependencies install: "pip install -r requirements.txt --use-mirrors" diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 9f661b0d..45cc48c8 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -1,8 +1,12 @@ import numpy as np import time -import h5py import re +try: + import h5py +except Exception, e: + print 'Warning: SimPEG table needs h5py to be installed.' + def natural_keys(text): ''' alist.sort(key=natural_keys) sorts in human order @@ -25,6 +29,52 @@ class SimPEGTable: self.inversions = hdf5InversionGroup(self,self.root.addGroup('inversions',soft=True)) + + def saveInversion(self, invObj, dataPath): + + invObj._invNode = self.inversions.addGroup('%d'%self.numChildren) + + # At the start of every iteration we will create a inversion iteration node. + def _doStartIteration_hdf5_inv(invObj): + invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) + invNodeIt.attrs['complete'] = False + invNodeIt.attrs['time'] = time.time() + invObj._invNodeIt = invNodeIt + invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) + + def _doEndIteration_hdf5_inv(invObj): + invNodeIt = invObj._invNodeIt + invNodeIt.attrs['time'] = time.time() - invNodeIt.attrs['time'] + invNodeIt.attrs['ctime'] = time.ctime() + invNodeIt.attrs['phi_d'] = invObj.phi_d + invNodeIt.attrs['phi_m'] = invObj.phi_m + invNodeIt.setArray('m', invObj.m) + invNodeIt.setArray('dpred', invObj.dpred) + invNodeIt.attrs['complete'] = True + invObj.hook(_doEndIteration_hdf5_inv, overwrite=True) + + def _doStartIteration_hdf5_opt(optObj): + optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) + optNodeIt.attrs['complete'] = False + optNodeIt.attrs['time'] = time.time() + optObj._optNodeIt = optNodeIt + + invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True) + + def _doEndIteration_hdf5_opt(optObj, xt): + optNodeIt = optObj._optNodeIt + optNodeIt.attrs['time'] = time.time() - optNodeIt.attrs['time'] + optNodeIt.attrs['ctime'] = time.ctime() + optNodeIt.setArray('m', xt) + optNodeIt.setArray('dpred', optObj.parent.dpred) + optNodeIt.setArray('searchDirection', optObj.searchDirection) + optNodeIt.attrs['complete'] = True + invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True) + + return invObj._invNode + + + class hdf5Group(object): """Has some low level support for wrapping the native HDF5-Group class""" @@ -93,54 +143,44 @@ class hdf5Group(object): child = self.childClass(self.T, child) return child + 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)>' % (self.__class__.__name__, self.path, self.numChildren, '' if self.numChildren == 1 else 's') + + class hdf5InversionGroup(hdf5Group): def __init__(self, T, groupNode): hdf5Group.__init__(self, T, groupNode) self.childClass = hdf5Inversion - def saveInversion(self, invObj, dataPath): - invObj._invNode = self.addGroup('%d'%self.numChildren) - - # At the start of every iteration we will create a inversion iteration node. - def _doStartIteration_hdf5_inv(invObj): - invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) - invNodeIt.attrs['complete'] = False - invNodeIt.attrs['time'] = time.time() - invObj._invNodeIt = invNodeIt - invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) - - def _doEndIteration_hdf5_inv(invObj): - invNodeIt = invObj._invNodeIt - invNodeIt.attrs['time'] = time.time() - invNodeIt.attrs['time'] - invNodeIt.attrs['ctime'] = time.ctime() - invNodeIt.attrs['phi_d'] = invObj.phi_d - invNodeIt.attrs['phi_m'] = invObj.phi_m - invNodeIt.setArray('m', invObj.m) - invNodeIt.setArray('dpred', invObj.dpred) - invNodeIt.attrs['complete'] = True - invObj.hook(_doEndIteration_hdf5_inv, overwrite=True) - - def _doStartIteration_hdf5_opt(optObj): - optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) - optNodeIt.attrs['complete'] = False - optNodeIt.attrs['time'] = time.time() - optObj._optNodeIt = optNodeIt - - invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True) - - def _doEndIteration_hdf5_opt(optObj, xt): - optNodeIt = optObj._optNodeIt - optNodeIt.attrs['time'] = time.time() - optNodeIt.attrs['time'] - optNodeIt.attrs['ctime'] = time.ctime() - optNodeIt.setArray('m', xt) - optNodeIt.setArray('dpred', optObj.parent.dpred) - optNodeIt.setArray('searchDirection', optObj.searchDirection) - optNodeIt.attrs['complete'] = True - invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True) - - return invObj._invNode - class hdf5Inversion(hdf5Group): def __init__(self, T, groupNode): From 18f9dc67d5d9bd8c6d90f8244d96795c83e07e5a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 14:17:56 -0800 Subject: [PATCH 19/35] bug fix. --- SimPEG/utils/Save.py | 2 +- requirements.txt | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 45cc48c8..4c646971 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -32,7 +32,7 @@ class SimPEGTable: def saveInversion(self, invObj, dataPath): - invObj._invNode = self.inversions.addGroup('%d'%self.numChildren) + invObj._invNode = self.inversions.addGroup('%d'%self.inversions.numChildren) # At the start of every iteration we will create a inversion iteration node. def _doStartIteration_hdf5_inv(invObj): diff --git a/requirements.txt b/requirements.txt index 77ac05b9..6b7df820 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,3 @@ numpy scipy ipython matplotlib -h5py From ccb4d4052de7b59ace47eeb8a10af81454a4a61a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 15:31:11 -0800 Subject: [PATCH 20/35] Change callHooks structure. --- SimPEG/inverse/Inversion.py | 42 ++++++++++----------------------- SimPEG/inverse/Optimize.py | 40 ++++++++++++------------------- SimPEG/tests/test_optimizers.py | 2 +- SimPEG/utils/__init__.py | 38 +++++++++++++++++++++++++---- 4 files changed, 62 insertions(+), 60 deletions(-) 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 6edacabd..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): diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/tests/test_optimizers.py index 654d1804..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.]) diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 622557d6..19aecec6 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -23,7 +23,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.') @@ -87,11 +87,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): - 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 class Counter(object): From 191d0e242c86b9c33705e76a2ee42ade121d32f3 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 16:10:31 -0800 Subject: [PATCH 21/35] Create a new inversion node every time inversion is ran. --- SimPEG/inverse/Optimize.py | 10 +++++++++- SimPEG/utils/Save.py | 18 ++++++++++++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 80a03632..41fd3c22 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -266,9 +266,17 @@ class Minimize(object): printStoppers(self, stoppers, pad='', stop=stop, done=done) + @callHooks('finish') def finish(self): - pass + """finish() + **finish** is called at the end of the optimization. + + :rtype: None + :return: None + + """ + pass def stoppingCriteria(self, inLS=False): if self._iter == 0: diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 4c646971..e7bec76d 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -29,10 +29,14 @@ class SimPEGTable: self.inversions = hdf5InversionGroup(self,self.root.addGroup('inversions',soft=True)) + def show(self): self.root.show() def saveInversion(self, invObj, dataPath): - invObj._invNode = self.inversions.addGroup('%d'%self.inversions.numChildren) + # Create a new inversion anytime this is run. + def _startup_hdf5_inv(invObj, m0): + invObj._invNode = self.inversions.addGroup('%d'%self.inversions.numChildren) + 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): @@ -53,6 +57,14 @@ class SimPEGTable: invNodeIt.attrs['complete'] = True invObj.hook(_doEndIteration_hdf5_inv, overwrite=True) + # Delete all iterates that did not finish. + def _finish_hdf5_inv(invObj): + for it in invObj._invNode: + if not it.attrs['complete']: + del self.f[it.path] + del invObj._invNode + invObj.hook(_finish_hdf5_inv, overwrite=True) + def _doStartIteration_hdf5_opt(optObj): optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) optNodeIt.attrs['complete'] = False @@ -71,7 +83,6 @@ class SimPEGTable: optNodeIt.attrs['complete'] = True invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True) - return invObj._invNode @@ -143,6 +154,9 @@ class hdf5Group(object): 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. From 8fb9f3683b509595d54c01bbf95878de5e65d2ac Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 17:52:40 -0800 Subject: [PATCH 22/35] Easy way to create a padded tensor mesh. --- SimPEG/mesh/TensorMesh.py | 6 ++- SimPEG/utils/__init__.py | 4 +- SimPEG/utils/lomutils.py | 23 ----------- SimPEG/utils/meshutils.py | 58 +++++++++++++++++++++++++++ docs/api_Utils.rst | 7 ++++ docs/examples/mesh/plot_TensorMesh.py | 21 ---------- 6 files changed, 73 insertions(+), 46 deletions(-) create mode 100644 SimPEG/utils/meshutils.py delete mode 100644 docs/examples/mesh/plot_TensorMesh.py diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 4c0841b8..65f94b3e 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -24,7 +24,11 @@ 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:: diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 19aecec6..4bb1e445 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,9 +3,11 @@ 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 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/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() From e25afacd1b65808e044ac260a5340bf3e9d0ec1f Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 17:55:09 -0800 Subject: [PATCH 23/35] refactor and generalize where things are being saved --- SimPEG/inverse/Inversion.py | 6 ++++++ SimPEG/inverse/Optimize.py | 11 +++++++++++ SimPEG/utils/Save.py | 35 ++++++++++++++++------------------- 3 files changed, 33 insertions(+), 19 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index a2ef2513..99e89e61 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -349,6 +349,12 @@ 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 diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 41fd3c22..3c4a6a5a 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -431,6 +431,17 @@ 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.setArray('m', self.xc) + group.setArray('dpred', self.parent.dpred) + + + class Remember(object): """ This mixin remembers all the things you tend to forget. diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index e7bec76d..a224f4c6 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -17,6 +17,15 @@ def natural_keys(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. @@ -41,20 +50,13 @@ class SimPEGTable: # At the start of every iteration we will create a inversion iteration node. def _doStartIteration_hdf5_inv(invObj): invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) - invNodeIt.attrs['complete'] = False - invNodeIt.attrs['time'] = time.time() + preIteration(invNodeIt) invObj._invNodeIt = invNodeIt invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) def _doEndIteration_hdf5_inv(invObj): - invNodeIt = invObj._invNodeIt - invNodeIt.attrs['time'] = time.time() - invNodeIt.attrs['time'] - invNodeIt.attrs['ctime'] = time.ctime() - invNodeIt.attrs['phi_d'] = invObj.phi_d - invNodeIt.attrs['phi_m'] = invObj.phi_m - invNodeIt.setArray('m', invObj.m) - invNodeIt.setArray('dpred', invObj.dpred) - invNodeIt.attrs['complete'] = True + invObj.save(invObj._invNodeIt) + postIteration(invObj._invNodeIt) invObj.hook(_doEndIteration_hdf5_inv, overwrite=True) # Delete all iterates that did not finish. @@ -63,24 +65,19 @@ class SimPEGTable: 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): optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) - optNodeIt.attrs['complete'] = False - optNodeIt.attrs['time'] = time.time() + preIteration(optNodeIt) optObj._optNodeIt = optNodeIt invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True) def _doEndIteration_hdf5_opt(optObj, xt): - optNodeIt = optObj._optNodeIt - optNodeIt.attrs['time'] = time.time() - optNodeIt.attrs['time'] - optNodeIt.attrs['ctime'] = time.ctime() - optNodeIt.setArray('m', xt) - optNodeIt.setArray('dpred', optObj.parent.dpred) - optNodeIt.setArray('searchDirection', optObj.searchDirection) - optNodeIt.attrs['complete'] = True + optObj.save(optObj._optNodeIt) + postIteration(optObj._optNodeIt) invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True) From f3c139ee9e72c03c146b8ea37610eb13e37fa4b1 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 5 Dec 2013 11:30:21 -0800 Subject: [PATCH 24/35] Warning on newton root finding. --- SimPEG/inverse/Optimize.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 3c4a6a5a..d247e324 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -834,7 +834,7 @@ class NewtonRoot(object): xt = x + dh 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) @@ -852,7 +852,10 @@ class NewtonRoot(object): 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 From 3221b3c7504df1ba613730b3923e0c8fc1de3ee3 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 5 Dec 2013 11:31:43 -0800 Subject: [PATCH 25/35] Updates to video function (skip frames) --- SimPEG/mesh/TensorView.py | 8 +++++--- SimPEG/utils/ModelBuilder.py | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) 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/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)) From 5e0fb8642d0e70595f4fa3d90b69887999d92e9f Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 11:23:01 -0800 Subject: [PATCH 26/35] MeatClasses to make every SimPEG object saveable to an hdf5 file. cleaned up imports in a lot of places. Made solver not copy matrices around for GS preconditioning. --- SimPEG/__init__.py | 2 + SimPEG/forward/LinearProblem.py | 3 + SimPEG/forward/Problem.py | 20 +++---- SimPEG/inverse/Inversion.py | 38 ++++++------ SimPEG/inverse/Optimize.py | 78 +++++++++++++----------- SimPEG/mesh/BaseMesh.py | 13 ++-- SimPEG/utils/Save.py | 102 +++++++++++++++++++++++++++++++- SimPEG/utils/Solver.py | 6 +- SimPEG/utils/__init__.py | 2 +- 9 files changed, 188 insertions(+), 76 deletions(-) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 1ea0601c..e8295092 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,3 +1,5 @@ +import numpy as np +import scipy.sparse as sp import utils from utils import Solver import mesh diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py index ef92957b..5a996c71 100644 --- a/SimPEG/forward/LinearProblem.py +++ b/SimPEG/forward/LinearProblem.py @@ -9,6 +9,9 @@ import matplotlib.pyplot as plt class LinearProblem(Problem): """docstring for LinearProblem""" + def __init__(self, *args, **kwargs): + Problem.__init__(self, *args, **kwargs) + def dpred(self, m, u=None): return self.G.dot(m) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 99747f7d..1f3d817c 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, np, sp norm = np.linalg.norm @@ -37,6 +35,8 @@ class Problem(object): to (locally) find how model parameters change the data, and optimize! """ + __metaclass__ = utils.Save.Savable + counter = None def __init__(self, mesh): @@ -85,7 +85,7 @@ class Problem(object): def dobs(self, value): self._dobs = value - @count + @utils.count def dpred(self, m, u=None): """ Predicted data. @@ -97,7 +97,7 @@ class Problem(object): u = self.field(m) return self.P*u - @count + @utils.count def dataResidual(self, m, u=None): """ :param numpy.array m: geophysical model @@ -117,7 +117,7 @@ class Problem(object): return self.dpred(m, u=u) - self.dobs - @timeIt + @utils.timeIt def J(self, m, v, u=None): """ :param numpy.array m: model @@ -147,7 +147,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 +161,7 @@ class Problem(object): raise NotImplementedError('Jt is not yet implemented.') - @timeIt + @utils.timeIt def J_approx(self, m, v, u=None): """ @@ -176,7 +176,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 @@ -241,7 +241,7 @@ class Problem(object): 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 + eps = np.linalg.norm(utils.mkvc(dobs),2)*1e-5 Wd = 1/(abs(dobs)*std+eps) return dobs, Wd diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 99e89e61..86a5dd14 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,7 +1,5 @@ -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 @@ -9,6 +7,8 @@ from SimPEG.inverse import IterationPrinters, StoppingCriteria class BaseInversion(object): """docstring for BaseInversion""" + __metaclass__ = utils.Save.Savable + maxIter = 1 #: Maximum number of iterations name = 'BaseInversion' @@ -17,11 +17,11 @@ 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) + utils.setKwargs(self, **kwargs) self.prob = prob self.reg = reg self.opt = opt @@ -46,7 +46,7 @@ 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 + eps = np.linalg.norm(utils.mkvc(self.prob.dobs),2)*1e-5 self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) return self._Wd @Wd.setter @@ -70,7 +70,7 @@ class BaseInversion(object): def phi_d_target(self, value): self._phi_d_target = value - @timeIt + @utils.timeIt def run(self, m0): """run(m0) @@ -89,7 +89,7 @@ class BaseInversion(object): return self.m - @callHooks('startup') + @utils.callHooks('startup') def startup(self, m0): """ **startup** is called at the start of any new run call. @@ -109,7 +109,7 @@ class BaseInversion(object): self.phi_d_last = np.nan self.phi_m_last = np.nan - @callHooks('doStartIteration') + @utils.callHooks('doStartIteration') def doStartIteration(self): """ **doStartIteration** is called at the end of each run iteration. @@ -120,7 +120,7 @@ class BaseInversion(object): self._beta = self.getBeta() - @callHooks('doEndIteration') + @utils.callHooks('doEndIteration') def doEndIteration(self): """ **doEndIteration** is called at the end of each run iteration. @@ -179,7 +179,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): @@ -189,7 +189,7 @@ class BaseInversion(object): """ printStoppers(self, self.stoppers) - @callHooks('finish') + @utils.callHooks('finish') def finish(self): """finish() @@ -197,7 +197,7 @@ class BaseInversion(object): """ pass - @timeIt + @utils.timeIt def evalFunction(self, m, return_g=True, return_H=True): """evalFunction(m, return_g=True, return_H=True) @@ -207,7 +207,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) @@ -237,7 +237,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) @@ -257,10 +257,10 @@ class BaseInversion(object): """ # 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 = utils.mkvc(R) return 0.5*np.vdot(R, R) - @timeIt + @utils.timeIt def dataObjDeriv(self, m, u=None): """dataObjDeriv(m, u=None) @@ -302,7 +302,7 @@ class BaseInversion(object): return dmisfit - @timeIt + @utils.timeIt def dataObj2Deriv(self, m, v, u=None): """dataObj2Deriv(m, v, u=None) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index d247e324..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) @@ -190,7 +189,7 @@ class Minimize(object): def parent(self, value): self._parent = value - @callHooks('startup') + @utils.callHooks('startup') def startup(self, x0): """ **startup** is called at the start of any new minimize call. @@ -215,8 +214,8 @@ class Minimize(object): self.f_last = np.nan self.x_last = x0 - @count - @callHooks('doStartIteration') + @utils.count + @utils.callHooks('doStartIteration') def doStartIteration(self): """doStartIteration() @@ -238,9 +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) - @callHooks('printIter') + @utils.callHooks('printIter') def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -250,7 +249,7 @@ class Minimize(object): """ 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): """ @@ -263,10 +262,10 @@ 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) - @callHooks('finish') + @utils.callHooks('finish') def finish(self): """finish() @@ -282,10 +281,10 @@ class Minimize(object): 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 - @callHooks('projection') + @utils.timeIt + @utils.callHooks('projection') def projection(self, p): """projection(p) @@ -299,7 +298,7 @@ class Minimize(object): """ return p - @timeIt + @utils.timeIt def findSearchDirection(self): """findSearchDirection() @@ -330,7 +329,7 @@ class Minimize(object): """ return -self.g - @count + @utils.count def scaleSearchDirection(self, p): """scaleSearchDirection(p) @@ -349,7 +348,7 @@ class Minimize(object): nameLS = "Armijo linesearch" #: The line-search name - @timeIt + @utils.timeIt def modifySearchDirection(self, p): """modifySearchDirection(p) @@ -387,7 +386,7 @@ class Minimize(object): return self._LS_xt, self._iterLS < self.maxIterLS - @count + @utils.count def modifySearchDirectionBreak(self, p): """modifySearchDirectionBreak(p) @@ -409,8 +408,8 @@ class Minimize(object): print 'The linesearch got broken. Boo.' return p, False - @count - @callHooks('doEndIteration') + @utils.count + @utils.callHooks('doEndIteration') def doEndIteration(self, xt): """doEndIteration(xt) @@ -437,6 +436,8 @@ class Minimize(object): 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) @@ -526,7 +527,7 @@ class ProjectedGradient(Minimize, Remember): self.aSet_prev = self.activeSet(x0) - @count + @utils.count def projection(self, x): """projection(x) @@ -535,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) @@ -544,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) @@ -553,7 +554,7 @@ class ProjectedGradient(Minimize, Remember): """ return np.logical_not(self.activeSet(x)) - @count + @utils.count def bindingSet(self, x): """bindingSet(x) @@ -566,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() @@ -611,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) @@ -646,6 +647,9 @@ class BFGS(Minimize, Remember): name = 'BFGS' nbfgs = 10 + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + @property def bfgsH0(self): """ @@ -711,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) @@ -758,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) @@ -768,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 @@ -801,7 +811,7 @@ class NewtonRoot(object): doLS = True def __init__(self, **kwargs): - setKwargs(self, **kwargs) + utils.setKwargs(self, **kwargs) def root(self, fun, x): """root(fun, x) @@ -875,7 +885,7 @@ if __name__ == '__main__': print 'test the newtonRoot finding.' - fun = lambda x, return_g=True: np.sin(x) if not return_g else ( 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=True).root(fun,x) print pnt diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 62ec4251..ced9bd59 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,8 @@ class BaseMesh(object): :param numpy.array,list x0: Origin of the mesh (dim, ) """ + __metaclass__ = utils.Save.Savable + def __init__(self, n, x0=None): # Check inputs @@ -78,7 +79,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 +92,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 +102,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 +309,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/utils/Save.py b/SimPEG/utils/Save.py index a224f4c6..318723ce 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -7,6 +7,9 @@ try: except Exception, e: print 'Warning: SimPEG table needs h5py to be installed.' + +SAVEABLES = {} + def natural_keys(text): ''' alist.sort(key=natural_keys) sorts in human order @@ -57,6 +60,7 @@ class SimPEGTable: 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. @@ -78,11 +82,11 @@ class SimPEGTable: 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""" @@ -183,7 +187,7 @@ class hdf5Group(object): return s def __str__(self): - return '<%s "%s" (%d member%s)>' % (self.__class__.__name__, self.path, self.numChildren, '' if self.numChildren == 1 else 's') + 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])) @@ -204,3 +208,97 @@ class hdf5InversionIteration(hdf5Group): def __init__(self, T, groupNode): hdf5Group.__init__(self, T, groupNode) self.parentClass = hdf5Inversion + + + +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): + """ + """ + assert type(obj.__class__) is Savable, 'Can only save objects that are Savable objects.' + + def doSave(grp, name, val): + if type(val.__class__) is Savable: + subgrp = grp.addGroup(name) + saveInitArgs(val, subgrp) + elif type(val) is np.ndarray: + grp.setArray(name, val) + elif type(val) in [list, tuple]: + # Split up, and save each element + for i, v in enumerate(val): + doSave(grp, name + '[%d]'%i, v) + else: + # just try saving it as an attr + grp.attrs[name] = val + + 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) + + +def loadSavable(node): + + 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=utils.Save.natural_keys) + kwargs.sort(key=utils.Save.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 '[' 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(load(val)) + 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 '[' 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] = load(val) + else: + key = name[7:] + KWARGS[key] = val + + cls = get(node, '__class__') + if cls in SAVEABLES: + return SAVEABLES[cls](*ARGS,**KWARGS) + else: + print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls + return (cls, ARGS, KWARGS) + 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 9d0946cf..f75fed3d 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -45,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) From 9e26543931b6a0e35f5698fd586370d51b0a556a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 14:19:37 -0800 Subject: [PATCH 27/35] fixes to the meta classes including soft linking inside the database. These are reflected in the loaded outputs which return the same object. --- SimPEG/inverse/Inversion.py | 2 +- SimPEG/mesh/BaseMesh.py | 1 - SimPEG/mesh/LogicallyOrthogonalMesh.py | 96 +++++++++++++------------ SimPEG/mesh/TensorMesh.py | 60 ++++++++-------- SimPEG/regularization/Regularization.py | 41 +++++------ SimPEG/utils/Save.py | 74 ++++++++++++++----- 6 files changed, 157 insertions(+), 117 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 86a5dd14..19280e0a 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -187,7 +187,7 @@ class BaseInversion(object): **printDone** is called at the end of the inversion routine. """ - printStoppers(self, self.stoppers) + utils.printStoppers(self, self.stoppers) @utils.callHooks('finish') def finish(self): diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index ced9bd59..6151d20b 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -11,7 +11,6 @@ class BaseMesh(object): :param numpy.array,list x0: Origin of the mesh (dim, ) """ - __metaclass__ = utils.Save.Savable def __init__(self, n, x0=None): 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 65f94b3e..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): """ @@ -35,6 +32,9 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): mesh = TensorMesh([10, 12, 15]) """ + + __metaclass__ = utils.Save.Savable + _meshType = 'TENSOR' def __init__(self, h_in, x0=None): @@ -52,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) @@ -170,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 @@ -181,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 @@ -192,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 @@ -203,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 @@ -214,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 @@ -225,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 @@ -236,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 @@ -247,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 @@ -262,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 @@ -289,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 @@ -311,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 @@ -410,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/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 41938152..5c50ff59 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -1,9 +1,21 @@ -from SimPEG.utils import sdiag, count, timeIt, setKwargs -import numpy as np +from SimPEG import utils, np class Regularization(object): """docstring for Regularization""" + __metaclass__ = utils.Save.Savable + + alpha_s = 1e-6 + alpha_x = 1.0 + alpha_y = 1.0 + alpha_z = 1.0 + + 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: @@ -16,43 +28,32 @@ class Regularization(object): @property def Ws(self): if getattr(self,'_Ws', None) is None: - self._Ws = sdiag(self.mesh.vol) + self._Ws = utils.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) + self._Wx = self.mesh.cellGradx*utils.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) + self._Wy = self.mesh.cellGrady*utils.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) + self._Wz = self.mesh.cellGradz*utils.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 + @utils.timeIt def modelObj(self, m): mresid = m - self.mref @@ -67,7 +68,7 @@ class Regularization(object): return mobj - @timeIt + @utils.timeIt def modelObjDeriv(self, m): """ @@ -103,7 +104,7 @@ class Regularization(object): return mobjDeriv - @timeIt + @utils.timeIt def modelObj2Deriv(self): mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 318723ce..d9b3dd58 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -5,7 +5,7 @@ import re try: import h5py except Exception, e: - print 'Warning: SimPEG table needs h5py to be installed.' + print 'Warning: SimPEG.utils.Save needs h5py to be installed.' SAVEABLES = {} @@ -47,7 +47,10 @@ class SimPEGTable: # Create a new inversion anytime this is run. def _startup_hdf5_inv(invObj, m0): - invObj._invNode = self.inversions.addGroup('%d'%self.inversions.numChildren) + node = self.inversions.addGroup('%d'%self.inversions.numChildren) + saveSavable(invObj,node.addGroup('rebuild')) + results = node.addGroup('results') + invObj._invNode = results invObj.hook(_startup_hdf5_inv, overwrite=True) # At the start of every iteration we will create a inversion iteration node. @@ -196,18 +199,25 @@ class hdf5InversionGroup(hdf5Group): 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 = hdf5InversionIteration + 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 = hdf5Inversion + self.parentClass = hdf5InversionResults @@ -225,38 +235,61 @@ class Savable(type): return newClass -def saveSavable(obj, group): +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: - subgrp = grp.addGroup(name) - saveInitArgs(val, subgrp) - elif type(val) is np.ndarray: - grp.setArray(name, val) + 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 - grp.attrs[name] = val + try: + grp.attrs[name] = val + except Exception, e: + print 'Warning: Could not save %s, problems may arise is 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): +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=utils.Save.natural_keys) - kwargs.sort(key=utils.Save.natural_keys) + args.sort(key=natural_keys) + kwargs.sort(key=natural_keys) def get(node,key): if key in node.children: return node[key] @@ -266,6 +299,7 @@ def loadSavable(node): 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 @@ -273,7 +307,7 @@ def loadSavable(node): else: ARGS[ind].append(val) elif issubclass(val.__class__,hdf5Group): - ARGS.append(load(val)) + ARGS.append(loadSavable(val,pointers=pointers)) else: ind = int(name[4:]) ARGS.append(val) @@ -282,6 +316,7 @@ def loadSavable(node): 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 @@ -290,15 +325,18 @@ def loadSavable(node): KWARGS[key].append(val) elif issubclass(val.__class__,hdf5Group): key = name[7:] - KWARGS[key] = load(val) + KWARGS[key] = loadSavable(val,pointers=pointers) else: key = name[7:] KWARGS[key] = val cls = get(node, '__class__') if cls in SAVEABLES: - return SAVEABLES[cls](*ARGS,**KWARGS) + out = SAVEABLES[cls](*ARGS, **KWARGS) + out._savable = node + pointers.append(out) # Because this is recursive. + return out else: print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls - return (cls, ARGS, KWARGS) + return (cls, ARGS, KWARGS, node) From b5eb72e097a1fa83828e3e67238bbf4d324fd077 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 14:27:59 -0800 Subject: [PATCH 28/35] Moved DC and Linear to the examples directory. --- .../{forward/DCProblem.py => examples/DC.py} | 0 .../LinearProblem.py => examples/Linear.py} | 18 +++++++----------- SimPEG/examples/__init__.py | 2 ++ SimPEG/forward/__init__.py | 2 -- SimPEG/tests/test_forward_DCproblem.py | 2 +- 5 files changed, 10 insertions(+), 14 deletions(-) rename SimPEG/{forward/DCProblem.py => examples/DC.py} (100%) rename SimPEG/{forward/LinearProblem.py => examples/Linear.py} (76%) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/examples/DC.py similarity index 100% rename from SimPEG/forward/DCProblem.py rename to SimPEG/examples/DC.py diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/examples/Linear.py similarity index 76% rename from SimPEG/forward/LinearProblem.py rename to SimPEG/examples/Linear.py index 5a996c71..b593edf8 100644 --- a/SimPEG/forward/LinearProblem.py +++ b/SimPEG/examples/Linear.py @@ -1,16 +1,12 @@ -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, regularization, np import matplotlib.pyplot as plt -class LinearProblem(Problem): +class LinearProblem(forward.Problem): """docstring for LinearProblem""" def __init__(self, *args, **kwargs): - Problem.__init__(self, *args, **kwargs) + forward.Problem.__init__(self, *args, **kwargs) def dpred(self, m, u=None): return self.G.dot(m) @@ -24,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) @@ -63,9 +59,9 @@ if __name__ == '__main__': prob, m_true = example(100) M = prob.mesh - reg = Regularization(M) - opt = InexactGaussNewton(maxIter=20) - inv = Inversion(prob,reg,opt,beta0=1e-4) + reg = regularization.Regularization(M) + opt = inverse.InexactGaussNewton(maxIter=20) + inv = inverse.Inversion(prob,reg,opt,beta0=1e-4) m0 = np.zeros_like(m_true) mrec = inv.run(m0) 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/__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/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index 847d1055..c81c885a 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -3,7 +3,7 @@ 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 09938ed7df6661c76dff2186a45975af24f8202c Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 16:53:07 -0800 Subject: [PATCH 29/35] SimPEG data object --- SimPEG/__init__.py | 5 +++-- SimPEG/data/__init__.py | 20 ++++++++++++++++++++ 2 files changed, 23 insertions(+), 2 deletions(-) create mode 100644 SimPEG/data/__init__.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index e8295092..f455e866 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -3,10 +3,11 @@ import scipy.sparse as sp import utils from utils import Solver import mesh -import inverse -import visualize +import data import forward import regularization +import inverse +import visualize 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 From 783cd405387b15ecb95f46462e2e7bbb41c5fadb Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 16:54:38 -0800 Subject: [PATCH 30/35] Bug fixes in Save, incorporation of data into the inversion object, taken out of problem. --- SimPEG/examples/Linear.py | 30 +++++++++----------------- SimPEG/forward/Problem.py | 43 +++++++++++-------------------------- SimPEG/inverse/Inversion.py | 26 +++++++++++----------- SimPEG/utils/Save.py | 28 ++++++++++++++---------- 4 files changed, 53 insertions(+), 74 deletions(-) diff --git a/SimPEG/examples/Linear.py b/SimPEG/examples/Linear.py index b593edf8..79bec85c 100644 --- a/SimPEG/examples/Linear.py +++ b/SimPEG/examples/Linear.py @@ -34,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.Regularization(M) opt = inverse.InexactGaussNewton(maxIter=20) - inv = inverse.Inversion(prob,reg,opt,beta0=1e-4) - m0 = np.zeros_like(m_true) + inv = inverse.Inversion(prob,reg,opt,data) + m0 = np.zeros_like(data.mtrue) mrec = inv.run(m0) @@ -72,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/forward/Problem.py b/SimPEG/forward/Problem.py index 1f3d817c..1412ff5b 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -1,4 +1,4 @@ -from SimPEG import utils, np, sp +from SimPEG import utils, data, np, sp norm = np.linalg.norm @@ -37,9 +37,11 @@ class Problem(object): __metaclass__ = utils.Save.Savable - counter = None + counter = None #: A SimPEG.utils.Counter object - def __init__(self, mesh): + + def __init__(self, mesh, *args, **kwargs): + utils.setKwargs(self, **kwargs) self.mesh = mesh @property @@ -65,26 +67,6 @@ 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 - @utils.count def dpred(self, m, u=None): """ @@ -98,7 +80,7 @@ class Problem(object): return self.P*u @utils.count - def dataResidual(self, m, u=None): + def dataResidual(self, m, data, u=None): """ :param numpy.array m: geophysical model :param numpy.array u: fields @@ -115,7 +97,7 @@ 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 @utils.timeIt def J(self, m, v, u=None): @@ -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(utils.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/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 19280e0a..926f719a 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -5,7 +5,8 @@ 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 @@ -20,11 +21,12 @@ class BaseInversion(object): 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): + 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(utils.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,7 +65,7 @@ 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 @@ -256,7 +258,7 @@ 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 = self.Wd*self.prob.dataResidual(m, self.data, u=u) R = utils.mkvc(R) return 0.5*np.vdot(R, R) @@ -296,7 +298,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) dmisfit = self.prob.Jt(m, self.Wd * R, u=u) @@ -341,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. @@ -360,8 +362,8 @@ 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) @@ -377,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/utils/Save.py b/SimPEG/utils/Save.py index d9b3dd58..70ce1d2e 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -50,14 +50,15 @@ class SimPEGTable: 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): - invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) - preIteration(invNodeIt) - invObj._invNodeIt = invNodeIt + invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) + preIteration(invObj._invNodeIt) invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) def _doEndIteration_hdf5_inv(invObj): @@ -68,6 +69,7 @@ class SimPEGTable: # 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] @@ -76,10 +78,8 @@ class SimPEGTable: invObj.hook(_finish_hdf5_inv, overwrite=True) def _doStartIteration_hdf5_opt(optObj): - optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) - preIteration(optNodeIt) - optObj._optNodeIt = optNodeIt - + 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): @@ -332,10 +332,16 @@ def loadSavable(node, pointers=None): cls = get(node, '__class__') if cls in SAVEABLES: - out = SAVEABLES[cls](*ARGS, **KWARGS) - out._savable = node - pointers.append(out) # Because this is recursive. - return out + 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) From 7c9254841135ba92adccac4fb15c9598cc9ce481 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 17:04:08 -0800 Subject: [PATCH 31/35] Fixed DC resistivity tests. --- SimPEG/tests/test_forward_DCproblem.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index c81c885a..36216bf7 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -44,23 +44,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) + 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 +67,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 From 32277eff949738cb7ec930125b27afc322116aa9 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 7 Dec 2013 00:01:25 -0800 Subject: [PATCH 32/35] Documentation for Regularization. --- SimPEG/regularization/Regularization.py | 131 ++++++++++++++++++++++-- SimPEG/utils/Save.py | 2 +- docs/api_Optimize.rst | 9 ++ docs/api_Problem.rst | 4 +- 4 files changed, 134 insertions(+), 12 deletions(-) diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 5c50ff59..a8406d9b 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -1,14 +1,97 @@ -from SimPEG import utils, np +from SimPEG import utils, np, sp class Regularization(object): - """docstring for Regularization""" + """**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 = 1e-6 - alpha_x = 1.0 - alpha_y = 1.0 - alpha_z = 1.0 + alpha_s = 1e-6 #: Smallness weight + alpha_x = 1.0 #: Weight for the first derivative in the x direction + alpha_y = 1.0 #: Weight for the first derivative in the y direction + alpha_z = 1.0 #: Weight for the first derivative in the z direction + alpha_xx = 0.0 #: Weight for the second derivative in the x direction + alpha_yy = 0.0 #: Weight for the second derivative in the y direction + alpha_zz = 0.0 #: Weight for the second derivative in the z direction counter = None @@ -34,21 +117,42 @@ class Regularization(object): @property def Wx(self): if getattr(self, '_Wx', None) is None: - self._Wx = self.mesh.cellGradx*utils.sdiag(self.mesh.vol) + Ave_x_vol = self.mesh.aveCC2F[:self.mesh.nFv[0],:]*self.mesh.vol + self._Wx = utils.sdiag(Ave_x_vol**0.5)*self.mesh.cellGradx return self._Wx @property def Wy(self): if getattr(self, '_Wy', None) is None: - self._Wy = self.mesh.cellGrady*utils.sdiag(self.mesh.vol) + Ave_y_vol = self.mesh.aveCC2F[self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2]),:]*self.mesh.vol + self._Wy = utils.sdiag(Ave_y_vol**0.5)*self.mesh.cellGrady return self._Wy @property def Wz(self): if getattr(self, '_Wz', None) is None: - self._Wz = self.mesh.cellGradz*utils.sdiag(self.mesh.vol) + Ave_z_vol = self.mesh.aveCC2F[np.sum(self.mesh.nFv[:2]):,:]*self.mesh.vol + self._Wz = utils.sdiag(Ave_z_vol**0.5)*self.mesh.cellGradz return self._Wz + @property + def Wxx(self): + if getattr(self, '_Wxx', None) is None: + self._Wxx = self.mesh.faceDivx*self.mesh.cellGradx*utils.sdiag(self.mesh.vol) + return self._Wxx + + @property + def Wyy(self): + if getattr(self, '_Wyy', None) is None: + self._Wyy = self.mesh.faceDivy*self.mesh.cellGrady*utils.sdiag(self.mesh.vol) + return self._Wyy + + @property + def Wzz(self): + if getattr(self, '_Wzz', None) is None: + self._Wzz = self.mesh.faceDivz*self.mesh.cellGradz*utils.sdiag(self.mesh.vol) + return self._Wzz + def pnorm(self, r): return 0.5*r.dot(r) @@ -60,11 +164,14 @@ class Regularization(object): mobj = self.alpha_s * self.pnorm( self.Ws * mresid ) mobj += self.alpha_x * self.pnorm( self.Wx * mresid ) + mobj += self.alpha_xx * self.pnorm( self.Wxx * mresid ) if self.mesh.dim > 1: mobj += self.alpha_y * self.pnorm( self.Wy * mresid ) + mobj += self.alpha_yy * self.pnorm( self.Wyy * mresid ) if self.mesh.dim > 2: mobj += self.alpha_z * self.pnorm( self.Wz * mresid ) + mobj += self.alpha_zz * self.pnorm( self.Wzz * mresid ) return mobj @@ -95,11 +202,14 @@ class Regularization(object): mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid) mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid) + mobjDeriv = mobjDeriv + self.alpha_xx * self.Wxx.T * ( self.Wxx * mresid) if self.mesh.dim > 1: mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid) + mobjDeriv = mobjDeriv + self.alpha_yy * self.Wyy.T * ( self.Wyy * mresid) if self.mesh.dim > 2: mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid) + mobjDeriv = mobjDeriv + self.alpha_zz * self.Wzz.T * ( self.Wzz * mresid) return mobjDeriv @@ -110,11 +220,14 @@ class Regularization(object): mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx + mobj2Deriv = mobj2Deriv + self.alpha_xx * self.Wxx.T * self.Wxx if self.mesh.dim > 1: mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy + mobj2Deriv = mobj2Deriv + self.alpha_yy * self.Wyy.T * self.Wyy if self.mesh.dim > 2: mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz + mobj2Deriv = mobj2Deriv + self.alpha_zz * self.Wzz.T * self.Wzz return mobj2Deriv diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 70ce1d2e..cad7cf58 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -43,7 +43,7 @@ class SimPEGTable: def show(self): self.root.show() - def saveInversion(self, invObj, dataPath): + def saveInversion(self, invObj): # Create a new inversion anytime this is run. def _startup_hdf5_inv(invObj, m0): diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst index 8ca74e27..15c31e6d 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.regularization.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: From 21d0a772bcdd2fa79d9d294fe0ab1d49cc21e6f5 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 18 Dec 2013 10:35:59 -0800 Subject: [PATCH 33/35] Moved Regularization to inverse package. --- README.md | 2 +- SimPEG/__init__.py | 1 - SimPEG/examples/DC.py | 3 +-- SimPEG/examples/Linear.py | 4 ++-- .../{regularization => inverse}/Regularization.py | 14 +++++++------- SimPEG/inverse/__init__.py | 1 + SimPEG/mesh/InnerProducts.py | 8 ++++---- SimPEG/regularization/__init__.py | 1 - SimPEG/tests/test_forward_DCproblem.py | 3 +-- SimPEG/tests/test_forward_problem.py | 10 ++++------ docs/api_Optimize.rst | 2 +- 11 files changed, 22 insertions(+), 27 deletions(-) rename SimPEG/{regularization => inverse}/Regularization.py (92%) delete mode 100644 SimPEG/regularization/__init__.py 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 f455e866..f421bc82 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -5,7 +5,6 @@ from utils import Solver import mesh import data import forward -import regularization import inverse import visualize import examples diff --git a/SimPEG/examples/DC.py b/SimPEG/examples/DC.py index 074cdb8b..ad98d6c2 100644 --- a/SimPEG/examples/DC.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/examples/Linear.py b/SimPEG/examples/Linear.py index 79bec85c..dbee7dff 100644 --- a/SimPEG/examples/Linear.py +++ b/SimPEG/examples/Linear.py @@ -1,4 +1,4 @@ -from SimPEG import mesh, forward, inverse, regularization, np +from SimPEG import mesh, forward, inverse, np import matplotlib.pyplot as plt @@ -51,7 +51,7 @@ if __name__ == '__main__': prob, data = example(100) M = prob.mesh - reg = regularization.Regularization(M) + reg = inverse.Regularization(M) opt = inverse.InexactGaussNewton(maxIter=20) inv = inverse.Inversion(prob,reg,opt,data) m0 = np.zeros_like(data.mtrue) diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/inverse/Regularization.py similarity index 92% rename from SimPEG/regularization/Regularization.py rename to SimPEG/inverse/Regularization.py index a8406d9b..30c4b86d 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/inverse/Regularization.py @@ -111,46 +111,46 @@ class Regularization(object): @property def Ws(self): if getattr(self,'_Ws', None) is None: - self._Ws = utils.sdiag(self.mesh.vol) + self._Ws = utils.sdiag(self.mesh.vol**0.5) return self._Ws @property def Wx(self): if getattr(self, '_Wx', None) is None: - Ave_x_vol = self.mesh.aveCC2F[:self.mesh.nFv[0],:]*self.mesh.vol + Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFv[0]].T*self.mesh.vol self._Wx = utils.sdiag(Ave_x_vol**0.5)*self.mesh.cellGradx return self._Wx @property def Wy(self): if getattr(self, '_Wy', None) is None: - Ave_y_vol = self.mesh.aveCC2F[self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2]),:]*self.mesh.vol + 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**0.5)*self.mesh.cellGrady return self._Wy @property def Wz(self): if getattr(self, '_Wz', None) is None: - Ave_z_vol = self.mesh.aveCC2F[np.sum(self.mesh.nFv[:2]):,:]*self.mesh.vol + Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.nFv[:2]):].T*self.mesh.vol self._Wz = utils.sdiag(Ave_z_vol**0.5)*self.mesh.cellGradz return self._Wz @property def Wxx(self): if getattr(self, '_Wxx', None) is None: - self._Wxx = self.mesh.faceDivx*self.mesh.cellGradx*utils.sdiag(self.mesh.vol) + self._Wxx = utils.sdiag(self.mesh.vol**0.5)*self.mesh.faceDivx*self.mesh.cellGradx return self._Wxx @property def Wyy(self): if getattr(self, '_Wyy', None) is None: - self._Wyy = self.mesh.faceDivy*self.mesh.cellGrady*utils.sdiag(self.mesh.vol) + self._Wyy = utils.sdiag(self.mesh.vol**0.5)*self.mesh.faceDivy*self.mesh.cellGrady return self._Wyy @property def Wzz(self): if getattr(self, '_Wzz', None) is None: - self._Wzz = self.mesh.faceDivz*self.mesh.cellGradz*utils.sdiag(self.mesh.vol) + self._Wzz = utils.sdiag(self.mesh.vol**0.5)*self.mesh.faceDivz*self.mesh.cellGradz return self._Wzz 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/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/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 36216bf7..90312fdf 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -6,7 +6,6 @@ from SimPEG.forward import Problem from SimPEG.examples.DC import * from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve -from SimPEG.regularization import Regularization from SimPEG import inverse @@ -48,7 +47,7 @@ class DCProblemTests(unittest.TestCase): # Now set up the problem to do some minimization 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, data, beta0=1e4) self.inv = inv 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/docs/api_Optimize.rst b/docs/api_Optimize.rst index 15c31e6d..a30e55f7 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Optimize.rst @@ -29,7 +29,7 @@ Beta Schedule Regularization ************** -.. automodule:: SimPEG.regularization.Regularization +.. automodule:: SimPEG.inverse.Regularization :show-inheritance: :members: :undoc-members: From 43b81e1396fa88d894f71b2182fabac33ddad29b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 18 Dec 2013 11:11:32 -0800 Subject: [PATCH 34/35] Simplifications to Regularization code. --- SimPEG/inverse/Regularization.py | 120 ++++++++++++++----------------- SimPEG/utils/__init__.py | 9 +++ 2 files changed, 62 insertions(+), 67 deletions(-) diff --git a/SimPEG/inverse/Regularization.py b/SimPEG/inverse/Regularization.py index 30c4b86d..8c6f19ef 100644 --- a/SimPEG/inverse/Regularization.py +++ b/SimPEG/inverse/Regularization.py @@ -85,13 +85,13 @@ class Regularization(object): __metaclass__ = utils.Save.Savable - alpha_s = 1e-6 #: Smallness weight - alpha_x = 1.0 #: Weight for the first derivative in the x direction - alpha_y = 1.0 #: Weight for the first derivative in the y direction - alpha_z = 1.0 #: Weight for the first derivative in the z direction - alpha_xx = 0.0 #: Weight for the second derivative in the x direction - alpha_yy = 0.0 #: Weight for the second derivative in the y direction - alpha_zz = 0.0 #: Weight for the second derivative in the z direction + 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 @@ -110,124 +110,110 @@ class Regularization(object): @property def Ws(self): + """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = utils.sdiag(self.mesh.vol**0.5) + 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**0.5)*self.mesh.cellGradx + 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**0.5)*self.mesh.cellGrady + 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**0.5)*self.mesh.cellGradz + 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**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + 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**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + 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**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + self._Wzz = utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz return self._Wzz - def pnorm(self, r): - return 0.5*r.dot(r) + @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): - mresid = m - self.mref - - mobj = self.alpha_s * self.pnorm( self.Ws * mresid ) - - mobj += self.alpha_x * self.pnorm( self.Wx * mresid ) - mobj += self.alpha_xx * self.pnorm( self.Wxx * mresid ) - - if self.mesh.dim > 1: - mobj += self.alpha_y * self.pnorm( self.Wy * mresid ) - mobj += self.alpha_yy * self.pnorm( self.Wyy * mresid ) - if self.mesh.dim > 2: - mobj += self.alpha_z * self.pnorm( self.Wz * mresid ) - mobj += self.alpha_zz * self.pnorm( self.Wzz * mresid ) - - return mobj + r = self.W * (m - self.mref) + return 0.5*r.dot(r) @utils.timeIt def modelObjDeriv(self, m): """ - In 1D: + The regularization is: .. 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 + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} - \\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}}) + So the derivative is straight forward: + .. math:: - \\frac{ \partial^2 m_{\\text{obj}} }{\partial m^2} = - \\alpha_s W_s^{\\top} W_s + - \\alpha_x W_x^{\\top} W_x + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - - 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) - mobjDeriv = mobjDeriv + self.alpha_xx * self.Wxx.T * ( self.Wxx * mresid) - - if self.mesh.dim > 1: - mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid) - mobjDeriv = mobjDeriv + self.alpha_yy * self.Wyy.T * ( self.Wyy * mresid) - if self.mesh.dim > 2: - mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid) - mobjDeriv = mobjDeriv + self.alpha_zz * self.Wzz.T * ( self.Wzz * mresid) - - return mobjDeriv - + return self.W.T * ( self.W * (m - self.mref) ) @utils.timeIt def modelObj2Deriv(self): + """ - mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws + The regularization is: - mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx - mobj2Deriv = mobj2Deriv + self.alpha_xx * self.Wxx.T * self.Wxx + .. math:: - if self.mesh.dim > 1: - mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy - mobj2Deriv = mobj2Deriv + self.alpha_yy * self.Wyy.T * self.Wyy - if self.mesh.dim > 2: - mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz - mobj2Deriv = mobj2Deriv + self.alpha_zz * self.Wzz.T * self.Wzz + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} - return mobj2Deriv + So the second derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W} + + """ + return self.W.T * self.W diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index f75fed3d..8d9b0504 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -124,6 +124,15 @@ def callHooks(match): 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): """ From 53233e25e59c894b90603b17cb0e09d0baa3526e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 18 Dec 2013 11:19:42 -0800 Subject: [PATCH 35/35] documentation updates. --- SimPEG/utils/Save.py | 2 +- docs/index.rst | 21 +++++++++++++++------ 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index cad7cf58..4ad550de 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -266,7 +266,7 @@ def saveSavable(obj, group, debug=False): try: grp.attrs[name] = val except Exception, e: - print 'Warning: Could not save %s, problems may arise is loading.' % name + print 'Warning: Could not save %s, problems may arise in loading.' % name group.attrs['__class__'] = obj.__class__.__name__ for arg in obj._kwargs_init: 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