Merge pull request #94 from simpeg/TreeMesh_accel

Tree mesh acceleration

Author: lheagy

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.2.1
+current_version = 0.3.0
 files = setup.py discretize/__init__.py docs/conf.py
 

.gitignore

@@ -42,6 +42,8 @@ docs/warnings.txt
 .DS_Store
 docs/auto_examples
 docs/modules
+tree_ext.cpp
+interputils_cython.c
 
 # Jupyter
 *.ipynb

.travis.yml

@@ -34,6 +34,7 @@ install:
   - conda install --quiet --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib pillow cython ipython;
   - conda install --quiet --yes -c menpo vtk
   - pip install -r requirements_dev.txt
+  - python --version
   - python setup.py install
   - make build
   - export MPLBACKEND="agg"

appveyor.yml

@@ -4,56 +4,53 @@ build: off
 
 environment:
   matrix:
-    - PYTHON: 2.7
-      MINICONDA: C:\Miniconda
-      PYTHON_ARCH: 32
-    - PYTHON: 3.5
-      MINICONDA: C:\Miniconda3
-      PYTHON_ARCH: 32
-    - PYTHON: 3.6
-      MINICONDA: C:\Miniconda3
-      PYTHON_ARCH: 32
-    - PYTHON: 2.7
-      MINICONDA: C:\Miniconda
-      PYTHON_ARCH: 64
-    - PYTHON: 3.5
-      MINICONDA: C:\Miniconda3
-      PYTHON_ARCH: 64
-    - PYTHON: 3.6
-      MINICONDA: C:\Miniconda3
-      PYTHON_ARCH: 64
+    - PYTHON_VERSION: 2.7
+      CONDA: C:\Miniconda
+
+    - PYTHON_VERSION: 3.5
+      CONDA: C:\Miniconda35
+
+    - PYTHON_VERSION: 3.6
+      CONDA: C:\Miniconda36
 
   pypi_password:
     secure: S2mZCxm5LzIF4KFDss/f7g==
 
+image: Visual Studio 2015
+
+platform:
+- x86
+- x64
+
 init:
-  - "ECHO %PYTHON% %MINICONDA%"
+  - ps: if ($Env:PLATFORM -eq "x64") { $Env:CONDA = "${Env:CONDA}-x64" }
+  - ps: Write-Host $Env:PYTHON_VERSION
+  - ps: Write-Host $Env:CONDA
+  - ps: Write-Host $Env:PLATFORM
 
 install:
-  - "set PATH=%MINICONDA%;%MINICONDA%\\Scripts;%PATH%"
+  # Load the conda root environment, configure and install some packages
+  - '"%CONDA%\Scripts\activate.bat"'
   - conda config --set always_yes yes --set changeps1 no
   - conda update -q conda
   - conda info -a
-  - "conda create -q -n test-environment python=%PYTHON% numpy scipy matplotlib cython ipython pillow wheel"
-  - activate test-environment
-  - conda install --quiet --yes -c menpo vtk
+  - conda install numpy scipy matplotlib cython ipython pillow wheel
+  - conda install -c menpo vtk
   - pip install -r requirements_dev.txt
   - python setup.py install
 
-  # Use the Agg backend in Matplotlib
-  - "set MPLBACKEND=Agg"
-
 test_script:
   - cd tests
-  - nosetests base -v -s
+  - nosetests -v -s base/test_interpolation.py base/test_tensor_io.py
+  - nosetests -v -s tree
   - cd ..
 
 after_test:
   # This step builds your wheels.
   # Again, you only need build.cmd if you're building C extensions for
   # 64-bit Python 3.3/3.4. And you need to use %PYTHON% to get the correct
   # interpreter
-  - python setup.py sdist bdist_wheel
+  - python setup.py bdist_wheel
 
 artifacts:
   # bdist_wheel puts your built wheel in the dist directory
@@ -65,6 +62,6 @@ deploy_script:
   - echo [pypi] >> %USERPROFILE%\\.pypirc
   - echo username=lheagy >> %USERPROFILE%\\.pypirc
   - echo password=%pypi_password% >> %USERPROFILE%\\.pypirc
-  - twine upload dist\* --skip-existing
+  # - twine upload dist\* --skip-existing
   # deploy on dev or master
