How to extract 6 neasrest neigbors for particles

Hello everyone.

I have a text file (kindly see the attached text file) containing neighbors of particles (ranging from 1 to 64). I want to extract 6 nearest neighbors of every particle in a seperate file.

Following is the sample from the output file. Column 1 and 2 show the particle number and column 3 show the distance.

1 19 0.874004

1 41 1.130905

1 42 1.131306

1 45 0.841402

1 50 1.393870

1 52 0.961277

1 61 0.802248

1 63 1.327221

1 64 1.408026

2 16 0.819803

2 33 1.153725

2 37 0.954796

2 39 1.439849

2 42 1.419793

2 58 1.140030

2 60 0.912613

2 63 1.454144

3 6 0.843783

3 7 1.071621

3 9 1.009899

3 10 1.254684

3 13 1.108856

3 14 1.015914

3 20 0.885181

3 35 1.263609

3 45 1.356776

3 51 1.268250

4 5 1.197906

4 40 1.238262

4 49 1.131004

4 53 1.207126

4 54 0.848770

4 57 1.278407

4 59 1.398415

4 62 0.782374

4 64 0.872153

The output which I want is

1 19 0.874004

1 41 1.130905

1 42 1.131306

1 45 0.841402

1 50 1.393870

1 52 0.961277

2 16 0.819803

2 33 1.153725

2 37 0.954796

2 39 1.439849

2 42 1.419793

2 58 1.140030

3 6 0.843783

3 7 1.071621

3 9 1.009899

3 10 1.254684

3 13 1.108856

3 14 1.015914

4 5 1.197906

4 40 1.238262

4 49 1.131004

4 53 1.207126

4 54 0.848770

4 57 1.278407

Any help will be appreciated.

Thanks in advance,

Abhinav

[Rik]

What have you tried already? This doesn't sound too hard with a loop. What issues did you have implementing it?

[Abhinav Srivastava]

Following is my code:

% Code for determing 6 nearest neighbors
clear 
clc
a=load('particle-distance.txt');
b=a;
count = 0;
natoms = 64;
for i = 1:natoms
    for j = 1:length(a)
        if a(j,1) == i
            count = count + 1;
            if count >= 6
               str = [int2str(a(j,1)),'    ',int2str(a(j,2)),'    ',num2str(a(j,3))] ;
               disp(str)
            end
          end
    end
end

Answer 1

[DGM]

This can probably be improved, but this was what I did:

A = [1 19 0.874004;
	1 41 1.130905;
	1 42 1.131306;
	1 45 0.841402;
	1 50 1.393870;
	1 52 0.961277;
	1 61 0.802248;
	1 63 1.327221;
	1 64 1.408026;
	2 16 0.819803;
	2 33 1.153725;
	2 37 0.954796;
	2 39 1.439849;
	2 42 1.419793;
	2 58 1.140030;
	2 60 0.912613;
	2 63 1.454144;
	3 6 0.843783;
	3 7 1.071621;
	3 9 1.009899;
	3 10 1.254684;
	3 13 1.108856;
	3 14 1.015914;
	3 20 0.885181;
	3 35 1.263609;
	3 45 1.356776;
	3 51 1.268250;
	4 5 1.197906;
	4 40 1.238262;
	4 49 1.131004;
	4 53 1.207126;
	4 54 0.848770;
	4 57 1.278407;
	4 59 1.398415;
	4 62 0.782374;
	4 64 0.872153]

% sort by distance
[~,dmap] = sort(A(:,3),'ascend');
A = A(dmap,:);

npoints = numel(unique(A(:,1)));
idx = [];
for p = 1:npoints
    % this will find at most 6 matches
    idx = [idx find(A(:,1)==p,6)];
end
A = A(idx,:)

gives

