knnsearch() fails when using kdtree NSMethod
Show older comments
Hi,
I don't think I understand the knnsearch() function. I read the documentation https://www.mathworks.com/help/stats/classification-using-nearest-neighbors.html#bsehylk explaining the difference between the exhaustive search method (which finds the distance between my query point and every single particle) and the kdtree method (which first groups my data into nodes in order to make the search easier and faster). However, for some strange reason, when I try to apply the kdtree method to my data it fails where I wouldn't expect it to.
Some background about the problem:
I simply want to be able to find the nearest neighbours of different particles within a lattice. In this example below I have a Face-Centered-Cubic (FCC) lattice and I want to get the indices of the nearest neighbours for a given query particle. In the bulk of a FCC lattice each particle should have 12 nearest neighbours, and each neighbour will have the same distance from the query particle.
When I compare the two methods, however, choosing the kdtree method returns an incorrect result where only 3 neighbours are returned. Meanwhile the exhaustive method gets it right with all 12 neighbours. I am not sure why one method fails while the other succeeds? Maybe it has something to do with the many particles tied for nearest? I thought it might be due to the overly large bucket size, but forcing that to 10 didn't improve anything. I would greatly like to reap the benefits of the kdtree search method because while this example lattice is a small 4 by 4 by 4, in general I want my code to be quick and adequate for examining far larger structures.
Green Dots represent the found neighbours. The Red Dot is the query particle.
Can anyone please explain why this is happening?
My code that actually displays the strange behaviour is very short, but I have three helper functions that I use regularly down below as well so you can reproduce what I've been looking at. I would be surprised if it were my helper functions causing issues.
%% Setting up the Problem:
unitCell.latticeVec = [1.2, 1.2, 0; 1.2, 0, 1.2; 0, 1.2, 1.2];
unitCell.basis = {0, 0, 0, "A"};
unitCell.species = {"A", "blue", 1};
trackList = generateLattice(unitCell, dimensions=[4, 4, 4]);
%% Displaying the Problem:
i = 23; %Arbitrarily chosen, but it should have 12 nearest neighbours.
%This wouldn't work if i were 1 or 64.
%I need to replace the query particle in my trackList with a particle at an
%arbitrarily far away distance in order to prevent knnsearch() from finding
%that the nearest neighbour to the query particle is just itself.
treeNeibs = knnsearch([trackList(1:i-1, 1:3); [inf, inf, inf]; trackList(i+1:size(trackList, 1), 1:3)], trackList(i, 1:3), 'IncludeTies', true, 'NSMethod','kdtree');
exhaustiveNeibs = knnsearch([trackList(1:i-1, 1:3); [inf, inf, inf]; trackList(i+1:size(trackList, 1), 1:3)], trackList(i, 1:3), 'IncludeTies', true, 'NSMethod','exhaustive');
%Display the tree method:
display3DParticles(trackList);
scatter3(trackList(i, 1), trackList(i, 2), trackList(i, 3), 180, "red", 'filled');
scatter3(trackList(cell2mat(treeNeibs), 1), trackList(cell2mat(treeNeibs), 2), trackList(cell2mat(treeNeibs), 3), 120, "green", 'filled');
title("Using Tree Method");
disp("treeNeibs: ");
disp(treeNeibs);
%Display the exhaustive method:
display3DParticles(trackList);
scatter3(trackList(i, 1), trackList(i, 2), trackList(i, 3), 180, "red", 'filled');
scatter3(trackList(cell2mat(exhaustiveNeibs), 1), trackList(cell2mat(exhaustiveNeibs), 2), trackList(cell2mat(exhaustiveNeibs), 3), 120, "green", 'filled');
title("Using Exhaustive Method");
disp("exhaustiveNeibs: ");
disp(exhaustiveNeibs);
%% Helper Functions:
%Some of these are written weirdly and inefficiently, but they're not the
%problem.
function display3DParticles(trackList, varargin)
%display3DParticles(trackList, unitCell=unitCell,
%faceResolution=faceResolution, opacity=opacity);
%In:
%trackList: This is a list containing [x, y, z, diameter, species] of
%each particle. All of the x, y, and z values should be positive.
%unitCell: a custom structure containing information about the unit cell
%of your lattice. Do note that the code will let you generate
%unphysical unit cells that make no sense. unitCells are of the form:
%latticeVectors (array), species (cellArray), basis (cellArray). Same
%form as created by generateTriangularLattice().
%fileName: What name to save this image under if you want to save it.
%saveImage: Whether you want to save the generated image as a file.
%resolution: How many extra pixels per original pixels do you want to create in order to
%make the image look not terrible. Maybe don't bother doing this for
%experimental images, but good to do for fabricated data.
%backroundColour: Whether or not to make the background white or black.
%Defaults to black.
p = inputParser;
addRequired(p, "trackList");
addOptional(p, "unitCell", @(x) isstruct(x));
addOptional(p, "fileName", "imageSpheres", @(x) isstring(x));
addOptional(p, "saveImage", false, @(x) isscalar(x) && islogical(x));
addOptional(p, "faceResolution", 25, @(x) isscalar(x) && isnumeric(x));
addOptional(p, "backgroundColour", "black", @(x) isstring(x));
addOptional(p, "opacity", 0.7, @(x) isscalar(x) && isnumeric(x));
parse(p, trackList, varargin{:});
fileName = p.Results.fileName;
saveImage = p.Results.saveImage;
faceResolution = p.Results.faceResolution;
backgroundColour = p.Results.backgroundColour;
opacity = p.Results.opacity;
%fileName = checkSafeSave(fileName);
%Maybe if unitCell unspecified make all particles blue or something.
if(ismember("unitCell", p.UsingDefaults))
maxSpecies = max(trackList(:, 4));
unitCell.species = cell(maxSpecies, 3);
unitCell.species{:, 2} = "blue";
else
unitCell = p.Results.unitCell;
end
%Get a colour map for each particle.
particleSpecies = trackList(:, 5);
colour = arrayfun(@(x) unitCell.species{x, 2}, particleSpecies);
cmap = validatecolor(colour, 'multiple'); %This will rewrite "blue" to [0, 0, 1] for example. Also works with hexcodes.
figure();
hold on;
sphereCell = cell(size(trackList, 1), 6);
[sphereCell{:, 6}] = deal(opacity);
% sphereCell{:, 6} = opacity;
% sphereCell{:, 1:3} = trackList(:, 1:3);
% sphereCell{:, 5} = trackList(:, 5);
% sphereCell{:, 4} = trackList(:, 4) / 2;
%Comma Separated List Assignment?
for i=1:size(trackList, 1)
sphereCell{i, 1} = trackList(i, 1);
sphereCell{i, 2} = trackList(i, 2);
sphereCell{i, 3} = trackList(i, 3);
sphereCell{i, 4} = trackList(i, 4) / 2;
sphereCell{i, 5} = cmap(i, 1:3);
end
viewSpheres(sphereCell, faceResolution=faceResolution);
end
function tracklist = generateLattice(unitCell, varargin)
%tracklist = generateLattice(unitCell, dimensions=dimensions);
%In:
%unitCell: a custom structure containing information about the unit cell
%of your lattice. Do note that the code will let you generate
%unphysical unit cells that make no sense. unitCells are of the form:
%latticeVec (array), species (cellArray), basis (cellArray).
%dimensions: How many unit cells you wish to tile your space with.
%If you want 2D set dimensions(3)=1.
%Out:
%tracklist is an array where the first column is x position, second
%column is y position, 3rd column is z position (assumed 0), 4th column
%is particle diameter, and 5th column is which numbered species the
%particle is. [x, y, z, diameter, species]
p = inputParser;
addRequired(p, "unitCell", @(x) isstruct(x));
addOptional(p, "dimensions", [10, 10, 1], @(x) isnumeric(x) && (size(x, 2) == 3));
parse(p, unitCell, varargin{:});
dimensions = p.Results.dimensions;
%If the 3rd lattice vector is left unspecified then put a dummy lattice
%vector in it's place.
if(size(unitCell.latticeVec, 1) == 2)
currentLatticeVectors = unitCell.latticeVec;
unitCell.latticeVec = vertcat(currentLatticeVectors, [0, 0, 0]);
end
numUnitCells = dimensions(1)*dimensions(2)*dimensions(3);
numParticlesInBasis = size(unitCell.basis, 1);
tracklist = zeros(numUnitCells*numParticlesInBasis, 5);
%Populate the tracklist.
for i=0:dimensions(1)-1
for j=0:dimensions(2)-1
for k=0:dimensions(3)-1
for m=1:numParticlesInBasis
index = i*dimensions(2)*dimensions(3)*numParticlesInBasis+j*dimensions(3)*numParticlesInBasis+k*numParticlesInBasis+m;
species = unitCell.basis{m, 4};
speciesIndex = find(vertcat(unitCell.species{:, 1}) == species);
speciesDiameter = unitCell.species{speciesIndex, 3};
tracklist(index, 4) = speciesDiameter;
tracklist(index, 5) = speciesIndex;
%Multiply our basis vector by the primitive lattice
%vectors to get the actual displacement.
actualVecFromBasis = cell2mat(unitCell.basis(m, 1:3)) * unitCell.latticeVec;
tracklist(index, 1:3) = unitCell.latticeVec(1, :)*i + unitCell.latticeVec(2, :)*j + unitCell.latticeVec(3, :)*k + actualVecFromBasis;
end
end
end
end
end
function viewSpheres(spheres, varargin)
%viewSpheres(spheres, faceResolution, opacityRatio);
%Spheres is an cellArray containing the sphere data we want to plot.
%It is of the form: {x, y, z, radius, color, opacity}.
%faceResolution shows how much detail you want on each sphere. Make
%this higher to better images.
%opacityRatio is an additional multiplier that can be applied to make
%everything more opaque.
p = inputParser;
addRequired(p, "spheres", @(x) iscell(x));
addOptional(p, "faceResolution", 25, @(x) isscalar(x) && isnumeric(x));
addOptional(p, "opacityRatio", 1, @(x) isscalar(x) && isnumeric(x));
parse(p, spheres, varargin{:});
faceResolution = p.Results.faceResolution;
opacityRatio = p.Results.opacityRatio;
[X, Y, Z] = sphere(faceResolution);
%Maybe some stuff with getting the right figure.
hold on;
for i=1:size(spheres, 1)
rad = spheres{i, 4};
surf(X*rad+spheres{i, 1}, Y*rad+spheres{i, 2}, Z*rad+spheres{i, 3}, FaceColor=spheres{i, 5}, FaceAlpha=spheres{i, 6}*opacityRatio, LineStyle="none");
end
axis equal;
%More configurations in order to make the viewing screen good enough.
end
Accepted Answer
This might be an explanation for the different results (from AI):
KD-Trees optimize lookups by comparing bounding boxes and splitting hyperplanes. When multiple data points reside at an identical mathematical distance from your query point, floating-point round-off errors can cause minor variations. [1, 2, 3]
- A point that you expect to have a distance of exactly 1.4142 might be computed as 1.414213562373095 along one tree branch and 1.414213562373096 along another.
- Because the distance is technically greater by a microscopic fraction, the KD-Tree algorithm classifies it as further away than the \(k\)-th smallest distance, failing to identify it as a tie.
Note: Octave returns all 12 nearest neighbours with the kdtree option.
You might want to report it as a potential bug under
But comparing floating point numbers for exact equality is always an errorprone task.
3 Comments
PS:
I found out that floating-point errors are really the cause of your problem.
If you call knnsearch as
[treeNeibs,D] = knnsearch([trackList(1:i-1, 1:3); [inf, inf, inf]; trackList(i+1:size(trackList, 1), 1:3)], trackList(i, 1:3), 'K',12,'IncludeTies', true, 'NSMethod','kdtree');
and compare the distances D, you will find that
isequal(D{1}(1),D{1}(2))
isequal(D{1}(1),D{1}(3))
returns true, but
isequal(D{1}(1),D{1}(4))
isequal(D{1}(1),D{1}(5))
isequal(D{1}(1),D{1}(6))
isequal(D{1}(1),D{1}(7))
isequal(D{1}(1),D{1}(8))
isequal(D{1}(1),D{1}(9))
isequal(D{1}(1),D{1}(10))
isequal(D{1}(1),D{1}(11))
isequal(D{1}(1),D{1}(12))
returns false.
%% Setting up the Problem:
unitCell.latticeVec = [1.2, 1.2, 0; 1.2, 0, 1.2; 0, 1.2, 1.2];
unitCell.basis = {0, 0, 0, "A"};
unitCell.species = {"A", "blue", 1};
trackList = generateLattice(unitCell, dimensions=[4, 4, 4]);
%% Displaying the Problem:
i = 23; %Arbitrarily chosen, but it should have 12 nearest neighbours.
%This wouldn't work if i were 1 or 64.
%I need to replace the query particle in my trackList with a particle at an
%arbitrarily far away distance in order to prevent knnsearch() from finding
%that the nearest neighbour to the query particle is just itself.
[treeNeibs,D] = knnsearch([trackList(1:i-1, 1:3); trackList(i+1:size(trackList, 1), 1:3)], trackList(i, 1:3), 'K',12,'IncludeTies', true, 'NSMethod','kdtree');
Di = typecast(cell2mat(D),'uint64');
dec2bin(Di,64)
function tracklist = generateLattice(unitCell, varargin)
%tracklist = generateLattice(unitCell, dimensions=dimensions);
%In:
%unitCell: a custom structure containing information about the unit cell
%of your lattice. Do note that the code will let you generate
%unphysical unit cells that make no sense. unitCells are of the form:
%latticeVec (array), species (cellArray), basis (cellArray).
%dimensions: How many unit cells you wish to tile your space with.
%If you want 2D set dimensions(3)=1.
%Out:
%tracklist is an array where the first column is x position, second
%column is y position, 3rd column is z position (assumed 0), 4th column
%is particle diameter, and 5th column is which numbered species the
%particle is. [x, y, z, diameter, species]
p = inputParser;
addRequired(p, "unitCell", @(x) isstruct(x));
addOptional(p, "dimensions", [10, 10, 1], @(x) isnumeric(x) && (size(x, 2) == 3));
parse(p, unitCell, varargin{:});
dimensions = p.Results.dimensions;
%If the 3rd lattice vector is left unspecified then put a dummy lattice
%vector in it's place.
if(size(unitCell.latticeVec, 1) == 2)
currentLatticeVectors = unitCell.latticeVec;
unitCell.latticeVec = vertcat(currentLatticeVectors, [0, 0, 0]);
end
numUnitCells = dimensions(1)*dimensions(2)*dimensions(3);
numParticlesInBasis = size(unitCell.basis, 1);
tracklist = zeros(numUnitCells*numParticlesInBasis, 5);
%Populate the tracklist.
for i=0:dimensions(1)-1
for j=0:dimensions(2)-1
for k=0:dimensions(3)-1
for m=1:numParticlesInBasis
index = i*dimensions(2)*dimensions(3)*numParticlesInBasis+j*dimensions(3)*numParticlesInBasis+k*numParticlesInBasis+m;
species = unitCell.basis{m, 4};
speciesIndex = find(vertcat(unitCell.species{:, 1}) == species);
speciesDiameter = unitCell.species{speciesIndex, 3};
tracklist(index, 4) = speciesDiameter;
tracklist(index, 5) = speciesIndex;
%Multiply our basis vector by the primitive lattice
%vectors to get the actual displacement.
actualVecFromBasis = cell2mat(unitCell.basis(m, 1:3)) * unitCell.latticeVec;
tracklist(index, 1:3) = unitCell.latticeVec(1, :)*i + unitCell.latticeVec(2, :)*j + unitCell.latticeVec(3, :)*k + actualVecFromBasis;
end
end
end
end
end
PPS:
Interestingly, the command
[exhaustiveNeibs,D] = knnsearch([trackList(1:i-1, 1:3); [inf, inf, inf]; trackList(i+1:size(trackList, 1), 1:3)], trackList(i, 1:3), 'IncludeTies', true, 'NSMethod','exhaustive');
Di = typecast(cell2mat(D),'uint64');
dec2bin(Di,64)
gives the same distance vector D as
[treeNeibs,D] = knnsearch([trackList(1:i-1, 1:3); trackList(i+1:size(trackList, 1), 1:3)], trackList(i, 1:3), 'K',12,'IncludeTies', true, 'NSMethod','kdtree');
Di = typecast(cell2mat(D),'uint64');
dec2bin(Di,64)
thus a vector where only the first three elements are equal in binary representation.
Thus maybe there is a check with "isapprox" for 'NSMethod','exhaustive' while there is a check with "isequal" for 'NSMethod','kdtree' internally in knnsearch. Only the developers know.
Ian
18 minutes ago
More Answers (0)
Categories
Find more on Mathematics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
