BUG: incorrect nearest neighbor search using workers=-1

[domist07]

Describe your issue.

If I using KDTree().query() with the argument workers and the values 1 or -1 I got different results. The result, using multiple workers, is wrong. The error only appears using many points in the query task (maybe you have to increase the number of points of coo2 in my code example). Using fewer points, all is fine.

Reproducing Code Example

import numpy as np
from scipy.spatial import KDTree

coo1 = np.random.rand(int(3e3), 2)
coo2 = np.random.rand(int(3e8), 2)

tree = KDTree(coo1)
_, idx_single = tree.query(coo2, k=1, workers=1)
idx_single_bin = np.bincount(idx_single)

_, idx_multi = tree.query(coo2, k=1, workers=-1)
idx_multi_bin = np.bincount(idx_multi)

if idx_single_bin[0] != idx_multi_bin[0]:
    raise ValueError("result is not equal")

Error message

Traceback (most recent call last):
  File "d:\studom\CodeProjects\voxel_compare\issue-workers.py", line 1
5, in <module>
    raise ValueError("result is not equal")
ValueError: result is not equal

SciPy/NumPy/Python version and system information

1.11.1 1.25.0 sys.version_info(major=3, minor=11, micro=4, releaselevel='final', serial=0)
lapack_armpl_info:
  NOT AVAILABLE
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
  NOT AVAILABLE
openblas_clapack_info:
  NOT AVAILABLE
flame_info:
  NOT AVAILABLE
accelerate_info:
  NOT AVAILABLE
atlas_3_10_threads_info:
  NOT AVAILABLE
atlas_3_10_info:
  NOT AVAILABLE
atlas_threads_info:
  NOT AVAILABLE
atlas_info:
  NOT AVAILABLE
lapack_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas']
    library_dirs = ['.../Miniconda3/envs/env\\Library\\lib']
    language = f77
blas_armpl_info:
  NOT AVAILABLE
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
atlas_3_10_blas_threads_info:
  NOT AVAILABLE
atlas_3_10_blas_info:
  NOT AVAILABLE
atlas_blas_threads_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
blas_info:
    libraries = ['cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['.../Miniconda3/envs/env\\Library\\lib']
    include_dirs = ['.../Miniconda3/envs/env\\Library\\include']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    libraries = ['cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['.../Miniconda3/envs/env\\Library\\lib']
    include_dirs = ['.../Miniconda3/envs/env\\Library\\include']
    language = f77
lapack_opt_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas', 'cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['.../Miniconda3/envs/env\\Library\\lib']
    language = f77
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    include_dirs = ['.../Miniconda3/envs/env\\Library\\include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2,AVX512F,AVX512CD,AVX512_SKX
    not found = AVX512_CLX,AVX512_CNL,AVX512_ICL

[j-bowhay]

cc @sturlamolden

[sturlamolden]

Several people have messed with the multithreaded code since I last layed my hands on it. It might be that someone has not understood it properly.

But before I start to debug, can you please make sure there are no ties in your data set? Ties in the context of a kd-tree are points with equal distance to the query point. In the presence of ties it might be that both queries in theory is correct, but if someone has made a race condition results cound be different. We should consider such a race condition a bug, but it is a different kind of bug than a simple logics error.

I am tempted to move the multithreaded code into C++ (threads are now in the standard, just to make it more temper proof).

[sturlamolden]

Also get rid of np.bincount. I don’t care for a wild goose chase that turns out to be an error in NumPy. Use e.g. np.array_equal – and also compare the single-threaded result with itself, as a negative control.

[domist07]

Thanks for the hints. I improved my Code Example. No similar distances are detected, but the error still exists:

import numpy as np
from scipy.spatial import KDTree

# initialize test data
coo1 = np.random.rand(int(3e3), 2)
coo2 = np.random.rand(int(3e8), 2)
# initialize KDTree
tree = KDTree(coo1)

# single core neighbor search
dist_single, idx_single = tree.query(coo2, k=1, workers=1)
# check, if singe core result itself is equal
if not np.array_equal(idx_single, idx_single):
    raise ValueError("single itself is not equal")

# multi core neighbor search
dist_multi, idx_multi = tree.query(coo2, k=1, workers=-1)
# check, if multi core result itself is equal
if not np.array_equal(idx_multi, idx_multi):
    raise ValueError("multi itself is not equal")

# check for equal distances
dist_single_u = np.unique(dist_single)
if not len(dist_single_u) == len(dist_single):
    print("similar distances detected")

# compare single vs multi core neighbor search
if not np.array_equal(idx_single, idx_multi):
    raise ValueError("results are not equal")

[sturlamolden]

Which of the tests fail?

[sturlamolden]

Your check for equal distances is wrong. You are taking the query results but you need to check all 3e3 x 3e8 distances against eachother. Basically you have ((3e3 x 3e8)**2 - 3e3x3e8)/2 pairs to compare. Compute all 3e3 x 3e8 distances and push them into a set to verify that they are all unique.

[domist07]

Which of the tests fail?

Only the last one, so I got the message: ValueError("results are not equal")

[domist07]

Your check for equal distances is wrong. You are taking the query results but you need to check all 3e3 x 3e8 distances against eachother. Basically you have ((3e3 x 3e8)**2 - 3e3x3e8)/2 pairs to compare. Compute all 3e3 x 3e8 distances and push them into a set to verify that they are all unique.

Did you agree with this code?

from scipy.spatial import distance

# check for equal distances
distances = distance.pdist(coo2)
equal_distances = np.isclose(distances, 0)
if np.any(equal_distances):
    print("similar distances detected")

But the problem is: I can not compute it, because it needs too much memory 😞 .

Did you have any idea, how to compute it, using less memory or do another check?

[sturlamolden]

Did you agree with this code?

from scipy.spatial import distance

# check for equal distances
distances = distance.pdist(coo2)
equal_distances = np.isclose(distances, 0)
if np.any(equal_distances):
    print("similar distances detected")

No, I do not.

But the problem is: I can not compute it, because it needs too much memory 😞 .

Did you have any idea, how to compute it, using less memory or do another check?

Memory is not the problem, a small loop in C will do this with almost zero memory overhead.

The real issue is time.