-  # - ps: If ($env:APPVEYOR_REPO_TAG -eq "true" ) { Invoke-Expression "twine upload dist/* --skip-existing" 2>$null } Else { write-output "Not on a tag, won't deploy to pypi"}
+  - ps: If ($env:APPVEYOR_REPO_TAG -eq "true" ) { Invoke-Expression "twine upload dist\*.whl --skip-existing" 2>$null } Else { write-output "Not on a tag, won't deploy to pypi"}

discretize/MeshIO.py

@@ -537,93 +537,19 @@ def writeUBC(mesh, fileName, models=None, directory='', comment_lines=''):
 
 class TreeMeshIO(object):
 
-    def writeUBC(mesh, fileName, models=None):
-        """Write UBC ocTree mesh and model files from a
-        octree mesh and model.
-
-        :param string fileName: File to write to
-        :param dict models: Models in a dict, where each key is the filename
-        """
-
-        # Calculate information to write in the file.
-        # Number of cells in the underlying mesh
-        nCunderMesh = np.array([h.size for h in mesh.h], dtype=np.int64)
-        # The top-south-west most corner of the mesh
-        tswCorn = mesh.x0 + np.array([0, 0, np.sum(mesh.h[2])])
-        # Smallest cell size
-        smallCell = np.array([h.min() for h in mesh.h])
-        # Number of cells
-        nrCells = mesh.nC
-
-        # Extract information about the cells.
-        cellPointers = np.array([c._pointer for c in mesh])
-        cellW = np.array([mesh._levelWidth(i) for i in cellPointers[:, -1]])
-        # Need to shift the pointers to work with UBC indexing
-        # UBC Octree indexes always the top-left-close (top-south-west) corner
-        # first and orders the cells in z(top-down), x, y vs x, y, z(bottom-up)
-        # Shift index up by 1
-        ubcCellPt = cellPointers[:, 0:-1].copy() + np.array([1., 1., 1.])
-        # Need re-index the z index to be from the top-left-close corner and to
-        # be from the global top.
-        ubcCellPt[:, 2] = (nCunderMesh[-1] + 2) - (ubcCellPt[:, 2] + cellW)
-
-        # Reorder the ubcCellPt
-        ubcReorder = np.argsort(
-            ubcCellPt.view(', '.join(3*['float'])),
-            axis=0,
-            order=['f2', 'f1', 'f0']
-        )[:, 0]
-        # Make a array with the pointers and the withs,
-        # that are order in the ubc ordering
-        indArr = np.concatenate(
-            (ubcCellPt[ubcReorder, :], cellW[ubcReorder].reshape((-1, 1))),
-            axis=1
-        )
-        # Write the UBC octree mesh file
-        head = (
-            '{:.0f} {:.0f} {:.0f}\n'.format(
-                nCunderMesh[0], nCunderMesh[1], nCunderMesh[2]
-            ) +
-            '{:.4f} {:.4f} {:.4f}\n'.format(
-                tswCorn[0], tswCorn[1], tswCorn[2]
-            ) +
-            '{:.3f} {:.3f} {:.3f}\n'.format(
-                smallCell[0], smallCell[1], smallCell[2]
-            ) +
-            '{:.0f}'.format(nrCells)
-        )
-        np.savetxt(fileName, indArr, fmt='%i', header=head, comments='')
-
-        # Print the models
-        # Assign the model('s) to the object
-        if models is not None:
-            for item in six.iteritems(models):
-                # Save the data
-                np.savetxt(item[0], item[1][ubcReorder], fmt='%3.5e')
-
     @classmethod
     def readUBC(TreeMesh, meshFile):
