Cannot find the closest point in 4 dimension using `stackBased::fcp` or `stackFree` algorithm

[ahadji05]

Hi, apologies for the big text that follows!

I am trying to implement a find-closest-point in 4 dimensions using this library in the following way.

Initially, I create randomly a number of points in 4D that I use as my search space and then define a query-point that I use to find the closest neighbor using a brute-force O(N) implementation that looks like this:

float sqrtdiff( float4 p, float x, float y, float z, float w ){
    return std::sqrt(std::pow(p.x - x,2) + std::pow(p.y - y,2) + std::pow(p.z - z,2) + std::pow(p.w - w,2));
}

void find_closest( int &closest, float4 p, int numData, const float * x, const float * y, const float * z, const float * w ){
    float Ebest = 1e6;
    for( int i(0); i < numData; ++i ){
        float E = sqrtdiff(p,x[i],y[i],z[i],w[i]);
        if( E < Ebest ){
            Ebest = E;
            closest = i;
        }
    }
}

Then, using the same search data and query point I perform the equivalent implementation using this library as follows:

Initially, I call the routine the builds the k-d tree structure: cukd::buildTree(data,numData);

and then I launch a very simple kernel with just one cuda thread to retrieve the closest neighbor of the same query-point: find_closest_kernel<<<1,1>>>(p, data, numData);

the cuda kernal is the following:

__global__ void find_closest_kernel(float4 queryPoint, float4 *data, int numData) {
    int IDofClosestDataPoint = cukd::stackFree::fcp(queryPoint,data,numData);
    printf("cuda-closest: %d\n",IDofClosestDataPoint);
}

For example, with the following 20 search-points and the query point -> x=10, y=20, z=10, w=20 ID DISTANCE X Y Z W

0: 25.7876 26, 8, -2, 9 1: 73.3689 21, 6, 39, -45 2: 51.5946 -30, 18, 33, 43 3: 37.2425 29, 45, 11, 40 4: 71.225 -28, -36, -8, 33 5: 40.3113 32, -3, -14, 26 6: 52.2398 -38, 15, -6, 8 7: 89.4763 25, -36, 27, -46 8: 92.2009 49, -38, 14, -40 9: 66.1967 -43, -11, -14, 14 10: 79.0759 36, 18, -37, -38 11: 43.4396 -7, 5, 47, 18 12: 88.0114 -43, -6, 4, -45 13: 31.289 -9, 15, 2, 43 14: 63.3167 50, 22, -21, -18 15: 75.1332 -6, -48, 16, -7 16: 69.0652 -19, -41, 2, 32 17: 104.8 29, -49, 41, -50 18: 20.445 30, 24, 9, 19 19: 79.1833 -43, 37, -44, 36

The result I get is: cuda-closest: 6

which is not what I expected since the closest point in search space is point with ID = 18

Did I use the library as I should, or do I miss something?

[ingowald]

It's absolutely possible that there's a bug in there, and I'd be happy to look.

BUT - one thing to check first: the k-d tree builder works by re-arranging the points in the array, so the point with index "6" after the buildtree is not necessarily the same as the point 6 in the input data - it may very well be the right closest point. Could you print/check the actual coordinates of the point found in the query?

[ahadji05]

YES, you are right, thanks 🙏 indeed if I print the coordinates as below the searching identifies the closet point!

__global__ void find_closest_kernel(float4 queryPoint, float4 *data, int numData) {
    int IDofClosestDataPoint = cukd::stackFree::fcp(queryPoint,data,numData);

    printf("cuda-closest: %d\n",IDofClosestDataPoint);

    printf("%f, %f, %f, %f\n",data[IDofClosestDataPoint].x, data[IDofClosestDataPoint].y, 
        data[IDofClosestDataPoint].z, data[IDofClosestDataPoint].w);
}

Is there any way to keep track of the rearranged indices so that based on the new index I can retrieve the old index somehow?

[ahadji05]

The reason that would be helpful is because the index that is stored in the initial array is needed to retrieve data from another dataset that is sorted in the same order as the array of coordinates!

I guess the way to do this is by implementing a custom data-type with a payLoad member variable.

[ahadji05]

OK, I managed to keep track of the original index via a custom data structure based on a payLoad parameter!

The code snippet below shows the working approach.

struct PointPlusIndex {
    float4 pos;
    int    index;
};

struct PointPlusIndex_traits: public cukd::default_data_traits<float4> {
    using point_t = float4;

    static inline __device__ __host__
    float4 get_point(const PointPlusIndex &data)
    { return data.pos; }

    static inline __device__ __host__
    float  get_dim(const PointPlusIndex &data)
    { return 4; }

    static inline __device__ __host__
    float  get_coord(const PointPlusIndex &data, int dim)
    { return cukd::get_coord(get_point(data),dim); }
};

[ingowald]

Nice; yes that's the way it's supposed to work!