[Yt-svn] yt-commit r1531 - in trunk: tests yt yt/_amr_utils yt/extensions/volume_rendering yt/lagos

mturk at wrangler.dreamhost.com mturk at wrangler.dreamhost.com
Fri Nov 20 19:54:18 PST 2009


Author: mturk
Date: Fri Nov 20 19:54:16 2009
New Revision: 1531
URL: http://yt.enzotools.org/changeset/1531

Log:
Moving all of the utility functions to amr_utils.  Eases compilation on Kraken
for CNL.  Thanks to Britton for a lot of this...



Added:
   trunk/yt/_amr_utils/
   trunk/yt/_amr_utils/CICDeposit.pyx
   trunk/yt/_amr_utils/DepthFirstOctree.pyx
   trunk/yt/_amr_utils/FixedInterpolator.c
   trunk/yt/_amr_utils/FixedInterpolator.h
   trunk/yt/_amr_utils/Interpolators.pyx
   trunk/yt/_amr_utils/PointsInVolume.pyx
   trunk/yt/_amr_utils/RayIntegrators.pyx
   trunk/yt/_amr_utils/VolumeIntegrator.pyx
   trunk/yt/amr_utils.pyx
Removed:
   trunk/yt/extensions/volume_rendering/FixedInterpolator.c
   trunk/yt/extensions/volume_rendering/FixedInterpolator.h
   trunk/yt/extensions/volume_rendering/VolumeIntegrator.c
   trunk/yt/extensions/volume_rendering/VolumeIntegrator.pyx
   trunk/yt/lagos/DepthFirstOctree.c
   trunk/yt/lagos/DepthFirstOctree.pyx
   trunk/yt/lagos/Interpolators.c
   trunk/yt/lagos/Interpolators.pyx
   trunk/yt/lagos/PointsInVolume.c
   trunk/yt/lagos/PointsInVolume.pyx
Modified:
   trunk/tests/test_lagos.py
   trunk/yt/extensions/volume_rendering/__init__.py
   trunk/yt/extensions/volume_rendering/setup.py
   trunk/yt/lagos/BaseDataTypes.py
   trunk/yt/lagos/HelperFunctions.py
   trunk/yt/lagos/__init__.py
   trunk/yt/lagos/setup.py
   trunk/yt/setup.py

Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py	(original)
+++ trunk/tests/test_lagos.py	Fri Nov 20 19:54:16 2009
@@ -454,7 +454,7 @@
 class TestSliceDataType(DataTypeTestingBase, LagosTestingBase, unittest.TestCase):
     def setUp(self):
         DataTypeTestingBase.setUp(self)
-        self.data = self.hierarchy.slice(0,0.5)
+        self.data = self.hierarchy.slice(0,0.5, center=[0.5, 0.5, 0.5])
 
 class TestCuttingPlane(DataTypeTestingBase, LagosTestingBase, unittest.TestCase):
     def setUp(self):

Added: trunk/yt/_amr_utils/CICDeposit.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/CICDeposit.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,80 @@
+"""
+Simle integrators for the radiative transfer equation
+
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: CASA/University of Colorado
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2008 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def CICDeposit_3(np.ndarray[np.float64_t, ndim=1] posx,
+                 np.ndarray[np.float64_t, ndim=1] posy,
+                 np.ndarray[np.float64_t, ndim=1] posz,
+                 np.ndarray[np.float32_t, ndim=1] mass,
+                 np.int64_t npositions,
+                 np.ndarray[np.float32_t, ndim=3] field,
+                 np.ndarray[np.float64_t, ndim=1] leftEdge,
+                 np.ndarray[np.int32_t, ndim=1] gridDimension,
+                 np.float64_t cellSize):
+
+    cdef int i1, j1, k1, n
+    cdef double xpos, ypos, zpos
+    cdef double fact, edge0, edge1, edge2
+    cdef double le0, le1, le2
+    cdef float dx, dy, dz, dx2, dy2, dz2
+
+    edge0 = (<float> gridDimension[0]) - 0.5001
+    edge1 = (<float> gridDimension[1]) - 0.5001
+    edge2 = (<float> gridDimension[2]) - 0.5001
+    fact = 1.0 / cellSize
+
+    le0 = leftEdge[0]
+    le1 = leftEdge[1]
+    le2 = leftEdge[2]
+
+    for n in range(npositions):
+
+        # Compute the position of the central cell
+        xpos = fmin(fmax((posx[n] - le0)*fact, 0.5001), edge0)
+        ypos = fmin(fmax((posy[n] - le1)*fact, 0.5001), edge1)
+        zpos = fmin(fmax((posz[n] - le2)*fact, 0.5001), edge2)
+
+        i1  = <int> (xpos + 0.5)
+        j1  = <int> (ypos + 0.5)
+        k1  = <int> (zpos + 0.5)
+
+        # Compute the weights
+        dx = (<float> i1) + 0.5 - xpos
+        dy = (<float> j1) + 0.5 - ypos
+        dz = (<float> k1) + 0.5 - zpos
+        dx2 =  1.0 - dx
+        dy2 =  1.0 - dy
+        dz2 =  1.0 - dz
+
+        # Interpolate from field into sumfield
+        field[i1-1,j1-1,k1-1] += mass[n] * dx  * dy  * dz
+        field[i1  ,j1-1,k1-1] += mass[n] * dx2 * dy  * dz
+        field[i1-1,j1  ,k1-1] += mass[n] * dx  * dy2 * dz
+        field[i1  ,j1  ,k1-1] += mass[n] * dx2 * dy2 * dz
+        field[i1-1,j1-1,k1  ] += mass[n] * dx  * dy  * dz2
+        field[i1  ,j1-1,k1  ] += mass[n] * dx2 * dy  * dz2
+        field[i1-1,j1  ,k1  ] += mass[n] * dx  * dy2 * dz2
+        field[i1  ,j1  ,k1  ] += mass[n] * dx2 * dy2 * dz2

Added: trunk/yt/_amr_utils/DepthFirstOctree.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/DepthFirstOctree.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,167 @@
+"""
+This is a recursive function to return a depth-first octree
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2008 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef class position:
+    cdef public int output_pos, refined_pos
+    def __cinit__(self):
+        self.output_pos = 0
+        self.refined_pos = 0
+
+cdef class OctreeGrid:
+    cdef public object child_indices, fields, left_edges, dimensions, dx
+    cdef public int level
+    def __cinit__(self,
+                  np.ndarray[np.int32_t, ndim=3] child_indices,
+                  np.ndarray[np.float64_t, ndim=4] fields,
+                  np.ndarray[np.float64_t, ndim=1] left_edges,
+                  np.ndarray[np.int32_t, ndim=1] dimensions,
+                  np.ndarray[np.float64_t, ndim=1] dx,
+                  int level):
+        self.child_indices = child_indices
+        self.fields = fields
+        self.left_edges = left_edges
+        self.dimensions = dimensions
+        self.dx = dx
+        self.level = level
+
+cdef class OctreeGridList:
+    cdef public object grids
+    def __cinit__(self, grids):
+        self.grids = grids
+
+    def __getitem__(self, int item):
+        return self.grids[item]
+
+ at cython.boundscheck(False)
+def RecurseOctreeDepthFirst(int i_i, int j_i, int k_i,
+                            int i_f, int j_f, int k_f,
+                            position curpos, int gi, 
+                            np.ndarray[np.float64_t, ndim=2] output,
+                            np.ndarray[np.int32_t, ndim=1] refined,
+                            OctreeGridList grids):
+    cdef int i, i_off, j, j_off, k, k_off, ci, fi
+    cdef int child_i, child_j, child_k
+    cdef OctreeGrid child_grid
+    cdef OctreeGrid grid = grids[gi-1]
+    cdef np.ndarray[np.int32_t, ndim=3] child_indices = grid.child_indices
+    cdef np.ndarray[np.int32_t, ndim=1] dimensions = grid.dimensions
+    cdef np.ndarray[np.float64_t, ndim=4] fields = grid.fields
+    cdef np.ndarray[np.float64_t, ndim=1] leftedges = grid.left_edges
+    cdef np.float64_t dx = grid.dx[0]
+    cdef np.float64_t child_dx
+    cdef np.ndarray[np.float64_t, ndim=1] child_leftedges
+    cdef np.float64_t cx, cy, cz
+    for k_off in range(k_f):
+        k = k_off + k_i
+        cz = (leftedges[2] + k*dx)
+        for j_off in range(j_f):
+            j = j_off + j_i
+            cy = (leftedges[1] + j*dx)
+            for i_off in range(i_f):
+                i = i_off + i_i
+                cx = (leftedges[0] + i*dx)
+                ci = grid.child_indices[i,j,k]
+                if ci == -1:
+                    for fi in range(fields.shape[0]):
+                        output[curpos.output_pos,fi] = fields[fi,i,j,k]
+                    refined[curpos.refined_pos] = 0
+                    curpos.output_pos += 1
+                    curpos.refined_pos += 1
+                else:
+                    refined[curpos.refined_pos] = 1
+                    curpos.refined_pos += 1
+                    child_grid = grids[ci-1]
+                    child_dx = child_grid.dx[0]
+                    child_leftedges = child_grid.left_edges
+                    child_i = int((cx - child_leftedges[0])/child_dx)
+                    child_j = int((cy - child_leftedges[1])/child_dx)
+                    child_k = int((cz - child_leftedges[2])/child_dx)
+                    s = RecurseOctreeDepthFirst(child_i, child_j, child_k, 2, 2, 2,
+                                        curpos, ci, output, refined, grids)
+    return s
+
+ at cython.boundscheck(False)
+def RecurseOctreeByLevels(int i_i, int j_i, int k_i,
+                          int i_f, int j_f, int k_f,
+                          np.ndarray[np.int32_t, ndim=1] curpos,
+                          int gi, 
+                          np.ndarray[np.float64_t, ndim=2] output,
+                          np.ndarray[np.int32_t, ndim=2] genealogy,
+                          np.ndarray[np.float64_t, ndim=2] corners,
+                          OctreeGridList grids):
+    cdef np.int32_t i, i_off, j, j_off, k, k_off, ci, fi
+    cdef int child_i, child_j, child_k
+    cdef OctreeGrid child_grid
+    cdef OctreeGrid grid = grids[gi-1]
+    cdef int level = grid.level
+    cdef np.ndarray[np.int32_t, ndim=3] child_indices = grid.child_indices
+    cdef np.ndarray[np.int32_t, ndim=1] dimensions = grid.dimensions
+    cdef np.ndarray[np.float64_t, ndim=4] fields = grid.fields
+    cdef np.ndarray[np.float64_t, ndim=1] leftedges = grid.left_edges
+    cdef np.float64_t dx = grid.dx[0]
+    cdef np.float64_t child_dx
+    cdef np.ndarray[np.float64_t, ndim=1] child_leftedges
+    cdef np.float64_t cx, cy, cz
+    cdef int cp
+    for i_off in range(i_f):
+        i = i_off + i_i
+        cx = (leftedges[0] + i*dx)
+        if i_f > 2: print k, cz
+        for j_off in range(j_f):
+            j = j_off + j_i
+            cy = (leftedges[1] + j*dx)
+            for k_off in range(k_f):
+                k = k_off + k_i
+                cz = (leftedges[2] + k*dx)
+                cp = curpos[level]
+                corners[cp, 0] = cx 
+                corners[cp, 1] = cy 
+                corners[cp, 2] = cz
+                genealogy[curpos[level], 2] = level
+                # always output data
+                for fi in range(fields.shape[0]):
+                    output[cp,fi] = fields[fi,i,j,k]
+                ci = child_indices[i,j,k]
+                if ci > -1:
+                    child_grid = grids[ci-1]
+                    child_dx = child_grid.dx[0]
+                    child_leftedges = child_grid.left_edges
+                    child_i = int((cx-child_leftedges[0])/child_dx)
+                    child_j = int((cy-child_leftedges[1])/child_dx)
+                    child_k = int((cz-child_leftedges[2])/child_dx)
+                    # set current child id to id of next cell to examine
+                    genealogy[cp, 0] = curpos[level+1] 
+                    # set next parent id to id of current cell
+                    genealogy[curpos[level+1]:curpos[level+1]+8, 1] = cp
+                    s = RecurseOctreeByLevels(child_i, child_j, child_k, 2, 2, 2,
+                                              curpos, ci, output, genealogy,
+                                              corners, grids)
+                curpos[level] += 1
+    return s
+

