Incorrect indexes or distances on ImmutableKdTree

[mattj23]

I'm not sure if I'm misunderstanding something about how Kiddo handles indices or distances, but I've seen some cases where the distances that Kiddo is reporting under the within function does not match the distance calculated manually for the indices kiddo returns.

For instance, the following code uses kiddo version 5.2.2 from crates.io. It has 109 points which come from an actual application of mine. The code creates a 2D ImmutableKdTree from these 109 points, then tests the point at index 4 against the tree.

The expected vector is the brute force results of checking the whole list against point 4. It expects that indices 1 through 7 are within range.

Kiddo returns indices 1 through 4, then 43, 44, and 49. The squared euclidean distances kiddo returns match the manually computed versions for indices 1-7, but not for 43, 44, and 49. The manually computed distances for 43, 44, and 49 are two orders of magnitude above the distance threshold used in the within function.

It looks to me like kiddo is returning indices in a different order than they are in the original list, which is making me wonder if I'm misunderstanding something about how the indices are stored.

use kiddo::immutable::float::kdtree;
use kiddo::SquaredEuclidean;

fn main() {
    let points = test_data();
    let tree = kdtree::ImmutableKdTree::<f64, usize, 2, 32>::new_from_slice(&points);
    let dist = 6.0;
    let dist_sq = dist * dist;
    let test_index = 4;

    let result = tree.within::<SquaredEuclidean>(&points[test_index], dist_sq);

    let mut expected = Vec::new();
    for (i, p) in points.iter().enumerate() {
        let d_sq = distance_squared(&points[test_index], p);
        if d_sq <= dist_sq {
            expected.push(i);
        }
    }

    println!("Expected: {:?}", expected);
    let mut found = result.iter().map(|r| r.item).collect::<Vec<_>>();
    found.sort();

    println!("Found: {:?}", found);

    for r in result {
        let neighbor = &points[r.item];
        let checked_squared = distance_squared(&points[test_index], neighbor);
        println!("Point i={}, dist^2 kiddo={}, manual={}", r.item, r.distance, checked_squared);
    }

}

fn distance_squared(a: &[f64; 2], b: &[f64; 2]) -> f64 {
    (a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2)
}

fn test_data() -> Vec<[f64; 2]> {
    vec![
        [3317.0756414929883, 811.7122408967787],
        [3318.251752812044, 810.4325124656744],
        [3319.4278641310993, 809.1527840345702],
        [3320.603975450155, 807.8730556034661],
        [3321.7800867692104, 806.5933271723619],
        [3322.956198088266, 805.3135987412577],
        [3324.1323094073214, 804.0338703101535],
        [3325.1233671193845, 802.5234025770475],
        [3326.1144248314476, 801.0129348439415],
        [3325.995136117253, 799.9751975016522],
        [3325.8758474030583, 798.9374601593631],
        [3326.3636855948803, 797.2830881719838],
        [3326.851523786702, 795.6287161846043],
        [3327.339361978524, 793.9743441972249],
        [3328.3787056579517, 793.9846642699185],
        [3329.4180493373797, 793.9949843426123],
        [3330.046449352118, 795.7532903690851],
        [3330.6748493668565, 797.511596395558],
        [3331.303249381595, 799.2699024220307],
        [3331.9316493963333, 801.0282084485035],
        [3333.407854588155, 801.0867536751596],
        [3334.8840597799767, 801.1452989018157],
        [3336.605294699901, 801.2866065912483],
        [3338.3265296198256, 801.4279142806807],
        [3339.4855974129955, 799.8317465248002],
        [3340.644665206165, 798.2355787689198],
        [3341.803732999335, 796.6394110130393],
        [3342.962800792505, 795.0432432571588],
        [3344.3555759433593, 795.3023618471735],
        [3345.748351094214, 795.5614804371883],
        [3347.2328465494534, 795.8232509905426],
        [3348.08670862355, 795.9153914507592],
        [3349.559617942115, 795.9708806897642],
        [3350.6376181983205, 795.9827334003453],
        [3351.715618454526, 795.9945861109263],
        [3352.94028646958, 795.9971006088998],
        [3354.1649544846337, 795.9996151068733],
        [3355.4140794822115, 795.9998075661966],
        [3356.663204479789, 796.00000002552],
        [3357.913204479789, 796.00000001276],
        [3359.163204479789, 796.0],
        [3360.413204479789, 796.0],
        [3361.663204479789, 796.0],
        [3362.913204479789, 796.0],
        [3364.163204479789, 796.0],
        [3365.413204479789, 796.0],
        [3366.663204479789, 796.0],
        [3367.913204479789, 796.0],
        [3369.163204479789, 796.0],
        [3370.413204479789, 796.0],
        [3371.663204479789, 796.0],
        [3372.913204479789, 795.4473296001761],
        [3374.163204479789, 794.8946592003522],
        [3375.413204479789, 793.4432227738059],
        [3376.663204479789, 791.9917863472598],
        [3377.913204479789, 790.5397362845654],
        [3379.163204479789, 789.0876862218711],
        [3380.413204479789, 788.3528340151082],
        [3381.663204479789, 787.6179818083453],
        [3382.913204479789, 787.7815929334186],
        [3384.163204479789, 787.9452040584916],
        [3385.413204479789, 788.1091338254666],
        [3386.663204479789, 788.2730635924416],
        [3387.913204479789, 789.2799232655191],
        [3389.163204479789, 790.2867829385964],
        [3390.413204479789, 791.4385968429092],
        [3391.663204479789, 792.590410747222],
        [3392.909947993963, 793.0832482743491],
        [3394.156691508137, 793.5760858014762],
        [3395.351695376479, 793.218240953876],
        [3396.5466992448205, 792.8603961062759],
        [3397.5341084798883, 792.3360608910593],
        [3398.5215177149557, 791.8117256758427],
        [3399.864388200669, 791.948360011792],
        [3400.362416082534, 793.5833457391662],
        [3400.8604439643987, 795.2183314665405],
        [3401.358471846263, 796.8533171939147],
        [3401.856499728128, 798.4883029212889],
        [3402.5676357016137, 797.3490994260089],
        [3403.278771675099, 796.2098959307291],
        [3403.8182218519833, 795.104044878583],
        [3404.357672028868, 793.9981938264369],
        [3404.1893224114624, 792.2313615329177],
        [3404.0209727940573, 790.4645292393985],
        [3403.852623176652, 788.6976969458792],
        [3403.684273559247, 786.93086465236],
        [3403.5159239418413, 785.1640323588408],
        [3404.614894836291, 784.8565517117287],
        [3405.713865730741, 784.5490710646166],
        [3407.460747589397, 785.3224266341614],
        [3409.2076294480535, 786.0957822037062],
        [3410.9545113067097, 786.869137773251],
        [3412.701393165366, 787.6424933427959],
        [3414.448275024022, 788.4158489123406],
        [3416.195156882678, 789.1892044818854],
        [3417.9420387413343, 789.9625600514303],
        [3419.6889205999905, 790.7359156209751],
        [3421.435802458647, 791.5092711905199],
        [3423.182684317303, 792.2826267600648],
        [3424.9295661759593, 793.0559823296096],
        [3426.6764480346155, 793.8293378991544],
        [3428.423329893272, 794.6026934686993],
        [3430.170211751928, 795.3760490382441],
        [3431.917093610584, 796.1494046077889],
        [3433.66397546924, 796.9227601773337],
        [3435.4108573278963, 797.6961157468785],
        [3437.1577391865526, 798.4694713164233],
        [3438.904621045209, 799.2428268859682],
        [3440.651502903865, 800.016182455513],
    ]
}