A =
            1           61      0.80225
            1           45       0.8414
            1           19        0.874
            1           52      0.96128
            1           41       1.1309
            1           42       1.1313
            2           16       0.8198
            2           60      0.91261
            2           37       0.9548
            2           58         1.14
            2           33       1.1537
            2           42       1.4198
            3            6      0.84378
            3           20      0.88518
            3            9       1.0099
            3           14       1.0159
            3            7       1.0716
            3           13       1.1089
            4           62      0.78237
            4           54      0.84877
            4           64      0.87215
            4           49        1.131
            4            5       1.1979
            4           53       1.2071

[Abhinav Srivastava]

Thank you so much.

I really appreciate your help.

Best,

Abhinav

[Rik]

You can improve performance by using logical indexing:

L = false(size(A,1));
for p = unique(A(:,1)).'
    % this will find at most 6 matches
    L(find(A(:,1)==p,6))=true;
end
A = A(L,:)

A = 24×3

1.0000 61.0000 0.8022 1.0000 45.0000 0.8414 1.0000 19.0000 0.8740 1.0000 52.0000 0.9613 1.0000 41.0000 1.1309 1.0000 42.0000 1.1313 2.0000 16.0000 0.8198 2.0000 60.0000 0.9126 2.0000 37.0000 0.9548 2.0000 58.0000 1.1400

You could probably avoid the double sorting by using ismember.

[DGM]

At least with the testing I did, the speed of the two approaches is the same (about 65us versus 70us for logical). That might change for larger datasets.

However, you're right in that I don't even need to do the second sort, because the find() routine will inherently sort them out. That would save about 9us.

[Rik]

Personally I try to write code like this, because then it doesn't matter much if it is used for a much larger dataset. Once you get in the habit, you will less often write very slow code (although there is probably a lot to gain in most of my code). Even code operating on smaller datasets can matter, if it is called often enough. 65us vs 70us isn't relevant, unless you need to run this several times a second, or once for every frame you're displaying.

[DGM]

