[Git][debian-gis-team/pykdtree][master] 6 commits: New upstream version 1.3.3+ds
Antonio Valentino
gitlab at salsa.debian.org
Tue Nov 10 07:51:35 GMT 2020
Antonio Valentino pushed to branch master at Debian GIS Project / pykdtree
Commits:
2dcb254f by Antonio Valentino at 2020-11-10T06:37:22+00:00
New upstream version 1.3.3+ds
- - - - -
d7fc9dea by Antonio Valentino at 2020-11-10T06:37:23+00:00
Update upstream source from tag 'upstream/1.3.3+ds'
Update to upstream version '1.3.3+ds'
with Debian dir 7017dedcc20debc602e4ea6ecc29ff5aa90afc4e
- - - - -
920ca3a5 by Antonio Valentino at 2020-11-10T07:43:03+01:00
New upstream release
- - - - -
0142f680 by Antonio Valentino at 2020-11-10T07:28:16+00:00
Ensure .c files generation
- - - - -
7e1213b5 by Antonio Valentino at 2020-11-10T07:35:39+00:00
Drop debian/source/lintian-overrides
- - - - -
d6b3c241 by Antonio Valentino at 2020-11-10T07:35:51+00:00
Set distribution to unstable
- - - - -
15 changed files:
- + .github/workflows/deploy-wheels.yml
- MANIFEST.in
- + Makefile
- README.rst
- debian/changelog
- debian/rules
- − debian/source/lintian-overrides
- − pykdtree/_kdtree_core.c
- pykdtree/_kdtree_core.c.mako
- − pykdtree/kdtree.c
- pykdtree/kdtree.pyx
- pykdtree/test_tree.py
- + scripts/build-manylinux-wheels.sh
- setup.cfg
- setup.py
Changes:
=====================================
.github/workflows/deploy-wheels.yml
=====================================
@@ -0,0 +1,76 @@
+---
+name: Deploy wheels
+
+on:
+ push:
+ release:
+ types:
+ - published
+jobs:
+ build:
+
+ runs-on: ${{ matrix.os }}
+ strategy:
+ fail-fast: false
+ matrix:
+ os: [windows-latest, macos-latest]
+ python-version: [3.6, 3.7, 3.8, 3.9]
+ include:
+ # Using pythons inside a docker image to provide all the Linux
+ # python-versions permutations.
+ - name: manylinux 64-bit
+ os: ubuntu-latest
+ python-version: 3.8
+ docker-image: manylinux1_x86_64
+ - name: manylinux 32-bit
+ os: ubuntu-latest
+ python-version: 3.8
+ docker-image: manylinux1_i686
+
+ steps:
+ - uses: actions/checkout at v2
+ - run: |
+ git fetch --prune --unshallow
+
+ - name: Set up Python ${{ matrix.python-version }}
+ uses: actions/setup-python at v1
+ with:
+ python-version: ${{ matrix.python-version }}
+
+ - name: Install dependencies
+ run: |
+ python -m pip install -U -q pip pytest wheel setuptools twine numpy
+
+ - name: Build and install macOS/Windows wheel
+ if: matrix.os != 'ubuntu-latest'
+ run: |
+ python setup.py -q bdist_wheel
+ pip install --find-links=./dist/ pykdtree
+
+ - name: Test
+ if: matrix.os != 'ubuntu-latest'
+ run: |
+ pytest -v --pyargs pykdtree
+
+ - name: Build Linux wheels inside docker
+ if: matrix.os == 'ubuntu-latest'
+ run: |
+ docker run \
+ -e PLAT=${{ matrix.docker-image }} \
+ -e USE_OMP=1 \
+ -v `pwd`:/io \
+ quay.io/pypa/${{ matrix.docker-image }} \
+ /io/scripts/build-manylinux-wheels.sh
+
+ - name: Upload wheel(s) as build artifacts
+ uses: actions/upload-artifact at v2
+ with:
+ name: wheels
+ path: dist/*.whl
+
+ - name: Publish package to PyPI
+ if: github.event.action == 'published'
+ env:
+ TWINE_USERNAME: ${{ secrets.pypi_username }}
+ TWINE_PASSWORD: ${{ secrets.pypi_password }}
+ run: twine upload --skip-existing dist/*.whl
=====================================
MANIFEST.in
=====================================
@@ -1,2 +1,4 @@
exclude pykdtree/render_template.py
include LICENSE.txt
+include README.rst
+include pykdtree/*.pyx
=====================================
Makefile
=====================================
@@ -0,0 +1,14 @@
+.PHONY: clean build
+
+build: pykdtree/kdtree.c pykdtree/_kdtree_core.c
+ python setup.py build
+
+pykdtree/kdtree.c: pykdtree/kdtree.pyx
+ cython pykdtree/kdtree.pyx
+
+pykdtree/_kdtree_core.c: pykdtree/_kdtree_core.c.mako
+ cd pykdtree && python render_template.py
+
+clean:
+ rm -rf build/
+ rm -f pykdtree/*.so
=====================================
README.rst
=====================================
@@ -20,14 +20,18 @@ Queries are optionally multithreaded using OpenMP.
Installation
------------
-Default build of pykdtree with OpenMP enabled queries using libgomp
+
+By default pykdtree is built with OpenMP enabled queries using libgomp except
+on OSX systems using the clang compiler (conda environments use a separate
+compiler).
.. code-block:: bash
$ cd <pykdtree_dir>
$ python setup.py install
-If it fails with undefined compiler flags or you want to use another OpenMP implementation please modify setup.py at the indicated point to match your system.
+If it fails with undefined compiler flags or you want to use another OpenMP
+implementation please modify setup.py at the indicated point to match your system.
Building without OpenMP support is controlled by the USE_OMP environment variable
@@ -43,6 +47,12 @@ Note evironment variables are by default not exported when using sudo so in this
$ USE_OMP=0 sudo -E python setup.py install
+Pykdtree can also be installed with conda via the conda-forge channel:
+
+.. code-block:: bash
+
+ $ conda install -c conda-forge pykdtree
+
Usage
-----
The usage of pykdtree is similar to scipy.spatial.cKDTree so for now refer to its documentation
@@ -125,6 +135,10 @@ turned off by adding the following to `appveyor.yml` in the
Changelog
---------
+v1.3.3 : Add compatibility to python 3.9
+
+v1.3.2 : Change OSX installation to not use OpenMP without conda interpreter
+
v1.3.1 : Fix masking in the "query" method introduced in 1.3.0
v1.3.0 : Keyword argument "mask" added to "query" method. OpenMP compilation now works for MS Visual Studio compiler
=====================================
debian/changelog
=====================================
@@ -1,14 +1,16 @@
-pykdtree (1.3.1+ds-2) UNRELEASED; urgency=medium
+pykdtree (1.3.3+ds-1) unstable; urgency=medium
[ Antonio Valentino ]
* Do not import autogenerated C sources.
* debian/rules:
- remove redundant call to cython3
+ - ensure .c source files generation
+ * Drop debian/source/lintian-overrides, no longer necessary.
[ Bas Couwenberg ]
* Update watch file to use repacksuffix.
- -- Antonio Valentino <antonio.valentino at tiscali.it> Wed, 14 Oct 2020 06:17:39 +0000
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Tue, 10 Nov 2020 07:35:44 +0000
pykdtree (1.3.1+ds-1) unstable; urgency=medium
=====================================
debian/rules
=====================================
@@ -13,12 +13,9 @@ export PYBUILD_NAME=pykdtree
dh $@ --with python3 --buildsystem=pybuild
override_dh_auto_configure:
- (cd $(CURDIR)/pykdtree && python3 render_template.py)
- if [ -f pykdtree/kdtree.c ] && [ ! -f pykdtree/kdtree.c.orig ]; then \
- mv pykdtree/kdtree.c pykdtree/kdtree.c.orig; \
- fi
- $(RM) pykdtree/kdtree.c
dh_auto_configure
+ cython3 pykdtree/kdtree.pyx
+ cd pykdtree && python3 render_template.py
override_dh_python3:
dh_python3 -ppython3-pykdtree
@@ -26,8 +23,5 @@ override_dh_python3:
override_dh_auto_clean:
dh_auto_clean
- if [ -f pykdtree/kdtree.c.orig ]; then \
- $(RM) pykdtree/kdtree.c; \
- mv pykdtree/kdtree.c.orig pykdtree/kdtree.c; \
- fi
$(RM) -r pykdtree.egg-info
+ $(RM) pykdtree/kdtree.c pykdtree/_kdtree_core.c
=====================================
debian/source/lintian-overrides deleted
=====================================
@@ -1,3 +0,0 @@
-# Source is not actually repacked, it swiched from PyPI tarballs to GitHub releases
-missing-explanation-for-repacked-upstream-tarball
-
=====================================
pykdtree/_kdtree_core.c deleted
=====================================
@@ -1,1417 +0,0 @@
-/*
-pykdtree, Fast kd-tree implementation with OpenMP-enabled queries
-
-Copyright (C) 2013 - present Esben S. Nielsen
-
-This program is free software: you can redistribute it and/or modify it under
-the terms of the GNU Lesser General Public License as published by the Free
-Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
-This program is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
-details.
-
-You should have received a copy of the GNU Lesser General Public License along
-with this program. If not, see <http://www.gnu.org/licenses/>.
-*/
-
-/*
-This kd-tree implementation is based on the scipy.spatial.cKDTree by
-Anne M. Archibald and libANN by David M. Mount and Sunil Arya.
-*/
-
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <stdint.h>
-#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; }
-
-#ifdef _MSC_VER
-#define restrict __restrict
-#endif
-
-
-typedef struct
-{
- float cut_val;
- int8_t cut_dim;
- uint32_t start_idx;
- uint32_t n;
- float cut_bounds_lv;
- float cut_bounds_hv;
- struct Node_float *left_child;
- struct Node_float *right_child;
-} Node_float;
-
-typedef struct
-{
- float *bbox;
- int8_t no_dims;
- uint32_t *pidx;
- struct Node_float *root;
-} Tree_float;
-
-
-typedef struct
-{
- double cut_val;
- int8_t cut_dim;
- uint32_t start_idx;
- uint32_t n;
- double cut_bounds_lv;
- double cut_bounds_hv;
- struct Node_double *left_child;
- struct Node_double *right_child;
-} Node_double;
-
-typedef struct
-{
- double *bbox;
- int8_t no_dims;
- uint32_t *pidx;
- struct Node_double *root;
-} Tree_double;
-
-
-
-void insert_point_float(uint32_t *closest_idx, float *closest_dist, uint32_t pidx, float cur_dist, uint32_t k);
-void get_bounding_box_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, float *bbox);
-int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim,
- float *cut_val, uint32_t *n_lo);
-Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
-Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, float *bbox);
-Node_float * create_node_float(uint32_t start_idx, uint32_t n, int is_leaf);
-void delete_subtree_float(Node_float *root);
-void delete_tree_float(Tree_float *tree);
-void print_tree_float(Node_float *root, int level);
-float calc_dist_float(float *point1_coord, float *point2_coord, int8_t no_dims);
-float get_cube_offset_float(int8_t dim, float *point_coord, float *bbox);
-float get_min_dist_float(float *point_coord, int8_t no_dims, float *bbox);
-void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
- uint32_t k, uint32_t *restrict closest_idx, float *restrict closest_dist);
-void search_leaf_float_mask(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
- uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, float *restrict closest_dist);
-void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord,
- float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint32_t * closest_idx, float *closest_dist);
-void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
- uint32_t num_points, uint32_t k, float distance_upper_bound,
- float eps, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists);
-
-
-void insert_point_double(uint32_t *closest_idx, double *closest_dist, uint32_t pidx, double cur_dist, uint32_t k);
-void get_bounding_box_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, double *bbox);
-int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim,
- double *cut_val, uint32_t *n_lo);
-Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
-Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, double *bbox);
-Node_double * create_node_double(uint32_t start_idx, uint32_t n, int is_leaf);
-void delete_subtree_double(Node_double *root);
-void delete_tree_double(Tree_double *tree);
-void print_tree_double(Node_double *root, int level);
-double calc_dist_double(double *point1_coord, double *point2_coord, int8_t no_dims);
-double get_cube_offset_double(int8_t dim, double *point_coord, double *bbox);
-double get_min_dist_double(double *point_coord, int8_t no_dims, double *bbox);
-void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
- uint32_t k, uint32_t *restrict closest_idx, double *restrict closest_dist);
-void search_leaf_double_mask(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
- uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, double *restrict closest_dist);
-void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord,
- double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint32_t * closest_idx, double *closest_dist);
-void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
- uint32_t num_points, uint32_t k, double distance_upper_bound,
- double eps, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists);
-
-
-
-/************************************************
-Insert point into priority queue
-Params:
- closest_idx : index queue
- closest_dist : distance queue
- pidx : permutation index of data points
- cur_dist : distance to point inserted
- k : number of neighbours
-************************************************/
-void insert_point_float(uint32_t *closest_idx, float *closest_dist, uint32_t pidx, float cur_dist, uint32_t k)
-{
- int i;
- for (i = k - 1; i > 0; i--)
- {
- if (closest_dist[i - 1] > cur_dist)
- {
- closest_dist[i] = closest_dist[i - 1];
- closest_idx[i] = closest_idx[i - 1];
- }
- else
- {
- break;
- }
- }
- closest_idx[i] = pidx;
- closest_dist[i] = cur_dist;
-}
-
-/************************************************
-Get the bounding box of a set of points
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims: number of dimensions
- n : number of points
- bbox : bounding box (return)
-************************************************/
-void get_bounding_box_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, float *bbox)
-{
- float cur;
- int8_t bbox_idx, i, j;
- uint32_t i2;
-
- /* Use first data point to initialize */
- for (i = 0; i < no_dims; i++)
- {
- bbox[2 * i] = bbox[2 * i + 1] = PA(0, i);
- }
-
- /* Update using rest of data points */
- for (i2 = 1; i2 < n; i2++)
- {
- for (j = 0; j < no_dims; j++)
- {
- bbox_idx = 2 * j;
- cur = PA(i2, j);
- if (cur < bbox[bbox_idx])
- {
- bbox[bbox_idx] = cur;
- }
- else if (cur > bbox[bbox_idx + 1])
- {
- bbox[bbox_idx + 1] = cur;
- }
- }
- }
-}
-
-/************************************************
-Partition a range of data points by manipulation the permutation index.
-The sliding midpoint rule is used for the partitioning.
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims: number of dimensions
- start_idx : index of first data point to use
- n : number of data points
- bbox : bounding box of data points
- cut_dim : dimension used for partition (return)
- cut_val : value of cutting point (return)
- n_lo : number of point below cutting plane (return)
-************************************************/
-int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim, float *cut_val, uint32_t *n_lo)
-{
- int8_t dim = 0, i;
- uint32_t p, q, i2;
- float size = 0, min_val, max_val, split, side_len, cur_val;
- uint32_t end_idx = start_idx + n - 1;
-
- /* Find largest bounding box side */
- for (i = 0; i < no_dims; i++)
- {
- side_len = bbox[2 * i + 1] - bbox[2 * i];
- if (side_len > size)
- {
- dim = i;
- size = side_len;
- }
- }
-
- min_val = bbox[2 * dim];
- max_val = bbox[2 * dim + 1];
-
- /* Check for zero length or inconsistent */
- if (min_val >= max_val)
- return 1;
-
- /* Use middle for splitting */
- split = (min_val + max_val) / 2;
-
- /* Partition all data points around middle */
- p = start_idx;
- q = end_idx;
- while (p <= q)
- {
- if (PA(p, dim) < split)
- {
- p++;
- }
- else if (PA(q, dim) >= split)
- {
- /* Guard for underflow */
- if (q > 0)
- {
- q--;
- }
- else
- {
- break;
- }
- }
- else
- {
- PASWAP(p, q);
- p++;
- q--;
- }
- }
-
- /* Check for empty splits */
- if (p == start_idx)
- {
- /* No points less than split.
- Split at lowest point instead.
- Minimum 1 point will be in lower box.
- */
-
- uint32_t j = start_idx;
- split = PA(j, dim);
- for (i2 = start_idx + 1; i2 <= end_idx; i2++)
- {
- /* Find lowest point */
- cur_val = PA(i2, dim);
- if (cur_val < split)
- {
- j = i2;
- split = cur_val;
- }
- }
- PASWAP(j, start_idx);
- p = start_idx + 1;
- }
- else if (p == end_idx + 1)
- {
- /* No points greater than split.
- Split at highest point instead.
- Minimum 1 point will be in higher box.
- */
-
- uint32_t j = end_idx;
- split = PA(j, dim);
- for (i2 = start_idx; i2 < end_idx; i2++)
- {
- /* Find highest point */
- cur_val = PA(i2, dim);
- if (cur_val > split)
- {
- j = i2;
- split = cur_val;
- }
- }
- PASWAP(j, end_idx);
- p = end_idx;
- }
-
- /* Set return values */
- *cut_dim = dim;
- *cut_val = split;
- *n_lo = p - start_idx;
- return 0;
-}
-
-/************************************************
-Construct a sub tree over a range of data points.
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims: number of dimensions
- start_idx : index of first data point to use
- n : number of data points
- bsp : number of points per leaf
- bbox : bounding box of set of data points
-************************************************/
-Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, float *bbox)
-{
- /* Create new node */
- int is_leaf = (n <= bsp);
- Node_float *root = create_node_float(start_idx, n, is_leaf);
- int rval;
- int8_t cut_dim;
- uint32_t n_lo;
- float cut_val, lv, hv;
- if (is_leaf)
- {
- /* Make leaf node */
- root->cut_dim = -1;
- }
- else
- {
- /* Make split node */
- /* Partition data set and set node info */
- rval = partition_float(pa, pidx, no_dims, start_idx, n, bbox, &cut_dim, &cut_val, &n_lo);
- if (rval == 1)
- {
- root->cut_dim = -1;
- return root;
- }
- root->cut_val = cut_val;
- root->cut_dim = cut_dim;
-
- /* Recurse on both subsets */
- lv = bbox[2 * cut_dim];
- hv = bbox[2 * cut_dim + 1];
-
- /* Set bounds for cut dimension */
- root->cut_bounds_lv = lv;
- root->cut_bounds_hv = hv;
-
- /* Update bounding box before call to lower subset and restore after */
- bbox[2 * cut_dim + 1] = cut_val;
- root->left_child = (struct Node_float *)construct_subtree_float(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_float *)construct_subtree_float(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
- bbox[2 * cut_dim] = lv;
- }
- return root;
-}
-
-/************************************************
-Construct a tree over data points.
-Params:
- pa : data points
- no_dims: number of dimensions
- n : number of data points
- bsp : number of points per leaf
-************************************************/
-Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp)
-{
- Tree_float *tree = (Tree_float *)malloc(sizeof(Tree_float));
- uint32_t i;
- uint32_t *pidx;
- float *bbox;
-
- tree->no_dims = no_dims;
-
- /* Initialize permutation array */
- pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
- for (i = 0; i < n; i++)
- {
- pidx[i] = i;
- }
-
- bbox = (float *)malloc(2 * sizeof(float) * no_dims);
- get_bounding_box_float(pa, pidx, no_dims, n, bbox);
- tree->bbox = bbox;
-
- /* Construct subtree on full dataset */
- tree->root = (struct Node_float *)construct_subtree_float(pa, pidx, no_dims, 0, n, bsp, bbox);
-
- tree->pidx = pidx;
- return tree;
-}
-
-/************************************************
-Create a tree node.
-Params:
- start_idx : index of first data point to use
- n : number of data points
-************************************************/
-Node_float* create_node_float(uint32_t start_idx, uint32_t n, int is_leaf)
-{
- Node_float *new_node;
- if (is_leaf)
- {
- /*
- Allocate only the part of the struct that will be used in a leaf node.
- 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_float *)malloc(sizeof(Node_float) - 2 * sizeof(Node_float *));
- }
- else
- {
- new_node = (Node_float *)malloc(sizeof(Node_float));
- }
- new_node->n = n;
- new_node->start_idx = start_idx;
- return new_node;
-}
-
-/************************************************
-Delete subtree
-Params:
- root : root node of subtree to delete
-************************************************/
-void delete_subtree_float(Node_float *root)
-{
- if (root->cut_dim != -1)
- {
- delete_subtree_float((Node_float *)root->left_child);
- delete_subtree_float((Node_float *)root->right_child);
- }
- free(root);
-}
-
-/************************************************
-Delete tree
-Params:
- tree : Tree struct of kd tree
-************************************************/
-void delete_tree_float(Tree_float *tree)
-{
- delete_subtree_float((Node_float *)tree->root);
- free(tree->bbox);
- free(tree->pidx);
- free(tree);
-}
-
-/************************************************
-Print
-************************************************/
-void print_tree_float(Node_float *root, int level)
-{
- int i;
- for (i = 0; i < level; i++)
- {
- printf(" ");
- }
- printf("(cut_val: %f, cut_dim: %i)\n", root->cut_val, root->cut_dim);
- if (root->cut_dim != -1)
- print_tree_float((Node_float *)root->left_child, level + 1);
- if (root->cut_dim != -1)
- print_tree_float((Node_float *)root->right_child, level + 1);
-}
-
-/************************************************
-Calculate squared cartesian distance between points
-Params:
- point1_coord : point 1
- point2_coord : point 2
-************************************************/
-float calc_dist_float(float *point1_coord, float *point2_coord, int8_t no_dims)
-{
- /* Calculate squared distance */
- float 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
-************************************************/
-float get_cube_offset_float(int8_t dim, float *point_coord, float *bbox)
-{
- float 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
-************************************************/
-float get_min_dist_float(float *point_coord, int8_t no_dims, float *bbox)
-{
- float cube_offset = 0, cube_offset_dim;
- int8_t i;
-
- for (i = 0; i < no_dims; i++)
- {
- cube_offset_dim = get_cube_offset_float(i, point_coord, bbox);
- cube_offset += cube_offset_dim * cube_offset_dim;
- }
-
- return cube_offset;
-}
-
-/************************************************
-Search a leaf node for closest point
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims : number of dimensions
- start_idx : index of first data point to use
- size : number of data points
- point_coord : query point
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
- uint32_t k, uint32_t *restrict closest_idx, float *restrict closest_dist)
-{
- float cur_dist;
- uint32_t i;
- /* Loop through all points in leaf */
- for (i = 0; i < n; i++)
- {
- /* Get distance to query point */
- cur_dist = calc_dist_float(&PA(start_idx + i, 0), point_coord, no_dims);
- /* Update closest info if new point is closest so far*/
- if (cur_dist < closest_dist[k - 1])
- {
- insert_point_float(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
- }
- }
-}
-
-
-/************************************************
-Search a leaf node for closest point with data point mask
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims : number of dimensions
- start_idx : index of first data point to use
- size : number of data points
- point_coord : query point
- mask : boolean array of invalid (True) and valid (False) data points
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_leaf_float_mask(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
- uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, float *restrict closest_dist)
-{
- float cur_dist;
- uint32_t i;
- /* Loop through all points in leaf */
- for (i = 0; i < n; i++)
- {
- /* Is this point masked out? */
- if (mask[pidx[start_idx + i]])
- {
- continue;
- }
- /* Get distance to query point */
- cur_dist = calc_dist_float(&PA(start_idx + i, 0), point_coord, no_dims);
- /* Update closest info if new point is closest so far*/
- if (cur_dist < closest_dist[k - 1])
- {
- insert_point_float(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
- }
- }
-}
-
-/************************************************
-Search subtree for nearest to query point
-Params:
- root : root node of subtree
- pa : data points
- pidx : permutation index of data points
- no_dims : number of dimensions
- point_coord : query point
- min_dist : minumum distance to nearest neighbour
- mask : boolean array of invalid (True) and valid (False) data points
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord,
- float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask,
- uint32_t *closest_idx, float *closest_dist)
-{
- int8_t dim;
- float dist_left, dist_right;
- float new_offset;
- float box_diff;
-
- /* Skip if distance bound exeeded */
- if (min_dist > distance_upper_bound)
- {
- return;
- }
-
- dim = root->cut_dim;
-
- /* Handle leaf node */
- if (dim == -1)
- {
- if (mask)
- {
- search_leaf_float_mask(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, mask, closest_idx, closest_dist);
- }
- else
- {
- search_leaf_float(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
- }
- return;
- }
-
- /* Get distance to cutting plane */
- new_offset = point_coord[dim] - root->cut_val;
-
- if (new_offset < 0)
- {
- /* Left of cutting plane */
- dist_left = min_dist;
- if (dist_left < closest_dist[k - 1] * eps_fac)
- {
- /* Search left subtree if minimum distance is below limit */
- search_splitnode_float((Node_float *)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.
- See Algorithms for Fast Vector Quantization
- Sunil Arya and David M. Mount. */
- box_diff = root->cut_bounds_lv - point_coord[dim];
- if (box_diff < 0)
- {
- box_diff = 0;
- }
- dist_right = min_dist - box_diff * box_diff + new_offset * new_offset;
- if (dist_right < closest_dist[k - 1] * eps_fac)
- {
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_float((Node_float *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
- }
- }
- else
- {
- /* Right of cutting plane */
- dist_right = min_dist;
- if (dist_right < closest_dist[k - 1] * eps_fac)
- {
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_float((Node_float *)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.
- See Algorithms for Fast Vector Quantization
- Sunil Arya and David M. Mount. */
- box_diff = point_coord[dim] - root->cut_bounds_hv;
- if (box_diff < 0)
- {
- box_diff = 0;
- }
- dist_left = min_dist - box_diff * box_diff + new_offset * new_offset;
- if (dist_left < closest_dist[k - 1] * eps_fac)
- {
- /* Search left subtree if minimum distance is below limit*/
- search_splitnode_float((Node_float *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
- }
- }
-}
-
-/************************************************
-Search for nearest neighbour for a set of query points
-Params:
- tree : Tree struct of kd tree
- pa : data points
- pidx : permutation index of data points
- point_coords : query points
- num_points : number of query points
- mask : boolean array of invalid (True) and valid (False) data points
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
- uint32_t num_points, uint32_t k, float distance_upper_bound,
- float eps, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists)
-{
- float min_dist;
- float eps_fac = 1 / ((1 + eps) * (1 + eps));
- int8_t no_dims = tree->no_dims;
- float *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_float *root = (Node_float *)tree->root;
-
- /* Queries are OpenMP enabled */
- #pragma omp parallel
- {
- /* The low chunk size is important to avoid L2 cache trashing
- for spatial coherent query datasets
- */
- #pragma omp for private(i, j) schedule(static, 100) nowait
- for (i = 0; i < local_num_points; i++)
- {
- for (j = 0; j < k; j++)
- {
- closest_idxs[i * k + j] = UINT32_MAX;
- closest_dists[i * k + j] = DBL_MAX;
- }
- min_dist = get_min_dist_float(point_coords + no_dims * i, no_dims, bbox);
- search_splitnode_float(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]);
- }
- }
-}
-
-/************************************************
-Insert point into priority queue
-Params:
- closest_idx : index queue
- closest_dist : distance queue
- pidx : permutation index of data points
- cur_dist : distance to point inserted
- k : number of neighbours
-************************************************/
-void insert_point_double(uint32_t *closest_idx, double *closest_dist, uint32_t pidx, double cur_dist, uint32_t k)
-{
- int i;
- for (i = k - 1; i > 0; i--)
- {
- if (closest_dist[i - 1] > cur_dist)
- {
- closest_dist[i] = closest_dist[i - 1];
- closest_idx[i] = closest_idx[i - 1];
- }
- else
- {
- break;
- }
- }
- closest_idx[i] = pidx;
- closest_dist[i] = cur_dist;
-}
-
-/************************************************
-Get the bounding box of a set of points
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims: number of dimensions
- n : number of points
- bbox : bounding box (return)
-************************************************/
-void get_bounding_box_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, double *bbox)
-{
- double cur;
- int8_t bbox_idx, i, j;
- uint32_t i2;
-
- /* Use first data point to initialize */
- for (i = 0; i < no_dims; i++)
- {
- bbox[2 * i] = bbox[2 * i + 1] = PA(0, i);
- }
-
- /* Update using rest of data points */
- for (i2 = 1; i2 < n; i2++)
- {
- for (j = 0; j < no_dims; j++)
- {
- bbox_idx = 2 * j;
- cur = PA(i2, j);
- if (cur < bbox[bbox_idx])
- {
- bbox[bbox_idx] = cur;
- }
- else if (cur > bbox[bbox_idx + 1])
- {
- bbox[bbox_idx + 1] = cur;
- }
- }
- }
-}
-
-/************************************************
-Partition a range of data points by manipulation the permutation index.
-The sliding midpoint rule is used for the partitioning.
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims: number of dimensions
- start_idx : index of first data point to use
- n : number of data points
- bbox : bounding box of data points
- cut_dim : dimension used for partition (return)
- cut_val : value of cutting point (return)
- n_lo : number of point below cutting plane (return)
-************************************************/
-int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim, double *cut_val, uint32_t *n_lo)
-{
- int8_t dim = 0, i;
- uint32_t p, q, i2;
- double size = 0, min_val, max_val, split, side_len, cur_val;
- uint32_t end_idx = start_idx + n - 1;
-
- /* Find largest bounding box side */
- for (i = 0; i < no_dims; i++)
- {
- side_len = bbox[2 * i + 1] - bbox[2 * i];
- if (side_len > size)
- {
- dim = i;
- size = side_len;
- }
- }
-
- min_val = bbox[2 * dim];
- max_val = bbox[2 * dim + 1];
-
- /* Check for zero length or inconsistent */
- if (min_val >= max_val)
- return 1;
-
- /* Use middle for splitting */
- split = (min_val + max_val) / 2;
-
- /* Partition all data points around middle */
- p = start_idx;
- q = end_idx;
- while (p <= q)
- {
- if (PA(p, dim) < split)
- {
- p++;
- }
- else if (PA(q, dim) >= split)
- {
- /* Guard for underflow */
- if (q > 0)
- {
- q--;
- }
- else
- {
- break;
- }
- }
- else
- {
- PASWAP(p, q);
- p++;
- q--;
- }
- }
-
- /* Check for empty splits */
- if (p == start_idx)
- {
- /* No points less than split.
- Split at lowest point instead.
- Minimum 1 point will be in lower box.
- */
-
- uint32_t j = start_idx;
- split = PA(j, dim);
- for (i2 = start_idx + 1; i2 <= end_idx; i2++)
- {
- /* Find lowest point */
- cur_val = PA(i2, dim);
- if (cur_val < split)
- {
- j = i2;
- split = cur_val;
- }
- }
- PASWAP(j, start_idx);
- p = start_idx + 1;
- }
- else if (p == end_idx + 1)
- {
- /* No points greater than split.
- Split at highest point instead.
- Minimum 1 point will be in higher box.
- */
-
- uint32_t j = end_idx;
- split = PA(j, dim);
- for (i2 = start_idx; i2 < end_idx; i2++)
- {
- /* Find highest point */
- cur_val = PA(i2, dim);
- if (cur_val > split)
- {
- j = i2;
- split = cur_val;
- }
- }
- PASWAP(j, end_idx);
- p = end_idx;
- }
-
- /* Set return values */
- *cut_dim = dim;
- *cut_val = split;
- *n_lo = p - start_idx;
- return 0;
-}
-
-/************************************************
-Construct a sub tree over a range of data points.
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims: number of dimensions
- start_idx : index of first data point to use
- n : number of data points
- bsp : number of points per leaf
- bbox : bounding box of set of data points
-************************************************/
-Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, double *bbox)
-{
- /* Create new node */
- int is_leaf = (n <= bsp);
- Node_double *root = create_node_double(start_idx, n, is_leaf);
- int rval;
- int8_t cut_dim;
- uint32_t n_lo;
- double cut_val, lv, hv;
- if (is_leaf)
- {
- /* Make leaf node */
- root->cut_dim = -1;
- }
- else
- {
- /* Make split node */
- /* Partition data set and set node info */
- rval = partition_double(pa, pidx, no_dims, start_idx, n, bbox, &cut_dim, &cut_val, &n_lo);
- if (rval == 1)
- {
- root->cut_dim = -1;
- return root;
- }
- root->cut_val = cut_val;
- root->cut_dim = cut_dim;
-
- /* Recurse on both subsets */
- lv = bbox[2 * cut_dim];
- hv = bbox[2 * cut_dim + 1];
-
- /* Set bounds for cut dimension */
- root->cut_bounds_lv = lv;
- root->cut_bounds_hv = hv;
-
- /* Update bounding box before call to lower subset and restore after */
- bbox[2 * cut_dim + 1] = cut_val;
- root->left_child = (struct Node_double *)construct_subtree_double(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_double *)construct_subtree_double(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
- bbox[2 * cut_dim] = lv;
- }
- return root;
-}
-
-/************************************************
-Construct a tree over data points.
-Params:
- pa : data points
- no_dims: number of dimensions
- n : number of data points
- bsp : number of points per leaf
-************************************************/
-Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp)
-{
- Tree_double *tree = (Tree_double *)malloc(sizeof(Tree_double));
- uint32_t i;
- uint32_t *pidx;
- double *bbox;
-
- tree->no_dims = no_dims;
-
- /* Initialize permutation array */
- pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
- for (i = 0; i < n; i++)
- {
- pidx[i] = i;
- }
-
- bbox = (double *)malloc(2 * sizeof(double) * no_dims);
- get_bounding_box_double(pa, pidx, no_dims, n, bbox);
- tree->bbox = bbox;
-
- /* Construct subtree on full dataset */
- tree->root = (struct Node_double *)construct_subtree_double(pa, pidx, no_dims, 0, n, bsp, bbox);
-
- tree->pidx = pidx;
- return tree;
-}
-
-/************************************************
-Create a tree node.
-Params:
- start_idx : index of first data point to use
- n : number of data points
-************************************************/
-Node_double* create_node_double(uint32_t start_idx, uint32_t n, int is_leaf)
-{
- Node_double *new_node;
- if (is_leaf)
- {
- /*
- Allocate only the part of the struct that will be used in a leaf node.
- 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_double *)malloc(sizeof(Node_double) - 2 * sizeof(Node_double *));
- }
- else
- {
- new_node = (Node_double *)malloc(sizeof(Node_double));
- }
- new_node->n = n;
- new_node->start_idx = start_idx;
- return new_node;
-}
-
-/************************************************
-Delete subtree
-Params:
- root : root node of subtree to delete
-************************************************/
-void delete_subtree_double(Node_double *root)
-{
- if (root->cut_dim != -1)
- {
- delete_subtree_double((Node_double *)root->left_child);
- delete_subtree_double((Node_double *)root->right_child);
- }
- free(root);
-}
-
-/************************************************
-Delete tree
-Params:
- tree : Tree struct of kd tree
-************************************************/
-void delete_tree_double(Tree_double *tree)
-{
- delete_subtree_double((Node_double *)tree->root);
- free(tree->bbox);
- free(tree->pidx);
- free(tree);
-}
-
-/************************************************
-Print
-************************************************/
-void print_tree_double(Node_double *root, int level)
-{
- int i;
- for (i = 0; i < level; i++)
- {
- printf(" ");
- }
- printf("(cut_val: %f, cut_dim: %i)\n", root->cut_val, root->cut_dim);
- if (root->cut_dim != -1)
- print_tree_double((Node_double *)root->left_child, level + 1);
- if (root->cut_dim != -1)
- print_tree_double((Node_double *)root->right_child, level + 1);
-}
-
-/************************************************
-Calculate squared cartesian distance between points
-Params:
- point1_coord : point 1
- point2_coord : point 2
-************************************************/
-double calc_dist_double(double *point1_coord, double *point2_coord, int8_t no_dims)
-{
- /* Calculate squared distance */
- double 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
-************************************************/
-double get_cube_offset_double(int8_t dim, double *point_coord, double *bbox)
-{
- double 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
-************************************************/
-double get_min_dist_double(double *point_coord, int8_t no_dims, double *bbox)
-{
- double cube_offset = 0, cube_offset_dim;
- int8_t i;
-
- for (i = 0; i < no_dims; i++)
- {
- cube_offset_dim = get_cube_offset_double(i, point_coord, bbox);
- cube_offset += cube_offset_dim * cube_offset_dim;
- }
-
- return cube_offset;
-}
-
-/************************************************
-Search a leaf node for closest point
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims : number of dimensions
- start_idx : index of first data point to use
- size : number of data points
- point_coord : query point
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
- uint32_t k, uint32_t *restrict closest_idx, double *restrict closest_dist)
-{
- double cur_dist;
- uint32_t i;
- /* Loop through all points in leaf */
- for (i = 0; i < n; i++)
- {
- /* Get distance to query point */
- cur_dist = calc_dist_double(&PA(start_idx + i, 0), point_coord, no_dims);
- /* Update closest info if new point is closest so far*/
- if (cur_dist < closest_dist[k - 1])
- {
- insert_point_double(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
- }
- }
-}
-
-
-/************************************************
-Search a leaf node for closest point with data point mask
-Params:
- pa : data points
- pidx : permutation index of data points
- no_dims : number of dimensions
- start_idx : index of first data point to use
- size : number of data points
- point_coord : query point
- mask : boolean array of invalid (True) and valid (False) data points
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_leaf_double_mask(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
- uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, double *restrict closest_dist)
-{
- double cur_dist;
- uint32_t i;
- /* Loop through all points in leaf */
- for (i = 0; i < n; i++)
- {
- /* Is this point masked out? */
- if (mask[pidx[start_idx + i]])
- {
- continue;
- }
- /* Get distance to query point */
- cur_dist = calc_dist_double(&PA(start_idx + i, 0), point_coord, no_dims);
- /* Update closest info if new point is closest so far*/
- if (cur_dist < closest_dist[k - 1])
- {
- insert_point_double(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
- }
- }
-}
-
-/************************************************
-Search subtree for nearest to query point
-Params:
- root : root node of subtree
- pa : data points
- pidx : permutation index of data points
- no_dims : number of dimensions
- point_coord : query point
- min_dist : minumum distance to nearest neighbour
- mask : boolean array of invalid (True) and valid (False) data points
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord,
- double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask,
- uint32_t *closest_idx, double *closest_dist)
-{
- int8_t dim;
- double dist_left, dist_right;
- double new_offset;
- double box_diff;
-
- /* Skip if distance bound exeeded */
- if (min_dist > distance_upper_bound)
- {
- return;
- }
-
- dim = root->cut_dim;
-
- /* Handle leaf node */
- if (dim == -1)
- {
- if (mask)
- {
- search_leaf_double_mask(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, mask, closest_idx, closest_dist);
- }
- else
- {
- search_leaf_double(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
- }
- return;
- }
-
- /* Get distance to cutting plane */
- new_offset = point_coord[dim] - root->cut_val;
-
- if (new_offset < 0)
- {
- /* Left of cutting plane */
- dist_left = min_dist;
- if (dist_left < closest_dist[k - 1] * eps_fac)
- {
- /* Search left subtree if minimum distance is below limit */
- search_splitnode_double((Node_double *)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.
- See Algorithms for Fast Vector Quantization
- Sunil Arya and David M. Mount. */
- box_diff = root->cut_bounds_lv - point_coord[dim];
- if (box_diff < 0)
- {
- box_diff = 0;
- }
- dist_right = min_dist - box_diff * box_diff + new_offset * new_offset;
- if (dist_right < closest_dist[k - 1] * eps_fac)
- {
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_double((Node_double *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
- }
- }
- else
- {
- /* Right of cutting plane */
- dist_right = min_dist;
- if (dist_right < closest_dist[k - 1] * eps_fac)
- {
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_double((Node_double *)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.
- See Algorithms for Fast Vector Quantization
- Sunil Arya and David M. Mount. */
- box_diff = point_coord[dim] - root->cut_bounds_hv;
- if (box_diff < 0)
- {
- box_diff = 0;
- }
- dist_left = min_dist - box_diff * box_diff + new_offset * new_offset;
- if (dist_left < closest_dist[k - 1] * eps_fac)
- {
- /* Search left subtree if minimum distance is below limit*/
- search_splitnode_double((Node_double *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
- }
- }
-}
-
-/************************************************
-Search for nearest neighbour for a set of query points
-Params:
- tree : Tree struct of kd tree
- pa : data points
- pidx : permutation index of data points
- point_coords : query points
- num_points : number of query points
- mask : boolean array of invalid (True) and valid (False) data points
- closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
-************************************************/
-void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
- uint32_t num_points, uint32_t k, double distance_upper_bound,
- double eps, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists)
-{
- double min_dist;
- double eps_fac = 1 / ((1 + eps) * (1 + eps));
- int8_t no_dims = tree->no_dims;
- double *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_double *root = (Node_double *)tree->root;
-
- /* Queries are OpenMP enabled */
- #pragma omp parallel
- {
- /* The low chunk size is important to avoid L2 cache trashing
- for spatial coherent query datasets
- */
- #pragma omp for private(i, j) schedule(static, 100) nowait
- for (i = 0; i < local_num_points; i++)
- {
- for (j = 0; j < k; j++)
- {
- closest_idxs[i * k + j] = UINT32_MAX;
- closest_dists[i * k + j] = DBL_MAX;
- }
- min_dist = get_min_dist_double(point_coords + no_dims * i, no_dims, bbox);
- search_splitnode_double(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]);
- }
- }
-}
=====================================
pykdtree/_kdtree_core.c.mako
=====================================
@@ -128,8 +128,8 @@ Params:
void get_bounding_box_${DTYPE}(${DTYPE} *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, ${DTYPE} *bbox)
{
${DTYPE} cur;
- int8_t bbox_idx, i, j;
- uint32_t i2;
+ int8_t i, j;
+ uint32_t bbox_idx, i2;
/* Use first data point to initialize */
for (i = 0; i < no_dims; i++)
=====================================
pykdtree/kdtree.c deleted
=====================================
The diff for this file was not included because it is too large.
=====================================
pykdtree/kdtree.pyx
=====================================
@@ -115,6 +115,8 @@ cdef class KDTree:
self.leafsize = <uint32_t>leafsize
if data_pts.ndim == 1:
self.ndim = 1
+ elif data_pts.shape[1] > 127:
+ raise ValueError('Max 127 dimensions allowed')
else:
self.ndim = <int8_t>data_pts.shape[1]
@@ -135,7 +137,7 @@ cdef class KDTree:
:Parameters:
query_pts : numpy array
- Query points with shape (n , dims)
+ Query points with shape (m, dims)
k : int
The number of nearest neighbours to return
eps : non-negative float
@@ -151,8 +153,8 @@ cdef class KDTree:
mask : numpy array, optional
Array of booleans where neighbors are considered invalid and
should not be returned. A mask value of True represents an
- invalid pixel. Mask should have shape (n,). By default all
- points are considered valid.
+ invalid pixel. Mask should have shape (n,) to match data points.
+ By default all points are considered valid.
"""
=====================================
pykdtree/test_tree.py
=====================================
@@ -351,3 +351,22 @@ def test3d_mask():
assert abs(dist[0] - 66759196.1053) < epsilon * dist[0]
assert abs(dist[1] - 3.) < epsilon * dist[1]
assert abs(dist[2] - 20001.) < epsilon * dist[2]
+
+def test128d_fail():
+ pts = 100
+ dims = 128
+ data_pts = np.arange(pts * dims).reshape(pts, dims)
+ try:
+ kdtree = KDTree(data_pts)
+ except ValueError as exc:
+ assert "Max 127 dimensions" in str(exc)
+ else:
+ raise Exception("Should not accept 129 dimensional data")
+
+def test127d_ok():
+ pts = 2
+ dims = 127
+ data_pts = np.arange(pts * dims).reshape(pts, dims)
+ kdtree = KDTree(data_pts)
+ dist, idx = kdtree.query(data_pts)
+ assert np.all(dist == 0)
=====================================
scripts/build-manylinux-wheels.sh
=====================================
@@ -0,0 +1,44 @@
+#!/bin/bash
+set -e -x
+
+# This is to be run by Docker inside a Docker image.
+# You can test it locally on a Linux machine by installing docker and running from this repo's root:
+# $ docker run -e PLAT=manylinux1_x86_64 -v `pwd`:/io quay.io/pypa/manylinux1_x86_64 /io/scripts/build-manylinux-wheels.sh
+
+# * The -e just defines an environment variable PLAT=[docker name] inside the
+# docker - auditwheel can't detect the docker name automatically.
+# * The -v gives a directory alias for passing files in and out of the docker
+# (/io is arbitrary). E.g the `setup.py` script would be accessed in the
+# docker via `/io/setup.py`.
+# * quay.io/pypa/manylinux1_x86_64 is the full docker image name. Docker
+# downloads it automatically.
+# * The last argument is a shell command that the Docker will execute.
+# Filenames must be from the Docker's perspective.
+
+# Wheels are initially generated as you would usually, but put in a temp
+# directory temp-wheels. The pip-cache is optional but can speed up local builds
+# having a real permanent pip-cache dir.
+mkdir -p /io/pip-cache
+mkdir -p /io/temp-wheels
+
+# Clean out any old existing wheels.
+find /io/temp-wheels/ -type f -delete
+
+# Iterate through available pythons.
+for PYBIN in /opt/python/cp3[5678]*/bin; do
+ "${PYBIN}/pip" install -q -U setuptools wheel nose --cache-dir /io/pip-cache
+ # Run the following in root of this repo.
+ (cd /io/ && USE_OMP=$USE_OMP "${PYBIN}/pip" install -q .)
+ (cd /io/ && USE_OMP=$USE_OMP "${PYBIN}/python" setup.py nosetests)
+ (cd /io/ && USE_OMP=$USE_OMP "${PYBIN}/python" setup.py -q bdist_wheel -d /io/temp-wheels)
+done
+
+"$PYBIN/pip" install -q auditwheel
+
+# Wheels aren't considered manylinux unless they have been through
+# auditwheel. Audited wheels go in /io/dist/.
+mkdir -p /io/dist/
+
+for whl in /io/temp-wheels/*.whl; do
+ auditwheel repair "$whl" --plat "$PLAT" -w /io/dist/
+done
=====================================
setup.cfg
=====================================
@@ -1,5 +1,5 @@
[bdist_rpm]
-requires=numpy
+requires=python3-numpy
release=1
=====================================
setup.py
=====================================
@@ -16,14 +16,33 @@
# with this program. If not, see <http://www.gnu.org/licenses/>.
import os
+import sys
from setuptools import setup, Extension
from setuptools.command.build_ext import build_ext
+
+def is_conda_interpreter():
+ """Is the running interpreter from Anaconda or miniconda?
+
+ See https://stackoverflow.com/a/21318941/433202
+
+ Examples::
+
+ 2.7.6 |Anaconda 1.8.0 (x86_64)| (default, Jan 10 2014, 11:23:15)
+ 2.7.6 |Continuum Analytics, Inc.| (default, Jan 10 2014, 11:23:15)
+ 3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 09:55:02)
+
+ """
+ return 'conda' in sys.version or 'Continuum' in sys.version
+
+
# Get OpenMP setting from environment
try:
use_omp = int(os.environ['USE_OMP'])
except KeyError:
- use_omp = True
+ # OpenMP is not supported with default clang
+ # Conda provides its own compiler which does support openmp
+ use_omp = 'darwin' not in sys.platform or is_conda_interpreter()
def set_builtin(name, value):
@@ -70,26 +89,31 @@ class build_ext_subclass(build_ext):
import numpy
self.include_dirs.append(numpy.get_include())
+with open('README.rst', 'r') as readme_file:
+ readme = readme_file.read()
setup(
name='pykdtree',
- version='1.3.1',
+ version='1.3.3',
+ url="https://github.com/storpipfugl/pykdtree",
description='Fast kd-tree implementation with OpenMP-enabled queries',
+ long_description=readme,
author='Esben S. Nielsen',
author_email='storpipfugl at gmail.com',
- packages = ['pykdtree'],
+ packages=['pykdtree'],
python_requires='>=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*',
install_requires=['numpy'],
setup_requires=['numpy'],
tests_require=['nose'],
zip_safe=False,
- test_suite = 'nose.collector',
- ext_modules = [Extension('pykdtree.kdtree',
- ['pykdtree/kdtree.c', 'pykdtree/_kdtree_core.c'])],
- cmdclass = {'build_ext': build_ext_subclass },
+ test_suite='nose.collector',
+ ext_modules=[Extension('pykdtree.kdtree',
+ ['pykdtree/kdtree.c', 'pykdtree/_kdtree_core.c'])],
+ cmdclass={'build_ext': build_ext_subclass},
classifiers=[
'Development Status :: 5 - Production/Stable',
- 'License :: OSI Approved :: GNU General Public License v3 (GPLv3)',
+ ('License :: OSI Approved :: '
+ 'GNU Lesser General Public License v3 (LGPLv3)'),
'Programming Language :: Python',
'Operating System :: OS Independent',
'Intended Audience :: Science/Research',
View it on GitLab: https://salsa.debian.org/debian-gis-team/pykdtree/-/compare/6d4081bf93b72963a73f8021dd2b591221f16e37...d6b3c24124cefeeb9a78fb82ee42aed5b9be0f33
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pykdtree/-/compare/6d4081bf93b72963a73f8021dd2b591221f16e37...d6b3c24124cefeeb9a78fb82ee42aed5b9be0f33
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/20201110/149cde28/attachment-0001.html>
More information about the Pkg-grass-devel
mailing list