[yt-svn] commit/yt: jzuhone: Merged in MatthewTurk/yt (pull request #2329)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Aug 24 11:15:58 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/6c9ed4a3752d/
Changeset: 6c9ed4a3752d
Branch: yt
User: jzuhone
Date: 2016-08-24 18:15:23+00:00
Summary: Merged in MatthewTurk/yt (pull request #2329)
Split out distance queue into its own Cython module
Affected #: 5 files
diff -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 -r 6c9ed4a3752d649d872cf47adffd419a3c8bf903 setup.py
--- a/setup.py
+++ b/setup.py
@@ -196,7 +196,7 @@
"particle_mesh_operations", "depth_first_octree", "fortran_reader",
"interpolators", "misc_utilities", "basic_octree", "image_utilities",
"points_in_volume", "quad_tree", "ray_integrators", "mesh_utilities",
- "amr_kdtools", "lenses",
+ "amr_kdtools", "lenses", "distance_queue"
]
for ext_name in lib_exts:
cython_extensions.append(
diff -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 -r 6c9ed4a3752d649d872cf47adffd419a3c8bf903 yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -23,15 +23,12 @@
from yt.utilities.lib.fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
from .particle_deposit cimport kernel_func, get_kernel_func, gind
+from yt.utilities.lib.distance_queue cimport NeighborList, Neighbor_compare, \
+ r2dist, DistanceQueue
cdef extern from "platform_dep.h":
void *alloca(int)
-cdef struct NeighborList
-cdef struct NeighborList:
- np.int64_t pn # Particle number
- np.float64_t r2 # radius**2
-
cdef class ParticleSmoothOperation:
# We assume each will allocate and define their own temporary storage
cdef kernel_func sph_kernel
@@ -39,10 +36,8 @@
cdef np.float64_t DW[3]
cdef int nfields
cdef int maxn
- cdef int curn
cdef bint periodicity[3]
# Note that we are preallocating here, so this is *not* threadsafe.
- cdef NeighborList *neighbors
cdef void (*pos_setup)(np.float64_t ipos[3], np.float64_t opos[3])
cdef void neighbor_process(self, int dim[3], np.float64_t left_edge[3],
np.float64_t dds[3], np.float64_t[:,:] ppos,
@@ -52,7 +47,7 @@
np.int64_t offset, np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
int *nsize, np.float64_t[:,:] oct_left_edges,
- np.float64_t[:,:] oct_dds)
+ np.float64_t[:,:] oct_dds, DistanceQueue dq)
cdef int neighbor_search(self, np.float64_t pos[3], OctreeContainer octree,
np.int64_t **nind, int *nsize,
np.int64_t nneighbors, np.int64_t domain_id,
@@ -65,10 +60,7 @@
np.int64_t offset,
np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
- int *nsize)
- cdef void neighbor_eval(self, np.int64_t pn, np.float64_t ppos[3],
- np.float64_t cpos[3])
- cdef void neighbor_reset(self)
+ int *nsize, DistanceQueue dq)
cdef void neighbor_find(self,
np.int64_t nneighbors,
np.int64_t *nind,
@@ -78,7 +70,7 @@
np.float64_t[:,:] ppos,
np.float64_t cpos[3],
np.float64_t[:,:] oct_left_edges,
- np.float64_t[:,:] oct_dds)
+ np.float64_t[:,:] oct_dds, DistanceQueue dq)
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **index_fields)
+ np.float64_t **index_fields, DistanceQueue dq)
diff -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 -r 6c9ed4a3752d649d872cf47adffd419a3c8bf903 yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -24,44 +24,6 @@
from oct_container cimport Oct, OctAllocationContainer, \
OctreeContainer, OctInfo
-cdef int Neighbor_compare(void *on1, void *on2) nogil:
- cdef NeighborList *n1
- cdef NeighborList *n2
- n1 = <NeighborList *> on1
- n2 = <NeighborList *> on2
- # Note that we set this up so that "greatest" evaluates to the *end* of the
- # list, so we can do standard radius comparisons.
- if n1.r2 < n2.r2:
- return -1
- elif n1.r2 == n2.r2:
- return 0
- else:
- return 1
-
- at cython.cdivision(True)
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.initializedcheck(False)
-cdef np.float64_t r2dist(np.float64_t ppos[3],
- np.float64_t cpos[3],
- np.float64_t DW[3],
- bint periodicity[3],
- np.float64_t max_dist2):
- cdef int i
- cdef np.float64_t r2, DR
- r2 = 0.0
- for i in range(3):
- DR = (ppos[i] - cpos[i])
- if not periodicity[i]:
- pass
- elif (DR > DW[i]/2.0):
- DR -= DW[i]
- elif (DR < -DW[i]/2.0):
- DR += DW[i]
- r2 += DR * DR
- if max_dist2 >= 0.0 and r2 > max_dist2:
- return -1.0
- return r2
cdef void spherical_coord_setup(np.float64_t ipos[3], np.float64_t opos[3]):
opos[0] = ipos[0] * sin(ipos[1]) * cos(ipos[2])
@@ -80,10 +42,6 @@
self.nvals = nvals
self.nfields = nfields
self.maxn = max_neighbors
-
- self.neighbors = <NeighborList *> malloc(
- sizeof(NeighborList) * self.maxn)
- self.neighbor_reset()
self.sph_kernel = get_kernel_func(kernel_name)
def initialize(self, *args):
@@ -247,6 +205,9 @@
cdef np.ndarray[np.uint8_t, ndim=1] visited
visited = np.zeros(mdom_ind.shape[0], dtype="uint8")
cdef int nproc = 0
+ # This should be thread-private if we ever go to OpenMP
+ cdef DistanceQueue dist_queue = DistanceQueue(self.maxn)
+ dist_queue._setup(self.DW, self.periodicity)
for i in range(oct_positions.shape[0]):
for j in range(3):
pos[j] = oct_positions[i, j]
@@ -260,7 +221,7 @@
dims, moi.left_edge, moi.dds, cart_positions, field_pointers, doff,
&nind, pind, pcount, offset, index_field_pointers,
particle_octree, domain_id, &nsize, oct_left_edges,
- oct_dds)
+ oct_dds, dist_queue)
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
if nind != NULL:
@@ -369,6 +330,9 @@
# refers to that oct's particles.
cdef int maxnei = 0
cdef int nproc = 0
+ # This should be thread-private if we ever go to OpenMP
+ cdef DistanceQueue dist_queue = DistanceQueue(self.maxn)
+ dist_queue._setup(self.DW, self.periodicity)
for i in range(doff.shape[0]):
if doff[i] < 0: continue
offset = pind[doff[i]]
@@ -380,7 +344,8 @@
pos[k] = positions[pind0, k]
self.neighbor_process_particle(pos, cart_positions, field_pointers,
doff, &nind, pind, pcount, pind0,
- NULL, particle_octree, domain_id, &nsize)
+ NULL, particle_octree, domain_id, &nsize,
+ dist_queue)
#print "VISITED", visited.sum(), visited.size,
#print 100.0*float(visited.sum())/visited.size
if nind != NULL:
@@ -463,55 +428,9 @@
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **ifields):
+ np.float64_t **ifields, DistanceQueue dq):
raise NotImplementedError
- cdef void neighbor_reset(self):
- self.curn = 0
- for i in range(self.maxn):
- self.neighbors[i].pn = -1
- self.neighbors[i].r2 = 1e300
-
- cdef void neighbor_eval(self, np.int64_t pn, np.float64_t ppos[3],
- np.float64_t cpos[3]):
- # Here's a python+numpy simulator of this:
- # http://paste.yt-project.org/show/5445/
- cdef int i, di
- cdef np.float64_t r2, r2_trunc
- if self.curn == self.maxn:
- # Truncate calculation if it's bigger than this in any dimension
- r2_trunc = self.neighbors[self.curn - 1].r2
- else:
- # Don't truncate our calculation
- r2_trunc = -1
- r2 = r2dist(ppos, cpos, self.DW, self.periodicity, r2_trunc)
- if r2 == -1:
- return
- if self.curn == 0:
- self.neighbors[0].r2 = r2
- self.neighbors[0].pn = pn
- self.curn += 1
- return
- # Now insert in a sorted way
- di = -1
- for i in range(self.curn - 1, -1, -1):
- # We are checking if i is less than us, to see if we should insert
- # to the right (i.e., i+1).
- if self.neighbors[i].r2 < r2:
- di = i
- break
- # The outermost one is already too small.
- if di == self.maxn - 1:
- return
- if (self.maxn - (di + 2)) > 0:
- memmove(<void *> (self.neighbors + di + 2),
- <void *> (self.neighbors + di + 1),
- sizeof(NeighborList) * (self.maxn - (di + 2)))
- self.neighbors[di + 1].r2 = r2
- self.neighbors[di + 1].pn = pn
- if self.curn < self.maxn:
- self.curn += 1
-
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -526,6 +445,7 @@
np.float64_t cpos[3],
np.float64_t[:,:] oct_left_edges,
np.float64_t[:,:] oct_dds,
+ DistanceQueue dq
):
# We are now given the number of neighbors, the indices into the
# domains for them, and the number of particles for each.
@@ -535,13 +455,13 @@
cdef np.float64_t ex[2]
cdef np.float64_t DR[2]
cdef np.float64_t cp, r2_trunc, r2, dist
- self.neighbor_reset()
+ dq.neighbor_reset()
for ni in range(nneighbors):
if nind[ni] == -1: continue
# terminate early if all 8 corners of oct are farther away than
# most distant currently known neighbor
- if oct_left_edges != None and self.curn == self.maxn:
- r2_trunc = self.neighbors[self.curn - 1].r2
+ if oct_left_edges != None and dq.curn == dq.maxn:
+ r2_trunc = dq.neighbors[dq.curn - 1].r2
# iterate over each dimension in the outer loop so we can
# consolidate temporary storage
# What this next bit does is figure out which component is the
@@ -577,7 +497,7 @@
pn = pinds[offset + i]
for j in range(3):
pos[j] = ppos[pn, j]
- self.neighbor_eval(pn, pos, cpos)
+ dq.neighbor_eval(pn, pos, cpos)
@cython.cdivision(True)
@cython.boundscheck(False)
@@ -592,7 +512,8 @@
np.float64_t **index_fields,
OctreeContainer octree, np.int64_t domain_id,
int *nsize, np.float64_t[:,:] oct_left_edges,
- np.float64_t[:,:] oct_dds):
+ np.float64_t[:,:] oct_dds,
+ DistanceQueue dq):
# Note that we assume that fields[0] == smoothing length in the native
# units supplied. We can now iterate over every cell in the block and
# every particle to find the nearest. We will use a priority heap.
@@ -610,17 +531,18 @@
nneighbors = self.neighbor_search(opos, octree,
nind, nsize, nneighbors, domain_id, &oct, 0)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts,
- pinds, ppos, opos, oct_left_edges, oct_dds)
+ pinds, ppos, opos, oct_left_edges,
+ oct_dds, dq)
# Now we have all our neighbors in our neighbor list.
- if self.curn <-1*self.maxn:
+ if dq.curn <-1*dq.maxn:
ntot = nntot = 0
for m in range(nneighbors):
if nind[0][m] < 0: continue
nntot += 1
ntot += pcounts[nind[0][m]]
- print "SOMETHING WRONG", self.curn, nneighbors, ntot, nntot
+ print "SOMETHING WRONG", dq.curn, nneighbors, ntot, nntot
self.process(offset, i, j, k, dim, opos, fields,
- index_fields)
+ index_fields, dq)
cpos[2] += dds[2]
cpos[1] += dds[1]
cpos[0] += dds[0]
@@ -637,7 +559,8 @@
np.int64_t offset,
np.float64_t **index_fields,
OctreeContainer octree,
- np.int64_t domain_id, int *nsize):
+ np.int64_t domain_id, int *nsize,
+ DistanceQueue dq):
# Note that we assume that fields[0] == smoothing length in the native
# units supplied. We can now iterate over every cell in the block and
# every particle to find the nearest. We will use a priority heap.
@@ -652,8 +575,8 @@
nneighbors = self.neighbor_search(opos, octree,
nind, nsize, nneighbors, domain_id, &oct, 0)
self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos,
- opos, None, None)
- self.process(offset, i, j, k, dim, opos, fields, index_fields)
+ opos, None, None, dq)
+ self.process(offset, i, j, k, dim, opos, fields, index_fields, dq)
cdef class VolumeWeightedSmooth(ParticleSmoothOperation):
# This smoothing function evaluates the field, *without* normalization, at
@@ -692,7 +615,7 @@
@cython.initializedcheck(False)
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **index_fields):
+ np.float64_t **index_fields, DistanceQueue dq):
# We have our i, j, k for our cell, as well as the cell position.
# We also have a list of neighboring particles with particle numbers.
cdef int n, fi
@@ -702,13 +625,13 @@
cdef np.int64_t pn
# We get back our mass
# rho_i = sum(j = 1 .. n) m_j * W_ij
- max_r = sqrt(self.neighbors[self.curn-1].r2)
+ max_r = sqrt(dq.neighbors[dq.curn-1].r2)
max_hsml = index_fields[0][gind(i,j,k,dim) + offset]
- for n in range(self.curn):
+ for n in range(dq.curn):
# No normalization for the moment.
# fields[0] is the smoothing length.
- r2 = self.neighbors[n].r2
- pn = self.neighbors[n].pn
+ r2 = dq.neighbors[n].r2
+ pn = dq.neighbors[n].pn
# Smoothing kernel weight function
mass = fields[0][pn]
hsml = fields[1][pn]
@@ -751,15 +674,15 @@
@cython.initializedcheck(False)
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **index_fields):
+ np.float64_t **index_fields, DistanceQueue dq):
# We have our i, j, k for our cell, as well as the cell position.
# We also have a list of neighboring particles with particle numbers.
cdef np.int64_t pn
# We get back our mass
# rho_i = sum(j = 1 .. n) m_j * W_ij
- pn = self.neighbors[0].pn
+ pn = dq.neighbors[0].pn
self.fp[gind(i,j,k,dim) + offset] = fields[0][pn]
- #self.fp[gind(i,j,k,dim) + offset] = self.neighbors[0].r2
+ #self.fp[gind(i,j,k,dim) + offset] = dq.neighbors[0].r2
return
nearest_smooth = NearestNeighborSmooth
@@ -785,18 +708,18 @@
@cython.initializedcheck(False)
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **index_fields):
+ np.float64_t **index_fields, DistanceQueue dq):
# We have our i, j, k for our cell, as well as the cell position.
# We also have a list of neighboring particles with particle numbers.
cdef np.int64_t pn, ni, di
cdef np.float64_t total_weight = 0.0, total_value = 0.0, r2, val, w
# We're going to do a very simple IDW average
- if self.neighbors[0].r2 == 0.0:
- pn = self.neighbors[0].pn
+ if dq.neighbors[0].r2 == 0.0:
+ pn = dq.neighbors[0].pn
self.fp[gind(i,j,k,dim) + offset] = fields[0][pn]
- for ni in range(self.curn):
- r2 = self.neighbors[ni].r2
- val = fields[0][self.neighbors[ni].pn]
+ for ni in range(dq.curn):
+ r2 = dq.neighbors[ni].r2
+ val = fields[0][dq.neighbors[ni].pn]
w = r2
for di in range(self.p2 - 1):
w *= r2
@@ -821,10 +744,10 @@
@cython.initializedcheck(False)
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **index_fields):
+ np.float64_t **index_fields, DistanceQueue dq):
cdef np.float64_t max_r
# We assume "offset" here is the particle index.
- max_r = sqrt(self.neighbors[self.curn-1].r2)
+ max_r = sqrt(dq.neighbors[dq.curn-1].r2)
fields[0][offset] = max_r
nth_neighbor_smooth = NthNeighborDistanceSmooth
@@ -842,16 +765,16 @@
@cython.initializedcheck(False)
cdef void process(self, np.int64_t offset, int i, int j, int k,
int dim[3], np.float64_t cpos[3], np.float64_t **fields,
- np.float64_t **index_fields):
+ np.float64_t **index_fields, DistanceQueue dq):
cdef np.float64_t r2, hsml, dens, mass, weight, lw
cdef int pn
# We assume "offset" here is the particle index.
- hsml = sqrt(self.neighbors[self.curn-1].r2)
+ hsml = sqrt(dq.neighbors[dq.curn-1].r2)
dens = 0.0
weight = 0.0
- for pn in range(self.curn):
- mass = fields[0][self.neighbors[pn].pn]
- r2 = self.neighbors[pn].r2
+ for pn in range(dq.curn):
+ mass = fields[0][dq.neighbors[pn].pn]
+ r2 = dq.neighbors[pn].r2
lw = self.sph_kernel(sqrt(r2) / hsml)
dens += mass * lw
weight = (4.0/3.0) * 3.1415926 * hsml**3
diff -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 -r 6c9ed4a3752d649d872cf47adffd419a3c8bf903 yt/utilities/lib/distance_queue.pxd
--- /dev/null
+++ b/yt/utilities/lib/distance_queue.pxd
@@ -0,0 +1,43 @@
+"""
+A queue for evaluating distances to discrete points
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+cimport cython
+cimport numpy as np
+import numpy as np
+from libc.stdlib cimport malloc, free
+from libc.string cimport memmove
+
+cdef struct NeighborList:
+ np.int64_t pn # Particle number
+ np.float64_t r2 # radius**2
+
+cdef int Neighbor_compare(void *on1, void *on2) nogil
+cdef np.float64_t r2dist(np.float64_t ppos[3],
+ np.float64_t cpos[3],
+ np.float64_t DW[3],
+ bint periodicity[3],
+ np.float64_t max_dist2)
+
+cdef class DistanceQueue:
+ cdef int maxn
+ cdef int curn
+ cdef np.float64_t DW[3]
+ cdef bint periodicity[3]
+ cdef NeighborList* neighbors # flat array
+ cdef void _setup(self, np.float64_t DW[3], bint periodicity[3])
+ cdef void neighbor_eval(self, np.int64_t pn, np.float64_t ppos[3],
+ np.float64_t cpos[3])
+ cdef void neighbor_reset(self)
diff -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 -r 6c9ed4a3752d649d872cf47adffd419a3c8bf903 yt/utilities/lib/distance_queue.pyx
--- /dev/null
+++ b/yt/utilities/lib/distance_queue.pyx
@@ -0,0 +1,149 @@
+"""
+Distance queue implementation
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+cimport numpy as np
+import numpy as np
+cimport cython
+
+cdef int Neighbor_compare(void *on1, void *on2) nogil:
+ cdef NeighborList *n1
+ cdef NeighborList *n2
+ n1 = <NeighborList *> on1
+ n2 = <NeighborList *> on2
+ # Note that we set this up so that "greatest" evaluates to the *end* of the
+ # list, so we can do standard radius comparisons.
+ if n1.r2 < n2.r2:
+ return -1
+ elif n1.r2 == n2.r2:
+ return 0
+ else:
+ return 1
+
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.initializedcheck(False)
+cdef np.float64_t r2dist(np.float64_t ppos[3],
+ np.float64_t cpos[3],
+ np.float64_t DW[3],
+ bint periodicity[3],
+ np.float64_t max_dist2):
+ cdef int i
+ cdef np.float64_t r2, DR
+ r2 = 0.0
+ for i in range(3):
+ DR = (ppos[i] - cpos[i])
+ if not periodicity[i]:
+ pass
+ elif (DR > DW[i]/2.0):
+ DR -= DW[i]
+ elif (DR < -DW[i]/2.0):
+ DR += DW[i]
+ r2 += DR * DR
+ if max_dist2 >= 0.0 and r2 > max_dist2:
+ return -1.0
+ return r2
+
+cdef class DistanceQueue:
+ """This is a distance queue object, designed to incrementally evaluate N
+ positions against a reference point. It is initialized with the maximum
+ number that are to be retained (i.e., maxn-nearest neighbors)."""
+ def __cinit__(self, int maxn):
+ cdef int i
+ self.maxn = maxn
+ self.curn = 0
+ self.neighbors = <NeighborList *> malloc(
+ sizeof(NeighborList) * self.maxn)
+ self.neighbor_reset()
+ for i in range(3):
+ self.DW[i] = 0
+ self.periodicity[i] = False
+
+ cdef void _setup(self, np.float64_t DW[3], bint periodicity[3]):
+ for i in range(3):
+ self.DW[i] = DW[i]
+ self.periodicity[i] = periodicity[i]
+
+ def setup(self, np.float64_t[:] DW, periodicity):
+ for i in range(3):
+ self.DW[i] = DW[i]
+ self.periodicity[i] = periodicity[i]
+
+ def __dealloc__(self):
+ free(self.neighbors)
+
+ cdef void neighbor_eval(self, np.int64_t pn, np.float64_t ppos[3],
+ np.float64_t cpos[3]):
+ # Here's a python+numpy simulator of this:
+ # http://paste.yt-project.org/show/5445/
+ cdef int i, di
+ cdef np.float64_t r2, r2_trunc
+ if self.curn == self.maxn:
+ # Truncate calculation if it's bigger than this in any dimension
+ r2_trunc = self.neighbors[self.curn - 1].r2
+ else:
+ # Don't truncate our calculation
+ r2_trunc = -1
+ r2 = r2dist(ppos, cpos, self.DW, self.periodicity, r2_trunc)
+ if r2 == -1:
+ return
+ if self.curn == 0:
+ self.neighbors[0].r2 = r2
+ self.neighbors[0].pn = pn
+ self.curn += 1
+ return
+ # Now insert in a sorted way
+ di = -1
+ for i in range(self.curn - 1, -1, -1):
+ # We are checking if i is less than us, to see if we should insert
+ # to the right (i.e., i+1).
+ if self.neighbors[i].r2 < r2:
+ di = i
+ break
+ # The outermost one is already too small.
+ if di == self.maxn - 1:
+ return
+ if (self.maxn - (di + 2)) > 0:
+ memmove(<void *> (self.neighbors + di + 2),
+ <void *> (self.neighbors + di + 1),
+ sizeof(NeighborList) * (self.maxn - (di + 2)))
+ self.neighbors[di + 1].r2 = r2
+ self.neighbors[di + 1].pn = pn
+ if self.curn < self.maxn:
+ self.curn += 1
+
+ cdef void neighbor_reset(self):
+ for i in range(self.maxn):
+ self.neighbors[i].r2 = 1e300
+ self.neighbors[i].pn = -1
+ self.curn = 0
+
+ def find_nearest(self, np.float64_t[:] center, np.float64_t[:,:] points):
+ """This function accepts a center and a set of [N,3] points, and it
+ will return the indices into the points array of the maxn closest
+ neighbors."""
+ cdef int i, j
+ cdef np.float64_t ppos[3], cpos[3]
+ self.neighbor_reset()
+ for i in range(3):
+ cpos[i] = center[i]
+ for j in range(points.shape[0]):
+ for i in range(3):
+ ppos[i] = points[j,i]
+ self.neighbor_eval(j, ppos, cpos)
+ rv = np.empty(self.curn, dtype="int64")
+ for i in range(self.curn):
+ rv[i] = self.neighbors[i].pn
+ return rv
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