[Git][debian-gis-team/pykdtree][master] 5 commits: New upstream version 1.4.1+ds
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Fri Jan 31 08:18:47 GMT 2025
Antonio Valentino pushed to branch master at Debian GIS Project / pykdtree
Commits:
beefc088 by Antonio Valentino at 2025-01-31T08:45:32+01:00
New upstream version 1.4.1+ds
- - - - -
feff7f82 by Antonio Valentino at 2025-01-31T08:45:33+01:00
Update upstream source from tag 'upstream/1.4.1+ds'
Update to upstream version '1.4.1+ds'
with Debian dir 20433656f19b6e61178725903a1fff7e8bc627c4
- - - - -
9b71b780 by Antonio Valentino at 2025-01-31T08:46:18+01:00
New upstream release
- - - - -
0b2fad5d by Antonio Valentino at 2025-01-31T08:51:21+01:00
Update dates in d/copyright
- - - - -
945155da by Antonio Valentino at 2025-01-31T09:03:28+01:00
Set distribution to unstable
- - - - -
9 changed files:
- .github/workflows/deploy-wheels.yml
- CHANGELOG.md
- debian/changelog
- debian/copyright
- pykdtree/__init__.py
- pykdtree/_kdtree_core.c.mako
- pykdtree/kdtree.pyx
- pykdtree/test_tree.py
- setup.py
Changes:
=====================================
.github/workflows/deploy-wheels.yml
=====================================
@@ -34,45 +34,43 @@ jobs:
matrix:
include:
- os: windows-2019
- cibw_archs: "AMD64 ARM64"
- artifact_name: "win"
- - os: macos-12
- cibw_archs: "x86_64 arm64"
- artifact_name: "mac"
- - os: "ubuntu-20.04"
+ cibw_archs: "AMD64"
+ - os: windows-2019
+ cibw_archs: "ARM64"
+ - os: macos-13
+ cibw_archs: "x86_64"
+ deploy_target: "13"
+ - os: macos-14
+ cibw_archs: "arm64"
+ deploy_target: "14"
+ - os: ubuntu-24.04-arm
cibw_archs: "aarch64"
- artifact_name: "ubuntu-aarch"
- - os: "ubuntu-20.04"
+ - os: ubuntu-22.04
cibw_archs: "x86_64"
- artifact_name: "ubuntu-x86_64"
steps:
- uses: actions/checkout at v4
- run: |
git fetch --prune --unshallow
- - name: Set up QEMU
- if: runner.os == 'Linux'
- uses: docker/setup-qemu-action at v3
- with:
- platforms: all
-
- name: Build wheels
- uses: pypa/cibuildwheel at v2.20.0
+ uses: pypa/cibuildwheel at v2.22.0
env:
CIBW_SKIP: "cp36-* cp37-* cp38-* pp* *i686 *-musllinux_aarch64"
CIBW_ARCHS: "${{ matrix.cibw_archs }}"
CIBW_TEST_COMMAND: "pytest -v --pyargs pykdtree"
CIBW_TEST_REQUIRES: "pytest"
- CIBW_TEST_SKIP: "*_arm64 *_universal2:arm64"
+ CIBW_TEST_SKIP: "*-win_arm64"
+ CIBW_BUILD_VERBOSITY: 1
# we use openmp (libomp) from homebrew which has a current limit of
- # macos 12 (Monterey): https://formulae.brew.sh/formula/libomp
- CIBW_ENVIRONMENT_MACOS: MACOSX_DEPLOYMENT_TARGET=12
+ # macos 13 (Ventura): https://formulae.brew.sh/formula/libomp
+ CIBW_ENVIRONMENT_MACOS: "MACOSX_DEPLOYMENT_TARGET=${{ matrix.deploy_target }}"
+ CIBW_BEFORE_BUILD_MACOS: "brew install libomp"
- name: Upload wheel(s) as build artifacts
uses: actions/upload-artifact at v4
with:
- name: wheels-${{ matrix.artifact_name }}
+ name: "wheels-${{ matrix.os }}-${{ matrix.cibw_archs }}"
path: ./wheelhouse/*.whl
upload_pypi:
@@ -84,29 +82,15 @@ jobs:
with:
name: sdist
path: dist
- - name: Download wheels artifact - win
- uses: actions/download-artifact at v4
- with:
- name: wheels-win
- path: dist
- - name: Download wheels artifact - mac
- uses: actions/download-artifact at v4
- with:
- name: wheels-mac
- path: dist
- - name: Download wheels artifact - ubuntu aarch
- uses: actions/download-artifact at v4
- with:
- name: wheels-ubuntu-aarch
- path: dist
- - name: Download wheels artifact - ubuntu x86_64
+ - name: Download wheels artifacts
uses: actions/download-artifact at v4
with:
- name: wheels-ubuntu-x86_64
+ pattern: wheels-*
+ merge-multiple: true
path: dist
- name: Publish package to PyPI
if: github.event.action == 'published'
- uses: pypa/gh-action-pypi-publish at v1.10.0
+ uses: pypa/gh-action-pypi-publish at v1.12.3
with:
user: ${{ secrets.pypi_username }}
password: ${{ secrets.pypi_password }}
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,42 @@
+## Version 1.4.1 (2025/01/29)
+
+### Issues Closed
+
+* [Issue 138](https://github.com/storpipfugl/pykdtree/issues/138) - Another integer overflow issue ([PR 139](https://github.com/storpipfugl/pykdtree/pull/139) by [@lkeegan](https://github.com/lkeegan))
+* [Issue 137](https://github.com/storpipfugl/pykdtree/issues/137) - Mako template needs to be re-rendered ([PR 139](https://github.com/storpipfugl/pykdtree/pull/139) by [@lkeegan](https://github.com/lkeegan))
+
+In this release 2 issues were closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 139](https://github.com/storpipfugl/pykdtree/pull/139) - Fix additional integer overflow issues ([138](https://github.com/storpipfugl/pykdtree/issues/138), [137](https://github.com/storpipfugl/pykdtree/issues/137))
+
+In this release 1 pull request was closed.
+
+
+## Version 1.4.0 (2025/01/27)
+
+### Issues Closed
+
+* [Issue 38](https://github.com/storpipfugl/pykdtree/issues/38) - pydktree breaks down for very large number of data points ([PR 135](https://github.com/storpipfugl/pykdtree/pull/135) by [@lkeegan](https://github.com/lkeegan))
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 126](https://github.com/storpipfugl/pykdtree/pull/126) - Fix the wrong syntax of pip install
+
+#### Features added
+
+* [PR 135](https://github.com/storpipfugl/pykdtree/pull/135) - Add automatic switching between 32-bit and 64-bit structures depending on number of points ([38](https://github.com/storpipfugl/pykdtree/issues/38), [38](https://github.com/storpipfugl/pykdtree/issues/38))
+
+In this release 2 pull requests were closed.
+
+
## Version 1.3.13 (2024/09/04)
### Issues Closed
=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+pykdtree (1.4.1+ds-1) unstable; urgency=medium
+
+ * New upstream release.
+ * Update dates in d/copyright.
+
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Fri, 31 Jan 2025 09:03:11 +0100
+
pykdtree (1.3.13+ds-1) unstable; urgency=medium
[ Bas Couwenberg ]
=====================================
debian/copyright
=====================================
@@ -25,7 +25,7 @@ License: LGPL-3+
Public License can be found in "/usr/share/common-licenses/LGPL-3".
Files: debian/*
-Copyright: 2013-2024, Antonio Valentino <antonio.valentino at tiscali.it>
+Copyright: 2013-2025, Antonio Valentino <antonio.valentino at tiscali.it>
License: GPL-3+
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
=====================================
pykdtree/__init__.py
=====================================
@@ -9,5 +9,5 @@ except ImportError as err:
"could not find it on your system. To enable better performance "
"OpenMP must be installed (ex. ``brew install omp`` on Mac with "
"HomeBrew). Otherwise, try installing Pykdtree from source (ex. "
- "``pip install --no-binary pykdtree --force-install pykdtree``)."
+ "``pip install --no-binary pykdtree --force-reinstall pykdtree``)."
) from err
=====================================
pykdtree/_kdtree_core.c.mako
=====================================
@@ -29,65 +29,148 @@ Anne M. Archibald and libANN by David M. Mount and Sunil Arya.
#include <float.h>
#define PA(i,d) (pa[no_dims * pidx[i] + d])
-#define PASWAP(a,b) { uint32_t tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
+% for ITYPE in ['int32_t', 'int64_t']:
+#define PASWAP_${ITYPE}(a,b) { u${ITYPE} tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
+% endfor
+
+#define IDX_MAX_int32_t UINT32_MAX
+#define IDX_MAX_int64_t UINT64_MAX
+#define DIST_MAX_float FLT_MAX
+#define DIST_MAX_double DBL_MAX
#ifdef _MSC_VER
#define restrict __restrict
#endif
% for DTYPE in ['float', 'double']:
+% for ITYPE in ['int32_t', 'int64_t']:
typedef struct
{
${DTYPE} cut_val;
int8_t cut_dim;
- uint32_t start_idx;
- uint32_t n;
+ u${ITYPE} start_idx;
+ u${ITYPE} n;
${DTYPE} cut_bounds_lv;
${DTYPE} cut_bounds_hv;
- struct Node_${DTYPE} *left_child;
- struct Node_${DTYPE} *right_child;
-} Node_${DTYPE};
+ struct Node_${DTYPE}_${ITYPE} *left_child;
+ struct Node_${DTYPE}_${ITYPE} *right_child;
+} Node_${DTYPE}_${ITYPE};
typedef struct
{
${DTYPE} *bbox;
int8_t no_dims;
- uint32_t *pidx;
- struct Node_${DTYPE} *root;
-} Tree_${DTYPE};
+ u${ITYPE} *pidx;
+ struct Node_${DTYPE}_${ITYPE} *root;
+} Tree_${DTYPE}_${ITYPE};
+% endfor
% endfor
% for DTYPE in ['float', 'double']:
-void insert_point_${DTYPE}(uint32_t *closest_idx, ${DTYPE} *closest_dist, uint32_t pidx, ${DTYPE} cur_dist, uint32_t k);
-void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, ${DTYPE} *bbox);
-int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *bbox, int8_t *cut_dim,
- ${DTYPE} *cut_val, uint32_t *n_lo);
-Tree_${DTYPE}* construct_tree_${DTYPE}(${DTYPE} *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
-Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, ${DTYPE} *bbox);
-Node_${DTYPE} * create_node_${DTYPE}(uint32_t start_idx, uint32_t n, int is_leaf);
-void delete_subtree_${DTYPE}(Node_${DTYPE} *root);
-void delete_tree_${DTYPE}(Tree_${DTYPE} *tree);
-void print_tree_${DTYPE}(Node_${DTYPE} *root, int level);
${DTYPE} calc_dist_${DTYPE}(${DTYPE} *point1_coord, ${DTYPE} *point2_coord, int8_t no_dims);
${DTYPE} get_cube_offset_${DTYPE}(int8_t dim, ${DTYPE} *point_coord, ${DTYPE} *bbox);
${DTYPE} get_min_dist_${DTYPE}(${DTYPE} *point_coord, int8_t no_dims, ${DTYPE} *bbox);
-void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
- uint32_t k, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist);
-void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
- uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist);
-void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, ${DTYPE} *point_coord,
- ${DTYPE} min_dist, uint32_t k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask, uint32_t * closest_idx, ${DTYPE} *closest_dist);
-void search_tree_${DTYPE}(Tree_${DTYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coords,
- uint32_t num_points, uint32_t k, ${DTYPE} distance_upper_bound,
- ${DTYPE} eps, uint8_t *mask, uint32_t *closest_idxs, ${DTYPE} *closest_dists);
+% for ITYPE in ['int32_t', 'int64_t']:
+
+void insert_point_${DTYPE}_${ITYPE}(u${ITYPE} *closest_idx, ${DTYPE} *closest_dist, u${ITYPE} pidx, ${DTYPE} cur_dist, u${ITYPE} k);
+void get_bounding_box_${DTYPE}_${ITYPE}(${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, u${ITYPE} n, ${DTYPE} *bbox);
+int partition_${DTYPE}_${ITYPE}(${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, ${DTYPE} *bbox, int8_t *cut_dim,
+ ${DTYPE} *cut_val, u${ITYPE} *n_lo);
+Tree_${DTYPE}_${ITYPE}* construct_tree_${DTYPE}_${ITYPE}(${DTYPE} *pa, int8_t no_dims, u${ITYPE} n, u${ITYPE} bsp);
+Node_${DTYPE}_${ITYPE}* construct_subtree_${DTYPE}_${ITYPE}(${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, u${ITYPE} bsp, ${DTYPE} *bbox);
+Node_${DTYPE}_${ITYPE} * create_node_${DTYPE}_${ITYPE}(u${ITYPE} start_idx, u${ITYPE} n, int is_leaf);
+void delete_subtree_${DTYPE}_${ITYPE}(Node_${DTYPE}_${ITYPE} *root);
+void delete_tree_${DTYPE}_${ITYPE}(Tree_${DTYPE}_${ITYPE} *tree);
+void print_tree_${DTYPE}_${ITYPE}(Node_${DTYPE}_${ITYPE} *root, int level);
+void search_leaf_${DTYPE}_${ITYPE}(${DTYPE} *restrict pa, u${ITYPE} *restrict pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, ${DTYPE} *restrict point_coord,
+ u${ITYPE} k, u${ITYPE} *restrict closest_idx, ${DTYPE} *restrict closest_dist);
+void search_leaf_${DTYPE}_${ITYPE}_mask(${DTYPE} *restrict pa, u${ITYPE} *restrict pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, ${DTYPE} *restrict point_coord,
+ u${ITYPE} k, uint8_t *restrict mask, u${ITYPE} *restrict closest_idx, ${DTYPE} *restrict closest_dist);
+void search_splitnode_${DTYPE}_${ITYPE}(Node_${DTYPE}_${ITYPE} *root, ${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, ${DTYPE} *point_coord,
+ ${DTYPE} min_dist, u${ITYPE} k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask, u${ITYPE} * closest_idx, ${DTYPE} *closest_dist);
+void search_tree_${DTYPE}_${ITYPE}(Tree_${DTYPE}_${ITYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coords,
+ u${ITYPE} num_points, u${ITYPE} k, ${DTYPE} distance_upper_bound,
+ ${DTYPE} eps, uint8_t *mask, u${ITYPE} *closest_idxs, ${DTYPE} *closest_dists);
+
+% endfor
% endfor
% for DTYPE in ['float', 'double']:
+/************************************************
+Calculate squared cartesian distance between points
+Params:
+ point1_coord : point 1
+ point2_coord : point 2
+************************************************/
+${DTYPE} calc_dist_${DTYPE}(${DTYPE} *point1_coord, ${DTYPE} *point2_coord, int8_t no_dims)
+{
+ /* Calculate squared distance */
+ ${DTYPE} dist = 0, dim_dist;
+ int8_t i;
+ for (i = 0; i < no_dims; i++)
+ {
+ dim_dist = point2_coord[i] - point1_coord[i];
+ dist += dim_dist * dim_dist;
+ }
+ return dist;
+}
+
+/************************************************
+Get squared distance from point to cube in specified dimension
+Params:
+ dim : dimension
+ point_coord : cartesian coordinates of point
+ bbox : cube
+************************************************/
+${DTYPE} get_cube_offset_${DTYPE}(int8_t dim, ${DTYPE} *point_coord, ${DTYPE} *bbox)
+{
+ ${DTYPE} dim_coord = point_coord[dim];
+
+ if (dim_coord < bbox[2 * dim])
+ {
+ /* Left of cube in dimension */
+ return dim_coord - bbox[2 * dim];
+ }
+ else if (dim_coord > bbox[2 * dim + 1])
+ {
+ /* Right of cube in dimension */
+ return dim_coord - bbox[2 * dim + 1];
+ }
+ else
+ {
+ /* Inside cube in dimension */
+ return 0.;
+ }
+}
+
+/************************************************
+Get minimum squared distance between point and cube.
+Params:
+ point_coord : cartesian coordinates of point
+ no_dims : number of dimensions
+ bbox : cube
+************************************************/
+${DTYPE} get_min_dist_${DTYPE}(${DTYPE} *point_coord, int8_t no_dims, ${DTYPE} *bbox)
+{
+ ${DTYPE} cube_offset = 0, cube_offset_dim;
+ int8_t i;
+
+ for (i = 0; i < no_dims; i++)
+ {
+ cube_offset_dim = get_cube_offset_${DTYPE}(i, point_coord, bbox);
+ cube_offset += cube_offset_dim * cube_offset_dim;
+ }
+
+ return cube_offset;
+}
+
+% for ITYPE in ['int32_t', 'int64_t']:
+
/************************************************
Insert point into priority queue
Params:
@@ -97,7 +180,7 @@ Params:
cur_dist : distance to point inserted
k : number of neighbours
************************************************/
-void insert_point_${DTYPE}(uint32_t *closest_idx, ${DTYPE} *closest_dist, uint32_t pidx, ${DTYPE} cur_dist, uint32_t k)
+void insert_point_${DTYPE}_${ITYPE}(u${ITYPE} *closest_idx, ${DTYPE} *closest_dist, u${ITYPE} pidx, ${DTYPE} cur_dist, u${ITYPE} k)
{
int i;
for (i = k - 1; i > 0; i--)
@@ -125,11 +208,11 @@ Params:
n : number of points
bbox : bounding box (return)
************************************************/
-void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, ${DTYPE} *bbox)
+void get_bounding_box_${DTYPE}_${ITYPE}(${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, u${ITYPE} n, ${DTYPE} *bbox)
{
${DTYPE} cur;
int8_t i, j;
- uint32_t bbox_idx, i2;
+ u${ITYPE} bbox_idx, i2;
/* Use first data point to initialize */
for (i = 0; i < no_dims; i++)
@@ -170,12 +253,12 @@ Params:
cut_val : value of cutting point (return)
n_lo : number of point below cutting plane (return)
************************************************/
-int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *bbox, int8_t *cut_dim, ${DTYPE} *cut_val, uint32_t *n_lo)
+int partition_${DTYPE}_${ITYPE}(${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, ${DTYPE} *bbox, int8_t *cut_dim, ${DTYPE} *cut_val, u${ITYPE} *n_lo)
{
int8_t dim = 0, i;
- uint32_t p, q, i2;
+ u${ITYPE} p, q, i2;
${DTYPE} size = 0, min_val, max_val, split, side_len, cur_val;
- uint32_t end_idx = start_idx + n - 1;
+ u${ITYPE} end_idx = start_idx + n - 1;
/* Find largest bounding box side */
for (i = 0; i < no_dims; i++)
@@ -221,7 +304,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
}
else
{
- PASWAP(p, q);
+ PASWAP_${ITYPE}(p, q);
p++;
q--;
}
@@ -235,7 +318,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
Minimum 1 point will be in lower box.
*/
- uint32_t j = start_idx;
+ u${ITYPE} j = start_idx;
split = PA(j, dim);
for (i2 = start_idx + 1; i2 <= end_idx; i2++)
{
@@ -247,7 +330,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
split = cur_val;
}
}
- PASWAP(j, start_idx);
+ PASWAP_${ITYPE}(j, start_idx);
p = start_idx + 1;
}
else if (p == end_idx + 1)
@@ -257,7 +340,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
Minimum 1 point will be in higher box.
*/
- uint32_t j = end_idx;
+ u${ITYPE} j = end_idx;
split = PA(j, dim);
for (i2 = start_idx; i2 < end_idx; i2++)
{
@@ -269,7 +352,7 @@ int partition_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t st
split = cur_val;
}
}
- PASWAP(j, end_idx);
+ PASWAP_${ITYPE}(j, end_idx);
p = end_idx;
}
@@ -291,14 +374,14 @@ Params:
bsp : number of points per leaf
bbox : bounding box of set of data points
************************************************/
-Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, ${DTYPE} *bbox)
+Node_${DTYPE}_${ITYPE}* construct_subtree_${DTYPE}_${ITYPE}(${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, u${ITYPE} bsp, ${DTYPE} *bbox)
{
/* Create new node */
int is_leaf = (n <= bsp);
- Node_${DTYPE} *root = create_node_${DTYPE}(start_idx, n, is_leaf);
+ Node_${DTYPE}_${ITYPE} *root = create_node_${DTYPE}_${ITYPE}(start_idx, n, is_leaf);
int rval;
int8_t cut_dim;
- uint32_t n_lo;
+ u${ITYPE} n_lo;
${DTYPE} cut_val, lv, hv;
if (is_leaf)
{
@@ -309,7 +392,7 @@ Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t n
{
/* Make split node */
/* Partition data set and set node info */
- rval = partition_${DTYPE}(pa, pidx, no_dims, start_idx, n, bbox, &cut_dim, &cut_val, &n_lo);
+ rval = partition_${DTYPE}_${ITYPE}(pa, pidx, no_dims, start_idx, n, bbox, &cut_dim, &cut_val, &n_lo);
if (rval == 1)
{
root->cut_dim = -1;
@@ -328,12 +411,12 @@ Node_${DTYPE}* construct_subtree_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t n
/* Update bounding box before call to lower subset and restore after */
bbox[2 * cut_dim + 1] = cut_val;
- root->left_child = (struct Node_${DTYPE} *)construct_subtree_${DTYPE}(pa, pidx, no_dims, start_idx, n_lo, bsp, bbox);
+ root->left_child = (struct Node_${DTYPE}_${ITYPE} *)construct_subtree_${DTYPE}_${ITYPE}(pa, pidx, no_dims, start_idx, n_lo, bsp, bbox);
bbox[2 * cut_dim + 1] = hv;
/* Update bounding box before call to higher subset and restore after */
bbox[2 * cut_dim] = cut_val;
- root->right_child = (struct Node_${DTYPE} *)construct_subtree_${DTYPE}(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
+ root->right_child = (struct Node_${DTYPE}_${ITYPE} *)construct_subtree_${DTYPE}_${ITYPE}(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
bbox[2 * cut_dim] = lv;
}
return root;
@@ -347,28 +430,28 @@ Params:
n : number of data points
bsp : number of points per leaf
************************************************/
-Tree_${DTYPE}* construct_tree_${DTYPE}(${DTYPE} *pa, int8_t no_dims, uint32_t n, uint32_t bsp)
+Tree_${DTYPE}_${ITYPE}* construct_tree_${DTYPE}_${ITYPE}(${DTYPE} *pa, int8_t no_dims, u${ITYPE} n, u${ITYPE} bsp)
{
- Tree_${DTYPE} *tree = (Tree_${DTYPE} *)malloc(sizeof(Tree_${DTYPE}));
- uint32_t i;
- uint32_t *pidx;
+ Tree_${DTYPE}_${ITYPE} *tree = (Tree_${DTYPE}_${ITYPE} *)malloc(sizeof(Tree_${DTYPE}_${ITYPE}));
+ u${ITYPE} i;
+ u${ITYPE} *pidx;
${DTYPE} *bbox;
tree->no_dims = no_dims;
/* Initialize permutation array */
- pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
+ pidx = (u${ITYPE} *)malloc(sizeof(u${ITYPE}) * n);
for (i = 0; i < n; i++)
{
pidx[i] = i;
}
bbox = (${DTYPE} *)malloc(2 * sizeof(${DTYPE}) * no_dims);
- get_bounding_box_${DTYPE}(pa, pidx, no_dims, n, bbox);
+ get_bounding_box_${DTYPE}_${ITYPE}(pa, pidx, no_dims, n, bbox);
tree->bbox = bbox;
/* Construct subtree on full dataset */
- tree->root = (struct Node_${DTYPE} *)construct_subtree_${DTYPE}(pa, pidx, no_dims, 0, n, bsp, bbox);
+ tree->root = (struct Node_${DTYPE}_${ITYPE} *)construct_subtree_${DTYPE}_${ITYPE}(pa, pidx, no_dims, 0, n, bsp, bbox);
tree->pidx = pidx;
return tree;
@@ -380,9 +463,9 @@ Params:
start_idx : index of first data point to use
n : number of data points
************************************************/
-Node_${DTYPE}* create_node_${DTYPE}(uint32_t start_idx, uint32_t n, int is_leaf)
+Node_${DTYPE}_${ITYPE}* create_node_${DTYPE}_${ITYPE}(u${ITYPE} start_idx, u${ITYPE} n, int is_leaf)
{
- Node_${DTYPE} *new_node;
+ Node_${DTYPE}_${ITYPE} *new_node;
if (is_leaf)
{
/*
@@ -390,11 +473,11 @@ Node_${DTYPE}* create_node_${DTYPE}(uint32_t start_idx, uint32_t n, int is_leaf)
This relies on the C99 specification of struct layout conservation and padding and
that dereferencing is never attempted for the node pointers in a leaf.
*/
- new_node = (Node_${DTYPE} *)malloc(sizeof(Node_${DTYPE}) - 2 * sizeof(Node_${DTYPE} *));
+ new_node = (Node_${DTYPE}_${ITYPE} *)malloc(sizeof(Node_${DTYPE}_${ITYPE}) - 2 * sizeof(Node_${DTYPE}_${ITYPE} *));
}
else
{
- new_node = (Node_${DTYPE} *)malloc(sizeof(Node_${DTYPE}));
+ new_node = (Node_${DTYPE}_${ITYPE} *)malloc(sizeof(Node_${DTYPE}_${ITYPE}));
}
new_node->n = n;
new_node->start_idx = start_idx;
@@ -406,12 +489,12 @@ Delete subtree
Params:
root : root node of subtree to delete
************************************************/
-void delete_subtree_${DTYPE}(Node_${DTYPE} *root)
+void delete_subtree_${DTYPE}_${ITYPE}(Node_${DTYPE}_${ITYPE} *root)
{
if (root->cut_dim != -1)
{
- delete_subtree_${DTYPE}((Node_${DTYPE} *)root->left_child);
- delete_subtree_${DTYPE}((Node_${DTYPE} *)root->right_child);
+ delete_subtree_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->left_child);
+ delete_subtree_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->right_child);
}
free(root);
}
@@ -421,9 +504,9 @@ Delete tree
Params:
tree : Tree struct of kd tree
************************************************/
-void delete_tree_${DTYPE}(Tree_${DTYPE} *tree)
+void delete_tree_${DTYPE}_${ITYPE}(Tree_${DTYPE}_${ITYPE} *tree)
{
- delete_subtree_${DTYPE}((Node_${DTYPE} *)tree->root);
+ delete_subtree_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)tree->root);
free(tree->bbox);
free(tree->pidx);
free(tree);
@@ -432,7 +515,7 @@ void delete_tree_${DTYPE}(Tree_${DTYPE} *tree)
/************************************************
Print
************************************************/
-void print_tree_${DTYPE}(Node_${DTYPE} *root, int level)
+void print_tree_${DTYPE}_${ITYPE}(Node_${DTYPE}_${ITYPE} *root, int level)
{
int i;
for (i = 0; i < level; i++)
@@ -441,77 +524,9 @@ void print_tree_${DTYPE}(Node_${DTYPE} *root, int level)
}
printf("(cut_val: %f, cut_dim: %i)\n", root->cut_val, root->cut_dim);
if (root->cut_dim != -1)
- print_tree_${DTYPE}((Node_${DTYPE} *)root->left_child, level + 1);
+ print_tree_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->left_child, level + 1);
if (root->cut_dim != -1)
- print_tree_${DTYPE}((Node_${DTYPE} *)root->right_child, level + 1);
-}
-
-/************************************************
-Calculate squared cartesian distance between points
-Params:
- point1_coord : point 1
- point2_coord : point 2
-************************************************/
-${DTYPE} calc_dist_${DTYPE}(${DTYPE} *point1_coord, ${DTYPE} *point2_coord, int8_t no_dims)
-{
- /* Calculate squared distance */
- ${DTYPE} dist = 0, dim_dist;
- int8_t i;
- for (i = 0; i < no_dims; i++)
- {
- dim_dist = point2_coord[i] - point1_coord[i];
- dist += dim_dist * dim_dist;
- }
- return dist;
-}
-
-/************************************************
-Get squared distance from point to cube in specified dimension
-Params:
- dim : dimension
- point_coord : cartesian coordinates of point
- bbox : cube
-************************************************/
-${DTYPE} get_cube_offset_${DTYPE}(int8_t dim, ${DTYPE} *point_coord, ${DTYPE} *bbox)
-{
- ${DTYPE} dim_coord = point_coord[dim];
-
- if (dim_coord < bbox[2 * dim])
- {
- /* Left of cube in dimension */
- return dim_coord - bbox[2 * dim];
- }
- else if (dim_coord > bbox[2 * dim + 1])
- {
- /* Right of cube in dimension */
- return dim_coord - bbox[2 * dim + 1];
- }
- else
- {
- /* Inside cube in dimension */
- return 0.;
- }
-}
-
-/************************************************
-Get minimum squared distance between point and cube.
-Params:
- point_coord : cartesian coordinates of point
- no_dims : number of dimensions
- bbox : cube
-************************************************/
-${DTYPE} get_min_dist_${DTYPE}(${DTYPE} *point_coord, int8_t no_dims, ${DTYPE} *bbox)
-{
- ${DTYPE} cube_offset = 0, cube_offset_dim;
- int8_t i;
-
- for (i = 0; i < no_dims; i++)
- {
- cube_offset_dim = get_cube_offset_${DTYPE}(i, point_coord, bbox);
- cube_offset += cube_offset_dim * cube_offset_dim;
- }
-
- return cube_offset;
+ print_tree_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->right_child, level + 1);
}
/************************************************
@@ -526,11 +541,11 @@ Params:
closest_idx : index of closest data point found (return)
closest_dist : distance to closest point (return)
************************************************/
-void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
- uint32_t k, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist)
+void search_leaf_${DTYPE}_${ITYPE}(${DTYPE} *restrict pa, u${ITYPE} *restrict pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, ${DTYPE} *restrict point_coord,
+ u${ITYPE} k, u${ITYPE} *restrict closest_idx, ${DTYPE} *restrict closest_dist)
{
${DTYPE} cur_dist;
- uint32_t i;
+ u${ITYPE} i;
/* Loop through all points in leaf */
for (i = 0; i < n; i++)
{
@@ -539,7 +554,7 @@ void search_leaf_${DTYPE}(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t
/* Update closest info if new point is closest so far*/
if (cur_dist < closest_dist[k - 1])
{
- insert_point_${DTYPE}(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
+ insert_point_${DTYPE}_${ITYPE}(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
}
}
}
@@ -558,11 +573,11 @@ Params:
closest_idx : index of closest data point found (return)
closest_dist : distance to closest point (return)
************************************************/
-void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, ${DTYPE} *restrict point_coord,
- uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, ${DTYPE} *restrict closest_dist)
+void search_leaf_${DTYPE}_${ITYPE}_mask(${DTYPE} *restrict pa, u${ITYPE} *restrict pidx, int8_t no_dims, u${ITYPE} start_idx, u${ITYPE} n, ${DTYPE} *restrict point_coord,
+ u${ITYPE} k, uint8_t *mask, u${ITYPE} *restrict closest_idx, ${DTYPE} *restrict closest_dist)
{
${DTYPE} cur_dist;
- uint32_t i;
+ u${ITYPE} i;
/* Loop through all points in leaf */
for (i = 0; i < n; i++)
{
@@ -576,7 +591,7 @@ void search_leaf_${DTYPE}_mask(${DTYPE} *restrict pa, uint32_t *restrict pidx, i
/* Update closest info if new point is closest so far*/
if (cur_dist < closest_dist[k - 1])
{
- insert_point_${DTYPE}(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
+ insert_point_${DTYPE}_${ITYPE}(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
}
}
}
@@ -594,9 +609,9 @@ Params:
closest_idx : index of closest data point found (return)
closest_dist : distance to closest point (return)
************************************************/
-void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, ${DTYPE} *point_coord,
- ${DTYPE} min_dist, uint32_t k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask,
- uint32_t *closest_idx, ${DTYPE} *closest_dist)
+void search_splitnode_${DTYPE}_${ITYPE}(Node_${DTYPE}_${ITYPE} *root, ${DTYPE} *pa, u${ITYPE} *pidx, int8_t no_dims, ${DTYPE} *point_coord,
+ ${DTYPE} min_dist, u${ITYPE} k, ${DTYPE} distance_upper_bound, ${DTYPE} eps_fac, uint8_t *mask,
+ u${ITYPE} *closest_idx, ${DTYPE} *closest_dist)
{
int8_t dim;
${DTYPE} dist_left, dist_right;
@@ -616,11 +631,11 @@ void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx
{
if (mask)
{
- search_leaf_${DTYPE}_mask(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, mask, closest_idx, closest_dist);
+ search_leaf_${DTYPE}_${ITYPE}_mask(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, mask, closest_idx, closest_dist);
}
else
{
- search_leaf_${DTYPE}(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
+ search_leaf_${DTYPE}_${ITYPE}(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
}
return;
}
@@ -635,7 +650,7 @@ void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx
if (dist_left < closest_dist[k - 1] * eps_fac)
{
/* Search left subtree if minimum distance is below limit */
- search_splitnode_${DTYPE}((Node_${DTYPE} *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
+ search_splitnode_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
/* Right of cutting plane. Update minimum distance.
@@ -650,7 +665,7 @@ void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx
if (dist_right < closest_dist[k - 1] * eps_fac)
{
/* Search right subtree if minimum distance is below limit*/
- search_splitnode_${DTYPE}((Node_${DTYPE} *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
+ search_splitnode_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
}
else
@@ -660,7 +675,7 @@ void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx
if (dist_right < closest_dist[k - 1] * eps_fac)
{
/* Search right subtree if minimum distance is below limit*/
- search_splitnode_${DTYPE}((Node_${DTYPE} *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
+ search_splitnode_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
/* Left of cutting plane. Update minimum distance.
@@ -675,7 +690,7 @@ void search_splitnode_${DTYPE}(Node_${DTYPE} *root, ${DTYPE} *pa, uint32_t *pidx
if (dist_left < closest_dist[k - 1] * eps_fac)
{
/* Search left subtree if minimum distance is below limit*/
- search_splitnode_${DTYPE}((Node_${DTYPE} *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
+ search_splitnode_${DTYPE}_${ITYPE}((Node_${DTYPE}_${ITYPE} *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
}
}
@@ -692,24 +707,20 @@ Params:
closest_idx : index of closest data point found (return)
closest_dist : distance to closest point (return)
************************************************/
-void search_tree_${DTYPE}(Tree_${DTYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coords,
- uint32_t num_points, uint32_t k, ${DTYPE} distance_upper_bound,
- ${DTYPE} eps, uint8_t *mask, uint32_t *closest_idxs, ${DTYPE} *closest_dists)
+void search_tree_${DTYPE}_${ITYPE}(Tree_${DTYPE}_${ITYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coords,
+ u${ITYPE} num_points, u${ITYPE} k, ${DTYPE} distance_upper_bound,
+ ${DTYPE} eps, uint8_t *mask, u${ITYPE} *closest_idxs, ${DTYPE} *closest_dists)
{
${DTYPE} min_dist;
${DTYPE} eps_fac = 1 / ((1 + eps) * (1 + eps));
int8_t no_dims = tree->no_dims;
${DTYPE} *bbox = tree->bbox;
- uint32_t *pidx = tree->pidx;
- uint32_t j = 0;
-#if defined(_MSC_VER) && defined(_OPENMP)
- int32_t i = 0;
- int32_t local_num_points = (int32_t) num_points;
-#else
- uint32_t i;
- uint32_t local_num_points = num_points;
-#endif
- Node_${DTYPE} *root = (Node_${DTYPE} *)tree->root;
+ u${ITYPE} *pidx = tree->pidx;
+ /* use 64-bit ints for indexing to avoid overflow, use signed ints to support all Openmp implementations */
+ int64_t i = 0;
+ int64_t j = 0;
+ int64_t local_num_points = (int64_t) num_points;
+ Node_${DTYPE}_${ITYPE} *root = (Node_${DTYPE}_${ITYPE} *)tree->root;
/* Queries are OpenMP enabled */
#pragma omp parallel
@@ -722,13 +733,14 @@ void search_tree_${DTYPE}(Tree_${DTYPE} *tree, ${DTYPE} *pa, ${DTYPE} *point_coo
{
for (j = 0; j < k; j++)
{
- closest_idxs[i * k + j] = UINT32_MAX;
- closest_dists[i * k + j] = DBL_MAX;
+ closest_idxs[i * k + j] = IDX_MAX_${ITYPE};
+ closest_dists[i * k + j] = DIST_MAX_${DTYPE};
}
min_dist = get_min_dist_${DTYPE}(point_coords + no_dims * i, no_dims, bbox);
- search_splitnode_${DTYPE}(root, pa, pidx, no_dims, point_coords + no_dims * i, min_dist,
+ search_splitnode_${DTYPE}_${ITYPE}(root, pa, pidx, no_dims, point_coords + no_dims * i, min_dist,
k, distance_upper_bound, eps_fac, mask, &closest_idxs[i * k], &closest_dists[i * k]);
}
}
}
% endfor
+% endfor
=====================================
pykdtree/kdtree.pyx
=====================================
@@ -17,51 +17,91 @@
import numpy as np
cimport numpy as np
-from libc.stdint cimport uint32_t, int8_t, uint8_t
+from libc.stdint cimport uint64_t, uint32_t, int8_t, uint8_t, UINT32_MAX
cimport cython
np.import_array()
# Node structure
-cdef struct node_float:
+cdef struct node_float_int32_t:
float cut_val
int8_t cut_dim
uint32_t start_idx
uint32_t n
float cut_bounds_lv
float cut_bounds_hv
- node_float *left_child
- node_float *right_child
+ node_float_int32_t *left_child
+ node_float_int32_t *right_child
-cdef struct tree_float:
+cdef struct tree_float_int32_t:
float *bbox
int8_t no_dims
uint32_t *pidx
- node_float *root
+ node_float_int32_t *root
-cdef struct node_double:
+cdef struct node_double_int32_t:
double cut_val
int8_t cut_dim
uint32_t start_idx
uint32_t n
double cut_bounds_lv
double cut_bounds_hv
- node_double *left_child
- node_double *right_child
+ node_double_int32_t *left_child
+ node_double_int32_t *right_child
-cdef struct tree_double:
+cdef struct tree_double_int32_t:
double *bbox
int8_t no_dims
uint32_t *pidx
- node_double *root
+ node_double_int32_t *root
-cdef extern tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp) nogil
-cdef extern void search_tree_float(tree_float *kdtree, float *pa, float *point_coords, uint32_t num_points, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists) nogil
-cdef extern void delete_tree_float(tree_float *kdtree)
+cdef struct node_float_int64_t:
+ float cut_val
+ int8_t cut_dim
+ uint64_t start_idx
+ uint64_t n
+ float cut_bounds_lv
+ float cut_bounds_hv
+ node_float_int64_t *left_child
+ node_float_int64_t *right_child
+
+cdef struct tree_float_int64_t:
+ float *bbox
+ int8_t no_dims
+ uint64_t *pidx
+ node_float_int64_t *root
-cdef extern tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp) nogil
-cdef extern void search_tree_double(tree_double *kdtree, double *pa, double *point_coords, uint32_t num_points, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists) nogil
-cdef extern void delete_tree_double(tree_double *kdtree)
+cdef struct node_double_int64_t:
+ double cut_val
+ int8_t cut_dim
+ uint64_t start_idx
+ uint64_t n
+ double cut_bounds_lv
+ double cut_bounds_hv
+ node_double_int64_t *left_child
+ node_double_int64_t *right_child
+
+cdef struct tree_double_int64_t:
+ double *bbox
+ int8_t no_dims
+ uint64_t *pidx
+ node_double_int64_t *root
+
+cdef extern tree_float_int32_t* construct_tree_float_int32_t(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp) nogil
+cdef extern void search_tree_float_int32_t(tree_float_int32_t *kdtree, float *pa, float *point_coords, uint32_t num_points, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists) nogil
+cdef extern void delete_tree_float_int32_t(tree_float_int32_t *kdtree)
+
+cdef extern tree_double_int32_t* construct_tree_double_int32_t(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp) nogil
+cdef extern void search_tree_double_int32_t(tree_double_int32_t *kdtree, double *pa, double *point_coords, uint32_t num_points, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists) nogil
+cdef extern void delete_tree_double_int32_t(tree_double_int32_t *kdtree)
+
+cdef extern tree_float_int64_t* construct_tree_float_int64_t(float *pa, int8_t no_dims, uint64_t n, uint64_t bsp) nogil
+cdef extern void search_tree_float_int64_t(tree_float_int64_t *kdtree, float *pa, float *point_coords, uint64_t num_points, uint64_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint64_t *closest_idxs, float *closest_dists) nogil
+cdef extern void delete_tree_float_int64_t(tree_float_int64_t *kdtree)
+
+cdef extern tree_double_int64_t* construct_tree_double_int64_t(double *pa, int8_t no_dims, uint64_t n, uint64_t bsp) nogil
+cdef extern void search_tree_double_int64_t(tree_double_int64_t *kdtree, double *pa, double *point_coords, uint64_t num_points, uint64_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint64_t *closest_idxs, double *closest_dists) nogil
+cdef extern void delete_tree_double_int64_t(tree_double_int64_t *kdtree)
cdef class KDTree:
"""kd-tree for fast nearest-neighbour lookup.
@@ -75,19 +115,24 @@ cdef class KDTree:
Maximum number of data points in tree leaf
"""
- cdef tree_float *_kdtree_float
- cdef tree_double *_kdtree_double
+ cdef tree_float_int32_t *_kdtree_float_int32_t
+ cdef tree_double_int32_t *_kdtree_double_int32_t
+ cdef tree_float_int64_t *_kdtree_float_int64_t
+ cdef tree_double_int64_t *_kdtree_double_int64_t
+ cdef readonly bint _use_int32_t
cdef readonly np.ndarray data_pts
cdef readonly np.ndarray data
cdef float *_data_pts_data_float
cdef double *_data_pts_data_double
- cdef readonly uint32_t n
+ cdef readonly uint64_t n
cdef readonly int8_t ndim
cdef readonly uint32_t leafsize
def __cinit__(KDTree self):
- self._kdtree_float = NULL
- self._kdtree_double = NULL
+ self._kdtree_float_int32_t = NULL
+ self._kdtree_double_int32_t = NULL
+ self._kdtree_float_int64_t = NULL
+ self._kdtree_double_int64_t = NULL
def __init__(KDTree self, np.ndarray data_pts not None, int leafsize=16):
@@ -116,7 +161,8 @@ cdef class KDTree:
self.data = self.data_pts
# Get tree info
- self.n = <uint32_t>data_pts.shape[0]
+ self.n = <uint64_t>data_pts.shape[0]
+ self._use_int32_t = self.n * data_pts.shape[1] < UINT32_MAX
self.leafsize = <uint32_t>leafsize
if data_pts.ndim == 1:
self.ndim = 1
@@ -127,13 +173,23 @@ cdef class KDTree:
# Release GIL and construct tree
if data_pts.dtype == np.float32:
- with nogil:
- self._kdtree_float = construct_tree_float(self._data_pts_data_float, self.ndim,
- self.n, self.leafsize)
+ if self._use_int32_t:
+ with nogil:
+ self._kdtree_float_int32_t = construct_tree_float_int32_t(self._data_pts_data_float, self.ndim,
+ <uint32_t>self.n, self.leafsize)
+ else:
+ with nogil:
+ self._kdtree_float_int64_t = construct_tree_float_int64_t(self._data_pts_data_float, self.ndim,
+ self.n, self.leafsize)
else:
- with nogil:
- self._kdtree_double = construct_tree_double(self._data_pts_data_double, self.ndim,
- self.n, self.leafsize)
+ if self._use_int32_t:
+ with nogil:
+ self._kdtree_double_int32_t = construct_tree_double_int32_t(self._data_pts_data_double, self.ndim,
+ <uint32_t>self.n, self.leafsize)
+ else:
+ with nogil:
+ self._kdtree_double_int64_t = construct_tree_double_int64_t(self._data_pts_data_double, self.ndim,
+ self.n, self.leafsize)
def query(KDTree self, np.ndarray query_pts not None, k=1, eps=0,
@@ -185,18 +241,27 @@ cdef class KDTree:
raise TypeError('Type mismatch. query points must be of type float32 when data points are of type float32')
# Get query info
- cdef uint32_t num_qpoints = query_pts.shape[0]
- cdef uint32_t num_n = k
- cdef np.ndarray[uint32_t, ndim=1] closest_idxs = np.empty(num_qpoints * k, dtype=np.uint32)
+ cdef uint64_t num_qpoints = query_pts.shape[0]
+ cdef uint64_t num_n = k
+ cdef np.ndarray[uint32_t, ndim=1] closest_idxs_int32_t
+ cdef np.ndarray[uint64_t, ndim=1] closest_idxs_int64_t
cdef np.ndarray[float, ndim=1] closest_dists_float
cdef np.ndarray[double, ndim=1] closest_dists_double
-
# Set up return arrays
- cdef uint32_t *closest_idxs_data = <uint32_t *>closest_idxs.data
+ cdef uint32_t *closest_idxs_data_int32_t
+ cdef uint64_t *closest_idxs_data_int64_t
cdef float *closest_dists_data_float
cdef double *closest_dists_data_double
-
+ if self._use_int32_t:
+ closest_idxs_int32_t = np.empty(num_qpoints * k, dtype=np.uint32)
+ closest_idxs = closest_idxs_int32_t
+ closest_idxs_data_int32_t = <uint32_t *>closest_idxs_int32_t.data
+ else:
+ closest_idxs_int64_t = np.empty(num_qpoints * k, dtype=np.uint64)
+ closest_idxs = closest_idxs_int64_t
+ closest_idxs_data_int64_t = <uint64_t *>closest_idxs_int64_t.data
+
# Get query points data
cdef np.ndarray[float, ndim=1] query_array_float
cdef np.ndarray[double, ndim=1] query_array_double
@@ -247,17 +312,28 @@ cdef class KDTree:
# Release GIL and query tree
if self.data_pts.dtype == np.float32:
- with nogil:
- search_tree_float(self._kdtree_float, self._data_pts_data_float,
- query_array_data_float, num_qpoints, num_n, dub_float, epsilon_float,
- query_mask_data, closest_idxs_data, closest_dists_data_float)
-
+ if self._use_int32_t:
+ with nogil:
+ search_tree_float_int32_t(self._kdtree_float_int32_t, self._data_pts_data_float,
+ query_array_data_float, <uint32_t>num_qpoints, <uint32_t>num_n, dub_float, epsilon_float,
+ query_mask_data, closest_idxs_data_int32_t, closest_dists_data_float)
+ else:
+ with nogil:
+ search_tree_float_int64_t(self._kdtree_float_int64_t, self._data_pts_data_float,
+ query_array_data_float, num_qpoints, num_n, dub_float, epsilon_float,
+ query_mask_data, closest_idxs_data_int64_t, closest_dists_data_float)
else:
- with nogil:
- search_tree_double(self._kdtree_double, self._data_pts_data_double,
- query_array_data_double, num_qpoints, num_n, dub_double, epsilon_double,
- query_mask_data, closest_idxs_data, closest_dists_data_double)
-
+ if self._use_int32_t:
+ with nogil:
+ search_tree_double_int32_t(self._kdtree_double_int32_t, self._data_pts_data_double,
+ query_array_data_double, <uint32_t>num_qpoints, <uint32_t>num_n, dub_double, epsilon_double,
+ query_mask_data, closest_idxs_data_int32_t, closest_dists_data_double)
+ else:
+ with nogil:
+ search_tree_double_int64_t(self._kdtree_double_int64_t, self._data_pts_data_double,
+ query_array_data_double, num_qpoints, num_n, dub_double, epsilon_double,
+ query_mask_data, closest_idxs_data_int64_t, closest_dists_data_double)
+
# Shape result
if k > 1:
closest_dists_res = closest_dists.reshape(num_qpoints, k)
@@ -281,7 +357,11 @@ cdef class KDTree:
return closest_dists_res, closest_idxs_res
def __dealloc__(KDTree self):
- if self._kdtree_float != NULL:
- delete_tree_float(self._kdtree_float)
- elif self._kdtree_double != NULL:
- delete_tree_double(self._kdtree_double)
+ if self._kdtree_float_int32_t != NULL:
+ delete_tree_float_int32_t(self._kdtree_float_int32_t)
+ elif self._kdtree_double_int32_t != NULL:
+ delete_tree_double_int32_t(self._kdtree_double_int32_t)
+ if self._kdtree_float_int64_t != NULL:
+ delete_tree_float_int64_t(self._kdtree_float_int64_t)
+ elif self._kdtree_double_int64_t != NULL:
+ delete_tree_double_int64_t(self._kdtree_double_int64_t)
=====================================
pykdtree/test_tree.py
=====================================
@@ -1,4 +1,5 @@
import numpy as np
+import pytest
from pykdtree.kdtree import KDTree
@@ -383,3 +384,47 @@ def test_empty_fail():
kdtree = KDTree(data_pts)
except ValueError as e:
assert 'non-empty' in str(e), str(e)
+
+ at pytest.mark.skip(reason="Requires ~50G RAM")
+def test_tree_n_lt_maxint32_n_query_k_gt_maxint32():
+ # n_points < UINT32_MAX but n_query * k > UINT32_MAX -> still uses 32-bit index
+ n_dim = 2
+ n_points = 2**20
+ n_query = 2**20 + 8
+ k = 2**12
+ data_pts = np.random.random((n_points, n_dim)).astype(np.float32)
+ query_pts = np.random.random((n_query, n_dim)).astype(np.float32)
+ data_pts[0] = query_pts[0]
+ data_pts[1533] = query_pts[15633]
+ data_pts[1048575] = query_pts[1048583]
+ kdtree = KDTree(data_pts)
+ dist, idx = kdtree.query(query_pts, k=k)
+ assert idx.shape == (n_query, k)
+ assert idx.dtype == np.uint32
+ assert idx[0][0] == 0
+ assert idx[15633][0] == 1533
+ assert idx[1048583][0] == 1048575
+ assert np.all(idx < data_pts.shape[0])
+ assert dist.shape == (n_query, k)
+ assert dist.dtype == np.float32
+
+ at pytest.mark.skip(reason="Requires ~50G RAM")
+def test_tree_n_points_n_dim_gt_maxint32():
+ # n_points < UINT32_MAX but n_points * n_dim > UINT32_MAX -> uses 64-bit index
+ n_dim = 2**6
+ n_points = 2**26 + 8
+ data_pts = np.random.random((n_points, n_dim)).astype(np.float32)
+ query_pts = np.random.random((3, n_dim)).astype(np.float32)
+ data_pts[0] = query_pts[0]
+ data_pts[874516] = query_pts[1]
+ data_pts[67108871] = query_pts[2]
+ kdtree = KDTree(data_pts)
+ dist, idx = kdtree.query(query_pts, k=4)
+ assert idx.shape == (3, 4)
+ assert idx.dtype == np.uint64
+ assert idx[0][0] == 0
+ assert idx[1][0] == 874516
+ assert idx[2][0] == 67108871
+ assert np.all(idx < data_pts.shape[0])
+ assert dist.shape == (3, 4)
+ assert dist.dtype == np.float32
=====================================
setup.py
=====================================
@@ -197,7 +197,7 @@ extensions = [
setup(
name='pykdtree',
- version='1.3.13',
+ version='1.4.1',
url="https://github.com/storpipfugl/pykdtree",
description='Fast kd-tree implementation with OpenMP-enabled queries',
long_description=readme,
View it on GitLab: https://salsa.debian.org/debian-gis-team/pykdtree/-/compare/0fcb65cedd95ca968646ee66c723f55541d9958e...945155da706be858f1bf88c0a70953d0043f4a5a
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pykdtree/-/compare/0fcb65cedd95ca968646ee66c723f55541d9958e...945155da706be858f1bf88c0a70953d0043f4a5a
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20250131/3b2cd21f/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list