This is an automated email from the git hooks/post-receive script.
a_valentino-guest pushed a commit to branch master
in repository pykdtree.
commit 4d17f6202805a0341d8c7b07916ffc7249eeeded
Author: Antonio Valentino <antonio.valent...@tiscali.it>
Date: Sat Jul 30 10:05:41 2016 +0000
Imported Upstream version 1.2.1
---
PKG-INFO | 2 +-
pykdtree.egg-info/PKG-INFO | 2 +-
pykdtree/_kdtree_core.c | 146 ++++++++++++++++++++++++++-------------------
pykdtree/test_tree.py | 21 +++++++
setup.py | 10 +++-
5 files changed, 114 insertions(+), 67 deletions(-)
diff --git a/PKG-INFO b/PKG-INFO
index 6fb385f..9b78cbe 100644
--- a/PKG-INFO
+++ b/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.1
Name: pykdtree
-Version: 1.1.1
+Version: 1.2.1
Summary: Fast kd-tree implementation with OpenMP-enabled queries
Home-page: UNKNOWN
Author: Esben S. Nielsen
diff --git a/pykdtree.egg-info/PKG-INFO b/pykdtree.egg-info/PKG-INFO
index 6fb385f..9b78cbe 100644
--- a/pykdtree.egg-info/PKG-INFO
+++ b/pykdtree.egg-info/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.1
Name: pykdtree
-Version: 1.1.1
+Version: 1.2.1
Summary: Fast kd-tree implementation with OpenMP-enabled queries
Home-page: UNKNOWN
Author: Esben S. Nielsen
diff --git a/pykdtree/_kdtree_core.c b/pykdtree/_kdtree_core.c
index 3cfbfdd..bde453f 100644
--- a/pykdtree/_kdtree_core.c
+++ b/pykdtree/_kdtree_core.c
@@ -31,6 +31,10 @@ Anne M. Archibald and libANN by David M. Mount and Sunil
Arya.
#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
{
@@ -160,21 +164,22 @@ Params:
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;
+ int8_t bbox_idx, i, j;
+ uint32_t i2;
/* Use first data point to initialize */
- for (int8_t i = 0; i < no_dims; i++)
+ for (i = 0; i < no_dims; i++)
{
bbox[2 * i] = bbox[2 * i + 1] = PA(0, i);
}
/* Update using rest of data points */
- for (uint32_t i = 1; i < n; i++)
+ for (i2 = 1; i2 < n; i2++)
{
- for (int8_t j = 0; j < no_dims; j++)
+ for (j = 0; j < no_dims; j++)
{
bbox_idx = 2 * j;
- cur = PA(i, j);
+ cur = PA(i2, j);
if (cur < bbox[bbox_idx])
{
bbox[bbox_idx] = cur;
@@ -203,13 +208,13 @@ Params:
************************************************/
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;
- uint32_t p, q;
+ 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 (int8_t i = 0; i < no_dims; i++)
+ for (i = 0; i < no_dims; i++)
{
side_len = bbox[2 * i + 1] - bbox[2 * i];
if (side_len > size)
@@ -268,13 +273,13 @@ int partition_float(float *pa, uint32_t *pidx, int8_t
no_dims, uint32_t start_id
uint32_t j = start_idx;
split = PA(j, dim);
- for (uint32_t i = start_idx + 1; i <= end_idx; i++)
+ for (i2 = start_idx + 1; i2 <= end_idx; i2++)
{
/* Find lowest point */
- cur_val = PA(i, dim);
+ cur_val = PA(i2, dim);
if (cur_val < split)
{
- j = i;
+ j = i2;
split = cur_val;
}
}
@@ -290,13 +295,13 @@ int partition_float(float *pa, uint32_t *pidx, int8_t
no_dims, uint32_t start_id
uint32_t j = end_idx;
split = PA(j, dim);
- for (uint32_t i = start_idx; i < end_idx; i++)
+ for (i2 = start_idx; i2 < end_idx; i2++)
{
/* Find highest point */
- cur_val = PA(i, dim);
+ cur_val = PA(i2, dim);
if (cur_val > split)
{
- j = i;
+ j = i2;
split = cur_val;
}
}
@@ -327,6 +332,10 @@ Node_float* construct_subtree_float(float *pa, uint32_t
*pidx, int8_t no_dims, u
/* 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 */
@@ -335,11 +344,6 @@ Node_float* construct_subtree_float(float *pa, uint32_t
*pidx, int8_t no_dims, u
else
{
/* Make split node */
- int rval;
- int8_t cut_dim;
- uint32_t n_lo;
- float cut_val;
-
/* 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)
@@ -351,8 +355,8 @@ Node_float* construct_subtree_float(float *pa, uint32_t
*pidx, int8_t no_dims, u
root->cut_dim = cut_dim;
/* Recurse on both subsets */
- float lv = bbox[2 * cut_dim];
- float hv = bbox[2 * cut_dim + 1];
+ lv = bbox[2 * cut_dim];
+ hv = bbox[2 * cut_dim + 1];
/* Set bounds for cut dimension */
root->cut_bounds_lv = lv;
@@ -382,16 +386,20 @@ Params:
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 */
- uint32_t *pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
- for (uint32_t i = 0; i < n; i++)
+ pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
+ for (i = 0; i < n; i++)
{
pidx[i] = i;
}
- float *bbox = (float *)malloc(2 * sizeof(float) * no_dims);
+ bbox = (float *)malloc(2 * sizeof(float) * no_dims);
get_bounding_box_float(pa, pidx, no_dims, n, bbox);
tree->bbox = bbox;
@@ -462,7 +470,8 @@ Print
************************************************/
void print_tree_float(Node_float *root, int level)
{
- for (int i = 0; i < level; i++)
+ int i;
+ for (i = 0; i < level; i++)
{
printf(" ");
}
@@ -483,7 +492,8 @@ float calc_dist_float(float *point1_coord, float
*point2_coord, int8_t no_dims)
{
/* Calculate squared distance */
float dist = 0, dim_dist;
- for (int8_t i = 0; i < no_dims; i++)
+ int8_t i;
+ for (i = 0; i < no_dims; i++)
{
dim_dist = point2_coord[i] - point1_coord[i];
dist += dim_dist * dim_dist;
@@ -529,8 +539,9 @@ Params:
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 (int8_t i = 0; i < no_dims; 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;
@@ -555,8 +566,9 @@ void search_leaf_float(float *restrict pa, uint32_t
*restrict pidx, int8_t no_di
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 (uint32_t i = 0; i < n; i++)
+ 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);
@@ -678,6 +690,7 @@ void search_tree_float(Tree_float *tree, float *pa, float
*point_coords,
int8_t no_dims = tree->no_dims;
float *bbox = tree->bbox;
uint32_t *pidx = tree->pidx;
+ uint32_t i, j;
Node_float *root = (Node_float *)tree->root;
/* Queries are OpenMP enabled */
@@ -686,10 +699,10 @@ void search_tree_float(Tree_float *tree, float *pa, float
*point_coords,
/* The low chunk size is important to avoid L2 cache trashing
for spatial coherent query datasets
*/
- #pragma omp for schedule(static, 100) nowait
- for (uint32_t i = 0; i < num_points; i++)
+ #pragma omp for private(i, j) schedule(static, 100) nowait
+ for (i = 0; i < num_points; i++)
{
- for (uint32_t j = 0; j < k; j++)
+ for (j = 0; j < k; j++)
{
closest_idxs[i * k + j] = UINT32_MAX;
closest_dists[i * k + j] = DBL_MAX;
@@ -741,21 +754,22 @@ Params:
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;
+ int8_t bbox_idx, i, j;
+ uint32_t i2;
/* Use first data point to initialize */
- for (int8_t i = 0; i < no_dims; i++)
+ for (i = 0; i < no_dims; i++)
{
bbox[2 * i] = bbox[2 * i + 1] = PA(0, i);
}
/* Update using rest of data points */
- for (uint32_t i = 1; i < n; i++)
+ for (i2 = 1; i2 < n; i2++)
{
- for (int8_t j = 0; j < no_dims; j++)
+ for (j = 0; j < no_dims; j++)
{
bbox_idx = 2 * j;
- cur = PA(i, j);
+ cur = PA(i2, j);
if (cur < bbox[bbox_idx])
{
bbox[bbox_idx] = cur;
@@ -784,13 +798,13 @@ Params:
************************************************/
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;
- uint32_t p, q;
+ 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 (int8_t i = 0; i < no_dims; i++)
+ for (i = 0; i < no_dims; i++)
{
side_len = bbox[2 * i + 1] - bbox[2 * i];
if (side_len > size)
@@ -849,13 +863,13 @@ int partition_double(double *pa, uint32_t *pidx, int8_t
no_dims, uint32_t start_
uint32_t j = start_idx;
split = PA(j, dim);
- for (uint32_t i = start_idx + 1; i <= end_idx; i++)
+ for (i2 = start_idx + 1; i2 <= end_idx; i2++)
{
/* Find lowest point */
- cur_val = PA(i, dim);
+ cur_val = PA(i2, dim);
if (cur_val < split)
{
- j = i;
+ j = i2;
split = cur_val;
}
}
@@ -871,13 +885,13 @@ int partition_double(double *pa, uint32_t *pidx, int8_t
no_dims, uint32_t start_
uint32_t j = end_idx;
split = PA(j, dim);
- for (uint32_t i = start_idx; i < end_idx; i++)
+ for (i2 = start_idx; i2 < end_idx; i2++)
{
/* Find highest point */
- cur_val = PA(i, dim);
+ cur_val = PA(i2, dim);
if (cur_val > split)
{
- j = i;
+ j = i2;
split = cur_val;
}
}
@@ -908,6 +922,10 @@ Node_double* construct_subtree_double(double *pa, uint32_t
*pidx, int8_t no_dims
/* 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 */
@@ -916,11 +934,6 @@ Node_double* construct_subtree_double(double *pa, uint32_t
*pidx, int8_t no_dims
else
{
/* Make split node */
- int rval;
- int8_t cut_dim;
- uint32_t n_lo;
- double cut_val;
-
/* 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)
@@ -932,8 +945,8 @@ Node_double* construct_subtree_double(double *pa, uint32_t
*pidx, int8_t no_dims
root->cut_dim = cut_dim;
/* Recurse on both subsets */
- double lv = bbox[2 * cut_dim];
- double hv = bbox[2 * cut_dim + 1];
+ lv = bbox[2 * cut_dim];
+ hv = bbox[2 * cut_dim + 1];
/* Set bounds for cut dimension */
root->cut_bounds_lv = lv;
@@ -963,16 +976,20 @@ Params:
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 */
- uint32_t *pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
- for (uint32_t i = 0; i < n; i++)
+ pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
+ for (i = 0; i < n; i++)
{
pidx[i] = i;
}
- double *bbox = (double *)malloc(2 * sizeof(double) * no_dims);
+ bbox = (double *)malloc(2 * sizeof(double) * no_dims);
get_bounding_box_double(pa, pidx, no_dims, n, bbox);
tree->bbox = bbox;
@@ -1043,7 +1060,8 @@ Print
************************************************/
void print_tree_double(Node_double *root, int level)
{
- for (int i = 0; i < level; i++)
+ int i;
+ for (i = 0; i < level; i++)
{
printf(" ");
}
@@ -1064,7 +1082,8 @@ double calc_dist_double(double *point1_coord, double
*point2_coord, int8_t no_di
{
/* Calculate squared distance */
double dist = 0, dim_dist;
- for (int8_t i = 0; i < no_dims; i++)
+ int8_t i;
+ for (i = 0; i < no_dims; i++)
{
dim_dist = point2_coord[i] - point1_coord[i];
dist += dim_dist * dim_dist;
@@ -1110,8 +1129,9 @@ Params:
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 (int8_t i = 0; i < no_dims; 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;
@@ -1136,8 +1156,9 @@ void search_leaf_double(double *restrict pa, uint32_t
*restrict pidx, int8_t no_
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 (uint32_t i = 0; i < n; i++)
+ 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);
@@ -1259,6 +1280,7 @@ void search_tree_double(Tree_double *tree, double *pa,
double *point_coords,
int8_t no_dims = tree->no_dims;
double *bbox = tree->bbox;
uint32_t *pidx = tree->pidx;
+ uint32_t i, j;
Node_double *root = (Node_double *)tree->root;
/* Queries are OpenMP enabled */
@@ -1267,10 +1289,10 @@ void search_tree_double(Tree_double *tree, double *pa,
double *point_coords,
/* The low chunk size is important to avoid L2 cache trashing
for spatial coherent query datasets
*/
- #pragma omp for schedule(static, 100) nowait
- for (uint32_t i = 0; i < num_points; i++)
+ #pragma omp for private(i, j) schedule(static, 100) nowait
+ for (i = 0; i < num_points; i++)
{
- for (uint32_t j = 0; j < k; j++)
+ for (j = 0; j < k; j++)
{
closest_idxs[i * k + j] = UINT32_MAX;
closest_dists[i * k + j] = DBL_MAX;
diff --git a/pykdtree/test_tree.py b/pykdtree/test_tree.py
index df89775..56f8497 100644
--- a/pykdtree/test_tree.py
+++ b/pykdtree/test_tree.py
@@ -269,6 +269,27 @@ def test3d_8n_ub_eps():
assert np.array_equal(idx, exp_idx)
assert np.allclose(dist, exp_dist)
+def test3d_large_query():
+ # Target idxs: 7, 93, 45
+ query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
+ [751763.125, -59925.969, 6326205.5],
+ [769957.188, -202418.125, 6321069.5]])
+
+ # Repeat the same points multiple times to get 60000 query points
+ n = 20000
+ query_pts = np.repeat(query_pts, n, axis=0)
+
+ kdtree = KDTree(data_pts_real)
+ dist, idx = kdtree.query(query_pts, sqr_dists=True)
+
+ epsilon = 1e-5
+ assert np.all(idx[:n] == 7)
+ assert np.all(idx[n:2*n] == 93)
+ assert np.all(idx[2*n:] == 45)
+ assert np.all(dist[:n] == 0)
+ assert np.all(abs(dist[n:2*n] - 3.) < epsilon * dist[n:2*n])
+ assert np.all(abs(dist[2*n:] - 20001.) < epsilon * dist[2*n:])
+
def test_scipy_comp():
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
diff --git a/setup.py b/setup.py
index 3706921..c0c0b22 100644
--- a/setup.py
+++ b/setup.py
@@ -37,10 +37,14 @@ class build_ext_subclass(build_ext):
else:
extra_compile_args = ['-std=c99', '-O3']
extra_link_args = []
+ elif comp == 'msvc':
+ extra_compile_args = ['/Ox']
+ extra_link_args = []
+ if use_omp:
+ extra_compile_args.append('/openmp')
else:
# Add support for more compilers here
- raise ValueError(('Compiler flags undefined for %s.',
- 'Please modify setup.py and add compiler flags')
+ raise ValueError('Compiler flags undefined for %s. Please modify
setup.py and add compiler flags'
% comp)
self.extensions[0].extra_compile_args = extra_compile_args
self.extensions[0].extra_link_args = extra_link_args
@@ -61,7 +65,7 @@ class build_ext_subclass(build_ext):
setup(
name='pykdtree',
- version='1.1.1',
+ version='1.2.1',
description='Fast kd-tree implementation with OpenMP-enabled queries',
author='Esben S. Nielsen',
author_email='storpipf...@gmail.com',
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-grass/pykdtree.git
_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