I've found a lot of cases, especially with relatively short vectors, where it seems the expense of growing the vector by variable-size chunks is less costly than repeatedly calculating the range of indices into a preallocated vector. I'm also am wary of the idea that logical indexing is generally faster than even seemingly wasteful approaches (a bit off-topic, but there's an example in there)

https://www.mathworks.com/matlabcentral/answers/409702-speed-of-masked-matrix-operations-in-single-vs-double

I know those are only observations on my part and would be terrible to presume as general rules -- but it's enough doubt that it does tend to make me less averse to doing things that aren't recommended, especially when the reward for optimization is low.

I try to test speed with presumed array sizes/types for my own work, but on the forum, half the time I manage to leave in vestigal artifacts of the fact that I started writing the code before I even decided how I was going to finish it (e.g. the redundant sort) -- so I don't know how much room I have to talk about optimal approaches to anything. I know I'm not exactly the most rigorous character around here, but I'm learning a lot of things since I started being active.

[Rik]

I don't pretend to be an expert, but I do know that using a coding style like I suggested is easier for Mathworks engineers to optimize than code that exploits behavior specific to the JIT in a particular release, or CPU, or OS, or whatever. So by using code like this, I can reasonably expect performance to improve more release to release than I otherwise would.

I don't know who the most rigorous person would be round here. I don't think I am. I view this forum as a challenge to solve puzzles, as well as help and educate people. Especially that last part is important in this context. If you teach good coding practices when it doesn't matter, people may be more likely to apply them when it does. Since you seem to be sticking around, I just wanted to make sure you are making a deliberate choice.

[DGM]

I suppose I never really thought about it that way. I probably am optimizing for my own very niche case (old releases on a specific os on old hardware) and have no idea what the typical performance is going to be. I do know that a lot of my optimizations in MIMT are really only useful in much older versions (R2009b) and are likely useless to 99.9% of users now.

Regarding my intent in being here, I only started answering questions to pass time while I was waiting for answers to my own. Now it's just sort of turned into an idle habit. I admit a large part of it is just entertainment and puzzles for me, though that's not all of it. I've been surprised to learn a lot of things I had never had reason to even try before, and it is nice to help people out (when it turns out that I can). I expect my opinion of my own habits will change and shape them moving forward. At least I'm writing (marginally) better formatted code since people have to read it now. I doubt that's the last of my cherished bad habits to suffer erosion.

[Rik]

Ever since I 'discovered' that R13 (v6.5) is very fast on some tasks (thanks to Jan), I made sure my code can run of that. I was already making sure my code runs on Octave so I would able to share things I wrote in a legal manner.

Since I already cover such a wide range, it wasn't a lot of extra effort to check how my code works on different releases. Feel free to check the documentation for my FEX submissions if you want to have a look where that obsession has led so far ;)

There are people on this forum who have been using Matlab before I was done with primary eductation. You can learn a lot if you're willing.

Regarding other people reading your code: you should treat your future self as a different person as well. You probably won't understand uncommented code you wrote 2 years ago. You should at least write comments like this.

[DGM]

I admit I was kind of impressed by that when I started using ifversion(). I'm way too used to expecting FEX things to be burdened with undocumented strict dependencies.

And yeah. I had to learn that if it's been long enough to forget how the code works, it's been long enough to forget the relevance of any amount of terse comments. As I get older, it's not just commenting code where I have to be thorough in my correspondence with future me.

[Abhinav Srivastava]

Hi,

I have successfully able to use your code for generating nearest six neighbors. How ever in my system, there are some particles which do not have exact six neigbors. So, I wanted to eleminate those lines from the text file.

Example,

1 9 1.484157

1 16 0.828089

1 21 1.481878

2 14 1.084812

2 19 0.212855

2 23 0.745430

2 27 0.709192

2 29 0.485757

2 36 1.351770

3 8 1.485509

3 10 0.789934

3 12 0.996034

3 22 1.273936

3 28 1.240151

3 31 1.209610

4 1 0.445649

4 5 0.491061

4 6 0.611696

4 9 1.460465

4 16 1.081862

4 24 1.222578

5 1 0.868254

5 4 0.491061

5 6 0.366953

5 9 1.269292

5 16 1.212136

5 18 1.255205

6 1 1.052711

6 4 0.611696

6 5 0.366953

6 18 1.467878

6 20 1.142525

6 34 0.708171

7 8 0.398548

7 10 1.086699

7 14 1.288270

7 32 0.816957

7 37 0.881323

8 3 1.485509

8 7 0.398548

8 10 1.228414

8 32 0.781834

8 37 0.906441

9 1 1.484157

9 4 1.460465

9 5 1.269292

9 16 0.938138

9 18 0.803028

9 24 1.018272

Here particle number 7, 8 ( first column) have 5 neighbors.

I want have a data file containing exactly 6 neighbors. It would be very nice if you can help.

Thanks in advance,

Abhinav

[DGM]

Sorry about that. The notifications on this site tend to show up out of order, so sometimes I completely miss followup questions if they don't show up until queue has moved several hours past their timestamp.

For example:

A = [1 19 0.874004;
    1 41 1.130905;
    1 50 1.393870;
    1 63 1.327221;
    1 64 1.408026;
    2 16 0.819803;
    2 33 1.153725;
    2 37 0.954796;
    2 39 1.439849;
    2 42 1.419793;
    2 58 1.140030;
    2 60 0.912613;
    2 63 1.454144;
    3 6 0.843783;
    3 7 1.071621;
    3 13 1.108856;
    3 35 1.263609;
    3 45 1.356776;
    4 5 1.197906;
    4 40 1.238262;
    4 49 1.131004;
    4 53 1.207126;
    4 54 0.848770;
    4 57 1.278407;
    4 59 1.398415;
    4 62 0.782374;
    4 64 0.872153]

n = 6; % minimum number of neighbors

[~,dmap] = sort(A(:,3),'ascend');
A = A(dmap,:);

npoints = numel(unique(A(:,1)));
idx = [];
for p = 1:npoints
    % this will find at most n matches
    % if there are fewer than n neighbors, skip it
    thisblock = find(A(:,1)==p,n);
    if numel(thisblock)>=n
        idx = [idx thisblock];
    end
end
A = A(idx,:)

[Abhinav Srivastava]

The above code worked.

Thank you so much for your help.