The output that I get on a windows machine, Rust edition 2024, rustc 1.87.0 is below. You can see that indices 1-4 match the manually computed distances, and 43, 44, and 49 do not.

Expected: [1, 2, 3, 4, 5, 6, 7]
Found: [1, 2, 3, 4, 43, 44, 49]
Point i=4, dist^2 kiddo=0, manual=0
Point i=3, dist^2 kiddo=3.0209426921869493, manual=3.0209426921869493
Point i=49, dist^2 kiddo=3.0209426921869493, manual=2477.398718831687
Point i=2, dist^2 kiddo=12.083770768747215, manual=12.083770768747215
Point i=43, dist^2 kiddo=12.083770768747797, manual=1804.1519531730094
Point i=1, dist^2 kiddo=27.18848422968167, manual=27.18848422968167
Point i=44, dist^2 kiddo=27.741809711405963, manual=1908.5472474494557

[JohnPiwinski]

I second this problem. I've been using Kiddo for a number theory research problem about counting lattice points in spherical caps, but the results seem to fluctuate by an error term proportional to the boundary.

My plan is to just write down my preliminary results and then rerun the analysis once this bug has been resolved.

I wish I were more of a programmer, because then I would offer to investigate this myself. Unfortunately, many of the library internals are quite complicated for me.

[cbueth]

Hey @mattj23, thanks for discussing the features of this package. I tried to reproduce your problem and ran the the code you shared. For me it returns the correct distances and incides with the version on master (2056051, 5.2.1). Switching to the 5.2.4 release branch I also get the same, correct outputs. Can you review if you still get the wrong outputs?

Expected: [1, 2, 3, 4, 5, 6, 7]
Found: [1, 2, 3, 4, 5, 6, 7]
Point i=4, dist^2 kiddo=0, manual=0
Point i=3, dist^2 kiddo=3.0209426921869493, manual=3.0209426921869493
Point i=5, dist^2 kiddo=3.0209426921869493, manual=3.0209426921869493
Point i=2, dist^2 kiddo=12.083770768747215, manual=12.083770768747215
Point i=6, dist^2 kiddo=12.083770768747797, manual=12.083770768747797
Point i=1, dist^2 kiddo=27.18848422968167, manual=27.18848422968167
Point i=7, dist^2 kiddo=27.741809711405963, manual=27.741809711405963

[cbueth]

[...] I've been using Kiddo for a number theory research problem about counting lattice points in spherical caps, but the results seem to fluctuate by an error term proportional to the boundary.

Hey @JohnPiwinski, you might be interested in the thoughts shared in #285. Did you implement a DistanceMetric for spherical caps? I think there are some grains of salt to watch out for. In your case it would be the great-circle distance as the shortest path along the surface of the sphere. Currently as per my impression, the kiddo crate is designed for Euclidean spaces (i.e. k-d trees in $\mathbb{R}^n$). Spherical caps live on the surface of a sphere (i.e. $\mathbb{S}^2$), so you’d need to map spherical coordinates to Cartesian coordinates (latitude/longitude to 3D points on a unit sphere). The grain of salt here: the dist1 method might be problematic for spherical distances, as the concept of distance along a single axis doesn’t directly translate to spherical geometry. I'd be happt yo hear what was your approach.