[yt-svn] commit/yt: atmyers: Merged in MatthewTurk/yt (pull request #1919)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Feb 3 09:09:14 PST 2016
1 new commit in yt:
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