-        """Read UBC 3D OcTree mesh and/or modelFiles
-
+        """Read UBC 3D OcTree mesh file
         Input:
         :param str meshFile: path to the UBC GIF OcTree mesh file to read
         :rtype: discretize.TreeMesh
         :return: The octree mesh
         """
-
-        # Read the file lines
         fileLines = np.genfromtxt(meshFile, dtype=str,
-            delimiter='\n', comments='!')
-        # Extract the data
-        nCunderMesh = np.array(fileLines[0].
-            split('!')[0].split(), dtype=int)
-        # I think this is the case?
-        # Format of file changed... First 3 values are the # of cells in the
-        # underlying mesh and remaining 6 values are padding for the core region.
-        nCunderMesh  = nCunderMesh[0:3]
-
-        if np.unique(nCunderMesh).size >1:
-            raise Exception('TreeMeshes have the same number of cell in all directions')
+                                  delimiter='\n', comments='!')
+        nCunderMesh = np.array(fileLines[0].split('!')[0].split(), dtype=int)
+        nCunderMesh = nCunderMesh[0:3]
+
         tswCorn = np.array(
             fileLines[1].split('!')[0].split(),
             dtype=float
@@ -632,44 +558,37 @@ def readUBC(TreeMesh, meshFile):
             fileLines[2].split('!')[0].split(),
             dtype=float
         )
-        nrCells = np.array(
-            fileLines[3].split('!')[0].split(),
-            dtype=float
-        )
         # Read the index array
-        indArr = np.genfromtxt((line.encode('utf8') for line in fileLines[4::]), dtype=np.int)
+        indArr = np.genfromtxt((line.encode('utf8') for line in fileLines[4::]),
+                               dtype=np.int)
 
-        # Calculate simpeg parameters
         h1, h2, h3 = [np.ones(nr)*sz for nr, sz in zip(nCunderMesh, smallCell)]
         x0 = tswCorn - np.array([0, 0, np.sum(h3)])
-        # Convert the index array to a points list that complies with TreeMesh
-        # Shift to start at 0
-        simpegCellPt = indArr[:, 0:-1].copy()
-        simpegCellPt[:, 2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:, 2] + indArr[:, 3])
-        # Need reindex the z index to be from the bottom-left-close corner
-        # and to be from the global bottom.
-        simpegCellPt = simpegCellPt - np.array([1., 1., 1.])
-
-        # Calculate the cell level
-        simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:, 3])
-        # Make a pointer matrix
-        simpegPointers = np.concatenate((simpegCellPt, simpegLevel.reshape((-1, 1))), axis=1)
-
-        # Make the tree mesh
+
+        ls = np.log2(nCunderMesh).astype(int)
+        if ls[0] == ls[1] and ls[1] == ls[2]:
+            max_level = ls[0]
+        else:
+            max_level = min(ls)+1
+
         mesh = TreeMesh([h1, h2, h3], x0=x0)
-        mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()])
 
-        # Figure out the reordering
-        mesh._simpegReorderUBC = np.argsort(
-            np.array([mesh._index(i) for i in simpegPointers.tolist()])
-        )
-        # mesh._simpegReorderUBC = np.argsort((np.array([[1, 1, 1, -1]])*simpegPointers).view(', '.join(4*['float'])), axis=0, order=['f3', 'f2', 'f1', 'f0'])[:, 0]
+        # Convert indArr to points in coordinates of underlying cpp tree
+        # indArr is ix, iy, iz(top-down) need it in ix, iy, iz (bottom-up)
 
+        levels = indArr[:, -1]
+        indArr = indArr[:, :-1]
+
+        indArr -= 1  # shift by 1....
+        indArr = 2*indArr + levels[:, None]  # get cell center index
+        indArr[:, 2] = 2*nCunderMesh[2] - indArr[:, 2]  # switch direction of iz
+        levels = max_level-np.log2(levels)  # calculate level
+
+        mesh.__setstate__((indArr, levels))
         return mesh
 
     def readModelUBC(mesh, fileName):
         """Read UBC OcTree model and get vector
-
         :param string fileName: path to the UBC GIF model file to read
         :rtype: numpy.ndarray
         :return: OcTree model
@@ -681,49 +600,91 @@ def readModelUBC(mesh, fileName):
                 out[f] = mesh.readModelUBC(f)
             return out
 
-        assert hasattr(mesh, '_simpegReorderUBC'), 'The file must have been loaded from a UBC format.'
         assert mesh.dim == 3
 
         modList = []
         modArr = np.loadtxt(fileName)
-        if len(modArr.shape) == 1:
-            modList.append(modArr[mesh._simpegReorderUBC])
-        else:
-            modList.append(modArr[mesh._simpegReorderUBC, :])
-        return modList
 
