[yt-svn] commit/yt: 7 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Feb 3 09:09:14 PST 2016
7 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/9b656920f242/
Changeset: 9b656920f242
Branch: yt
User: MatthewTurk
Date: 2015-12-28 17:56:50+00:00
Summary: Adding an indirection and implementing field pointers.
Affected #: 4 files
diff -r d9502d84802972ba3e644d2d5c47b87f73a6bd08 -r 9b656920f242a70db05ae9d89ca4bb5a1f87b723 yt/frontends/artio/_artio_caller.pyx
--- a/yt/frontends/artio/_artio_caller.pyx
+++ b/yt/frontends/artio/_artio_caller.pyx
@@ -17,6 +17,7 @@
from libc.stdlib cimport malloc, free
from libc.string cimport memcpy
import data_structures
+from yt.utilities.lib.misc_utilities import OnceIndirect
cdef extern from "platform_dep.h":
void *alloca(int)
@@ -1556,6 +1557,9 @@
if fields is None:
fields = []
nf = len(fields)
+ cdef np.float64_t[::cython.view.indirect, ::1] field_pointers
+ if nf > 0: field_pointers = OnceIndirect(fields)
+ cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
cdef np.ndarray[np.uint8_t, ndim=1, cast=True] mask
mask = self.mask(selector, -1)
cdef np.ndarray[np.int64_t, ndim=1] domain_ind
@@ -1570,17 +1574,9 @@
continue
domain_ind[sfc - self.sfc_start] = j
j += 1
- cdef np.float64_t **field_pointers
- cdef np.float64_t *field_vals
cdef np.float64_t pos[3]
cdef np.float64_t left_edge[3]
cdef int coords[3]
- cdef np.ndarray[np.float64_t, ndim=1] tarr
- field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf)
- field_vals = <np.float64_t*>alloca(sizeof(np.float64_t) * nf)
- for i in range(nf):
- tarr = fields[i]
- field_pointers[i] = <np.float64_t *> tarr.data
cdef int dims[3]
dims[0] = dims[1] = dims[2] = 1
cdef np.int64_t offset, moff
diff -r d9502d84802972ba3e644d2d5c47b87f73a6bd08 -r 9b656920f242a70db05ae9d89ca4bb5a1f87b723 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -138,5 +138,5 @@
cdef public int update_values
cdef void process(self, int dim[3], np.float64_t left_edge[3],
np.float64_t dds[3], np.int64_t offset,
- np.float64_t ppos[3], np.float64_t *fields,
+ np.float64_t ppos[3], np.float64_t[:] fields,
np.int64_t domain_ind)
diff -r d9502d84802972ba3e644d2d5c47b87f73a6bd08 -r 9b656920f242a70db05ae9d89ca4bb5a1f87b723 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -19,10 +19,21 @@
from libc.stdlib cimport malloc, free
cimport cython
from libc.math cimport sqrt
-
+from cpython cimport PyObject
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, \
OctreeContainer, OctInfo
+from cpython.array cimport array, clone
+from cython.view cimport memoryview as cymemview
+from yt.utilities.lib.misc_utilities import OnceIndirect
+
+cdef prepend_axes(np.ndarray arr, int naxes):
+ if arr.ndim == naxes:
+ return arr
+ # Avoid copies
+ arr2 = arr.view()
+ arr2.shape = (1,) * (naxes - arr2.ndim) + arr2.shape
+ return arr
cdef class ParticleDepositOperation:
def __init__(self, nvals, kernel_name):
@@ -47,15 +58,10 @@
if fields is None:
fields = []
nf = len(fields)
- cdef np.float64_t **field_pointers
- cdef np.float64_t *field_vals
+ cdef np.float64_t[::cython.view.indirect, ::1] field_pointers
+ if nf > 0: field_pointers = OnceIndirect(fields)
cdef np.float64_t pos[3]
- cdef np.ndarray[np.float64_t, ndim=1] tarr
- field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf)
- field_vals = <np.float64_t*>alloca(sizeof(np.float64_t) * nf)
- for i in range(nf):
- tarr = fields[i]
- field_pointers[i] = <np.float64_t *> tarr.data
+ cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
cdef int dims[3]
dims[0] = dims[1] = dims[2] = (1 << octree.oref)
cdef int nz = dims[0] * dims[1] * dims[2]
@@ -67,7 +73,7 @@
for i in range(positions.shape[0]):
# We should check if particle remains inside the Oct here
for j in range(nf):
- field_vals[j] = field_pointers[j][i]
+ field_vals[j] = field_pointers[j,i]
for j in range(3):
pos[j] = positions[i, j]
# This line should be modified to have it return the index into an
@@ -108,16 +114,11 @@
if fields is None:
fields = []
nf = len(fields)
- cdef np.float64_t **field_pointers
- cdef np.float64_t *field_vals
+ cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
+ cdef np.float64_t[::cython.view.indirect, ::1] field_pointers
+ if nf > 0: field_pointers = OnceIndirect(fields)
cdef np.float64_t pos[3]
- cdef np.ndarray[np.float64_t, ndim=1] tarr
- field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf)
- field_vals = <np.float64_t*>alloca(sizeof(np.float64_t) * nf)
cdef np.int64_t gid = getattr(gobj, "id", -1)
- for i in range(nf):
- tarr = fields[i]
- field_pointers[i] = <np.float64_t *> tarr.data
cdef np.float64_t dds[3]
cdef np.float64_t left_edge[3]
cdef int dims[3]
@@ -128,7 +129,7 @@
for i in range(positions.shape[0]):
# Now we process
for j in range(nf):
- field_vals[j] = field_pointers[j][i]
+ field_vals[j] = field_pointers[j,i]
for j in range(3):
pos[j] = positions[i, j]
self.process(dims, left_edge, dds, 0, pos, field_vals, gid)
@@ -138,19 +139,16 @@
cdef void process(self, int dim[3], np.float64_t left_edge[3],
np.float64_t dds[3], np.int64_t offset,
- np.float64_t ppos[3], np.float64_t *fields,
+ np.float64_t ppos[3], np.float64_t[:] fields,
np.int64_t domain_ind):
raise NotImplementedError
cdef class CountParticles(ParticleDepositOperation):
- cdef np.int64_t *count # float, for ease
- cdef public object ocount
+ cdef np.int64_t[:,:,:,:] count
def initialize(self):
# Create a numpy array accessible to python
- self.ocount = np.zeros(self.nvals, dtype="int64", order='F')
- cdef np.ndarray arr = self.ocount
- # alias the C-view for use in cython
- self.count = <np.int64_t*> arr.data
+ self.count = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -158,7 +156,7 @@
np.float64_t dds[3],
np.int64_t offset, # offset into IO field
np.float64_t ppos[3], # this particle's position
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
# here we do our thing; this is the kernel
@@ -166,10 +164,12 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- self.count[gind(ii[0], ii[1], ii[2], dim) + offset] += 1
+ self.count[offset, ii[0], ii[1], ii[2]] += 1
def finalize(self):
- return self.ocount.astype('f8')
+ arr = np.asarray(self.count)
+ arr.shape = self.nvals
+ return arr.astype("float64")
deposit_count = CountParticles
@@ -196,7 +196,7 @@
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
@@ -256,7 +256,7 @@
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
@@ -304,7 +304,7 @@
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
@@ -348,7 +348,7 @@
np.float64_t dds[3],
np.int64_t offset, # offset into IO field
np.float64_t ppos[3], # this particle's position
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
@@ -403,7 +403,7 @@
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
@@ -431,7 +431,7 @@
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
fields[0] = domain_ind
@@ -462,7 +462,7 @@
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
# This one is a bit slow. Every grid cell is going to be iterated
diff -r d9502d84802972ba3e644d2d5c47b87f73a6bd08 -r 9b656920f242a70db05ae9d89ca4bb5a1f87b723 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -22,6 +22,13 @@
from fp_utils cimport fmin, fmax, i64min, i64max
from yt.geometry.selection_routines cimport _ensure_code
+from libc.stdlib cimport malloc, free
+from libc.string cimport strcmp
+
+from cython.view cimport memoryview
+from cpython cimport buffer
+
+
cdef extern from "stdlib.h":
# NOTE that size_t might not be int
void *alloca(int)
@@ -957,3 +964,128 @@
dest[i,j,k] += dsp * (overlap[0]*overlap[1]*overlap[2])
else:
dest[i,j,k] = dsp
+
+# The OnceIndirect code is from:
+# http://stackoverflow.com/questions/10465091/assembling-a-cython-memoryview-from-numpy-arrays/12991519#12991519
+# This is under the CC-BY-SA license.
+
+cdef class OnceIndirect:
+ cdef object _objects
+ cdef void** buf
+ cdef int ndim
+ cdef int n_rows
+ cdef int buf_len
+ cdef Py_ssize_t* shape
+ cdef Py_ssize_t* strides
+ cdef Py_ssize_t* suboffsets
+ cdef Py_ssize_t itemsize
+ cdef bytes format
+ cdef int is_readonly
+
+ def __cinit__(self, object rows, want_writable=True, want_format=True, allow_indirect=False):
+ """
+ Set want_writable to False if you don't want writable data. (This may
+ prevent copies.)
+ Set want_format to False if your input doesn't support PyBUF_FORMAT (unlikely)
+ Set allow_indirect to True if you are ok with the memoryview being indirect
+ in dimensions other than the first. (This may prevent copies.)
+
+ An example usage:
+
+ cdef double[::cython.view.indirect, ::1] vectors =
+ OnceIndirect([object.vector for object in objects])
+ """
+ demand = buffer.PyBUF_INDIRECT if allow_indirect else buffer.PyBUF_STRIDES
+ if want_writable:
+ demand |= buffer.PyBUF_WRITABLE
+ if want_format:
+ demand |= buffer.PyBUF_FORMAT
+ self._objects = [memoryview(row, demand) for row in rows]
+ self.n_rows = len(self._objects)
+ self.buf_len = sizeof(void*) * self.n_rows
+ self.buf = <void**>malloc(self.buf_len)
+ self.ndim = 1 + self._objects[0].ndim
+ self.shape = <Py_ssize_t*>malloc(sizeof(Py_ssize_t) * self.ndim)
+ self.strides = <Py_ssize_t*>malloc(sizeof(Py_ssize_t) * self.ndim)
+ self.suboffsets = <Py_ssize_t*>malloc(sizeof(Py_ssize_t) * self.ndim)
+
+ cdef memoryview example_obj = self._objects[0]
+ self.itemsize = example_obj.itemsize
+
+ if want_format:
+ self.format = example_obj.view.format
+ else:
+ self.format = None
+ self.is_readonly |= example_obj.view.readonly
+
+ for dim in range(self.ndim):
+ if dim == 0:
+ self.shape[dim] = self.n_rows
+ self.strides[dim] = sizeof(void*)
+ self.suboffsets[dim] = 0
+ else:
+ self.shape[dim] = example_obj.view.shape[dim - 1]
+ self.strides[dim] = example_obj.view.strides[dim - 1]
+ if example_obj.view.suboffsets == NULL:
+ self.suboffsets[dim] = -1
+ else:
+ self.suboffsets[dim] = example_obj.suboffsets[dim - 1]
+
+ cdef memoryview obj
+ cdef int i = 0
+ for obj in self._objects:
+ assert_similar(example_obj, obj)
+ self.buf[i] = obj.view.buf
+ i += 1
+
+ def __getbuffer__(self, Py_buffer* buff, int flags):
+ if (flags & buffer.PyBUF_INDIRECT) != buffer.PyBUF_INDIRECT:
+ raise Exception("don't want to copy data")
+ if flags & buffer.PyBUF_WRITABLE and self.is_readonly:
+ raise Exception("couldn't provide writable, you should have demanded it earlier")
+ if flags & buffer.PyBUF_FORMAT:
+ if self.format is None:
+ raise Exception("couldn't provide format, you should have demanded it earlier")
+ buff.format = self.format
+ else:
+ buff.format = NULL
+
+ buff.buf = <void*>self.buf
+ buff.obj = self
+ buff.len = self.buf_len
+ buff.readonly = self.is_readonly
+ buff.ndim = self.ndim
+ buff.shape = self.shape
+ buff.strides = self.strides
+ buff.suboffsets = self.suboffsets
+ buff.itemsize = self.itemsize
+ buff.internal = NULL
+
+ def __dealloc__(self):
+ free(self.buf)
+ free(self.shape)
+ free(self.strides)
+ free(self.suboffsets)
+
+cdef int assert_similar(memoryview left_, memoryview right_) except -1:
+ cdef Py_buffer left = left_.view
+ cdef Py_buffer right = right_.view
+ assert left.ndim == right.ndim
+ cdef int i
+ for i in range(left.ndim):
+ assert left.shape[i] == right.shape[i], (left_.shape, right_.shape)
+ assert left.strides[i] == right.strides[i], (left_.strides, right_.strides)
+
+ if left.suboffsets == NULL:
+ assert right.suboffsets == NULL, (left_.suboffsets, right_.suboffsets)
+ else:
+ for i in range(left.ndim):
+ assert left.suboffsets[i] == right.suboffsets[i], (left_.suboffsets, right_.suboffsets)
+
+ if left.format == NULL:
+ assert right.format == NULL, (bytes(left.format), bytes(right.format))
+ else:
+ #alternatively, compare as Python strings:
+ #assert bytes(left.format) == bytes(right.format)
+ assert strcmp(left.format, right.format) == 0, (bytes(left.format), bytes(right.format))
+ return 0
https://bitbucket.org/yt_analysis/yt/commits/7062322f053f/
Changeset: 7062322f053f
Branch: yt
User: MatthewTurk
Date: 2015-12-28 21:06:44+00:00
Summary: Convert remaining deposit operations.
Affected #: 1 file
diff -r 9b656920f242a70db05ae9d89ca4bb5a1f87b723 -r 7062322f053f259a3f87579ad08a43070425ec5b yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -177,18 +177,14 @@
# Note that this does nothing at the edges. So it will give a poor
# estimate there, and since Octrees are mostly edges, this will be a very
# poor SPH kernel.
- cdef np.float64_t *data
- cdef public object odata
- cdef np.float64_t *temp
- cdef public object otemp
+ cdef np.float64_t[:,:,:,:] data
+ cdef np.float64_t[:,:,:,:] temp
def initialize(self):
- self.odata = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.odata
- self.data = <np.float64_t*> arr.data
- self.otemp = np.zeros(self.nvals, dtype="float64", order='F')
- arr = self.otemp
- self.temp = <np.float64_t*> arr.data
+ self.data = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ self.temp = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -228,14 +224,14 @@
dist = idist[0] + idist[1] + idist[2]
# Calculate distance in multiples of the smoothing length
dist = sqrt(dist) / fields[0]
- self.temp[gind(i,j,k,dim) + offset] = self.sph_kernel(dist)
- kernel_sum += self.temp[gind(i,j,k,dim) + offset]
+ self.temp[offset,i,j,k] = self.sph_kernel(dist)
+ kernel_sum += self.temp[offset,i,j,k]
# Having found the kernel, deposit accordingly into gdata
for i from ib0[0] <= i <= ib1[0]:
for j from ib0[1] <= j <= ib1[1]:
for k from ib0[2] <= k <= ib1[2]:
- dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
- self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
+ dist = self.temp[offset,i,j,k] / kernel_sum
+ self.data[offset,i,j,k] += fields[1] * dist
def finalize(self):
return self.odata
@@ -243,12 +239,10 @@
deposit_simple_smooth = SimpleSmooth
cdef class SumParticleField(ParticleDepositOperation):
- cdef np.float64_t *sum
- cdef public object osum
+ cdef np.float64_t[:,:,:,:] sum
def initialize(self):
- self.osum = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.osum
- self.sum = <np.float64_t*> arr.data
+ self.temp = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -263,11 +257,13 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.sum[gind(ii[0], ii[1], ii[2], dim) + offset] += fields[0]
+ self.sum[offset, ii[0], ii[1], ii[2]] += fields[0]
return
def finalize(self):
- return self.osum
+ sum = self.sum
+ sum.shape = self.nvals
+ return sum
deposit_sum = SumParticleField
@@ -275,28 +271,20 @@
# Thanks to Britton and MJ Turk for the link
# to a single-pass STD
# http://www.cs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
- cdef np.float64_t *mk
- cdef np.float64_t *qk
- cdef np.float64_t *i
- cdef public object omk
- cdef public object oqk
- cdef public object oi
+ cdef np.float64_t[:,:,:,:] mk
+ cdef np.float64_t[:,:,:,:] qk
+ cdef np.float64_t[:,:,:,:] i
def initialize(self):
# we do this in a single pass, but need two scalar
# per cell, M_k, and Q_k and also the number of particles
# deposited into each one
# the M_k term
- self.omk= np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray omkarr= self.omk
- self.mk= <np.float64_t*> omkarr.data
- # the Q_k term
- self.oqk= np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray oqkarr= self.oqk
- self.qk= <np.float64_t*> oqkarr.data
- # particle count
- self.oi = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray oiarr = self.oi
- self.i = <np.float64_t*> oiarr.data
+ self.mk = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ self.qk = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ self.i = prepend_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -312,35 +300,36 @@
cdef float k, mk, qk
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- cell_index = gind(ii[0], ii[1], ii[2], dim) + offset
- k = self.i[cell_index]
- mk = self.mk[cell_index]
- qk = self.qk[cell_index]
+ k = self.i[offset, ii[0], ii[1], ii[2]]
+ mk = self.mk[offset, ii[0], ii[1], ii[2]]
+ qk = self.qk[offset, ii[0], ii[1], ii[2]]
#print k, mk, qk, cell_index
if k == 0.0:
# Initialize cell values
- self.mk[cell_index] = fields[0]
+ self.mk[offset, ii[0], ii[1], ii[2]] = fields[0]
else:
- self.mk[cell_index] = mk + (fields[0] - mk) / k
- self.qk[cell_index] = qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
- self.i[cell_index] += 1
+ self.mk[offset, ii[0], ii[1], ii[2]] = mk + (fields[0] - mk) / k
+ self.qk[offset, ii[0], ii[1], ii[2]] = \
+ qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
+ self.i[offset, ii[0], ii[1], ii[2]] += 1
def finalize(self):
# This is the standard variance
# if we want sample variance divide by (self.oi - 1.0)
- std2 = self.oqk / self.oi
- std2[self.oi == 0.0] = 0.0
+ i = np.asarray(self.i)
+ std2 = np.asarray(self.qk) / i
+ std2[i == 0.0] = 0.0
+ std2.shape = self.nvals
return np.sqrt(std2)
deposit_std = StdParticleField
cdef class CICDeposit(ParticleDepositOperation):
- cdef np.float64_t *field
+ cdef np.float64_t[:,:,:,:] field
cdef public object ofield
def initialize(self):
- self.ofield = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.ofield
- self.field = <np.float64_t *> arr.data
+ self.field = prepend_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -373,29 +362,26 @@
for i in range(2):
for j in range(2):
for k in range(2):
- ii = gind(ind[0] - i, ind[1] - j, ind[2] - k, dim) + offset
- self.field[ii] += fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
+ self.field[offset, ind[0] - i, ind[1] - j, ind[2] - k] += \
+ fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
def finalize(self):
- return self.ofield
+ rv = np.asarray(self.field)
+ rv.shape = self.nvals
+ return rv
deposit_cic = CICDeposit
cdef class WeightedMeanParticleField(ParticleDepositOperation):
# Deposit both mass * field and mass into two scalars
# then in finalize divide mass * field / mass
- cdef np.float64_t *wf
- cdef public object owf
- cdef np.float64_t *w
- cdef public object ow
+ cdef np.float64_t[:,:,:,:] wf
+ cdef np.float64_t[:,:,:,:] w
def initialize(self):
- self.owf = np.zeros(self.nvals, dtype='float64', order='F')
- cdef np.ndarray wfarr = self.owf
- self.wf = <np.float64_t*> wfarr.data
-
- self.ow = np.zeros(self.nvals, dtype='float64', order='F')
- cdef np.ndarray warr = self.ow
- self.w = <np.float64_t*> warr.data
+ self.wf = prepend_axes(
+ np.zeros(self.nvals, dtype='float64', order='F'), 4)
+ self.w = prepend_axes(
+ np.zeros(self.nvals, dtype='float64', order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -410,11 +396,15 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.w[ gind(ii[0], ii[1], ii[2], dim) + offset] += fields[1]
- self.wf[gind(ii[0], ii[1], ii[2], dim) + offset] += fields[0] * fields[1]
+ self.w[offset, ii[0], ii[1], ii[2]] += fields[1]
+ self.wf[offset, ii[0], ii[1], ii[2]] += fields[0] * fields[1]
def finalize(self):
- return self.owf / self.ow
+ wf = np.asarray(self.owf)
+ w = np.asarray(self.ow)
+ rv = wf / w
+ rv.shape = self.nvals
+ return rv
deposit_weighted_mean = WeightedMeanParticleField
@@ -442,19 +432,14 @@
deposit_mesh_id = MeshIdentifier
cdef class NNParticleField(ParticleDepositOperation):
- cdef np.float64_t *nnfield
- cdef np.float64_t *distfield
- cdef public object onnfield
- cdef public object odistfield
+ cdef np.float64_t[:,:,:,:] nnfield
+ cdef np.float64_t[:,:,:,:] distfield
def initialize(self):
- self.onnfield = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.onnfield
- self.nnfield = <np.float64_t*> arr.data
-
- self.odistfield = np.zeros(self.nvals, dtype="float64", order='F')
- self.odistfield[:] = np.inf
- arr = self.odistfield
- self.distfield = <np.float64_t*> arr.data
+ self.nnfield = prepend_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.distfield = prepend_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.distfield[:] = np.inf
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -469,7 +454,6 @@
# over, and we're going to deposit particles in it.
cdef int i, j, k
cdef int ii[3]
- cdef np.int64_t ggind
cdef np.float64_t r2
cdef np.float64_t gpos[3]
gpos[0] = left_edge[0] + 0.5 * dds[0]
@@ -478,19 +462,20 @@
for j in range(dim[1]):
gpos[2] = left_edge[2] + 0.5 * dds[2]
for k in range(dim[2]):
- ggind = gind(i, j, k, dim) + offset
r2 = ((ppos[0] - gpos[0])*(ppos[0] - gpos[0]) +
(ppos[1] - gpos[1])*(ppos[1] - gpos[1]) +
(ppos[2] - gpos[2])*(ppos[2] - gpos[2]))
- if r2 < self.distfield[ggind]:
- self.distfield[ggind] = r2
- self.nnfield[ggind] = fields[0]
+ if r2 < self.distfield[offset,i,j,k]:
+ self.distfield[offset,i,j,k] = r2
+ self.nnfield[offset,i,j,k] = fields[0]
gpos[2] += dds[2]
gpos[1] += dds[1]
gpos[0] += dds[0]
return
def finalize(self):
- return self.onnfield
+ nn = np.asarray(self.nnfield)
+ nn.shape = self.nvals
+ return nn
deposit_nearest = NNParticleField
https://bitbucket.org/yt_analysis/yt/commits/2738eb221cb9/
Changeset: 2738eb221cb9
Branch: yt
User: MatthewTurk
Date: 2015-12-28 21:56:16+00:00
Summary: Reset usage of offset, fix array types.
Affected #: 3 files
diff -r 7062322f053f259a3f87579ad08a43070425ec5b -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -702,11 +702,11 @@
if cls is None:
raise YTParticleDepositionNotImplemented(method)
# We allocate number of zones, not number of octs
- op = cls(self.ActiveDimensions.prod(), kernel_name)
+ op = cls(self.ActiveDimensions, kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
- return vals.reshape(self.ActiveDimensions, order="C")
+ return vals.copy(order="C")
def write_to_gdf(self, gdf_path, fields, nprocs=1, field_units=None,
**kwargs):
diff -r 7062322f053f259a3f87579ad08a43070425ec5b -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -327,12 +327,12 @@
if cls is None:
raise YTParticleDepositionNotImplemented(method)
# We allocate number of zones, not number of octs
- op = cls(self.ActiveDimensions.prod(), kernel_name)
+ op = cls(tuple(self.ActiveDimensions), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
if vals is None: return
- return vals.reshape(self.ActiveDimensions, order="C")
+ return vals.copy(order="C")
def select_blocks(self, selector):
mask = self._get_selector_mask(selector)
diff -r 7062322f053f259a3f87579ad08a43070425ec5b -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -27,13 +27,13 @@
from cython.view cimport memoryview as cymemview
from yt.utilities.lib.misc_utilities import OnceIndirect
-cdef prepend_axes(np.ndarray arr, int naxes):
+cdef append_axes(np.ndarray arr, int naxes):
if arr.ndim == naxes:
return arr
# Avoid copies
arr2 = arr.view()
- arr2.shape = (1,) * (naxes - arr2.ndim) + arr2.shape
- return arr
+ arr2.shape = arr2.shape + (1,) * (naxes - arr2.ndim)
+ return arr2
cdef class ParticleDepositOperation:
def __init__(self, nvals, kernel_name):
@@ -96,7 +96,7 @@
if oct == NULL or (domain_id > 0 and oct.domain != domain_id):
continue
# Note that this has to be our local index, not our in-file index.
- offset = dom_ind[oct.domain_ind - moff] * nz
+ offset = dom_ind[oct.domain_ind - moff]
if offset < 0: continue
# Check that we found the oct ...
self.process(dims, oi.left_edge, oi.dds,
@@ -147,7 +147,7 @@
cdef np.int64_t[:,:,:,:] count
def initialize(self):
# Create a numpy array accessible to python
- self.count = prepend_axes(
+ self.count = append_axes(
np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
@@ -164,7 +164,7 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- self.count[offset, ii[0], ii[1], ii[2]] += 1
+ self.count[ii[0], ii[1], ii[2], offset] += 1
def finalize(self):
arr = np.asarray(self.count)
@@ -181,10 +181,10 @@
cdef np.float64_t[:,:,:,:] temp
def initialize(self):
- self.data = prepend_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
- self.temp = prepend_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ self.data = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.temp = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -224,14 +224,14 @@
dist = idist[0] + idist[1] + idist[2]
# Calculate distance in multiples of the smoothing length
dist = sqrt(dist) / fields[0]
- self.temp[offset,i,j,k] = self.sph_kernel(dist)
- kernel_sum += self.temp[offset,i,j,k]
+ self.temp[i,j,k,offset] = self.sph_kernel(dist)
+ kernel_sum += self.temp[i,j,k,offset]
# Having found the kernel, deposit accordingly into gdata
for i from ib0[0] <= i <= ib1[0]:
for j from ib0[1] <= j <= ib1[1]:
for k from ib0[2] <= k <= ib1[2]:
- dist = self.temp[offset,i,j,k] / kernel_sum
- self.data[offset,i,j,k] += fields[1] * dist
+ dist = self.temp[i,j,k,offset] / kernel_sum
+ self.data[i,j,k,offset] += fields[1] * dist
def finalize(self):
return self.odata
@@ -241,8 +241,8 @@
cdef class SumParticleField(ParticleDepositOperation):
cdef np.float64_t[:,:,:,:] sum
def initialize(self):
- self.temp = prepend_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ self.sum = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -257,11 +257,11 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.sum[offset, ii[0], ii[1], ii[2]] += fields[0]
+ self.sum[ii[0], ii[1], ii[2], offset] += fields[0]
return
def finalize(self):
- sum = self.sum
+ sum = np.asarray(self.sum)
sum.shape = self.nvals
return sum
@@ -279,12 +279,12 @@
# per cell, M_k, and Q_k and also the number of particles
# deposited into each one
# the M_k term
- self.mk = prepend_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
- self.qk = prepend_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
- self.i = prepend_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ self.mk = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.qk = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.i = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
cdef void process(self, int dim[3],
@@ -300,18 +300,18 @@
cdef float k, mk, qk
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- k = self.i[offset, ii[0], ii[1], ii[2]]
- mk = self.mk[offset, ii[0], ii[1], ii[2]]
- qk = self.qk[offset, ii[0], ii[1], ii[2]]
+ k = self.i[ii[0], ii[1], ii[2], offset]
+ mk = self.mk[ii[0], ii[1], ii[2], offset]
+ qk = self.qk[ii[0], ii[1], ii[2], offset]
#print k, mk, qk, cell_index
if k == 0.0:
# Initialize cell values
- self.mk[offset, ii[0], ii[1], ii[2]] = fields[0]
+ self.mk[ii[0], ii[1], ii[2], offset] = fields[0]
else:
- self.mk[offset, ii[0], ii[1], ii[2]] = mk + (fields[0] - mk) / k
- self.qk[offset, ii[0], ii[1], ii[2]] = \
+ self.mk[ii[0], ii[1], ii[2], offset] = mk + (fields[0] - mk) / k
+ self.qk[ii[0], ii[1], ii[2], offset] = \
qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
- self.i[offset, ii[0], ii[1], ii[2]] += 1
+ self.i[ii[0], ii[1], ii[2], offset] += 1
def finalize(self):
# This is the standard variance
@@ -328,7 +328,7 @@
cdef np.float64_t[:,:,:,:] field
cdef public object ofield
def initialize(self):
- self.field = prepend_axes(
+ self.field = append_axes(
np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
@@ -362,7 +362,7 @@
for i in range(2):
for j in range(2):
for k in range(2):
- self.field[offset, ind[0] - i, ind[1] - j, ind[2] - k] += \
+ self.field[ind[0] - i, ind[1] - j, ind[2] - k, offset] += \
fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
def finalize(self):
@@ -378,9 +378,9 @@
cdef np.float64_t[:,:,:,:] wf
cdef np.float64_t[:,:,:,:] w
def initialize(self):
- self.wf = prepend_axes(
+ self.wf = append_axes(
np.zeros(self.nvals, dtype='float64', order='F'), 4)
- self.w = prepend_axes(
+ self.w = append_axes(
np.zeros(self.nvals, dtype='float64', order='F'), 4)
@cython.cdivision(True)
@@ -396,8 +396,8 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.w[offset, ii[0], ii[1], ii[2]] += fields[1]
- self.wf[offset, ii[0], ii[1], ii[2]] += fields[0] * fields[1]
+ self.w[ii[0], ii[1], ii[2], offset] += fields[1]
+ self.wf[ii[0], ii[1], ii[2], offset] += fields[0] * fields[1]
def finalize(self):
wf = np.asarray(self.owf)
@@ -435,9 +435,9 @@
cdef np.float64_t[:,:,:,:] nnfield
cdef np.float64_t[:,:,:,:] distfield
def initialize(self):
- self.nnfield = prepend_axes(
+ self.nnfield = append_axes(
np.zeros(self.nvals, dtype="float64", order='F'), 4)
- self.distfield = prepend_axes(
+ self.distfield = append_axes(
np.zeros(self.nvals, dtype="float64", order='F'), 4)
self.distfield[:] = np.inf
@@ -465,9 +465,9 @@
r2 = ((ppos[0] - gpos[0])*(ppos[0] - gpos[0]) +
(ppos[1] - gpos[1])*(ppos[1] - gpos[1]) +
(ppos[2] - gpos[2])*(ppos[2] - gpos[2]))
- if r2 < self.distfield[offset,i,j,k]:
- self.distfield[offset,i,j,k] = r2
- self.nnfield[offset,i,j,k] = fields[0]
+ if r2 < self.distfield[i,j,k,offset]:
+ self.distfield[i,j,k,offset] = r2
+ self.nnfield[i,j,k,offset] = fields[0]
gpos[2] += dds[2]
gpos[1] += dds[1]
gpos[0] += dds[0]
https://bitbucket.org/yt_analysis/yt/commits/1b0f177a9c83/
Changeset: 1b0f177a9c83
Branch: yt
User: MatthewTurk
Date: 2016-02-02 19:44:27+00:00
Summary: Merging from upstream
Affected #: 209 files
diff -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 doc/extensions/notebook_sphinxext.py
--- a/doc/extensions/notebook_sphinxext.py
+++ /dev/null
@@ -1,241 +0,0 @@
-import errno
-import os
-import shutil
-import string
-import re
-import tempfile
-import uuid
-from sphinx.util.compat import Directive
-from docutils import nodes
-from docutils.parsers.rst import directives
-from IPython.config import Config
-from IPython.nbconvert import html, python
-from IPython.nbformat import current as nbformat
-from runipy.notebook_runner import NotebookRunner, NotebookError
-
-class NotebookDirective(Directive):
- """Insert an evaluated notebook into a document
-
- This uses runipy and nbconvert to transform a path to an unevaluated notebook
- into html suitable for embedding in a Sphinx document.
- """
- required_arguments = 1
- optional_arguments = 1
- option_spec = {'skip_exceptions': directives.flag}
- final_argument_whitespace = True
-
- def run(self): # check if there are spaces in the notebook name
- nb_path = self.arguments[0]
- if ' ' in nb_path: raise ValueError(
- "Due to issues with docutils stripping spaces from links, white "
- "space is not allowed in notebook filenames '{0}'".format(nb_path))
- # check if raw html is supported
- if not self.state.document.settings.raw_enabled:
- raise self.warning('"%s" directive disabled.' % self.name)
-
- cwd = os.getcwd()
- tmpdir = tempfile.mkdtemp()
- os.chdir(tmpdir)
-
- # get path to notebook
- nb_filename = self.arguments[0]
- nb_basename = os.path.basename(nb_filename)
- rst_file = self.state_machine.document.attributes['source']
- rst_dir = os.path.abspath(os.path.dirname(rst_file))
- nb_abs_path = os.path.abspath(os.path.join(rst_dir, nb_filename))
-
- # Move files around.
- rel_dir = os.path.relpath(rst_dir, setup.confdir)
- dest_dir = os.path.join(setup.app.builder.outdir, rel_dir)
- dest_path = os.path.join(dest_dir, nb_basename)
-
- image_dir, image_rel_dir = make_image_dir(setup, rst_dir)
-
- # Ensure desination build directory exists
- thread_safe_mkdir(os.path.dirname(dest_path))
-
- # Copy unevaluated notebook
- shutil.copyfile(nb_abs_path, dest_path)
-
- # Construct paths to versions getting copied over
- dest_path_eval = string.replace(dest_path, '.ipynb', '_evaluated.ipynb')
- dest_path_script = string.replace(dest_path, '.ipynb', '.py')
- rel_path_eval = string.replace(nb_basename, '.ipynb', '_evaluated.ipynb')
- rel_path_script = string.replace(nb_basename, '.ipynb', '.py')
-
- # Create python script vesion
- script_text = nb_to_python(nb_abs_path)
- f = open(dest_path_script, 'w')
- f.write(script_text.encode('utf8'))
- f.close()
-
- skip_exceptions = 'skip_exceptions' in self.options
-
- ret = evaluate_notebook(
- nb_abs_path, dest_path_eval, skip_exceptions=skip_exceptions)
-
- try:
- evaluated_text, resources = ret
- evaluated_text = write_notebook_output(
- resources, image_dir, image_rel_dir, evaluated_text)
- except ValueError:
- # This happens when a notebook raises an unhandled exception
- evaluated_text = ret
-
- # Create link to notebook and script files
- link_rst = "(" + \
- formatted_link(nb_basename) + "; " + \
- formatted_link(rel_path_eval) + "; " + \
- formatted_link(rel_path_script) + \
- ")"
-
- self.state_machine.insert_input([link_rst], rst_file)
-
- # create notebook node
- attributes = {'format': 'html', 'source': 'nb_path'}
- nb_node = notebook_node('', evaluated_text, **attributes)
- (nb_node.source, nb_node.line) = \
- self.state_machine.get_source_and_line(self.lineno)
-
- # add dependency
- self.state.document.settings.record_dependencies.add(nb_abs_path)
-
- # clean up
- os.chdir(cwd)
- shutil.rmtree(tmpdir, True)
-
- return [nb_node]
-
-
-class notebook_node(nodes.raw):
- pass
-
-def nb_to_python(nb_path):
- """convert notebook to python script"""
- exporter = python.PythonExporter()
- output, resources = exporter.from_filename(nb_path)
- return output
-
-def nb_to_html(nb_path):
- """convert notebook to html"""
- c = Config({'ExtractOutputPreprocessor':{'enabled':True}})
-
- exporter = html.HTMLExporter(template_file='full', config=c)
- notebook = nbformat.read(open(nb_path), 'json')
- output, resources = exporter.from_notebook_node(notebook)
- header = output.split('<head>', 1)[1].split('</head>',1)[0]
- body = output.split('<body>', 1)[1].split('</body>',1)[0]
-
- # http://imgur.com/eR9bMRH
- header = header.replace('<style', '<style scoped="scoped"')
- header = header.replace('body {\n overflow: visible;\n padding: 8px;\n}\n',
- '')
- header = header.replace("code,pre{", "code{")
-
- # Filter out styles that conflict with the sphinx theme.
- filter_strings = [
- 'navbar',
- 'body{',
- 'alert{',
- 'uneditable-input{',
- 'collapse{',
- ]
-
- filter_strings.extend(['h%s{' % (i+1) for i in range(6)])
-
- line_begin = [
- 'pre{',
- 'p{margin'
- ]
-
- filterfunc = lambda x: not any([s in x for s in filter_strings])
- header_lines = filter(filterfunc, header.split('\n'))
-
- filterfunc = lambda x: not any([x.startswith(s) for s in line_begin])
- header_lines = filter(filterfunc, header_lines)
-
- header = '\n'.join(header_lines)
-
- # concatenate raw html lines
- lines = ['<div class="ipynotebook">']
- lines.append(header)
- lines.append(body)
- lines.append('</div>')
- return '\n'.join(lines), resources
-
-def evaluate_notebook(nb_path, dest_path=None, skip_exceptions=False):
- # Create evaluated version and save it to the dest path.
- notebook = nbformat.read(open(nb_path), 'json')
- nb_runner = NotebookRunner(notebook, pylab=False)
- try:
- nb_runner.run_notebook(skip_exceptions=skip_exceptions)
- except NotebookError as e:
- print('')
- print(e)
- # Return the traceback, filtering out ANSI color codes.
- # http://stackoverflow.com/questions/13506033/filtering-out-ansi-escape-sequences
- return "Notebook conversion failed with the " \
- "following traceback: \n%s" % \
- re.sub(r'\\033[\[\]]([0-9]{1,2}([;@][0-9]{0,2})*)*[mKP]?', '',
- str(e))
-
- if dest_path is None:
- dest_path = 'temp_evaluated.ipynb'
- nbformat.write(nb_runner.nb, open(dest_path, 'w'), 'json')
- ret = nb_to_html(dest_path)
- if dest_path is 'temp_evaluated.ipynb':
- os.remove(dest_path)
- return ret
-
-def formatted_link(path):
- return "`%s <%s>`__" % (os.path.basename(path), path)
-
-def visit_notebook_node(self, node):
- self.visit_raw(node)
-
-def depart_notebook_node(self, node):
- self.depart_raw(node)
-
-def setup(app):
- setup.app = app
- setup.config = app.config
- setup.confdir = app.confdir
-
- app.add_node(notebook_node,
- html=(visit_notebook_node, depart_notebook_node))
-
- app.add_directive('notebook', NotebookDirective)
-
- retdict = dict(
- version='0.1',
- parallel_read_safe=True,
- parallel_write_safe=True
- )
-
- return retdict
-
-def make_image_dir(setup, rst_dir):
- image_dir = setup.app.builder.outdir + os.path.sep + '_images'
- rel_dir = os.path.relpath(setup.confdir, rst_dir)
- image_rel_dir = rel_dir + os.path.sep + '_images'
- thread_safe_mkdir(image_dir)
- return image_dir, image_rel_dir
-
-def write_notebook_output(resources, image_dir, image_rel_dir, evaluated_text):
- my_uuid = uuid.uuid4().hex
-
- for output in resources['outputs']:
- new_name = image_dir + os.path.sep + my_uuid + output
- new_relative_name = image_rel_dir + os.path.sep + my_uuid + output
- evaluated_text = evaluated_text.replace(output, new_relative_name)
- with open(new_name, 'wb') as f:
- f.write(resources['outputs'][output])
- return evaluated_text
-
-def thread_safe_mkdir(dirname):
- try:
- os.makedirs(dirname)
- except OSError as e:
- if e.errno != errno.EEXIST:
- raise
- pass
diff -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 doc/extensions/notebookcell_sphinxext.py
--- a/doc/extensions/notebookcell_sphinxext.py
+++ /dev/null
@@ -1,87 +0,0 @@
-import os
-import shutil
-import io
-import tempfile
-from sphinx.util.compat import Directive
-from docutils.parsers.rst import directives
-from IPython.nbformat import current
-from notebook_sphinxext import \
- notebook_node, visit_notebook_node, depart_notebook_node, \
- evaluate_notebook, make_image_dir, write_notebook_output
-
-
-class NotebookCellDirective(Directive):
- """Insert an evaluated notebook cell into a document
-
- This uses runipy and nbconvert to transform an inline python
- script into html suitable for embedding in a Sphinx document.
- """
- required_arguments = 0
- optional_arguments = 1
- has_content = True
- option_spec = {'skip_exceptions': directives.flag}
-
- def run(self):
- # check if raw html is supported
- if not self.state.document.settings.raw_enabled:
- raise self.warning('"%s" directive disabled.' % self.name)
-
- cwd = os.getcwd()
- tmpdir = tempfile.mkdtemp()
- os.chdir(tmpdir)
-
- rst_file = self.state_machine.document.attributes['source']
- rst_dir = os.path.abspath(os.path.dirname(rst_file))
-
- image_dir, image_rel_dir = make_image_dir(setup, rst_dir)
-
- # Construct notebook from cell content
- content = "\n".join(self.content)
- with open("temp.py", "w") as f:
- f.write(content)
-
- convert_to_ipynb('temp.py', 'temp.ipynb')
-
- skip_exceptions = 'skip_exceptions' in self.options
-
- evaluated_text, resources = evaluate_notebook(
- 'temp.ipynb', skip_exceptions=skip_exceptions)
-
- evaluated_text = write_notebook_output(
- resources, image_dir, image_rel_dir, evaluated_text)
-
- # create notebook node
- attributes = {'format': 'html', 'source': 'nb_path'}
- nb_node = notebook_node('', evaluated_text, **attributes)
- (nb_node.source, nb_node.line) = \
- self.state_machine.get_source_and_line(self.lineno)
-
- # clean up
- os.chdir(cwd)
- shutil.rmtree(tmpdir, True)
-
- return [nb_node]
-
-def setup(app):
- setup.app = app
- setup.config = app.config
- setup.confdir = app.confdir
-
- app.add_node(notebook_node,
- html=(visit_notebook_node, depart_notebook_node))
-
- app.add_directive('notebook-cell', NotebookCellDirective)
-
- retdict = dict(
- version='0.1',
- parallel_read_safe=True,
- parallel_write_safe=True
- )
-
- return retdict
-
-def convert_to_ipynb(py_file, ipynb_file):
- with io.open(py_file, 'r', encoding='utf-8') as f:
- notebook = current.reads(f.read(), format='py')
- with io.open(ipynb_file, 'w', encoding='utf-8') as f:
- current.write(notebook, f, format='ipynb')
diff -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 doc/extensions/numpydocmod/__init__.py
--- a/doc/extensions/numpydocmod/__init__.py
+++ /dev/null
@@ -1,1 +0,0 @@
-from numpydoc import setup
diff -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 doc/extensions/numpydocmod/comment_eater.py
--- a/doc/extensions/numpydocmod/comment_eater.py
+++ /dev/null
@@ -1,158 +0,0 @@
-from cStringIO import StringIO
-import compiler
-import inspect
-import textwrap
-import tokenize
-
-from compiler_unparse import unparse
-
-
-class Comment(object):
- """ A comment block.
- """
- is_comment = True
- def __init__(self, start_lineno, end_lineno, text):
- # int : The first line number in the block. 1-indexed.
- self.start_lineno = start_lineno
- # int : The last line number. Inclusive!
- self.end_lineno = end_lineno
- # str : The text block including '#' character but not any leading spaces.
- self.text = text
-
- def add(self, string, start, end, line):
- """ Add a new comment line.
- """
- self.start_lineno = min(self.start_lineno, start[0])
- self.end_lineno = max(self.end_lineno, end[0])
- self.text += string
-
- def __repr__(self):
- return '%s(%r, %r, %r)' % (self.__class__.__name__, self.start_lineno,
- self.end_lineno, self.text)
-
-
-class NonComment(object):
- """ A non-comment block of code.
- """
- is_comment = False
- def __init__(self, start_lineno, end_lineno):
- self.start_lineno = start_lineno
- self.end_lineno = end_lineno
-
- def add(self, string, start, end, line):
- """ Add lines to the block.
- """
- if string.strip():
- # Only add if not entirely whitespace.
- self.start_lineno = min(self.start_lineno, start[0])
- self.end_lineno = max(self.end_lineno, end[0])
-
- def __repr__(self):
- return '%s(%r, %r)' % (self.__class__.__name__, self.start_lineno,
- self.end_lineno)
-
-
-class CommentBlocker(object):
- """ Pull out contiguous comment blocks.
- """
- def __init__(self):
- # Start with a dummy.
- self.current_block = NonComment(0, 0)
-
- # All of the blocks seen so far.
- self.blocks = []
-
- # The index mapping lines of code to their associated comment blocks.
- self.index = {}
-
- def process_file(self, file):
- """ Process a file object.
- """
- for token in tokenize.generate_tokens(file.next):
- self.process_token(*token)
- self.make_index()
-
- def process_token(self, kind, string, start, end, line):
- """ Process a single token.
- """
- if self.current_block.is_comment:
- if kind == tokenize.COMMENT:
- self.current_block.add(string, start, end, line)
- else:
- self.new_noncomment(start[0], end[0])
- else:
- if kind == tokenize.COMMENT:
- self.new_comment(string, start, end, line)
- else:
- self.current_block.add(string, start, end, line)
-
- def new_noncomment(self, start_lineno, end_lineno):
- """ We are transitioning from a noncomment to a comment.
- """
- block = NonComment(start_lineno, end_lineno)
- self.blocks.append(block)
- self.current_block = block
-
- def new_comment(self, string, start, end, line):
- """ Possibly add a new comment.
-
- Only adds a new comment if this comment is the only thing on the line.
- Otherwise, it extends the noncomment block.
- """
- prefix = line[:start[1]]
- if prefix.strip():
- # Oops! Trailing comment, not a comment block.
- self.current_block.add(string, start, end, line)
- else:
- # A comment block.
- block = Comment(start[0], end[0], string)
- self.blocks.append(block)
- self.current_block = block
-
- def make_index(self):
- """ Make the index mapping lines of actual code to their associated
- prefix comments.
- """
- for prev, block in zip(self.blocks[:-1], self.blocks[1:]):
- if not block.is_comment:
- self.index[block.start_lineno] = prev
-
- def search_for_comment(self, lineno, default=None):
- """ Find the comment block just before the given line number.
-
- Returns None (or the specified default) if there is no such block.
- """
- if not self.index:
- self.make_index()
- block = self.index.get(lineno, None)
- text = getattr(block, 'text', default)
- return text
-
-
-def strip_comment_marker(text):
- """ Strip # markers at the front of a block of comment text.
- """
- lines = []
- for line in text.splitlines():
- lines.append(line.lstrip('#'))
- text = textwrap.dedent('\n'.join(lines))
- return text
-
-
-def get_class_traits(klass):
- """ Yield all of the documentation for trait definitions on a class object.
- """
- # FIXME: gracefully handle errors here or in the caller?
- source = inspect.getsource(klass)
- cb = CommentBlocker()
- cb.process_file(StringIO(source))
- mod_ast = compiler.parse(source)
- class_ast = mod_ast.node.nodes[0]
- for node in class_ast.code.nodes:
- # FIXME: handle other kinds of assignments?
- if isinstance(node, compiler.ast.Assign):
- name = node.nodes[0].name
- rhs = unparse(node.expr).strip()
- doc = strip_comment_marker(cb.search_for_comment(node.lineno, default=''))
- yield name, rhs, doc
-
diff -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 doc/extensions/numpydocmod/compiler_unparse.py
--- a/doc/extensions/numpydocmod/compiler_unparse.py
+++ /dev/null
@@ -1,860 +0,0 @@
-""" Turn compiler.ast structures back into executable python code.
-
- The unparse method takes a compiler.ast tree and transforms it back into
- valid python code. It is incomplete and currently only works for
- import statements, function calls, function definitions, assignments, and
- basic expressions.
-
- Inspired by python-2.5-svn/Demo/parser/unparse.py
-
- fixme: We may want to move to using _ast trees because the compiler for
- them is about 6 times faster than compiler.compile.
-"""
-
-import sys
-import cStringIO
-from compiler.ast import Const, Name, Tuple, Div, Mul, Sub, Add
-
-def unparse(ast, single_line_functions=False):
- s = cStringIO.StringIO()
- UnparseCompilerAst(ast, s, single_line_functions)
- return s.getvalue().lstrip()
-
-op_precedence = { 'compiler.ast.Power':3, 'compiler.ast.Mul':2, 'compiler.ast.Div':2,
- 'compiler.ast.Add':1, 'compiler.ast.Sub':1 }
-
-class UnparseCompilerAst:
- """ Methods in this class recursively traverse an AST and
- output source code for the abstract syntax; original formatting
- is disregarged.
- """
-
- #########################################################################
- # object interface.
- #########################################################################
-
- def __init__(self, tree, file = sys.stdout, single_line_functions=False):
- """ Unparser(tree, file=sys.stdout) -> None.
-
- Print the source for tree to file.
- """
- self.f = file
- self._single_func = single_line_functions
- self._do_indent = True
- self._indent = 0
- self._dispatch(tree)
- self._write("\n")
- self.f.flush()
-
- #########################################################################
- # Unparser private interface.
- #########################################################################
-
- ### format, output, and dispatch methods ################################
-
- def _fill(self, text = ""):
- "Indent a piece of text, according to the current indentation level"
- if self._do_indent:
- self._write("\n"+" "*self._indent + text)
- else:
- self._write(text)
-
- def _write(self, text):
- "Append a piece of text to the current line."
- self.f.write(text)
-
- def _enter(self):
- "Print ':', and increase the indentation."
- self._write(": ")
- self._indent += 1
-
- def _leave(self):
- "Decrease the indentation level."
- self._indent -= 1
-
- def _dispatch(self, tree):
- "_dispatcher function, _dispatching tree type T to method _T."
- if isinstance(tree, list):
- for t in tree:
- self._dispatch(t)
- return
- meth = getattr(self, "_"+tree.__class__.__name__)
- if tree.__class__.__name__ == 'NoneType' and not self._do_indent:
- return
- meth(tree)
-
-
- #########################################################################
- # compiler.ast unparsing methods.
- #
- # There should be one method per concrete grammar type. They are
- # organized in alphabetical order.
- #########################################################################
-
- def _Add(self, t):
- self.__binary_op(t, '+')
-
- def _And(self, t):
- self._write(" (")
- for i, node in enumerate(t.nodes):
- self._dispatch(node)
- if i != len(t.nodes)-1:
- self._write(") and (")
- self._write(")")
-
- def _AssAttr(self, t):
- """ Handle assigning an attribute of an object
- """
- self._dispatch(t.expr)
- self._write('.'+t.attrname)
-
- def _Assign(self, t):
- """ Expression Assignment such as "a = 1".
-
- This only handles assignment in expressions. Keyword assignment
- is handled separately.
- """
- self._fill()
- for target in t.nodes:
- self._dispatch(target)
- self._write(" = ")
- self._dispatch(t.expr)
- if not self._do_indent:
- self._write('; ')
-
- def _AssName(self, t):
- """ Name on left hand side of expression.
-
- Treat just like a name on the right side of an expression.
- """
- self._Name(t)
-
- def _AssTuple(self, t):
- """ Tuple on left hand side of an expression.
- """
-
- # _write each elements, separated by a comma.
- for element in t.nodes[:-1]:
- self._dispatch(element)
- self._write(", ")
-
- # Handle the last one without writing comma
- last_element = t.nodes[-1]
- self._dispatch(last_element)
-
- def _AugAssign(self, t):
- """ +=,-=,*=,/=,**=, etc. operations
- """
-
- self._fill()
- self._dispatch(t.node)
- self._write(' '+t.op+' ')
- self._dispatch(t.expr)
- if not self._do_indent:
- self._write(';')
-
- def _Bitand(self, t):
- """ Bit and operation.
- """
-
- for i, node in enumerate(t.nodes):
- self._write("(")
- self._dispatch(node)
- self._write(")")
- if i != len(t.nodes)-1:
- self._write(" & ")
-
- def _Bitor(self, t):
- """ Bit or operation
- """
-
- for i, node in enumerate(t.nodes):
- self._write("(")
- self._dispatch(node)
- self._write(")")
- if i != len(t.nodes)-1:
- self._write(" | ")
-
- def _CallFunc(self, t):
- """ Function call.
- """
- self._dispatch(t.node)
- self._write("(")
- comma = False
- for e in t.args:
- if comma: self._write(", ")
- else: comma = True
- self._dispatch(e)
- if t.star_args:
- if comma: self._write(", ")
- else: comma = True
- self._write("*")
- self._dispatch(t.star_args)
- if t.dstar_args:
- if comma: self._write(", ")
- else: comma = True
- self._write("**")
- self._dispatch(t.dstar_args)
- self._write(")")
-
- def _Compare(self, t):
- self._dispatch(t.expr)
- for op, expr in t.ops:
- self._write(" " + op + " ")
- self._dispatch(expr)
-
- def _Const(self, t):
- """ A constant value such as an integer value, 3, or a string, "hello".
- """
- self._dispatch(t.value)
-
- def _Decorators(self, t):
- """ Handle function decorators (eg. @has_units)
- """
- for node in t.nodes:
- self._dispatch(node)
-
- def _Dict(self, t):
- self._write("{")
- for i, (k, v) in enumerate(t.items):
- self._dispatch(k)
- self._write(": ")
- self._dispatch(v)
- if i < len(t.items)-1:
- self._write(", ")
- self._write("}")
-
- def _Discard(self, t):
- """ Node for when return value is ignored such as in "foo(a)".
- """
- self._fill()
- self._dispatch(t.expr)
-
- def _Div(self, t):
- self.__binary_op(t, '/')
-
- def _Ellipsis(self, t):
- self._write("...")
-
- def _From(self, t):
- """ Handle "from xyz import foo, bar as baz".
- """
- # fixme: Are From and ImportFrom handled differently?
- self._fill("from ")
- self._write(t.modname)
- self._write(" import ")
- for i, (name,asname) in enumerate(t.names):
- if i != 0:
- self._write(", ")
- self._write(name)
- if asname is not None:
- self._write(" as "+asname)
-
- def _Function(self, t):
- """ Handle function definitions
- """
- if t.decorators is not None:
- self._fill("@")
- self._dispatch(t.decorators)
- self._fill("def "+t.name + "(")
- defaults = [None] * (len(t.argnames) - len(t.defaults)) + list(t.defaults)
- for i, arg in enumerate(zip(t.argnames, defaults)):
- self._write(arg[0])
- if arg[1] is not None:
- self._write('=')
- self._dispatch(arg[1])
- if i < len(t.argnames)-1:
- self._write(', ')
- self._write(")")
- if self._single_func:
- self._do_indent = False
- self._enter()
- self._dispatch(t.code)
- self._leave()
- self._do_indent = True
-
- def _Getattr(self, t):
- """ Handle getting an attribute of an object
- """
- if isinstance(t.expr, (Div, Mul, Sub, Add)):
- self._write('(')
- self._dispatch(t.expr)
- self._write(')')
- else:
- self._dispatch(t.expr)
-
- self._write('.'+t.attrname)
-
- def _If(self, t):
- self._fill()
-
- for i, (compare,code) in enumerate(t.tests):
- if i == 0:
- self._write("if ")
- else:
- self._write("elif ")
- self._dispatch(compare)
- self._enter()
- self._fill()
- self._dispatch(code)
- self._leave()
- self._write("\n")
-
- if t.else_ is not None:
- self._write("else")
- self._enter()
- self._fill()
- self._dispatch(t.else_)
- self._leave()
- self._write("\n")
-
- def _IfExp(self, t):
- self._dispatch(t.then)
- self._write(" if ")
- self._dispatch(t.test)
-
- if t.else_ is not None:
- self._write(" else (")
- self._dispatch(t.else_)
- self._write(")")
-
- def _Import(self, t):
- """ Handle "import xyz.foo".
- """
- self._fill("import ")
-
- for i, (name,asname) in enumerate(t.names):
- if i != 0:
- self._write(", ")
- self._write(name)
- if asname is not None:
- self._write(" as "+asname)
-
- def _Keyword(self, t):
- """ Keyword value assignment within function calls and definitions.
- """
- self._write(t.name)
- self._write("=")
- self._dispatch(t.expr)
-
- def _List(self, t):
- self._write("[")
- for i,node in enumerate(t.nodes):
- self._dispatch(node)
- if i < len(t.nodes)-1:
- self._write(", ")
- self._write("]")
-
- def _Module(self, t):
- if t.doc is not None:
- self._dispatch(t.doc)
- self._dispatch(t.node)
-
- def _Mul(self, t):
- self.__binary_op(t, '*')
-
- def _Name(self, t):
- self._write(t.name)
-
- def _NoneType(self, t):
- self._write("None")
-
- def _Not(self, t):
- self._write('not (')
- self._dispatch(t.expr)
- self._write(')')
-
- def _Or(self, t):
- self._write(" (")
- for i, node in enumerate(t.nodes):
- self._dispatch(node)
- if i != len(t.nodes)-1:
- self._write(") or (")
- self._write(")")
-
- def _Pass(self, t):
- self._write("pass\n")
-
- def _Printnl(self, t):
- self._fill("print ")
- if t.dest:
- self._write(">> ")
- self._dispatch(t.dest)
- self._write(", ")
- comma = False
- for node in t.nodes:
- if comma: self._write(', ')
- else: comma = True
- self._dispatch(node)
-
- def _Power(self, t):
- self.__binary_op(t, '**')
-
- def _Return(self, t):
- self._fill("return ")
- if t.value:
- if isinstance(t.value, Tuple):
- text = ', '.join([ name.name for name in t.value.asList() ])
- self._write(text)
- else:
- self._dispatch(t.value)
- if not self._do_indent:
- self._write('; ')
-
- def _Slice(self, t):
- self._dispatch(t.expr)
- self._write("[")
- if t.lower:
- self._dispatch(t.lower)
- self._write(":")
- if t.upper:
- self._dispatch(t.upper)
- #if t.step:
- # self._write(":")
- # self._dispatch(t.step)
- self._write("]")
-
- def _Sliceobj(self, t):
- for i, node in enumerate(t.nodes):
- if i != 0:
- self._write(":")
- if not (isinstance(node, Const) and node.value is None):
- self._dispatch(node)
-
- def _Stmt(self, tree):
- for node in tree.nodes:
- self._dispatch(node)
-
- def _Sub(self, t):
- self.__binary_op(t, '-')
-
- def _Subscript(self, t):
- self._dispatch(t.expr)
- self._write("[")
- for i, value in enumerate(t.subs):
- if i != 0:
- self._write(",")
- self._dispatch(value)
- self._write("]")
-
- def _TryExcept(self, t):
- self._fill("try")
- self._enter()
- self._dispatch(t.body)
- self._leave()
-
- for handler in t.handlers:
- self._fill('except ')
- self._dispatch(handler[0])
- if handler[1] is not None:
- self._write(', ')
- self._dispatch(handler[1])
- self._enter()
- self._dispatch(handler[2])
- self._leave()
-
- if t.else_:
- self._fill("else")
- self._enter()
- self._dispatch(t.else_)
- self._leave()
-
- def _Tuple(self, t):
-
- if not t.nodes:
- # Empty tuple.
- self._write("()")
- else:
- self._write("(")
-
- # _write each elements, separated by a comma.
- for element in t.nodes[:-1]:
- self._dispatch(element)
- self._write(", ")
-
- # Handle the last one without writing comma
- last_element = t.nodes[-1]
- self._dispatch(last_element)
-
- self._write(")")
-
- def _UnaryAdd(self, t):
- self._write("+")
- self._dispatch(t.expr)
-
- def _UnarySub(self, t):
- self._write("-")
- self._dispatch(t.expr)
-
- def _With(self, t):
- self._fill('with ')
- self._dispatch(t.expr)
- if t.vars:
- self._write(' as ')
- self._dispatch(t.vars.name)
- self._enter()
- self._dispatch(t.body)
- self._leave()
- self._write('\n')
-
- def _int(self, t):
- self._write(repr(t))
-
- def __binary_op(self, t, symbol):
- # Check if parenthesis are needed on left side and then dispatch
- has_paren = False
- left_class = str(t.left.__class__)
- if (left_class in op_precedence.keys() and
- op_precedence[left_class] < op_precedence[str(t.__class__)]):
- has_paren = True
- if has_paren:
- self._write('(')
- self._dispatch(t.left)
- if has_paren:
- self._write(')')
- # Write the appropriate symbol for operator
- self._write(symbol)
- # Check if parenthesis are needed on the right side and then dispatch
- has_paren = False
- right_class = str(t.right.__class__)
- if (right_class in op_precedence.keys() and
- op_precedence[right_class] < op_precedence[str(t.__class__)]):
- has_paren = True
- if has_paren:
- self._write('(')
- self._dispatch(t.right)
- if has_paren:
- self._write(')')
-
- def _float(self, t):
- # if t is 0.1, str(t)->'0.1' while repr(t)->'0.1000000000001'
- # We prefer str here.
- self._write(str(t))
-
- def _str(self, t):
- self._write(repr(t))
-
- def _tuple(self, t):
- self._write(str(t))
-
- #########################################################################
- # These are the methods from the _ast modules unparse.
- #
- # As our needs to handle more advanced code increase, we may want to
- # modify some of the methods below so that they work for compiler.ast.
- #########################################################################
-
-# # stmt
-# def _Expr(self, tree):
-# self._fill()
-# self._dispatch(tree.value)
-#
-# def _Import(self, t):
-# self._fill("import ")
-# first = True
-# for a in t.names:
-# if first:
-# first = False
-# else:
-# self._write(", ")
-# self._write(a.name)
-# if a.asname:
-# self._write(" as "+a.asname)
-#
-## def _ImportFrom(self, t):
-## self._fill("from ")
-## self._write(t.module)
-## self._write(" import ")
-## for i, a in enumerate(t.names):
-## if i == 0:
-## self._write(", ")
-## self._write(a.name)
-## if a.asname:
-## self._write(" as "+a.asname)
-## # XXX(jpe) what is level for?
-##
-#
-# def _Break(self, t):
-# self._fill("break")
-#
-# def _Continue(self, t):
-# self._fill("continue")
-#
-# def _Delete(self, t):
-# self._fill("del ")
-# self._dispatch(t.targets)
-#
-# def _Assert(self, t):
-# self._fill("assert ")
-# self._dispatch(t.test)
-# if t.msg:
-# self._write(", ")
-# self._dispatch(t.msg)
-#
-# def _Exec(self, t):
-# self._fill("exec ")
-# self._dispatch(t.body)
-# if t.globals:
-# self._write(" in ")
-# self._dispatch(t.globals)
-# if t.locals:
-# self._write(", ")
-# self._dispatch(t.locals)
-#
-# def _Print(self, t):
-# self._fill("print ")
-# do_comma = False
-# if t.dest:
-# self._write(">>")
-# self._dispatch(t.dest)
-# do_comma = True
-# for e in t.values:
-# if do_comma:self._write(", ")
-# else:do_comma=True
-# self._dispatch(e)
-# if not t.nl:
-# self._write(",")
-#
-# def _Global(self, t):
-# self._fill("global")
-# for i, n in enumerate(t.names):
-# if i != 0:
-# self._write(",")
-# self._write(" " + n)
-#
-# def _Yield(self, t):
-# self._fill("yield")
-# if t.value:
-# self._write(" (")
-# self._dispatch(t.value)
-# self._write(")")
-#
-# def _Raise(self, t):
-# self._fill('raise ')
-# if t.type:
-# self._dispatch(t.type)
-# if t.inst:
-# self._write(", ")
-# self._dispatch(t.inst)
-# if t.tback:
-# self._write(", ")
-# self._dispatch(t.tback)
-#
-#
-# def _TryFinally(self, t):
-# self._fill("try")
-# self._enter()
-# self._dispatch(t.body)
-# self._leave()
-#
-# self._fill("finally")
-# self._enter()
-# self._dispatch(t.finalbody)
-# self._leave()
-#
-# def _excepthandler(self, t):
-# self._fill("except ")
-# if t.type:
-# self._dispatch(t.type)
-# if t.name:
-# self._write(", ")
-# self._dispatch(t.name)
-# self._enter()
-# self._dispatch(t.body)
-# self._leave()
-#
-# def _ClassDef(self, t):
-# self._write("\n")
-# self._fill("class "+t.name)
-# if t.bases:
-# self._write("(")
-# for a in t.bases:
-# self._dispatch(a)
-# self._write(", ")
-# self._write(")")
-# self._enter()
-# self._dispatch(t.body)
-# self._leave()
-#
-# def _FunctionDef(self, t):
-# self._write("\n")
-# for deco in t.decorators:
-# self._fill("@")
-# self._dispatch(deco)
-# self._fill("def "+t.name + "(")
-# self._dispatch(t.args)
-# self._write(")")
-# self._enter()
-# self._dispatch(t.body)
-# self._leave()
-#
-# def _For(self, t):
-# self._fill("for ")
-# self._dispatch(t.target)
-# self._write(" in ")
-# self._dispatch(t.iter)
-# self._enter()
-# self._dispatch(t.body)
-# self._leave()
-# if t.orelse:
-# self._fill("else")
-# self._enter()
-# self._dispatch(t.orelse)
-# self._leave
-#
-# def _While(self, t):
-# self._fill("while ")
-# self._dispatch(t.test)
-# self._enter()
-# self._dispatch(t.body)
-# self._leave()
-# if t.orelse:
-# self._fill("else")
-# self._enter()
-# self._dispatch(t.orelse)
-# self._leave
-#
-# # expr
-# def _Str(self, tree):
-# self._write(repr(tree.s))
-##
-# def _Repr(self, t):
-# self._write("`")
-# self._dispatch(t.value)
-# self._write("`")
-#
-# def _Num(self, t):
-# self._write(repr(t.n))
-#
-# def _ListComp(self, t):
-# self._write("[")
-# self._dispatch(t.elt)
-# for gen in t.generators:
-# self._dispatch(gen)
-# self._write("]")
-#
-# def _GeneratorExp(self, t):
-# self._write("(")
-# self._dispatch(t.elt)
-# for gen in t.generators:
-# self._dispatch(gen)
-# self._write(")")
-#
-# def _comprehension(self, t):
-# self._write(" for ")
-# self._dispatch(t.target)
-# self._write(" in ")
-# self._dispatch(t.iter)
-# for if_clause in t.ifs:
-# self._write(" if ")
-# self._dispatch(if_clause)
-#
-# def _IfExp(self, t):
-# self._dispatch(t.body)
-# self._write(" if ")
-# self._dispatch(t.test)
-# if t.orelse:
-# self._write(" else ")
-# self._dispatch(t.orelse)
-#
-# unop = {"Invert":"~", "Not": "not", "UAdd":"+", "USub":"-"}
-# def _UnaryOp(self, t):
-# self._write(self.unop[t.op.__class__.__name__])
-# self._write("(")
-# self._dispatch(t.operand)
-# self._write(")")
-#
-# binop = { "Add":"+", "Sub":"-", "Mult":"*", "Div":"/", "Mod":"%",
-# "LShift":">>", "RShift":"<<", "BitOr":"|", "BitXor":"^", "BitAnd":"&",
-# "FloorDiv":"//", "Pow": "**"}
-# def _BinOp(self, t):
-# self._write("(")
-# self._dispatch(t.left)
-# self._write(")" + self.binop[t.op.__class__.__name__] + "(")
-# self._dispatch(t.right)
-# self._write(")")
-#
-# boolops = {_ast.And: 'and', _ast.Or: 'or'}
-# def _BoolOp(self, t):
-# self._write("(")
-# self._dispatch(t.values[0])
-# for v in t.values[1:]:
-# self._write(" %s " % self.boolops[t.op.__class__])
-# self._dispatch(v)
-# self._write(")")
-#
-# def _Attribute(self,t):
-# self._dispatch(t.value)
-# self._write(".")
-# self._write(t.attr)
-#
-## def _Call(self, t):
-## self._dispatch(t.func)
-## self._write("(")
-## comma = False
-## for e in t.args:
-## if comma: self._write(", ")
-## else: comma = True
-## self._dispatch(e)
-## for e in t.keywords:
-## if comma: self._write(", ")
-## else: comma = True
-## self._dispatch(e)
-## if t.starargs:
-## if comma: self._write(", ")
-## else: comma = True
-## self._write("*")
-## self._dispatch(t.starargs)
-## if t.kwargs:
-## if comma: self._write(", ")
-## else: comma = True
-## self._write("**")
-## self._dispatch(t.kwargs)
-## self._write(")")
-#
-# # slice
-# def _Index(self, t):
-# self._dispatch(t.value)
-#
-# def _ExtSlice(self, t):
-# for i, d in enumerate(t.dims):
-# if i != 0:
-# self._write(': ')
-# self._dispatch(d)
-#
-# # others
-# def _arguments(self, t):
-# first = True
-# nonDef = len(t.args)-len(t.defaults)
-# for a in t.args[0:nonDef]:
-# if first:first = False
-# else: self._write(", ")
-# self._dispatch(a)
-# for a,d in zip(t.args[nonDef:], t.defaults):
-# if first:first = False
-# else: self._write(", ")
-# self._dispatch(a),
-# self._write("=")
-# self._dispatch(d)
-# if t.vararg:
-# if first:first = False
-# else: self._write(", ")
-# self._write("*"+t.vararg)
-# if t.kwarg:
-# if first:first = False
-# else: self._write(", ")
-# self._write("**"+t.kwarg)
-#
-## def _keyword(self, t):
-## self._write(t.arg)
-## self._write("=")
-## self._dispatch(t.value)
-#
-# def _Lambda(self, t):
-# self._write("lambda ")
-# self._dispatch(t.args)
-# self._write(": ")
-# self._dispatch(t.body)
-
-
-
diff -r 2738eb221cb9fa51dbeaf7498c3385b730f091ce -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 doc/extensions/numpydocmod/docscrape.py
--- a/doc/extensions/numpydocmod/docscrape.py
+++ /dev/null
@@ -1,500 +0,0 @@
-"""Extract reference documentation from the NumPy source tree.
-
-"""
-
-import inspect
-import textwrap
-import re
-import pydoc
-from StringIO import StringIO
-from warnings import warn
-
-class Reader(object):
- """A line-based string reader.
-
- """
- def __init__(self, data):
- """
- Parameters
- ----------
- data : str
- String with lines separated by '\n'.
-
- """
- if isinstance(data,list):
- self._str = data
- else:
- self._str = data.split('\n') # store string as list of lines
-
- self.reset()
-
- def __getitem__(self, n):
- return self._str[n]
-
- def reset(self):
- self._l = 0 # current line nr
-
- def read(self):
- if not self.eof():
- out = self[self._l]
- self._l += 1
- return out
- else:
- return ''
-
- def seek_next_non_empty_line(self):
- for l in self[self._l:]:
- if l.strip():
- break
- else:
- self._l += 1
-
- def eof(self):
- return self._l >= len(self._str)
-
- def read_to_condition(self, condition_func):
- start = self._l
- for line in self[start:]:
- if condition_func(line):
- return self[start:self._l]
- self._l += 1
- if self.eof():
- return self[start:self._l+1]
- return []
-
- def read_to_next_empty_line(self):
- self.seek_next_non_empty_line()
- def is_empty(line):
- return not line.strip()
- return self.read_to_condition(is_empty)
-
- def read_to_next_unindented_line(self):
- def is_unindented(line):
- return (line.strip() and (len(line.lstrip()) == len(line)))
- return self.read_to_condition(is_unindented)
-
- def peek(self,n=0):
- if self._l + n < len(self._str):
- return self[self._l + n]
- else:
- return ''
-
- def is_empty(self):
- return not ''.join(self._str).strip()
-
-
-class NumpyDocString(object):
- def __init__(self, docstring, config={}):
- docstring = textwrap.dedent(docstring).split('\n')
-
- self._doc = Reader(docstring)
- self._parsed_data = {
- 'Signature': '',
- 'Summary': [''],
- 'Extended Summary': [],
- 'Parameters': [],
- 'Returns': [],
- 'Raises': [],
- 'Warns': [],
- 'Other Parameters': [],
- 'Attributes': [],
- 'Methods': [],
- 'See Also': [],
- 'Notes': [],
- 'Warnings': [],
- 'References': '',
- 'Examples': '',
- 'index': {}
- }
-
- self._parse()
-
- def __getitem__(self,key):
- return self._parsed_data[key]
-
- def __setitem__(self,key,val):
- if not self._parsed_data.has_key(key):
- warn("Unknown section %s" % key)
- else:
- self._parsed_data[key] = val
-
- def _is_at_section(self):
- self._doc.seek_next_non_empty_line()
-
- if self._doc.eof():
- return False
-
- l1 = self._doc.peek().strip() # e.g. Parameters
-
- if l1.startswith('.. index::'):
- return True
-
- l2 = self._doc.peek(1).strip() # ---------- or ==========
- return l2.startswith('-'*len(l1)) or l2.startswith('='*len(l1))
-
- def _strip(self,doc):
- i = 0
- j = 0
- for i,line in enumerate(doc):
- if line.strip(): break
-
- for j,line in enumerate(doc[::-1]):
- if line.strip(): break
-
- return doc[i:len(doc)-j]
-
- def _read_to_next_section(self):
- section = self._doc.read_to_next_empty_line()
-
- while not self._is_at_section() and not self._doc.eof():
- if not self._doc.peek(-1).strip(): # previous line was empty
- section += ['']
-
- section += self._doc.read_to_next_empty_line()
-
- return section
-
- def _read_sections(self):
- while not self._doc.eof():
- data = self._read_to_next_section()
- name = data[0].strip()
-
- if name.startswith('..'): # index section
- yield name, data[1:]
- elif len(data) < 2:
- yield StopIteration
- else:
- yield name, self._strip(data[2:])
-
- def _parse_param_list(self,content):
- r = Reader(content)
- params = []
- while not r.eof():
- header = r.read().strip()
- if ' : ' in header:
- arg_name, arg_type = header.split(' : ')[:2]
- else:
- arg_name, arg_type = header, ''
-
- desc = r.read_to_next_unindented_line()
- desc = dedent_lines(desc)
-
- params.append((arg_name,arg_type,desc))
-
- return params
-
-
- _name_rgx = re.compile(r"^\s*(:(?P<role>\w+):`(?P<name>[a-zA-Z0-9_.-]+)`|"
- r" (?P<name2>[a-zA-Z0-9_.-]+))\s*", re.X)
- def _parse_see_also(self, content):
- """
- func_name : Descriptive text
- continued text
- another_func_name : Descriptive text
- func_name1, func_name2, :meth:`func_name`, func_name3
-
- """
- items = []
-
- def parse_item_name(text):
- """Match ':role:`name`' or 'name'"""
- m = self._name_rgx.match(text)
- if m:
- g = m.groups()
- if g[1] is None:
- return g[3], None
- else:
- return g[2], g[1]
- raise ValueError("%s is not a item name" % text)
-
- def push_item(name, rest):
- if not name:
- return
- name, role = parse_item_name(name)
- items.append((name, list(rest), role))
- del rest[:]
-
- current_func = None
- rest = []
-
- for line in content:
- if not line.strip(): continue
-
- m = self._name_rgx.match(line)
- if m and line[m.end():].strip().startswith(':'):
- push_item(current_func, rest)
- current_func, line = line[:m.end()], line[m.end():]
- rest = [line.split(':', 1)[1].strip()]
- if not rest[0]:
- rest = []
- elif not line.startswith(' '):
- push_item(current_func, rest)
- current_func = None
- if ',' in line:
- for func in line.split(','):
- if func.strip():
- push_item(func, [])
- elif line.strip():
- current_func = line
- elif current_func is not None:
- rest.append(line.strip())
- push_item(current_func, rest)
- return items
-
- def _parse_index(self, section, content):
- """
- .. index: default
- :refguide: something, else, and more
-
- """
- def strip_each_in(lst):
- return [s.strip() for s in lst]
-
- out = {}
- section = section.split('::')
- if len(section) > 1:
- out['default'] = strip_each_in(section[1].split(','))[0]
- for line in content:
- line = line.split(':')
- if len(line) > 2:
- out[line[1]] = strip_each_in(line[2].split(','))
- return out
-
- def _parse_summary(self):
- """Grab signature (if given) and summary"""
- if self._is_at_section():
- return
-
- summary = self._doc.read_to_next_empty_line()
- summary_str = " ".join([s.strip() for s in summary]).strip()
- if re.compile('^([\w., ]+=)?\s*[\w\.]+\(.*\)$').match(summary_str):
- self['Signature'] = summary_str
- if not self._is_at_section():
- self['Summary'] = self._doc.read_to_next_empty_line()
- else:
- self['Summary'] = summary
-
- if not self._is_at_section():
- self['Extended Summary'] = self._read_to_next_section()
-
- def _parse(self):
- self._doc.reset()
- self._parse_summary()
-
- for (section,content) in self._read_sections():
- if not section.startswith('..'):
- section = ' '.join([s.capitalize() for s in section.split(' ')])
- if section in ('Parameters', 'Returns', 'Raises', 'Warns',
- 'Other Parameters', 'Attributes', 'Methods'):
- self[section] = self._parse_param_list(content)
- elif section.startswith('.. index::'):
- self['index'] = self._parse_index(section, content)
- elif section == 'See Also':
- self['See Also'] = self._parse_see_also(content)
- else:
- self[section] = content
-
- # string conversion routines
-
- def _str_header(self, name, symbol='-'):
- return [name, len(name)*symbol]
-
- def _str_indent(self, doc, indent=4):
- out = []
- for line in doc:
- out += [' '*indent + line]
- return out
-
- def _str_signature(self):
- if self['Signature']:
- return [self['Signature'].replace('*','\*')] + ['']
- else:
- return ['']
-
- def _str_summary(self):
- if self['Summary']:
- return self['Summary'] + ['']
- else:
- return []
-
- def _str_extended_summary(self):
- if self['Extended Summary']:
- return self['Extended Summary'] + ['']
- else:
- return []
-
- def _str_param_list(self, name):
- out = []
- if self[name]:
- out += self._str_header(name)
- for param,param_type,desc in self[name]:
- out += ['%s : %s' % (param, param_type)]
- out += self._str_indent(desc)
- out += ['']
- return out
-
- def _str_section(self, name):
- out = []
- if self[name]:
- out += self._str_header(name)
- out += self[name]
- out += ['']
- return out
-
- def _str_see_also(self, func_role):
- if not self['See Also']: return []
- out = []
- out += self._str_header("See Also")
- last_had_desc = True
- for func, desc, role in self['See Also']:
- if role:
- link = ':%s:`%s`' % (role, func)
- elif func_role:
- link = ':%s:`%s`' % (func_role, func)
- else:
- link = "`%s`_" % func
- if desc or last_had_desc:
- out += ['']
- out += [link]
- else:
- out[-1] += ", %s" % link
- if desc:
- out += self._str_indent([' '.join(desc)])
- last_had_desc = True
- else:
- last_had_desc = False
- out += ['']
- return out
-
- def _str_index(self):
- idx = self['index']
- out = []
- out += ['.. index:: %s' % idx.get('default','')]
- for section, references in idx.iteritems():
- if section == 'default':
- continue
- out += [' :%s: %s' % (section, ', '.join(references))]
- return out
-
- def __str__(self, func_role=''):
- out = []
- out += self._str_signature()
- out += self._str_summary()
- out += self._str_extended_summary()
- for param_list in ('Parameters', 'Returns', 'Other Parameters',
- 'Raises', 'Warns'):
- out += self._str_param_list(param_list)
- out += self._str_section('Warnings')
- out += self._str_see_also(func_role)
- for s in ('Notes','References','Examples'):
- out += self._str_section(s)
- for param_list in ('Attributes', 'Methods'):
- out += self._str_param_list(param_list)
- out += self._str_index()
- return '\n'.join(out)
-
-
-def indent(str,indent=4):
- indent_str = ' '*indent
- if str is None:
- return indent_str
- lines = str.split('\n')
- return '\n'.join(indent_str + l for l in lines)
-
-def dedent_lines(lines):
- """Deindent a list of lines maximally"""
- return textwrap.dedent("\n".join(lines)).split("\n")
-
-def header(text, style='-'):
- return text + '\n' + style*len(text) + '\n'
-
-
-class FunctionDoc(NumpyDocString):
- def __init__(self, func, role='func', doc=None, config={}):
- self._f = func
- self._role = role # e.g. "func" or "meth"
-
- if doc is None:
- if func is None:
- raise ValueError("No function or docstring given")
- doc = inspect.getdoc(func) or ''
- NumpyDocString.__init__(self, doc)
-
- if not self['Signature'] and func is not None:
- func, func_name = self.get_func()
- try:
- # try to read signature
- argspec = inspect.getargspec(func)
- argspec = inspect.formatargspec(*argspec)
- argspec = argspec.replace('*','\*')
- signature = '%s%s' % (func_name, argspec)
- except TypeError, e:
- signature = '%s()' % func_name
- self['Signature'] = signature
-
- def get_func(self):
- func_name = getattr(self._f, '__name__', self.__class__.__name__)
- if inspect.isclass(self._f):
- func = getattr(self._f, '__call__', self._f.__init__)
- else:
- func = self._f
- return func, func_name
-
- def __str__(self):
- out = ''
-
- func, func_name = self.get_func()
- signature = self['Signature'].replace('*', '\*')
-
- roles = {'func': 'function',
- 'meth': 'method'}
-
- if self._role:
- if not roles.has_key(self._role):
- print("Warning: invalid role %s" % self._role)
- out += '.. %s:: %s\n \n\n' % (roles.get(self._role,''),
- func_name)
-
- out += super(FunctionDoc, self).__str__(func_role=self._role)
- return out
-
-
-class ClassDoc(NumpyDocString):
- def __init__(self, cls, doc=None, modulename='', func_doc=FunctionDoc,
- config={}):
- if not inspect.isclass(cls) and cls is not None:
- raise ValueError("Expected a class or None, but got %r" % cls)
- self._cls = cls
-
- if modulename and not modulename.endswith('.'):
- modulename += '.'
- self._mod = modulename
-
- if doc is None:
- if cls is None:
- raise ValueError("No class or documentation string given")
- doc = pydoc.getdoc(cls)
-
- NumpyDocString.__init__(self, doc)
-
- if config.get('show_class_members', True):
- if not self['Methods']:
- self['Methods'] = [(name, '', '')
- for name in sorted(self.methods)]
- if not self['Attributes']:
- self['Attributes'] = [(name, '', '')
- for name in sorted(self.properties)]
-
- @property
- def methods(self):
- if self._cls is None:
- return []
- return [name for name,func in inspect.getmembers(self._cls)
- if not name.startswith('_') and callable(func)]
-
- @property
- def properties(self):
- if self._cls is None:
- return []
- return [name for name,func in inspect.getmembers(self._cls)
- if not name.startswith('_') and func is None]
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/d6b4e59ee4c6/
Changeset: d6b4e59ee4c6
Branch: yt
User: MatthewTurk
Date: 2016-02-02 20:14:57+00:00
Summary: Switch to F-ordering for the memoryviews.
Affected #: 1 file
diff -r 1b0f177a9c83133fc1fadc6bb9515af53d16d2a8 -r d6b4e59ee4c608f82ccabef58b1c929e6b15bcd0 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -164,7 +164,7 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- self.count[ii[0], ii[1], ii[2], offset] += 1
+ self.count[ii[2], ii[1], ii[0], offset] += 1
def finalize(self):
arr = np.asarray(self.count)
@@ -224,14 +224,14 @@
dist = idist[0] + idist[1] + idist[2]
# Calculate distance in multiples of the smoothing length
dist = sqrt(dist) / fields[0]
- self.temp[i,j,k,offset] = self.sph_kernel(dist)
- kernel_sum += self.temp[i,j,k,offset]
+ self.temp[k,j,i,offset] = self.sph_kernel(dist)
+ kernel_sum += self.temp[k,j,i,offset]
# Having found the kernel, deposit accordingly into gdata
for i from ib0[0] <= i <= ib1[0]:
for j from ib0[1] <= j <= ib1[1]:
for k from ib0[2] <= k <= ib1[2]:
- dist = self.temp[i,j,k,offset] / kernel_sum
- self.data[i,j,k,offset] += fields[1] * dist
+ dist = self.temp[k,j,i,offset] / kernel_sum
+ self.data[k,j,i,offset] += fields[1] * dist
def finalize(self):
return self.odata
@@ -257,7 +257,7 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.sum[ii[0], ii[1], ii[2], offset] += fields[0]
+ self.sum[ii[2], ii[1], ii[0], offset] += fields[0]
return
def finalize(self):
@@ -300,18 +300,18 @@
cdef float k, mk, qk
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- k = self.i[ii[0], ii[1], ii[2], offset]
- mk = self.mk[ii[0], ii[1], ii[2], offset]
- qk = self.qk[ii[0], ii[1], ii[2], offset]
+ k = self.i[ii[2], ii[1], ii[0], offset]
+ mk = self.mk[ii[2], ii[1], ii[0], offset]
+ qk = self.qk[ii[2], ii[1], ii[0], offset]
#print k, mk, qk, cell_index
if k == 0.0:
# Initialize cell values
- self.mk[ii[0], ii[1], ii[2], offset] = fields[0]
+ self.mk[ii[2], ii[1], ii[0], offset] = fields[0]
else:
- self.mk[ii[0], ii[1], ii[2], offset] = mk + (fields[0] - mk) / k
- self.qk[ii[0], ii[1], ii[2], offset] = \
+ self.mk[ii[2], ii[1], ii[0], offset] = mk + (fields[0] - mk) / k
+ self.qk[ii[2], ii[1], ii[0], offset] = \
qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
- self.i[ii[0], ii[1], ii[2], offset] += 1
+ self.i[ii[2], ii[1], ii[0], offset] += 1
def finalize(self):
# This is the standard variance
@@ -362,7 +362,7 @@
for i in range(2):
for j in range(2):
for k in range(2):
- self.field[ind[0] - i, ind[1] - j, ind[2] - k, offset] += \
+ self.field[ind[2] - k, ind[1] - j, ind[0] - i, offset] += \
fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
def finalize(self):
@@ -396,8 +396,8 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.w[ii[0], ii[1], ii[2], offset] += fields[1]
- self.wf[ii[0], ii[1], ii[2], offset] += fields[0] * fields[1]
+ self.w[ii[2], ii[1], ii[0], offset] += fields[1]
+ self.wf[ii[2], ii[1], ii[0], offset] += fields[0] * fields[1]
def finalize(self):
wf = np.asarray(self.owf)
@@ -465,9 +465,9 @@
r2 = ((ppos[0] - gpos[0])*(ppos[0] - gpos[0]) +
(ppos[1] - gpos[1])*(ppos[1] - gpos[1]) +
(ppos[2] - gpos[2])*(ppos[2] - gpos[2]))
- if r2 < self.distfield[i,j,k,offset]:
- self.distfield[i,j,k,offset] = r2
- self.nnfield[i,j,k,offset] = fields[0]
+ if r2 < self.distfield[k,j,i,offset]:
+ self.distfield[k,j,i,offset] = r2
+ self.nnfield[k,j,i,offset] = fields[0]
gpos[2] += dds[2]
gpos[1] += dds[1]
gpos[0] += dds[0]
https://bitbucket.org/yt_analysis/yt/commits/87aefc91df0d/
Changeset: 87aefc91df0d
Branch: yt
User: MatthewTurk
Date: 2016-02-02 22:45:10+00:00
Summary: Add decorators to turn off safety checks.
Affected #: 1 file
diff -r d6b4e59ee4c608f82ccabef58b1c929e6b15bcd0 -r 87aefc91df0d4b2164bb80b27f294c5342210d7d yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -49,6 +49,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
+ @cython.initializedcheck(False)
def process_octree(self, OctreeContainer octree,
np.ndarray[np.int64_t, ndim=1] dom_ind,
np.ndarray[np.float64_t, ndim=2] positions,
@@ -107,6 +108,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
+ @cython.initializedcheck(False)
def process_grid(self, gobj,
np.ndarray[np.float64_t, ndim=2] positions,
fields = None):
@@ -151,6 +153,9 @@
np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -187,6 +192,9 @@
np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -245,6 +253,9 @@
np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -287,6 +298,9 @@
np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -332,6 +346,9 @@
np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -384,6 +401,9 @@
np.zeros(self.nvals, dtype='float64', order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -416,6 +436,9 @@
self.update_values = 1
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
@@ -442,6 +465,9 @@
self.distfield[:] = np.inf
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
https://bitbucket.org/yt_analysis/yt/commits/187f84cc0441/
Changeset: 187f84cc0441
Branch: yt
User: atmyers
Date: 2016-02-03 17:09:06+00:00
Summary: Merged in MatthewTurk/yt (pull request #1919)
Convert particle deposit to memoryviews
Affected #: 6 files
diff -r 45dc609be9e893c48cdfc2e568eb26ca0a0f8cbe -r 187f84cc044157da920460ec8423f13d74659652 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -702,11 +702,11 @@
if cls is None:
raise YTParticleDepositionNotImplemented(method)
# We allocate number of zones, not number of octs
- op = cls(self.ActiveDimensions.prod(), kernel_name)
+ op = cls(self.ActiveDimensions, kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
- return vals.reshape(self.ActiveDimensions, order="C")
+ return vals.copy(order="C")
def write_to_gdf(self, gdf_path, fields, nprocs=1, field_units=None,
**kwargs):
diff -r 45dc609be9e893c48cdfc2e568eb26ca0a0f8cbe -r 187f84cc044157da920460ec8423f13d74659652 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -327,12 +327,12 @@
if cls is None:
raise YTParticleDepositionNotImplemented(method)
# We allocate number of zones, not number of octs
- op = cls(self.ActiveDimensions.prod(), kernel_name)
+ op = cls(tuple(self.ActiveDimensions), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
if vals is None: return
- return vals.reshape(self.ActiveDimensions, order="C")
+ return vals.copy(order="C")
def select_blocks(self, selector):
mask = self._get_selector_mask(selector)
diff -r 45dc609be9e893c48cdfc2e568eb26ca0a0f8cbe -r 187f84cc044157da920460ec8423f13d74659652 yt/frontends/artio/_artio_caller.pyx
--- a/yt/frontends/artio/_artio_caller.pyx
+++ b/yt/frontends/artio/_artio_caller.pyx
@@ -15,6 +15,7 @@
from libc.stdlib cimport malloc, free
from libc.string cimport memcpy
import data_structures
+from yt.utilities.lib.misc_utilities import OnceIndirect
cdef extern from "platform_dep.h":
void *alloca(int)
@@ -1554,6 +1555,9 @@
if fields is None:
fields = []
nf = len(fields)
+ cdef np.float64_t[::cython.view.indirect, ::1] field_pointers
+ if nf > 0: field_pointers = OnceIndirect(fields)
+ cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
cdef np.ndarray[np.uint8_t, ndim=1, cast=True] mask
mask = self.mask(selector, -1)
cdef np.ndarray[np.int64_t, ndim=1] domain_ind
@@ -1568,17 +1572,9 @@
continue
domain_ind[sfc - self.sfc_start] = j
j += 1
- cdef np.float64_t **field_pointers
- cdef np.float64_t *field_vals
cdef np.float64_t pos[3]
cdef np.float64_t left_edge[3]
cdef int coords[3]
- cdef np.ndarray[np.float64_t, ndim=1] tarr
- field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf)
- field_vals = <np.float64_t*>alloca(sizeof(np.float64_t) * nf)
- for i in range(nf):
- tarr = fields[i]
- field_pointers[i] = <np.float64_t *> tarr.data
cdef int dims[3]
dims[0] = dims[1] = dims[2] = 1
cdef np.int64_t offset, moff
diff -r 45dc609be9e893c48cdfc2e568eb26ca0a0f8cbe -r 187f84cc044157da920460ec8423f13d74659652 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -138,5 +138,5 @@
cdef public int update_values
cdef void process(self, int dim[3], np.float64_t left_edge[3],
np.float64_t dds[3], np.int64_t offset,
- np.float64_t ppos[3], np.float64_t *fields,
+ np.float64_t ppos[3], np.float64_t[:] fields,
np.int64_t domain_ind)
diff -r 45dc609be9e893c48cdfc2e568eb26ca0a0f8cbe -r 187f84cc044157da920460ec8423f13d74659652 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -19,10 +19,21 @@
from libc.stdlib cimport malloc, free
cimport cython
from libc.math cimport sqrt
-
+from cpython cimport PyObject
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, \
OctreeContainer, OctInfo
+from cpython.array cimport array, clone
+from cython.view cimport memoryview as cymemview
+from yt.utilities.lib.misc_utilities import OnceIndirect
+
+cdef append_axes(np.ndarray arr, int naxes):
+ if arr.ndim == naxes:
+ return arr
+ # Avoid copies
+ arr2 = arr.view()
+ arr2.shape = arr2.shape + (1,) * (naxes - arr2.ndim)
+ return arr2
cdef class ParticleDepositOperation:
def __init__(self, nvals, kernel_name):
@@ -38,6 +49,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
+ @cython.initializedcheck(False)
def process_octree(self, OctreeContainer octree,
np.ndarray[np.int64_t, ndim=1] dom_ind,
np.ndarray[np.float64_t, ndim=2] positions,
@@ -47,15 +59,10 @@
if fields is None:
fields = []
nf = len(fields)
- cdef np.float64_t **field_pointers
- cdef np.float64_t *field_vals
+ cdef np.float64_t[::cython.view.indirect, ::1] field_pointers
+ if nf > 0: field_pointers = OnceIndirect(fields)
cdef np.float64_t pos[3]
- cdef np.ndarray[np.float64_t, ndim=1] tarr
- field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf)
- field_vals = <np.float64_t*>alloca(sizeof(np.float64_t) * nf)
- for i in range(nf):
- tarr = fields[i]
- field_pointers[i] = <np.float64_t *> tarr.data
+ cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
cdef int dims[3]
dims[0] = dims[1] = dims[2] = (1 << octree.oref)
cdef int nz = dims[0] * dims[1] * dims[2]
@@ -67,7 +74,7 @@
for i in range(positions.shape[0]):
# We should check if particle remains inside the Oct here
for j in range(nf):
- field_vals[j] = field_pointers[j][i]
+ field_vals[j] = field_pointers[j,i]
for j in range(3):
pos[j] = positions[i, j]
# This line should be modified to have it return the index into an
@@ -90,7 +97,7 @@
if oct == NULL or (domain_id > 0 and oct.domain != domain_id):
continue
# Note that this has to be our local index, not our in-file index.
- offset = dom_ind[oct.domain_ind - moff] * nz
+ offset = dom_ind[oct.domain_ind - moff]
if offset < 0: continue
# Check that we found the oct ...
self.process(dims, oi.left_edge, oi.dds,
@@ -101,6 +108,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
+ @cython.initializedcheck(False)
def process_grid(self, gobj,
np.ndarray[np.float64_t, ndim=2] positions,
fields = None):
@@ -108,16 +116,11 @@
if fields is None:
fields = []
nf = len(fields)
- cdef np.float64_t **field_pointers
- cdef np.float64_t *field_vals
+ cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
+ cdef np.float64_t[::cython.view.indirect, ::1] field_pointers
+ if nf > 0: field_pointers = OnceIndirect(fields)
cdef np.float64_t pos[3]
- cdef np.ndarray[np.float64_t, ndim=1] tarr
- field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf)
- field_vals = <np.float64_t*>alloca(sizeof(np.float64_t) * nf)
cdef np.int64_t gid = getattr(gobj, "id", -1)
- for i in range(nf):
- tarr = fields[i]
- field_pointers[i] = <np.float64_t *> tarr.data
cdef np.float64_t dds[3]
cdef np.float64_t left_edge[3]
cdef int dims[3]
@@ -128,7 +131,7 @@
for i in range(positions.shape[0]):
# Now we process
for j in range(nf):
- field_vals[j] = field_pointers[j][i]
+ field_vals[j] = field_pointers[j,i]
for j in range(3):
pos[j] = positions[i, j]
self.process(dims, left_edge, dds, 0, pos, field_vals, gid)
@@ -138,27 +141,27 @@
cdef void process(self, int dim[3], np.float64_t left_edge[3],
np.float64_t dds[3], np.int64_t offset,
- np.float64_t ppos[3], np.float64_t *fields,
+ np.float64_t ppos[3], np.float64_t[:] fields,
np.int64_t domain_ind):
raise NotImplementedError
cdef class CountParticles(ParticleDepositOperation):
- cdef np.int64_t *count # float, for ease
- cdef public object ocount
+ cdef np.int64_t[:,:,:,:] count
def initialize(self):
# Create a numpy array accessible to python
- self.ocount = np.zeros(self.nvals, dtype="int64", order='F')
- cdef np.ndarray arr = self.ocount
- # alias the C-view for use in cython
- self.count = <np.int64_t*> arr.data
+ self.count = append_axes(
+ np.zeros(self.nvals, dtype="int64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset, # offset into IO field
np.float64_t ppos[3], # this particle's position
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
# here we do our thing; this is the kernel
@@ -166,10 +169,12 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- self.count[gind(ii[0], ii[1], ii[2], dim) + offset] += 1
+ self.count[ii[2], ii[1], ii[0], offset] += 1
def finalize(self):
- return self.ocount.astype('f8')
+ arr = np.asarray(self.count)
+ arr.shape = self.nvals
+ return arr.astype("float64")
deposit_count = CountParticles
@@ -177,26 +182,25 @@
# Note that this does nothing at the edges. So it will give a poor
# estimate there, and since Octrees are mostly edges, this will be a very
# poor SPH kernel.
- cdef np.float64_t *data
- cdef public object odata
- cdef np.float64_t *temp
- cdef public object otemp
+ cdef np.float64_t[:,:,:,:] data
+ cdef np.float64_t[:,:,:,:] temp
def initialize(self):
- self.odata = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.odata
- self.data = <np.float64_t*> arr.data
- self.otemp = np.zeros(self.nvals, dtype="float64", order='F')
- arr = self.otemp
- self.temp = <np.float64_t*> arr.data
+ self.data = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.temp = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
@@ -228,14 +232,14 @@
dist = idist[0] + idist[1] + idist[2]
# Calculate distance in multiples of the smoothing length
dist = sqrt(dist) / fields[0]
- self.temp[gind(i,j,k,dim) + offset] = self.sph_kernel(dist)
- kernel_sum += self.temp[gind(i,j,k,dim) + offset]
+ self.temp[k,j,i,offset] = self.sph_kernel(dist)
+ kernel_sum += self.temp[k,j,i,offset]
# Having found the kernel, deposit accordingly into gdata
for i from ib0[0] <= i <= ib1[0]:
for j from ib0[1] <= j <= ib1[1]:
for k from ib0[2] <= k <= ib1[2]:
- dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
- self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
+ dist = self.temp[k,j,i,offset] / kernel_sum
+ self.data[k,j,i,offset] += fields[1] * dist
def finalize(self):
return self.odata
@@ -243,31 +247,34 @@
deposit_simple_smooth = SimpleSmooth
cdef class SumParticleField(ParticleDepositOperation):
- cdef np.float64_t *sum
- cdef public object osum
+ cdef np.float64_t[:,:,:,:] sum
def initialize(self):
- self.osum = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.osum
- self.sum = <np.float64_t*> arr.data
+ self.sum = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.sum[gind(ii[0], ii[1], ii[2], dim) + offset] += fields[0]
+ self.sum[ii[2], ii[1], ii[0], offset] += fields[0]
return
def finalize(self):
- return self.osum
+ sum = np.asarray(self.sum)
+ sum.shape = self.nvals
+ return sum
deposit_sum = SumParticleField
@@ -275,36 +282,31 @@
# Thanks to Britton and MJ Turk for the link
# to a single-pass STD
# http://www.cs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
- cdef np.float64_t *mk
- cdef np.float64_t *qk
- cdef np.float64_t *i
- cdef public object omk
- cdef public object oqk
- cdef public object oi
+ cdef np.float64_t[:,:,:,:] mk
+ cdef np.float64_t[:,:,:,:] qk
+ cdef np.float64_t[:,:,:,:] i
def initialize(self):
# we do this in a single pass, but need two scalar
# per cell, M_k, and Q_k and also the number of particles
# deposited into each one
# the M_k term
- self.omk= np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray omkarr= self.omk
- self.mk= <np.float64_t*> omkarr.data
- # the Q_k term
- self.oqk= np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray oqkarr= self.oqk
- self.qk= <np.float64_t*> oqkarr.data
- # particle count
- self.oi = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray oiarr = self.oi
- self.i = <np.float64_t*> oiarr.data
+ self.mk = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.qk = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.i = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
@@ -312,43 +314,47 @@
cdef float k, mk, qk
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- cell_index = gind(ii[0], ii[1], ii[2], dim) + offset
- k = self.i[cell_index]
- mk = self.mk[cell_index]
- qk = self.qk[cell_index]
+ k = self.i[ii[2], ii[1], ii[0], offset]
+ mk = self.mk[ii[2], ii[1], ii[0], offset]
+ qk = self.qk[ii[2], ii[1], ii[0], offset]
#print k, mk, qk, cell_index
if k == 0.0:
# Initialize cell values
- self.mk[cell_index] = fields[0]
+ self.mk[ii[2], ii[1], ii[0], offset] = fields[0]
else:
- self.mk[cell_index] = mk + (fields[0] - mk) / k
- self.qk[cell_index] = qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
- self.i[cell_index] += 1
+ self.mk[ii[2], ii[1], ii[0], offset] = mk + (fields[0] - mk) / k
+ self.qk[ii[2], ii[1], ii[0], offset] = \
+ qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
+ self.i[ii[2], ii[1], ii[0], offset] += 1
def finalize(self):
# This is the standard variance
# if we want sample variance divide by (self.oi - 1.0)
- std2 = self.oqk / self.oi
- std2[self.oi == 0.0] = 0.0
+ i = np.asarray(self.i)
+ std2 = np.asarray(self.qk) / i
+ std2[i == 0.0] = 0.0
+ std2.shape = self.nvals
return np.sqrt(std2)
deposit_std = StdParticleField
cdef class CICDeposit(ParticleDepositOperation):
- cdef np.float64_t *field
+ cdef np.float64_t[:,:,:,:] field
cdef public object ofield
def initialize(self):
- self.ofield = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.ofield
- self.field = <np.float64_t *> arr.data
+ self.field = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset, # offset into IO field
np.float64_t ppos[3], # this particle's position
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
@@ -373,48 +379,52 @@
for i in range(2):
for j in range(2):
for k in range(2):
- ii = gind(ind[0] - i, ind[1] - j, ind[2] - k, dim) + offset
- self.field[ii] += fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
+ self.field[ind[2] - k, ind[1] - j, ind[0] - i, offset] += \
+ fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
def finalize(self):
- return self.ofield
+ rv = np.asarray(self.field)
+ rv.shape = self.nvals
+ return rv
deposit_cic = CICDeposit
cdef class WeightedMeanParticleField(ParticleDepositOperation):
# Deposit both mass * field and mass into two scalars
# then in finalize divide mass * field / mass
- cdef np.float64_t *wf
- cdef public object owf
- cdef np.float64_t *w
- cdef public object ow
+ cdef np.float64_t[:,:,:,:] wf
+ cdef np.float64_t[:,:,:,:] w
def initialize(self):
- self.owf = np.zeros(self.nvals, dtype='float64', order='F')
- cdef np.ndarray wfarr = self.owf
- self.wf = <np.float64_t*> wfarr.data
-
- self.ow = np.zeros(self.nvals, dtype='float64', order='F')
- cdef np.ndarray warr = self.ow
- self.w = <np.float64_t*> warr.data
+ self.wf = append_axes(
+ np.zeros(self.nvals, dtype='float64', order='F'), 4)
+ self.w = append_axes(
+ np.zeros(self.nvals, dtype='float64', order='F'), 4)
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
cdef int ii[3]
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i]) / dds[i])
- self.w[ gind(ii[0], ii[1], ii[2], dim) + offset] += fields[1]
- self.wf[gind(ii[0], ii[1], ii[2], dim) + offset] += fields[0] * fields[1]
+ self.w[ii[2], ii[1], ii[0], offset] += fields[1]
+ self.wf[ii[2], ii[1], ii[0], offset] += fields[0] * fields[1]
def finalize(self):
- return self.owf / self.ow
+ wf = np.asarray(self.owf)
+ w = np.asarray(self.ow)
+ rv = wf / w
+ rv.shape = self.nvals
+ return rv
deposit_weighted_mean = WeightedMeanParticleField
@@ -426,12 +436,15 @@
self.update_values = 1
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
fields[0] = domain_ind
@@ -442,34 +455,31 @@
deposit_mesh_id = MeshIdentifier
cdef class NNParticleField(ParticleDepositOperation):
- cdef np.float64_t *nnfield
- cdef np.float64_t *distfield
- cdef public object onnfield
- cdef public object odistfield
+ cdef np.float64_t[:,:,:,:] nnfield
+ cdef np.float64_t[:,:,:,:] distfield
def initialize(self):
- self.onnfield = np.zeros(self.nvals, dtype="float64", order='F')
- cdef np.ndarray arr = self.onnfield
- self.nnfield = <np.float64_t*> arr.data
-
- self.odistfield = np.zeros(self.nvals, dtype="float64", order='F')
- self.odistfield[:] = np.inf
- arr = self.odistfield
- self.distfield = <np.float64_t*> arr.data
+ self.nnfield = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.distfield = append_axes(
+ np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.distfield[:] = np.inf
@cython.cdivision(True)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.initializedcheck(False)
cdef void process(self, int dim[3],
np.float64_t left_edge[3],
np.float64_t dds[3],
np.int64_t offset,
np.float64_t ppos[3],
- np.float64_t *fields,
+ np.float64_t[:] fields,
np.int64_t domain_ind
):
# This one is a bit slow. Every grid cell is going to be iterated
# over, and we're going to deposit particles in it.
cdef int i, j, k
cdef int ii[3]
- cdef np.int64_t ggind
cdef np.float64_t r2
cdef np.float64_t gpos[3]
gpos[0] = left_edge[0] + 0.5 * dds[0]
@@ -478,19 +488,20 @@
for j in range(dim[1]):
gpos[2] = left_edge[2] + 0.5 * dds[2]
for k in range(dim[2]):
- ggind = gind(i, j, k, dim) + offset
r2 = ((ppos[0] - gpos[0])*(ppos[0] - gpos[0]) +
(ppos[1] - gpos[1])*(ppos[1] - gpos[1]) +
(ppos[2] - gpos[2])*(ppos[2] - gpos[2]))
- if r2 < self.distfield[ggind]:
- self.distfield[ggind] = r2
- self.nnfield[ggind] = fields[0]
+ if r2 < self.distfield[k,j,i,offset]:
+ self.distfield[k,j,i,offset] = r2
+ self.nnfield[k,j,i,offset] = fields[0]
gpos[2] += dds[2]
gpos[1] += dds[1]
gpos[0] += dds[0]
return
def finalize(self):
- return self.onnfield
+ nn = np.asarray(self.nnfield)
+ nn.shape = self.nvals
+ return nn
deposit_nearest = NNParticleField
diff -r 45dc609be9e893c48cdfc2e568eb26ca0a0f8cbe -r 187f84cc044157da920460ec8423f13d74659652 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -22,6 +22,13 @@
from fp_utils cimport fmin, fmax, i64min, i64max
from yt.geometry.selection_routines cimport _ensure_code
+from libc.stdlib cimport malloc, free
+from libc.string cimport strcmp
+
+from cython.view cimport memoryview
+from cpython cimport buffer
+
+
cdef extern from "stdlib.h":
# NOTE that size_t might not be int
void *alloca(int)
@@ -966,3 +973,128 @@
dest[i,j,k] += dsp * (overlap[0]*overlap[1]*overlap[2])
else:
dest[i,j,k] = dsp
+
+# The OnceIndirect code is from:
+# http://stackoverflow.com/questions/10465091/assembling-a-cython-memoryview-from-numpy-arrays/12991519#12991519
+# This is under the CC-BY-SA license.
+
+cdef class OnceIndirect:
+ cdef object _objects
+ cdef void** buf
+ cdef int ndim
+ cdef int n_rows
+ cdef int buf_len
+ cdef Py_ssize_t* shape
+ cdef Py_ssize_t* strides
+ cdef Py_ssize_t* suboffsets
+ cdef Py_ssize_t itemsize
+ cdef bytes format
+ cdef int is_readonly
+
+ def __cinit__(self, object rows, want_writable=True, want_format=True, allow_indirect=False):
+ """
+ Set want_writable to False if you don't want writable data. (This may
+ prevent copies.)
+ Set want_format to False if your input doesn't support PyBUF_FORMAT (unlikely)
+ Set allow_indirect to True if you are ok with the memoryview being indirect
+ in dimensions other than the first. (This may prevent copies.)
+
+ An example usage:
+
+ cdef double[::cython.view.indirect, ::1] vectors =
+ OnceIndirect([object.vector for object in objects])
+ """
+ demand = buffer.PyBUF_INDIRECT if allow_indirect else buffer.PyBUF_STRIDES
+ if want_writable:
+ demand |= buffer.PyBUF_WRITABLE
+ if want_format:
+ demand |= buffer.PyBUF_FORMAT
+ self._objects = [memoryview(row, demand) for row in rows]
+ self.n_rows = len(self._objects)
+ self.buf_len = sizeof(void*) * self.n_rows
+ self.buf = <void**>malloc(self.buf_len)
+ self.ndim = 1 + self._objects[0].ndim
+ self.shape = <Py_ssize_t*>malloc(sizeof(Py_ssize_t) * self.ndim)
+ self.strides = <Py_ssize_t*>malloc(sizeof(Py_ssize_t) * self.ndim)
+ self.suboffsets = <Py_ssize_t*>malloc(sizeof(Py_ssize_t) * self.ndim)
+
+ cdef memoryview example_obj = self._objects[0]
+ self.itemsize = example_obj.itemsize
+
+ if want_format:
+ self.format = example_obj.view.format
+ else:
+ self.format = None
+ self.is_readonly |= example_obj.view.readonly
+
+ for dim in range(self.ndim):
+ if dim == 0:
+ self.shape[dim] = self.n_rows
+ self.strides[dim] = sizeof(void*)
+ self.suboffsets[dim] = 0
+ else:
+ self.shape[dim] = example_obj.view.shape[dim - 1]
+ self.strides[dim] = example_obj.view.strides[dim - 1]
+ if example_obj.view.suboffsets == NULL:
+ self.suboffsets[dim] = -1
+ else:
+ self.suboffsets[dim] = example_obj.suboffsets[dim - 1]
+
+ cdef memoryview obj
+ cdef int i = 0
+ for obj in self._objects:
+ assert_similar(example_obj, obj)
+ self.buf[i] = obj.view.buf
+ i += 1
+
+ def __getbuffer__(self, Py_buffer* buff, int flags):
+ if (flags & buffer.PyBUF_INDIRECT) != buffer.PyBUF_INDIRECT:
+ raise Exception("don't want to copy data")
+ if flags & buffer.PyBUF_WRITABLE and self.is_readonly:
+ raise Exception("couldn't provide writable, you should have demanded it earlier")
+ if flags & buffer.PyBUF_FORMAT:
+ if self.format is None:
+ raise Exception("couldn't provide format, you should have demanded it earlier")
+ buff.format = self.format
+ else:
+ buff.format = NULL
+
+ buff.buf = <void*>self.buf
+ buff.obj = self
+ buff.len = self.buf_len
+ buff.readonly = self.is_readonly
+ buff.ndim = self.ndim
+ buff.shape = self.shape
+ buff.strides = self.strides
+ buff.suboffsets = self.suboffsets
+ buff.itemsize = self.itemsize
+ buff.internal = NULL
+
+ def __dealloc__(self):
+ free(self.buf)
+ free(self.shape)
+ free(self.strides)
+ free(self.suboffsets)
+
+cdef int assert_similar(memoryview left_, memoryview right_) except -1:
+ cdef Py_buffer left = left_.view
+ cdef Py_buffer right = right_.view
+ assert left.ndim == right.ndim
+ cdef int i
+ for i in range(left.ndim):
+ assert left.shape[i] == right.shape[i], (left_.shape, right_.shape)
+ assert left.strides[i] == right.strides[i], (left_.strides, right_.strides)
+
+ if left.suboffsets == NULL:
+ assert right.suboffsets == NULL, (left_.suboffsets, right_.suboffsets)
+ else:
+ for i in range(left.ndim):
+ assert left.suboffsets[i] == right.suboffsets[i], (left_.suboffsets, right_.suboffsets)
+
+ if left.format == NULL:
+ assert right.format == NULL, (bytes(left.format), bytes(right.format))
+ else:
+ #alternatively, compare as Python strings:
+ #assert bytes(left.format) == bytes(right.format)
+ assert strcmp(left.format, right.format) == 0, (bytes(left.format), bytes(right.format))
+ return 0
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list