Added: trunk/yt/_amr_utils/FixedInterpolator.c
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/FixedInterpolator.c	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,74 @@
+/************************************************************************
+* Copyright (C) 2009 Matthew Turk.  All Rights Reserved.
+*
+* This file is part of yt.
+*
+* yt is free software; you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation; either version 3 of the License, or
+* (at your option) any later version.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+************************************************************************/
+
+
+//
+// A small, tiny, itty bitty module for computation-intensive interpolation
+// that I can't seem to make fast in Cython
+//
+
+#include "FixedInterpolator.h"
+
+#define VINDEX(A,B,C) data[((((A)+ci[0])*(ds[1]+1)+((B)+ci[1]))*(ds[2]+1)+ci[2]+(C))]
+//  (((C*ds[1])+B)*ds[0]+A)
+
+npy_float64 fast_interpolate(int *ds, int *ci, npy_float64 *dp,
+                             npy_float64 *data)
+{
+    int i;
+    npy_float64 dv, dm[3];
+    for(i=0;i<3;i++)dm[i] = (1.0 - dp[i]);
+    dv  = 0.0;
+    dv += VINDEX(0,0,0) * (dm[0]*dm[1]*dm[2]);
+    dv += VINDEX(0,0,1) * (dm[0]*dm[1]*dp[2]);
+    dv += VINDEX(0,1,0) * (dm[0]*dp[1]*dm[2]);
+    dv += VINDEX(0,1,1) * (dm[0]*dp[1]*dp[2]);
+    dv += VINDEX(1,0,0) * (dp[0]*dm[1]*dm[2]);
+    dv += VINDEX(1,0,1) * (dp[0]*dm[1]*dp[2]);
+    dv += VINDEX(1,1,0) * (dp[0]*dp[1]*dm[2]);
+    dv += VINDEX(1,1,1) * (dp[0]*dp[1]*dp[2]);
+    /*assert(dv < -20);*/
+    return dv;
+}
+
+npy_float64 trilinear_interpolate(int *ds, int *ci, npy_float64 *dp,
+				  npy_float64 *data)
+{
+    /* dims is one less than the dimensions of the array */
+    int i;
+    npy_float64 dm[3], vz[4];
+  //dp is the distance to the plane.  dm is val, dp = 1-val
+    for(i=0;i<3;i++)dm[i] = (1.0 - dp[i]);
+    
+  //First interpolate in z
+    vz[0] = dm[2]*VINDEX(0,0,0) + dp[2]*VINDEX(0,0,1);
+    vz[1] = dm[2]*VINDEX(0,1,0) + dp[2]*VINDEX(0,1,1);
+    vz[2] = dm[2]*VINDEX(1,0,0) + dp[2]*VINDEX(1,0,1);
+    vz[3] = dm[2]*VINDEX(1,1,0) + dp[2]*VINDEX(1,1,1);
+
+  //Then in y
+    vz[0] = dm[1]*vz[0] + dp[1]*vz[1];
+    vz[1] = dm[1]*vz[2] + dp[1]*vz[3];
+
+  //Then in x
+    vz[0] = dm[0]*vz[0] + dp[0]*vz[1];
+    /*assert(dv < -20);*/
+    return vz[0];
+}

Added: trunk/yt/_amr_utils/FixedInterpolator.h
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/FixedInterpolator.h	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,40 @@
+/************************************************************************
+* Copyright (C) 2009 Matthew Turk.  All Rights Reserved.
+*
+* This file is part of yt.
+*
+* yt is free software; you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation; either version 3 of the License, or
+* (at your option) any later version.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+************************************************************************/
+
+
+//
+// A small, tiny, itty bitty module for computation-intensive interpolation
+// that I can't seem to make fast in Cython
+//
+
+#include "Python.h"
+
+#include <stdio.h>
+#include <math.h>
+#include <signal.h>
+#include <ctype.h>
+
+#include "numpy/ndarrayobject.h"
+
+npy_float64 fast_interpolate(int ds[3], int ci[3], npy_float64 dp[3],
+                             npy_float64 *data);
+
+npy_float64 trilinear_interpolate(int ds[3], int ci[3], npy_float64 dp[3],
+				  npy_float64 *data);

