[yt-svn] commit/yt: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Feb 3 10:22:16 PST 2016
2 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/ca2c63baa630/
Changeset: ca2c63baa630
Branch: yt
User: ngoldbaum
Date: 2016-02-03 18:19:33+00:00
Summary: Backed out changeset be1961831c22
Affected #: 1 file
diff -r be1961831c2209fce82d3770b7748634bffc2581 -r ca2c63baa6304474084019f7d4f728f1a712484e 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/41a14ee1c248/
Changeset: 41a14ee1c248
Branch: yt
User: ngoldbaum
Date: 2016-02-03 18:21:33+00:00
Summary: Back out PR 1919
Affected #: 6 files
diff -r ca2c63baa6304474084019f7d4f728f1a712484e -r 41a14ee1c248f09c88debe90b57c43fc769a407f 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, kernel_name)
+ op = cls(self.ActiveDimensions.prod(), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
- return vals.copy(order="C")
+ return vals.reshape(self.ActiveDimensions, order="C")
def write_to_gdf(self, gdf_path, fields, nprocs=1, field_units=None,
**kwargs):
diff -r ca2c63baa6304474084019f7d4f728f1a712484e -r 41a14ee1c248f09c88debe90b57c43fc769a407f 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(tuple(self.ActiveDimensions), kernel_name)
+ op = cls(self.ActiveDimensions.prod(), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
if vals is None: return
- return vals.copy(order="C")
+ return vals.reshape(self.ActiveDimensions, order="C")
def select_blocks(self, selector):
mask = self._get_selector_mask(selector)
diff -r ca2c63baa6304474084019f7d4f728f1a712484e -r 41a14ee1c248f09c88debe90b57c43fc769a407f yt/frontends/artio/_artio_caller.pyx
--- a/yt/frontends/artio/_artio_caller.pyx
+++ b/yt/frontends/artio/_artio_caller.pyx
@@ -15,7 +15,6 @@
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)
@@ -1555,9 +1554,6 @@
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
@@ -1572,9 +1568,17 @@
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 ca2c63baa6304474084019f7d4f728f1a712484e -r 41a14ee1c248f09c88debe90b57c43fc769a407f 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 ca2c63baa6304474084019f7d4f728f1a712484e -r 41a14ee1c248f09c88debe90b57c43fc769a407f yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -19,21 +19,10 @@
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):
@@ -49,7 +38,6 @@
@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,
@@ -59,10 +47,15 @@
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_pointers
+ cdef np.float64_t *field_vals
cdef np.float64_t pos[3]
- cdef np.float64_t[:] field_vals = np.empty(nf, dtype="float64")
+ 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 << octree.oref)
cdef int nz = dims[0] * dims[1] * dims[2]
@@ -74,7 +67,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
@@ -97,7 +90,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]
+ offset = dom_ind[oct.domain_ind - moff] * nz
if offset < 0: continue
# Check that we found the oct ...
self.process(dims, oi.left_edge, oi.dds,
@@ -108,7 +101,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- @cython.initializedcheck(False)
def process_grid(self, gobj,
np.ndarray[np.float64_t, ndim=2] positions,
fields = None):
@@ -116,11 +108,16 @@
if fields is None:
fields = []
nf = len(fields)
- 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 **field_pointers
+ cdef np.float64_t *field_vals
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]
@@ -131,7 +128,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)
@@ -141,27 +138,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
+ cdef np.int64_t *count # float, for ease
+ cdef public object ocount
def initialize(self):
# Create a numpy array accessible to python
- self.count = append_axes(
- np.zeros(self.nvals, dtype="int64", order='F'), 4)
+ 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
@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
@@ -169,12 +166,10 @@
cdef int i
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- self.count[ii[2], ii[1], ii[0], offset] += 1
+ self.count[gind(ii[0], ii[1], ii[2], dim) + offset] += 1
def finalize(self):
- arr = np.asarray(self.count)
- arr.shape = self.nvals
- return arr.astype("float64")
+ return self.ocount.astype('f8')
deposit_count = CountParticles
@@ -182,25 +177,26 @@
# 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 np.float64_t[:,:,:,:] temp
+ cdef np.float64_t *data
+ cdef public object odata
+ cdef np.float64_t *temp
+ cdef public object otemp
def initialize(self):
- 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)
+ 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
@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]
@@ -232,14 +228,14 @@
dist = idist[0] + idist[1] + idist[2]
# Calculate distance in multiples of the smoothing length
dist = sqrt(dist) / fields[0]
- self.temp[k,j,i,offset] = self.sph_kernel(dist)
- kernel_sum += self.temp[k,j,i,offset]
+ self.temp[gind(i,j,k,dim) + offset] = self.sph_kernel(dist)
+ kernel_sum += self.temp[gind(i,j,k,dim) + 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[k,j,i,offset] / kernel_sum
- self.data[k,j,i,offset] += fields[1] * dist
+ dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
+ self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
def finalize(self):
return self.odata
@@ -247,34 +243,31 @@
deposit_simple_smooth = SimpleSmooth
cdef class SumParticleField(ParticleDepositOperation):
- cdef np.float64_t[:,:,:,:] sum
+ cdef np.float64_t *sum
+ cdef public object osum
def initialize(self):
- self.sum = append_axes(
- np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.osum = np.zeros(self.nvals, dtype="float64", order='F')
+ cdef np.ndarray arr = self.osum
+ self.sum = <np.float64_t*> arr.data
@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[ii[2], ii[1], ii[0], offset] += fields[0]
+ self.sum[gind(ii[0], ii[1], ii[2], dim) + offset] += fields[0]
return
def finalize(self):
- sum = np.asarray(self.sum)
- sum.shape = self.nvals
- return sum
+ return self.osum
deposit_sum = SumParticleField
@@ -282,31 +275,36 @@
# 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 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
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.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)
+ 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
@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]
@@ -314,47 +312,43 @@
cdef float k, mk, qk
for i in range(3):
ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
- 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]
+ 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]
#print k, mk, qk, cell_index
if k == 0.0:
# Initialize cell values
- self.mk[ii[2], ii[1], ii[0], offset] = fields[0]
+ self.mk[cell_index] = fields[0]
else:
- 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
+ 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
def finalize(self):
# This is the standard variance
# if we want sample variance divide by (self.oi - 1.0)
- i = np.asarray(self.i)
- std2 = np.asarray(self.qk) / i
- std2[i == 0.0] = 0.0
- std2.shape = self.nvals
+ std2 = self.oqk / self.oi
+ std2[self.oi == 0.0] = 0.0
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.field = append_axes(
- np.zeros(self.nvals, dtype="float64", order='F'), 4)
+ self.ofield = np.zeros(self.nvals, dtype="float64", order='F')
+ cdef np.ndarray arr = self.ofield
+ self.field = <np.float64_t *> arr.data
@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
):
@@ -379,52 +373,48 @@
for i in range(2):
for j in range(2):
for k in range(2):
- self.field[ind[2] - k, ind[1] - j, ind[0] - i, offset] += \
- fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
+ 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]
def finalize(self):
- rv = np.asarray(self.field)
- rv.shape = self.nvals
- return rv
+ return self.ofield
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 np.float64_t[:,:,:,:] w
+ cdef np.float64_t *wf
+ cdef public object owf
+ cdef np.float64_t *w
+ cdef public object ow
def initialize(self):
- 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)
+ 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
@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[ii[2], ii[1], ii[0], offset] += fields[1]
- self.wf[ii[2], ii[1], ii[0], offset] += fields[0] * fields[1]
+ 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]
def finalize(self):
- wf = np.asarray(self.owf)
- w = np.asarray(self.ow)
- rv = wf / w
- rv.shape = self.nvals
- return rv
+ return self.owf / self.ow
deposit_weighted_mean = WeightedMeanParticleField
@@ -436,15 +426,12 @@
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
@@ -455,31 +442,34 @@
deposit_mesh_id = MeshIdentifier
cdef class NNParticleField(ParticleDepositOperation):
- cdef np.float64_t[:,:,:,:] nnfield
- cdef np.float64_t[:,:,:,:] distfield
+ cdef np.float64_t *nnfield
+ cdef np.float64_t *distfield
+ cdef public object onnfield
+ cdef public object odistfield
def initialize(self):
- 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
+ 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
@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]
@@ -488,20 +478,19 @@
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[k,j,i,offset]:
- self.distfield[k,j,i,offset] = r2
- self.nnfield[k,j,i,offset] = fields[0]
+ if r2 < self.distfield[ggind]:
+ self.distfield[ggind] = r2
+ self.nnfield[ggind] = fields[0]
gpos[2] += dds[2]
gpos[1] += dds[1]
gpos[0] += dds[0]
return
def finalize(self):
- nn = np.asarray(self.nnfield)
- nn.shape = self.nvals
- return nn
+ return self.onnfield
deposit_nearest = NNParticleField
diff -r ca2c63baa6304474084019f7d4f728f1a712484e -r 41a14ee1c248f09c88debe90b57c43fc769a407f yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -22,13 +22,6 @@
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)
@@ -973,128 +966,3 @@
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