-    def writeVTK(mesh, fileName, models=None):
+        ubc_order = mesh._ubc_order
+        # order_ubc will re-order from treemesh ordering to UBC ordering
+        # need the opposite operation
+        un_order = np.empty_like(ubc_order)
+        un_order[ubc_order] = np.arange(len(ubc_order))
+
+        model = modArr[un_order].copy() # ensure a contiguous array
+        return model
+
+    def writeUBC(mesh, fileName, models=None):
+        """Write UBC ocTree mesh and model files from a
+        octree mesh and model.
+        :param string fileName: File to write to
+        :param dict models: Models in a dict, where each key is the filename
+        """
+        uniform_hs = np.array([np.allclose(h, h[0]) for h in mesh.h])
+        if np.any(~uniform_hs):
+            raise Exception('UBC form does not support variable cell widths')
+        nCunderMesh = np.array([h.size for h in mesh.h], dtype=np.int64)
+        tswCorn = mesh.x0 + np.array([0, 0, np.sum(mesh.h[2])])
+        smallCell = np.array([h[0] for h in mesh.h])
+        nrCells = mesh.nC
+
+        indArr, levels = mesh._ubc_indArr
+        ubc_order = mesh._ubc_order
+
+        indArr = indArr[ubc_order]
+        levels = levels[ubc_order]
+
+        # Write the UBC octree mesh file
+        head = (
+            '{:.0f} {:.0f} {:.0f}\n'.format(
+                nCunderMesh[0], nCunderMesh[1], nCunderMesh[2]
+            ) +
+            '{:.4f} {:.4f} {:.4f}\n'.format(
+                tswCorn[0], tswCorn[1], tswCorn[2]
+            ) +
+            '{:.3f} {:.3f} {:.3f}\n'.format(
+                smallCell[0], smallCell[1], smallCell[2]
+            ) +
+            '{:.0f}'.format(nrCells)
+        )
+        np.savetxt(fileName, np.c_[indArr, levels], fmt='%i', header=head, comments='')
+
+        # Print the models
+        # Assign the model('s) to the object
+        if models is not None:
+            for item in six.iteritems(models):
+                # Save the data
+                np.savetxt(item[0], item[1][ubc_order], fmt='%3.5e')
+
+
+    def writeVTK(self, fileName, models=None):
         """Function to write a VTU file from a TreeMesh and model."""
         import vtk
         from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION
-        from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray
-
-        if str(type(mesh)).split()[-1][1:-2] not in 'discretize.TreeMesh.TreeMesh':
-            raise IOError('mesh is not a TreeMesh.')
+        from vtk.util.numpy_support import numpy_to_vtk
 
         # Make the data parts for the vtu object
         # Points
-        mesh.number()
-        ptsMat = mesh._gridN + mesh.x0
+        ptsMat = np.vstack((self.gridN, self.gridhN))
 
         vtkPts = vtk.vtkPoints()
         vtkPts.SetData(numpy_to_vtk(ptsMat, deep=True))
         # Cells
-        cellConn = np.array([c.nodes for c in mesh], dtype=np.int64)
+        cellArray = [c for c in self]
+        cellConn = np.array([cell.nodes for cell in cellArray])
 
-        cellsMat = np.concatenate((np.ones((cellConn.shape[0], 1), dtype=np.int64)*cellConn.shape[1], cellConn), axis=1).ravel()
+        cellsMat = np.concatenate((np.ones((cellConn.shape[0], 1))*cellConn.shape[1], cellConn), axis=1).ravel()
         cellsArr = vtk.vtkCellArray()
         cellsArr.SetNumberOfCells(cellConn.shape[0])
-        cellsArr.SetCells(cellConn.shape[0], numpy_to_vtkIdTypeArray(cellsMat, deep =True))
+        cellsArr.SetCells(cellConn.shape[0], numpy_to_vtk(cellsMat, deep=True, array_type=vtk.VTK_ID_TYPE))
 
         # Make the object
         vtuObj = vtk.vtkUnstructuredGrid()
         vtuObj.SetPoints(vtkPts)
         vtuObj.SetCells(vtk.VTK_VOXEL, cellsArr)
         # Add the level of refinement as a cell array
-        cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3, 2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())])
-        uniqueLevel, indLevel = np.unique(np.prod(cellSides, axis=1), return_inverse=True)
-        refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel, deep=1)
+        cell_levels = np.array([cell._level for cell in cellArray])
+        refineLevelArr = numpy_to_vtk(cell_levels, deep=1)
         refineLevelArr.SetName('octreeLevel')
         vtuObj.GetCellData().AddArray(refineLevelArr)
         # Assign the model('s) to the object