Added: trunk/yt/_amr_utils/Interpolators.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/Interpolators.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,116 @@
+"""
+Simple interpolators
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2008 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+ at cython.boundscheck(False)
+def UnilinearlyInterpolate(np.ndarray[np.float64_t, ndim=1] table,
+                           np.ndarray[np.float64_t, ndim=1] x_vals,
+                           np.ndarray[np.float64_t, ndim=1] x_bins,
+                           np.ndarray[np.int32_t, ndim=1] x_is,
+                           np.ndarray[np.float64_t, ndim=1] output):
+    cdef double x, xp, xm
+    cdef int i, x_i, y_i
+    for i in range(x_vals.shape[0]):
+        x_i = x_is[i]
+        x = x_vals[i]
+        dx_inv = 1.0 / (x_bins[x_i+1] - x_bins[x_i])
+        xp = (x - x_bins[x_i]) * dx_inv
+        xm = (x_bins[x_i+1] - x) * dx_inv
+        output[i]  = table[x_i  ] * (xm) \
+                   + table[x_i+1] * (xp)
+
+ at cython.boundscheck(False)
+def BilinearlyInterpolate(np.ndarray[np.float64_t, ndim=2] table,
+                          np.ndarray[np.float64_t, ndim=1] x_vals,
+                          np.ndarray[np.float64_t, ndim=1] y_vals,
+                          np.ndarray[np.float64_t, ndim=1] x_bins,
+                          np.ndarray[np.float64_t, ndim=1] y_bins,
+                          np.ndarray[np.int32_t, ndim=1] x_is,
+                          np.ndarray[np.int32_t, ndim=1] y_is,
+                          np.ndarray[np.float64_t, ndim=1] output):
+    cdef double x, xp, xm
+    cdef double y, yp, ym
+    cdef double dx_inv, dy_inv
+    cdef int i, x_i, y_i
+    for i in range(x_vals.shape[0]):
+        x_i = x_is[i]
+        y_i = y_is[i]
+        x = x_vals[i]
+        y = y_vals[i]
+        dx_inv = 1.0 / (x_bins[x_i+1] - x_bins[x_i])
+        dy_inv = 1.0 / (y_bins[y_i+1] - y_bins[y_i])
+        xp = (x - x_bins[x_i]) * dx_inv
+        yp = (y - y_bins[y_i]) * dy_inv
+        xm = (x_bins[x_i+1] - x) * dx_inv
+        ym = (y_bins[y_i+1] - y) * dy_inv
+        output[i]  = table[x_i  , y_i  ] * (xm*ym) \
+                   + table[x_i+1, y_i  ] * (xp*ym) \
+                   + table[x_i  , y_i+1] * (xm*yp) \
+                   + table[x_i+1, y_i+1] * (xp*yp)
+
+ at cython.boundscheck(False)
+def TrilinearlyInterpolate(np.ndarray[np.float64_t, ndim=3] table,
+                           np.ndarray[np.float64_t, ndim=1] x_vals,
+                           np.ndarray[np.float64_t, ndim=1] y_vals,
+                           np.ndarray[np.float64_t, ndim=1] z_vals,
+                           np.ndarray[np.float64_t, ndim=1] x_bins,
+                           np.ndarray[np.float64_t, ndim=1] y_bins,
+                           np.ndarray[np.float64_t, ndim=1] z_bins,
+                           np.ndarray[np.int_t, ndim=1] x_is,
+                           np.ndarray[np.int_t, ndim=1] y_is,
+                           np.ndarray[np.int_t, ndim=1] z_is,
+                           np.ndarray[np.float64_t, ndim=1] output):
+    cdef double x, xp, xm
+    cdef double y, yp, ym
+    cdef double z, zp, zm
+    cdef double dx_inv, dy_inv, dz_inv
+    cdef int i, x_i, y_i, z_i
+    for i in range(x_vals.shape[0]):
+        x_i = x_is[i]
+        y_i = y_is[i]
+        z_i = z_is[i]
+        x = x_vals[i]
+        y = y_vals[i]
+        z = z_vals[i]
+        dx_inv = 1.0 / (x_bins[x_i+1] - x_bins[x_i])
+        dy_inv = 1.0 / (y_bins[y_i+1] - y_bins[y_i])
+        dz_inv = 1.0 / (z_bins[z_i+1] - z_bins[z_i])
+        xp = (x - x_bins[x_i]) * dx_inv
+        yp = (y - y_bins[y_i]) * dy_inv
+        zp = (z - z_bins[z_i]) * dz_inv
+        xm = (x_bins[x_i+1] - x) * dx_inv
+        ym = (y_bins[y_i+1] - y) * dy_inv
+        zm = (z_bins[z_i+1] - z) * dz_inv
+        output[i]  = table[x_i  ,y_i  ,z_i  ] * (xm*ym*zm) \
+                   + table[x_i+1,y_i  ,z_i  ] * (xp*ym*zm) \
+                   + table[x_i  ,y_i+1,z_i  ] * (xm*yp*zm) \
+                   + table[x_i  ,y_i  ,z_i+1] * (xm*ym*zp) \
+                   + table[x_i+1,y_i  ,z_i+1] * (xp*ym*zp) \
+                   + table[x_i  ,y_i+1,z_i+1] * (xm*yp*zp) \
+                   + table[x_i+1,y_i+1,z_i  ] * (xp*yp*zm) \
+                   + table[x_i+1,y_i+1,z_i+1] * (xp*yp*zp)

Added: trunk/yt/_amr_utils/PointsInVolume.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/PointsInVolume.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,73 @@
+"""
+Checks for points contained in a volume
+
+Author: John Wise <jwise77 at gmail.com>
+Affiliation: Princeton
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2009 John.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def PointsInVolume(np.ndarray[np.float64_t, ndim=2] points,
+                   np.ndarray[np.int8_t, ndim=1] pmask,  # pixel mask
+                   np.ndarray[np.float64_t, ndim=1] left_edge,
+                   np.ndarray[np.float64_t, ndim=1] right_edge,
+                   np.ndarray[np.int32_t, ndim=3] mask,
+                   float dx):
+    cdef np.ndarray[np.int8_t, ndim=1] \
+         valid = np.zeros(points.shape[0], dtype='int8')
+    cdef int i, dim, count
+    cdef int ex
+    cdef double dx_inv
+    cdef unsigned int idx[3]
+    count = 0
+    dx_inv = 1.0 / dx
+    for i in xrange(points.shape[0]):
+        if pmask[i] == 0:
+            continue
+        ex = 1
+        for dim in xrange(3):
+            if points[i,dim] < left_edge[dim] or points[i,dim] > right_edge[dim]:
+                valid[i] = ex = 0
+                break
+        if ex == 1:
+            for dim in xrange(3):
+                idx[dim] = <unsigned int> \
+                           ((points[i,dim] - left_edge[dim]) * dx_inv)
+            if mask[idx[0], idx[1], idx[2]] == 1:
+                valid[i] = 1
+                count += 1
+
+    cdef np.ndarray[np.int32_t, ndim=1] result = np.empty(count, dtype='int32')
+    count = 0
+    for i in xrange(points.shape[0]):
+        if valid[i] == 1 and pmask[i] == 1:
+            result[count] = i
+            count += 1
+        
+    return result
+
+                
+                 

Added: trunk/yt/_amr_utils/RayIntegrators.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/RayIntegrators.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,374 @@
+"""
+Simle integrators for the radiative transfer equation
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2008 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from stdlib cimport malloc, free, abs
+
+cdef extern from "math.h":
+    double exp(double x)
+    float expf(float x)
+
+ at cython.boundscheck(False)
+def Transfer3D(np.ndarray[np.float64_t, ndim=3] i_s,
+               np.ndarray[np.float64_t, ndim=4] o_s,
+               np.ndarray[np.float64_t, ndim=4] e,
+               np.ndarray[np.float64_t, ndim=4] a,
+               int imin, int imax, int jmin, int jmax,
+               int kmin, int kmax, int istride, int jstride,
+               np.float64_t dx):
+    """
+    This function accepts an incoming slab (*i_s*), a buffer
+    for an outgoing set of values at every point in the grid (*o_s*),
+    an emission array (*e*), an absorption array (*a*), and dimensions of
+    the grid (*imin*, *imax*, *jmin*, *jmax*, *kmin*, *kmax*) as well
+    as strides in the *i* and *j* directions, and a *dx* of the grid being
+    integrated.
+    """
+    cdef int i, ii
+    cdef int j, jj
+    cdef int k, kk
+    cdef int n, nn
+    nn = o_s.shape[3] # This might be slow
+    cdef np.float64_t *temp = <np.float64_t *>malloc(sizeof(np.float64_t) * nn)
+    for i in range((imax-imin)*istride):
+        ii = i + imin*istride
+        for j in range((jmax-jmin)*jstride):
+            jj = j + jmin*jstride
+            # Not sure about the ordering of the loops here
+            for n in range(nn):
+                temp[n] = i_s[ii,jj,n] 
+            for k in range(kmax-kmin):
+                kk = k + kmin#*kstride, which doesn't make any sense
+                for n in range(nn):
+                    o_s[i,j,k,n] = temp[n] + dx*(e[i,j,k,n] - temp[n]*a[i,j,k,n])
+                    temp[n] = o_s[i,j,k,n]
+            for n in range(nn):
+                i_s[ii,jj,n] = temp[n]
+    free(temp)
+
+ at cython.boundscheck(False)
+def TransferShells(np.ndarray[np.float64_t, ndim=3] i_s,
+                   np.ndarray[np.float64_t, ndim=3] data,
+                   np.ndarray[np.float64_t, ndim=2] shells):
+    """
+    This function accepts an incoming slab (*i_s*), a buffer of *data*,
+    and a list of shells specified as [ (value, tolerance, r, g, b), ... ].
+    """
+    cdef int i, ii
+    cdef int j, jj
+    cdef int k, kk
+    cdef int n, nn
+    cdef np.float64_t dist
+    ii = data.shape[0]
+    jj = data.shape[1]
+    kk = data.shape[2]
+    nn = shells.shape[0]
+    cdef float rgba[4]
+    cdef float alpha
+    for i in range(ii):
+        for j in range(jj):
+            # Not sure about the ordering of the loops here
+            for k in range(kk):
+                for n in range(nn):
+                    dist = shells[n, 0] - data[i,j,k]
+                    if dist < 0: dist *= -1.0
+                    if dist < shells[n,1]:
+                        dist = exp(-dist/8.0)
+                        rgba[0] = shells[n,2]
+                        rgba[1] = shells[n,3]
+                        rgba[2] = shells[n,4]
+                        rgba[3] = shells[n,5]
+                        alpha = i_s[i,j,3]
+                        dist *= dist # This might improve appearance
+                        i_s[i,j,0] += (1.0 - alpha)*rgba[0]*dist*rgba[3]
+                        i_s[i,j,1] += (1.0 - alpha)*rgba[1]*dist*rgba[3]
+                        i_s[i,j,2] += (1.0 - alpha)*rgba[2]*dist*rgba[3]
+                        i_s[i,j,3] += (1.0 - alpha)*rgba[3]*dist*rgba[3]
+                        break
+
+ at cython.boundscheck(False)
+def Transfer1D(float i_s,
+               np.ndarray[np.float_t, ndim=1] o_s,
+               np.ndarray[np.float_t, ndim=1] e,
+               np.ndarray[np.float_t, ndim=1] a,
+               np.ndarray[np.float_t, ndim=1] dx,
+               int imin, int imax):
+    cdef int i
+    for i in range(imin, imax):
+        o_s[i] = i_s + dx[i]*(e[i] - i_s*a[i])
+        i_s = o_s[i]
+    return i_s
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def VoxelTraversal(np.ndarray[np.int_t, ndim=3] grid_mask,
+                   np.ndarray[np.float64_t, ndim=3] grid_t,
+                   np.ndarray[np.float64_t, ndim=3] grid_dt,
+                   np.ndarray[np.float64_t, ndim=1] left_edge,
+                   np.ndarray[np.float64_t, ndim=1] right_edge,
+                   np.ndarray[np.float64_t, ndim=1] dx,
+                   np.ndarray[np.float64_t, ndim=1] u,
+                   np.ndarray[np.float64_t, ndim=1] v):
+    # We're roughly following Amanatides & Woo
+    # Find the first place the ray hits the grid on its path
+    # Do left edge then right edge in each dim
+    cdef int i, x, y
+    cdef np.float64_t tl, tr, intersect_t, enter_t, exit_t, dt_tolerance
+    cdef np.ndarray[np.int64_t,   ndim=1] step = np.empty((3,), dtype=np.int64)
+    cdef np.ndarray[np.int64_t,   ndim=1] cur_ind = np.empty((3,), dtype=np.int64)
+    cdef np.ndarray[np.float64_t, ndim=1] tdelta = np.empty((3,), dtype=np.float64)
+    cdef np.ndarray[np.float64_t, ndim=1] tmax = np.empty((3,), dtype=np.float64)
+    cdef np.ndarray[np.float64_t, ndim=1] intersect = np.empty((3,), dtype=np.float64)
+    intersect_t = 1
+    dt_tolerance = 1e-6
+    # recall p = v * t + u
+    #  where p is position, v is our vector, u is the start point
+    for i in range(3):
+        # As long as we're iterating, set some other stuff, too
+        if(v[i] < 0): step[i] = -1
+        else: step[i] = 1
+        x = (i+1)%3
+        y = (i+2)%3
+        tl = (left_edge[i] - u[i])/v[i]
+        tr = (right_edge[i] - u[i])/v[i]
+        if (left_edge[x] <= (u[x] + tl*v[x]) <= right_edge[x]) and \
+           (left_edge[y] <= (u[y] + tl*v[y]) <= right_edge[y]) and \
+           (0.0 <= tl < intersect_t):
+            intersect_t = tl
+        if (left_edge[x] <= (u[x] + tr*v[x]) <= right_edge[x]) and \
+           (left_edge[y] <= (u[y] + tr*v[y]) <= right_edge[y]) and \
+           (0.0 <= tr < intersect_t):
+            intersect_t = tr
+    # if fully enclosed
+    if (left_edge[0] <= u[0] <= right_edge[0]) and \
+       (left_edge[1] <= u[1] <= right_edge[1]) and \
+       (left_edge[2] <= u[2] <= right_edge[2]):
+        intersect_t = 0.0
+    if not (0 <= intersect_t <= 1): return
+    # Now get the indices of the intersection
+    intersect = u + intersect_t * v
+    for i in range(3):
+        cur_ind[i] = np.floor((intersect[i] + 1e-8*dx[i] - left_edge[i])/dx[i])
+        tmax[i] = (((cur_ind[i]+step[i])*dx[i])+left_edge[i]-u[i])/v[i]
+        if cur_ind[i] == grid_mask.shape[i] and step[i] < 0:
+            cur_ind[i] = grid_mask.shape[i] - 1
+        if step[i] > 0: tmax[i] = (((cur_ind[i]+1)*dx[i])+left_edge[i]-u[i])/v[i]
+        if step[i] < 0: tmax[i] = (((cur_ind[i]+0)*dx[i])+left_edge[i]-u[i])/v[i]
+        tdelta[i] = (dx[i]/v[i])
+        if tdelta[i] < 0: tdelta[i] *= -1
+    # The variable intersect contains the point we first pierce the grid
+    enter_t = intersect_t
+    while 1:
+        if (not (0 <= cur_ind[0] < grid_mask.shape[0])) or \
+           (not (0 <= cur_ind[1] < grid_mask.shape[1])) or \
+           (not (0 <= cur_ind[2] < grid_mask.shape[2])):
+            break
+        # Note that we are calculating t on the fly, but we get *negative* t
+        # values from what they should be.
+        # If we've reached t = 1, we are done.
+        grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1
+        if (tmax[0] > 1.0) and (tmax[1] > 1.0) and (tmax[2] > 1.0):
+            grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = 1.0 - enter_t
+            break
+        if tmax[0] < tmax[1]:
+            if tmax[0] < tmax[2]:
+                grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+                grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[0] - enter_t
+                enter_t = tmax[0]
+                tmax[0] += tdelta[0]
+                cur_ind[0] += step[0]
+            else:
+                grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+                grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t
+                enter_t = tmax[2]
+                tmax[2] += tdelta[2]
+                cur_ind[2] += step[2]
+        else:
+            if tmax[1] < tmax[2]:
+                grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+                grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[1] - enter_t
+                enter_t = tmax[1]
+                tmax[1] += tdelta[1]
+                cur_ind[1] += step[1]
+            else:
+                grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+                grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t
+                enter_t = tmax[2]
+                tmax[2] += tdelta[2]
+                cur_ind[2] += step[2]
+    return
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def PlaneVoxelIntegration(np.ndarray[np.float64_t, ndim=1] left_edge,
+                          np.ndarray[np.float64_t, ndim=1] right_edge,
+                          np.ndarray[np.float64_t, ndim=1] dx,
+                          np.ndarray[np.float64_t, ndim=2] ug,
+                          np.ndarray[np.float64_t, ndim=1] v,
+                          np.ndarray[np.float64_t, ndim=2] image,
+                          np.ndarray[np.float64_t, ndim=3] data,
+                          np.ndarray[np.float64_t, ndim=2] shells):
+    # We're roughly following Amanatides & Woo on a ray-by-ray basis
+    # Note that for now it's just shells, but this can and should be
+    # generalized to transfer functions
+    cdef int i, x, y, vi
+    intersect_t = 1
+    dt_tolerance = 1e-6
+    cdef int nv = ug.shape[0]
+    cdef int nshells = shells.shape[0]
+    cdef np.ndarray[np.float64_t, ndim=1] u = np.empty((3,), dtype=np.float64)
+    # Copy things into temporary location for passing between functions
+    for vi in range(nv):
+        for i in range(3): u[i] = ug[vi, i]
+        integrate_ray(u, v, left_edge, right_edge, dx, 
+                      nshells, vi, data, shells, image)
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def integrate_ray(np.ndarray[np.float64_t, ndim=1] u,
+                  np.ndarray[np.float64_t, ndim=1] v,
+                  np.ndarray[np.float64_t, ndim=1] left_edge,
+                  np.ndarray[np.float64_t, ndim=1] right_edge,
+                  np.ndarray[np.float64_t, ndim=1] dx,
+                  int nshells, int ind,
+                  np.ndarray[np.float64_t, ndim=3] data,
+                  np.ndarray[np.float64_t, ndim=2] shells,
+                  np.ndarray[np.float64_t, ndim=2] image):
+    cdef int step[3], x, y, i, n
+    cdef np.float64_t intersect_t = 1
+    cdef np.float64_t dt_tolerance = 1e-6
+    cdef np.float64_t tl, tr, enter_t, exit_t
+    cdef np.int64_t cur_ind[3]
+    cdef np.float64_t tdelta[3]
+    cdef np.float64_t tmax[3]
+    cdef np.float64_t intersect[3]
+    cdef np.float64_t dt, dv
+    cdef np.float64_t dist, alpha
+    cdef np.float64_t one = 1.0
+    cdef int dims[3]
+    cdef np.float64_t rgba[4], temp_x, temp_y
+    for i in range(3):
+        # As long as we're iterating, set some other stuff, too
+        dims[i] = data.shape[i]
+        if(v[i] < 0): step[i] = -1
+        else: step[i] = 1
+        x = (i+1)%3
+        y = (i+2)%3
+        tl = (left_edge[i] - u[i])/v[i]
+        tr = (right_edge[i] - u[i])/v[i]
+        temp_x = (u[x] + tl*v[x])
+        temp_y = (u[y] + tl*v[y])
+        if (left_edge[x] <= temp_x) and (temp_x <= right_edge[x]) and \
+           (left_edge[y] <= temp_y) and (temp_y <= right_edge[y]) and \
+           (0.0 <= tl) and (tl < intersect_t):
+            intersect_t = tl
+        temp_x = (u[x] + tr*v[x])
+        temp_y = (u[y] + tr*v[y])
+        if (left_edge[x] <= temp_x) and (temp_x <= right_edge[x]) and \
+           (left_edge[y] <= temp_y) and (temp_y <= right_edge[y]) and \
+           (0.0 <= tr) and (tr < intersect_t):
+            intersect_t = tr
+    # if fully enclosed
+    if (left_edge[0] <= u[0] <= right_edge[0]) and \
+       (left_edge[1] <= u[1] <= right_edge[1]) and \
+       (left_edge[2] <= u[2] <= right_edge[2]):
+        intersect_t = 0.0
+    if not (0 <= intersect_t <= 1):
+        #print "Returning: intersect_t ==", intersect_t
+        return
+    # Now get the indices of the intersection
+    for i in range(3): intersect[i] = u[i] + intersect_t * v[i]
+    cdef int ncells = 0
+    for i in range(3):
+        cur_ind[i] = np.floor((intersect[i] + 1e-8*dx[i] - left_edge[i])/dx[i])
+        tmax[i] = (((cur_ind[i]+step[i])*dx[i])+left_edge[i]-u[i])/v[i]
+        if cur_ind[i] == dims[i] and step[i] < 0:
+            cur_ind[i] = dims[i] - 1
+        if step[i] > 0: tmax[i] = (((cur_ind[i]+1)*dx[i])+left_edge[i]-u[i])/v[i]
+        if step[i] < 0: tmax[i] = (((cur_ind[i]+0)*dx[i])+left_edge[i]-u[i])/v[i]
+        tdelta[i] = (dx[i]/v[i])
+        if tdelta[i] < 0: tdelta[i] *= -1
+    # The variable intersect contains the point we first pierce the grid
+    enter_t = intersect_t
+    if (not (0 <= cur_ind[0] < dims[0])) or \
+       (not (0 <= cur_ind[1] < dims[1])) or \
+       (not (0 <= cur_ind[2] < dims[2])):
+        #print "Returning: cur_ind", cur_ind[0], cur_ind[1], cur_ind[2]
+        #print "  dims:     ", dims[0], dims[1], dims[2]
+        #print "  intersect:",  intersect[0], intersect[1], intersect[2]
+        #print "  intersect:", intersect_t
+        #print "  u        :", u[0], u[1], u[2]
+        #
+        return
+    #print cur_ind[0], dims[0], cur_ind[1], dims[1], cur_ind[2], dims[2]
+    dv = data[cur_ind[0], cur_ind[1], cur_ind[2]]
+    #dt = 1e300
+    while 1:
+        if image[ind,3] >= 1.0: break
+        if (not (0 <= cur_ind[0] < dims[0])) or \
+           (not (0 <= cur_ind[1] < dims[1])) or \
+           (not (0 <= cur_ind[2] < dims[2])):
+            break
+        # Do our transfer here
+        for n in range(nshells):
+            dist = shells[n, 0] - dv
+            if dist < shells[n,1]:
+                dist = exp(-dist/8.0)
+                alpha = (1.0 - shells[n,5])*shells[n,5]#*dt
+                image[ind,0] += alpha*shells[n,2]*dist
+                image[ind,1] += alpha*shells[n,3]*dist
+                image[ind,2] += alpha*shells[n,4]*dist
+                image[ind,3] += alpha*shells[n,5]*dist
+                #image[ind,i] += rgba[i]*dist*rgba[3]/dt
+                #print rgba[i], image[ind,i], dist, dt
+                break
+        if (tmax[0] > 1.0) and (tmax[1] > 1.0) and (tmax[2] > 1.0):
+            dt = 1.0 - enter_t
+            break
+        if tmax[0] < tmax[1]:
+            if tmax[0] < tmax[2]:
+                dt = tmax[0] - enter_t
+                enter_t = tmax[0]
+                tmax[0] += tdelta[0]
+                cur_ind[0] += step[0]
+            else:
+                dt = tmax[2] - enter_t
+                enter_t = tmax[2]
+                tmax[2] += tdelta[2]
+                cur_ind[2] += step[2]
+        else:
+            if tmax[1] < tmax[2]:
+                dt = tmax[1] - enter_t
+                enter_t = tmax[1]
+                tmax[1] += tdelta[1]
+                cur_ind[1] += step[1]
+            else:
+                dt = tmax[2] - enter_t
+                enter_t = tmax[2]
+                tmax[2] += tdelta[2]
+                cur_ind[2] += step[2]
+        dv = data[cur_ind[0], cur_ind[1], cur_ind[2]]

Added: trunk/yt/_amr_utils/VolumeIntegrator.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/VolumeIntegrator.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,384 @@
+"""
+Simple integrators for the radiative transfer equation
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2009 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from stdlib cimport malloc, free, abs
+
+cdef inline int imax(int i0, int i1):
+    if i0 > i1: return i0
+    return i1
+
+cdef inline np.float64_t fmax(np.float64_t f0, np.float64_t f1):
+    if f0 > f1: return f0
+    return f1
+
+cdef inline int imin(int i0, int i1):
+    if i0 < i1: return i0
+    return i1
+
+cdef inline np.float64_t fmin(np.float64_t f0, np.float64_t f1):
+    if f0 < f1: return f0
+    return f1
+
+cdef inline int iclip(int i, int a, int b):
+    return imin(imax(i, a), b)
+
+cdef inline np.float64_t fclip(np.float64_t f,
+                      np.float64_t a, np.float64_t b):
+    return fmin(fmax(f, a), b)
+
+cdef extern from "math.h":
+    double exp(double x)
+    float expf(float x)
+    double floor(double x)
+    double ceil(double x)
+    double fmod(double x, double y)
+    double log2(double x)
+
+cdef extern from "FixedInterpolator.h":
+    np.float64_t fast_interpolate(int *ds, int *ci, np.float64_t *dp,
+                                  np.float64_t *data)
+cdef extern from "FixedInterpolator.h":
+    np.float64_t trilinear_interpolate(int *ds, int *ci, np.float64_t *dp,
+                                       np.float64_t *data)
+
+cdef class VectorPlane
+
+cdef class TransferFunctionProxy:
+    cdef np.float64_t x_bounds[2]
+    cdef np.float64_t *vs[4]
+    cdef int nbins
+    cdef np.float64_t dbin
+    cdef public object tf_obj
+    def __cinit__(self, tf_obj):
+        self.tf_obj = tf_obj
+        cdef np.ndarray[np.float64_t, ndim=1] temp
+        temp = tf_obj.red.y
+        self.vs[0] = <np.float64_t *> temp.data
+        temp = tf_obj.green.y
+        self.vs[1] = <np.float64_t *> temp.data
+        temp = tf_obj.blue.y
+        self.vs[2] = <np.float64_t *> temp.data
+        temp = tf_obj.alpha.y
+        self.vs[3] = <np.float64_t *> temp.data
+        self.x_bounds[0] = tf_obj.x_bounds[0]
+        self.x_bounds[1] = tf_obj.x_bounds[1]
+        self.nbins = tf_obj.nbins
+        self.dbin = (self.x_bounds[1] - self.x_bounds[0])/self.nbins
+
+    cdef void eval_transfer(self, np.float64_t dt, np.float64_t dv,
+                                    np.float64_t *rgba):
+        cdef int i
+        cdef int bin_id
+        cdef np.float64_t tf, trgba[4], bv, dx, dy, dd,ta
+        dx = self.dbin
+
+        # get source alpha first
+        # First locate our points
+        bin_id = iclip(<int> floor((dv - self.x_bounds[0]) / dx),
+                        0, self.nbins-2)
+            # Recall that linear interpolation is y0 + (x-x0) * dx/dy
+        bv = self.vs[3][bin_id] # This is x0
+        dy = self.vs[3][bin_id+1]-bv # dy
+        dd = dv-(self.x_bounds[0] + bin_id * dx) # x - x0
+            # This is our final value for transfer function on the entering face
+        tf = bv+dd*(dy/dx) 
+        ta = tf  # Store the source alpha
+        for i in range(3):
+            # Recall that linear interpolation is y0 + (x-x0) * dx/dy
+            bv = self.vs[i][bin_id] # This is x0
+            dy = self.vs[i][bin_id+1]-bv # dy
+            dd = dv-(self.x_bounds[0] + bin_id * dx) # x - x0
+            # This is our final value for transfer function on the entering face
+            tf = bv+dd*(dy/dx) 
+            # alpha blending
+            rgba[i] += (1. - rgba[3])*ta*tf*dt
+        #update alpha
+        rgba[3] += (1. - rgba[3])*ta*dt
+        
+        # We should really do some alpha blending.
+        # Front to back blending is defined as:
+        #  dst.rgb = dst.rgb + (1 - dst.a) * src.a * src.rgb
+        #  dst.a   = dst.a   + (1 - dst.a) * src.a     
+
+cdef class VectorPlane:
+    cdef public object avp_pos, avp_dir, acenter, aimage
+    cdef np.float64_t *vp_pos, *vp_dir, *center, *image,
+    cdef np.float64_t pdx, pdy, bounds[4]
+    cdef int nv
+    cdef public object ax_vec, ay_vec
+    cdef np.float64_t *x_vec, *y_vec
+
+    def __cinit__(self, 
+                  np.ndarray[np.float64_t, ndim=3] vp_pos,
+                  np.ndarray[np.float64_t, ndim=1] vp_dir,
+                  np.ndarray[np.float64_t, ndim=1] center,
+                  bounds,
+                  np.ndarray[np.float64_t, ndim=3] image,
+                  np.ndarray[np.float64_t, ndim=1] x_vec,
+                  np.ndarray[np.float64_t, ndim=1] y_vec):
+        cdef int i, j
+        self.avp_pos = vp_pos
+        self.avp_dir = vp_dir
+        self.acenter = center
+        self.aimage = image
+        self.ax_vec = x_vec
+        self.ay_vec = y_vec
+        self.vp_pos = <np.float64_t *> vp_pos.data
+        self.vp_dir = <np.float64_t *> vp_dir.data
+        self.center = <np.float64_t *> center.data
+        self.image = <np.float64_t *> image.data
+        self.x_vec = <np.float64_t *> x_vec.data
+        self.y_vec = <np.float64_t *> y_vec.data
+        self.nv = vp_pos.shape[0]
+        for i in range(4): self.bounds[i] = bounds[i]
+        self.pdx = (self.bounds[1] - self.bounds[0])/self.nv
+        self.pdy = (self.bounds[3] - self.bounds[2])/self.nv
+
+    cdef void get_start_stop(self, np.float64_t *ex, int *rv):
+        # Extrema need to be re-centered
+        cdef np.float64_t cx, cy
+        cx = cy = 0.0
+        for i in range(3):
+            cx += self.center[i] * self.x_vec[i]
+            cy += self.center[i] * self.y_vec[i]
+        rv[0] = <int> floor((ex[0] - cx - self.bounds[0])/self.pdx)
+        rv[1] = rv[0] + <int> ceil((ex[1] - ex[0])/self.pdx)
+        rv[2] = <int> floor((ex[2] - cy - self.bounds[2])/self.pdy)
+        rv[3] = rv[2] + <int> ceil((ex[3] - ex[2])/self.pdy)
+
+    cdef inline void copy_into(self, np.float64_t *fv, np.float64_t *tv,
+                        int i, int j, int nk):
+        # We know the first two dimensions of our from-vector, and our
+        # to-vector is flat and 'ni' long
+        cdef int k
+        for k in range(nk):
+            tv[k] = fv[(((k*self.nv)+j)*self.nv+i)]
+
+    cdef inline void copy_back(self, np.float64_t *fv, np.float64_t *tv,
+                        int i, int j, int nk):
+        cdef int k
+        for k in range(nk):
+            tv[(((k*self.nv)+j)*self.nv+i)] = fv[k]
+
+cdef class PartitionedGrid:
+    cdef public object my_data
+    cdef public object LeftEdge
+    cdef public object RightEdge
+    cdef np.float64_t *data
+    cdef np.float64_t left_edge[3]
+    cdef np.float64_t right_edge[3]
+    cdef np.float64_t dds[3]
+    cdef public np.float64_t min_dds
+    cdef int ns
+    cdef int dims[3]
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def __cinit__(self,
+                  np.ndarray[np.float64_t, ndim=3] data,
+                  np.ndarray[np.float64_t, ndim=1] left_edge,
+                  np.ndarray[np.float64_t, ndim=1] right_edge,
+                  np.ndarray[np.int64_t, ndim=1] dims):
+        # The data is likely brought in via a slice, so we copy it
+        cdef int i, j, k, size
+        self.LeftEdge = left_edge
+        self.RightEdge = right_edge
+        for i in range(3):
+            self.left_edge[i] = left_edge[i]
+            self.right_edge[i] = right_edge[i]
+            self.dims[i] = dims[i]
+            self.dds[i] = (self.right_edge[i] - self.left_edge[i])/dims[i]
+        self.my_data = data
+        self.data = <np.float64_t*> data.data
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def cast_plane(self, TransferFunctionProxy tf, VectorPlane vp):
+        # This routine will iterate over all of the vectors and cast each in
+        # turn.  Might benefit from a more sophisticated intersection check,
+        # like http://courses.csusm.edu/cs697exz/ray_box.htm
+        cdef int vi, vj, hit, i, ni, nj, nn
+        self.ns = 5 #* (1 + <int> log2(self.dds[0] / self.min_dds))
+        cdef int iter[4]
+        cdef np.float64_t v_pos[3], v_dir[3], rgba[4], extrema[4]
+        self.calculate_extent(vp, extrema)
+        vp.get_start_stop(extrema, iter)
+        for i in range(4): iter[i] = iclip(iter[i], 0, vp.nv)
+        hit = 0
+        for vj in range(iter[0], iter[1]):
+            for vi in range(iter[2], iter[3]):
+                vp.copy_into(vp.vp_pos, v_pos, vi, vj, 3)
+                vp.copy_into(vp.image, rgba, vi, vj, 4)
+                self.integrate_ray(v_pos, vp.vp_dir, rgba, tf)
+                vp.copy_back(rgba, vp.image, vi, vj, 4)
+        return hit
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef void calculate_extent(self, VectorPlane vp,
+                               np.float64_t extrema[4]):
+        # We do this for all eight corners
+        cdef np.float64_t *edges[2], temp
+        edges[0] = self.left_edge
+        edges[1] = self.right_edge
+        extrema[0] = extrema[2] = 1e300; extrema[1] = extrema[3] = -1e300
+        cdef int i, j, k
+        for i in range(2):
+            for j in range(2):
+                for k in range(2):
+                    # This should rotate it into the vector plane
+                    temp  = edges[i][0] * vp.x_vec[0]
+                    temp += edges[j][1] * vp.x_vec[1]
+                    temp += edges[k][2] * vp.x_vec[2]
+                    if temp < extrema[0]: extrema[0] = temp
+                    if temp > extrema[1]: extrema[1] = temp
+                    temp  = edges[i][0] * vp.y_vec[0]
+                    temp += edges[j][1] * vp.y_vec[1]
+                    temp += edges[k][2] * vp.y_vec[2]
+                    if temp < extrema[2]: extrema[2] = temp
+                    if temp > extrema[3]: extrema[3] = temp
+        #print extrema[0], extrema[1], extrema[2], extrema[3]
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef int integrate_ray(self, np.float64_t v_pos[3],
+                                 np.float64_t v_dir[3],
+                                 np.float64_t rgba[4],
+                                 TransferFunctionProxy tf):
+        cdef int cur_ind[3], step[3], x, y, i, n, flat_ind, hit, direction
+        cdef np.float64_t intersect_t = 1.0
+        cdef np.float64_t intersect[3], tmax[3], tdelta[3]
+        cdef np.float64_t enter_t, dist, alpha, dt
+        cdef np.float64_t tr, tl, temp_x, temp_y, dv
+        for i in range(3):
+            if (v_dir[i] < 0):
+                step[i] = -1
+            else:
+                step[i] = 1
+            x = (i+1) % 3
+            y = (i+2) % 3
+            tl = (self.left_edge[i] - v_pos[i])/v_dir[i]
+            tr = (self.right_edge[i] - v_pos[i])/v_dir[i]
+            temp_x = (v_pos[x] + tl*v_dir[x])
+            temp_y = (v_pos[y] + tl*v_dir[y])
+            if self.left_edge[x] <= temp_x and temp_x <= self.right_edge[x] and \
+               self.left_edge[y] <= temp_y and temp_y <= self.right_edge[y] and \
+               0.0 <= tl and tl < intersect_t:
+                direction = i
+                intersect_t = tl
+            temp_x = (v_pos[x] + tr*v_dir[x])
+            temp_y = (v_pos[y] + tr*v_dir[y])
+            if self.left_edge[x] <= temp_x and temp_x <= self.right_edge[x] and \
+               self.left_edge[y] <= temp_y and temp_y <= self.right_edge[y] and \
+               0.0 <= tr and tr < intersect_t:
+                direction = i
+                intersect_t = tr
+        if self.left_edge[0] <= v_pos[0] and v_pos[0] <= self.right_edge[0] and \
+           self.left_edge[1] <= v_pos[1] and v_pos[1] <= self.right_edge[1] and \
+           self.left_edge[2] <= v_pos[2] and v_pos[2] <= self.right_edge[2]:
+            intersect_t = 0.0
+        if not ((0.0 <= intersect_t) and (intersect_t < 1.0)): return 0
+        for i in range(3):
+            intersect[i] = v_pos[i] + intersect_t * v_dir[i]
+            cur_ind[i] = <int> floor((intersect[i] +
+                                      step[i]*1e-8*self.dds[i] -
+                                      self.left_edge[i])/self.dds[i])
+            tmax[i] = (((cur_ind[i]+step[i])*self.dds[i])+
+                        self.left_edge[i]-v_pos[i])/v_dir[i]
+            if cur_ind[i] == self.dims[i] and step[i] < 0:
+                cur_ind[i] = self.dims[i] - 1
+            if cur_ind[i] < 0 or cur_ind[i] >= self.dims[i]: return 0
+            if step[i] > 0:
+                tmax[i] = (((cur_ind[i]+1)*self.dds[i])
+                            +self.left_edge[i]-v_pos[i])/v_dir[i]
+            if step[i] < 0:
+                tmax[i] = (((cur_ind[i]+0)*self.dds[i])
+                            +self.left_edge[i]-v_pos[i])/v_dir[i]
+            tdelta[i] = (self.dds[i]/v_dir[i])
+            if tdelta[i] < 0: tdelta[i] *= -1
+        # We have to jumpstart our calculation
+        enter_t = intersect_t
+        while 1:
+            # dims here is one less than the dimensions of the data,
+            # but we are tracing on the grid, not on the data...
+            if (not (0 <= cur_ind[0] < self.dims[0])) or \
+               (not (0 <= cur_ind[1] < self.dims[1])) or \
+               (not (0 <= cur_ind[2] < self.dims[2])):
+                break
+            hit += 1
+            if tmax[0] < tmax[1]:
+                if tmax[0] < tmax[2]:
+                    self.sample_values(v_pos, v_dir, enter_t, tmax[0], cur_ind,
+                                       rgba, tf)
+                    cur_ind[0] += step[0]
+                    dt = fmin(tmax[0], 1.0) - enter_t
+                    enter_t = tmax[0]
+                    tmax[0] += tdelta[0]
+                else:
+                    self.sample_values(v_pos, v_dir, enter_t, tmax[2], cur_ind,
+                                       rgba, tf)
+                    cur_ind[2] += step[2]
+                    dt = fmin(tmax[2], 1.0) - enter_t
+                    enter_t = tmax[2]
+                    tmax[2] += tdelta[2]
+            else:
+                if tmax[1] < tmax[2]:
+                    self.sample_values(v_pos, v_dir, enter_t, tmax[1], cur_ind,
+                                       rgba, tf)
+                    cur_ind[1] += step[1]
+                    dt = fmin(tmax[1], 1.0) - enter_t
+                    enter_t = tmax[1]
+                    tmax[1] += tdelta[1]
+                else:
+                    self.sample_values(v_pos, v_dir, enter_t, tmax[2], cur_ind,
+                                       rgba, tf)
+                    cur_ind[2] += step[2]
+                    dt = fmin(tmax[2], 1.0) - enter_t
+                    enter_t = tmax[2]
+                    tmax[2] += tdelta[2]
+            if enter_t > 1.0: break
+        return hit
+
+    cdef void sample_values(self,
+                            np.float64_t v_pos[3],
+                            np.float64_t v_dir[3],
+                            np.float64_t enter_t,
+                            np.float64_t exit_t,
+                            int ci[3],
+                            np.float64_t *rgba,
+                            TransferFunctionProxy tf):
+        cdef np.float64_t cp[3], dp[3], temp, dt, t, dv
+        cdef int dti, i
+        dt = (exit_t - enter_t) / (self.ns-1) # five samples, so divide by four
+        for dti in range(self.ns - 1):
+            t = enter_t + dt * dti
+            for i in range(3):
+                cp[i] = v_pos[i] + t * v_dir[i]
+                dp[i] = fclip(fmod(cp[i], self.dds[i])/self.dds[i], 0, 1.0)
+            dv = trilinear_interpolate(self.dims, ci, dp, self.data)
+            tf.eval_transfer(dt, dv, rgba)

Added: trunk/yt/amr_utils.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/amr_utils.pyx	Fri Nov 20 19:54:16 2009
@@ -0,0 +1,43 @@
+"""
+Container file to hold all our Cython routines.  This is to avoid problems with
+static linking.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2008 Matthew Turk.  All Rights Reserved.
+
+  This file is part of yt.
+
+  yt is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+#cython embedsignature=True
+#cython cdivision=True
+#cython always_allow_keywords=True
+
+# Set up some imports
+import numpy as np
+cimport numpy as np
+cimport cython
+
+# We include all of our files
+
+include "_amr_utils/DepthFirstOctree.pyx"
+include "_amr_utils/Interpolators.pyx"
+include "_amr_utils/PointsInVolume.pyx"
+include "_amr_utils/RayIntegrators.pyx"
+include "_amr_utils/VolumeIntegrator.pyx"
+include "_amr_utils/CICDeposit.pyx"
\ No newline at end of file

Modified: trunk/yt/extensions/volume_rendering/__init__.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/__init__.py	(original)
+++ trunk/yt/extensions/volume_rendering/__init__.py	Fri Nov 20 19:54:16 2009
@@ -26,7 +26,7 @@
 import numpy as na
 
 from TransferFunction import TransferFunction, ColorTransferFunction
-from VolumeIntegrator import PartitionedGrid, VectorPlane, \
+from yt.amr_utils import PartitionedGrid, VectorPlane, \
                              TransferFunctionProxy
 from grid_partitioner import partition_all_grids, partition_grid, \
                              export_partitioned_grids, \

Modified: trunk/yt/extensions/volume_rendering/setup.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/setup.py	(original)
+++ trunk/yt/extensions/volume_rendering/setup.py	Fri Nov 20 19:54:16 2009
@@ -4,12 +4,11 @@
 
 import os.path
 
+os.system("cython -a yt/extensions/volume_rendering/VolumeIntegrator.pyx")
+
 def configuration(parent_package='',top_path=None):
     from numpy.distutils.misc_util import Configuration
     config = Configuration('volume_rendering',parent_package,top_path)
     config.make_config_py() # installs __config__.py
     config.make_svn_version_py()
-    config.add_extension("VolumeIntegrator",
-        ["yt/extensions/volume_rendering/VolumeIntegrator.c",
-         "yt/extensions/volume_rendering/FixedInterpolator.c",])
     return config

Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py	(original)
+++ trunk/yt/lagos/BaseDataTypes.py	Fri Nov 20 19:54:16 2009
@@ -28,7 +28,6 @@
 data_object_registry = {}
 
 from yt.lagos import *
-import PointsInVolume as PV
 
 def restore_grid_state(func):
     """
@@ -1168,7 +1167,7 @@
 
     def _get_point_indices(self, grid):
         if self._pixelmask.max() == 0: return []
-        k = PV.PointsInVolume(self._coord, self._pixelmask,
+        k = amr_utils.PointsInVolume(self._coord, self._pixelmask,
                               grid.LeftEdge, grid.RightEdge,
                               grid.child_mask, just_one(grid['dx']))
         return k

Modified: trunk/yt/lagos/HelperFunctions.py
==============================================================================
--- trunk/yt/lagos/HelperFunctions.py	(original)
+++ trunk/yt/lagos/HelperFunctions.py	Fri Nov 20 19:54:16 2009
@@ -25,7 +25,6 @@
 """
 
 from yt.lagos import *
-import Interpolators as IT
 
 class UnilinearFieldInterpolator:
     def __init__(self, table, boundaries, field_names):
@@ -46,7 +45,7 @@
             raise ValueError
 
         my_vals = na.zeros(x_vals.shape, dtype='float64')
-        IT.UnilinearlyInterpolate(self.table, x_vals, self.x_bins, x_i, my_vals)
+        amr_utils.UnilinearlyInterpolate(self.table, x_vals, self.x_bins, x_i, my_vals)
         return my_vals
 
 class BilinearFieldInterpolator:
@@ -77,7 +76,7 @@
                 y_i = na.minimum(na.maximum(y_i,0), len(self.y_bins)-2)
 
         my_vals = na.zeros(x_vals.shape, dtype='float64')
-        IT.BilinearlyInterpolate(self.table,
+        amr_utils.BilinearlyInterpolate(self.table,
                                  x_vals, y_vals, self.x_bins, self.y_bins,
                                  x_i, y_i, my_vals)
         return my_vals.reshape(orig_shape)
@@ -115,7 +114,7 @@
                 z_i = na.minimum(na.maximum(z_i,0), len(self.z_bins)-2)
 
         my_vals = na.zeros(x_vals.shape, dtype='float64')
-        IT.TrilinearlyInterpolate(self.table,
+        amr_utils.TrilinearlyInterpolate(self.table,
                                  x_vals, y_vals, z_vals,
                                  self.x_bins, self.y_bins, self.z_bins,
                                  x_i, y_i, z_i, my_vals)

Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py	(original)
+++ trunk/yt/lagos/__init__.py	Fri Nov 20 19:54:16 2009
@@ -55,6 +55,8 @@
 
 import time
 
+import yt.amr_utils as amr_utils
+
 from Cosmology import *
 from EnzoCosmology import *
 

Modified: trunk/yt/lagos/setup.py
==============================================================================
--- trunk/yt/lagos/setup.py	(original)
+++ trunk/yt/lagos/setup.py	Fri Nov 20 19:54:16 2009
@@ -20,12 +20,9 @@
     config.make_config_py() # installs __config__.py
     config.make_svn_version_py()
     config.add_extension("PointCombine", "yt/lagos/PointCombine.c", libraries=["m"])
-    config.add_extension("RTIntegrator", "yt/lagos/RTIntegrator.c", libraries=["m"])
-    config.add_extension("Interpolators", "yt/lagos/Interpolators.c")
-    config.add_extension("PointsInVolume", "yt/lagos/PointsInVolume.c")
-    #config.add_extension("DepthFirstOctree", "yt/lagos/DepthFirstOctree.c")
     config.add_subpackage("hop")
     config.add_subpackage("fof")
+    config.add_subpackage("parallelHOP")
     H5dir = check_for_hdf5()
     if H5dir is not None:
         include_dirs=[os.path.join(H5dir,"include")]
@@ -35,6 +32,6 @@
                              libraries=["m","hdf5"],
                              library_dirs=library_dirs, include_dirs=include_dirs)
     # Uncomment the next two lines if you want particle_density support
-    #config.add_extension("cic_deposit", ["yt/lagos/enzo_routines/cic_deposit.pyf",
-    #                                     "yt/lagos/enzo_routines/cic_deposit.f"])
+    config.add_extension("cic_deposit", ["yt/lagos/enzo_routines/cic_deposit.pyf",
+                                         "yt/lagos/enzo_routines/cic_deposit.f"])
     return config

Modified: trunk/yt/setup.py
==============================================================================
--- trunk/yt/setup.py	(original)
+++ trunk/yt/setup.py	Fri Nov 20 19:54:16 2009
@@ -9,6 +9,10 @@
     config.add_subpackage('fido')
     config.add_subpackage('reason')
     config.add_subpackage('extensions')
+    config.add_extension("amr_utils", 
+        ["yt/amr_utils.c", "yt/_amr_utils/FixedInterpolator.c"], 
+        include_dirs=["yt/_amr_utils/"],
+        libraries=["m"])
     config.make_config_py()
     config.make_svn_version_py()
     return config



More information about the yt-svn mailing list