I was interested if somebody has any experience with quick finding closest neighbours for given points. Lets say I have list of N points with x, y and z coordinates and want to find the closest neighbours for all points. The simplest way would be: for(i=1;N;i++){ for(j=1;N;i++){ r(i,j)=(x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2 if(r(i,j)<rmin){ rmin=r(i,j) x1=x(i) x2=x(j) y1=y(i) ..... } } } and so on... But the problem here is that time is increasing as O(N^2). I have found lots of talks over the web, that time can be reduced to O(N log N) what can be done by using various trees, sorting, etc. but I can't find anywhere programming code in any languages how this think can be treated. If anyone find or know something about this (not just Csharp but even in other languages, like C++, Fortran, etc.), it would be a great pleasure!
If I get you correctly, you have a group of points of which you want to determine which two are closest to each other? If so, I would have a global list of the points. I would have each point access this list and compare points to itself. I would then report back each of these distances to a manager. The manager will then sort thru to find the smallest distance and which two points it involved. Your approach seems kind if cluttered.
Hello I think you could save some O by not checking every point against every point. I.e. if you have checked the first point against every other point there is no reason to check any more points against that first point. Then you have checked the second point against all other points and there is no reason to ..... etc. I changed your code. I did not check the O or the code. Code (csharp): for(i=1;N;i++){ for(j=i;N;j++){ // Changed 1 to i here r(i,j)=(x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2 if(r(i,j)<rmin){ rmin=r(i,j) x1=x(i) x2=x(j) y1=y(i) ..... } } } Patrik Edit: I did some calculations. I find that it cuts the number of calculations in half but is still O(N^2). Edit Again: I see now that you wanted closest neighbour for each particle. A little different than the problem I solved. I still have some cool code in my next post that might at least give something to work off.
I was a little disappointed so I started matlab and invented something. Here is the code for the O(N^2) algorithm followed by my invented algorithm. For N = 100000 the O(N^2) code takes 162 seconds and the improved code takes 0.2 seconds. Code (csharp): clear n = 100000; x = rand(n,3); minD2 = 1; k = 0; tic for i = 1:n for j = i+1:n d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < minD2 minD2 = d; r = [i,j]; end k = k + 1; end end toc minD2 r k minD2 = 1; k = 0; tic x(1:n,4) = 1:n; x = sortrows(x); x(n+1,1) = 2; for i = 1:n j = i + 1; while (x(i,1)-x(j,1))^2 < minD2 d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < minD2 minD2 = d; r = [x(i,4),x(j,4)]; end k = k + 1; j = j + 1; end end toc minD2 r k The output was Code (csharp): >> NearestNeighbour Elapsed time is 161.684106 seconds. minD2 = 9.3041e-08 r = 60246 96101 k = 5.0000e+09 Elapsed time is 0.198043 seconds. minD2 = 9.3041e-08 r = 60246 96101 k = 3517366 It depends on you having an efficient sortrows algorithm. I did some numeric calculations on N and time and found something close to O(N*logN). Pretty cool. Patrik Edit: This code finds the shortest distance between any two particles. In my next post I find all closest neighbours.
This will find all closest neighbours of all particles. For N = 100000 the O(N^2) code is 293 seconds and improved code is 13 seconds. I shall do my own stuff now Code (csharp): clear n = 10000; x = rand(n,3); k = 0; tic for i = 1:n minD2 = 2; for j = 1:n if i ~= j d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < minD2 minD2 = d; r(i,1) = j; r(i,2) = minD2; end k = k + 1; end end end toc r1 = r; k clear r k = 0; tic x(1:n,4) = 1:n; x = sortrows(x); x = [-4 0 0 0; x]; x(n+2,1) = 4; for i = 2:n+1 j = i + 1; minD2 = 3; while (x(i,1)-x(j,1))^2 < minD2 d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < minD2 minD2 = d; r(x(i,4),1) = x(j,4); r(x(i,4),2) = minD2; end k = k + 1; j = j + 1; end j = i - 1; while (x(i,1)-x(j,1))^2 < minD2 d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < minD2 minD2 = d; r(x(i,4),1) = x(j,4); r(x(i,4),2) = minD2; end k = k + 1; j = j - 1; end end toc r2 = r; k % [r1 r2]
Thanks very much Kavorka for sharing your knowledge! It looks great. By the way, I was wondering if you have complete code of this script (with include, use and other statements - if there are such ones), which I could quickly copy and run to test with other numbers (lets say N=7000,8000,9000,10000,11000,12000) to check time increment profile over N? And which compiler do you use - that I could easily check it? I am usually using heapsort, which I tested many times and it works as O(N log N). Examples can be used from here: http://rosettacode.org/wiki/Sorting_algorithms/Heapsort#C.23 This sort is actually unstable, but it means that it can be not sorting values which are exactly the same. But as my coordinates are always a bit different, heapsort seems to be the best way. Which sorting is used by this sortrows(x) - I see it should be sorting by all three coordinates, right?
Hello No it sorts only on the first coordinate. That is the complete script. Matlab is an interpreted language so there is no compiler. The sortrows() is some built in function and I have no idea of its algorithm. You told me I could use any language so I used the one I always use for calculations. Matlab has extended support for manipulating vectors like I do on rows 27-30. I am sure you can figure out a way to port it to your language of choise. I don't have to declare variables in matlab, they are created by assigning to them.
Oh, ok, I will try to translate it to my language and will check out how sortrows() actually works. So if this is correct the example of coordinates: 3 8 1 1 4 9 2 1 3 should be sorted like that: 1 4 9 2 1 3 3 8 1 is that correct?
well, problem seems to be almost solved. By the way, as we are still fresh memory on this question, what you think would be the most efficient way in this code to find m-nearest neighbours (for example 2, 3-nearest neighbours for each point)?
I just translated this code to Fortran 90 and plotted time(N) law for Kavorka algorithm. program main integer :: nn,n,m,k double precision, dimension), allocatable :: x,y,z,r integer, dimension), allocatable :: index,ir double precision :: d, minD2, minD2a double precision :: t1, t2 integer :: useMod useMod = 1 nn=200 open(1,file="rez.dat",status="replace") open(2,file="rez2.dat",status="replace") write(1,*) '#index ','jj ','x ','y ','z ','x2 ','y2 ','z2 ','ir ','r ' write(2,*) '#n ','t ' do jj=1,nn print *, jj n = 100*jj print *, '1' allocate(index(n+2)) allocate(x(n+2)) allocate(y(n+2)) allocate(z(n+2)) allocate(ir(n+2)) allocate(r(n+2)) print *, '2' call cpu_time(t1) index(1)=1 x(1)=-10.0 y(1)=0.0 z(1)=0.0 index(n+2)=n+2 x(n+2)=10.0 y(n+2)=0.0 z(n+2)=0.0 print *, '3' do i =2,n+1 index(i)=i x(i)=rand() y(i)=rand() z(i)=rand() end do print *, '4' if(useMod.EQ.0) then do i =2,n+1 minD2a = 10.0 do j =2,n+1 if(i.NE.j) then da = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 if(da.LT.minD2a) then minD2a = da ir(i) = j r(i) = minD2a end if end if end do end do end if print *, '5' if(useMod.EQ.1) then ! call sort(n,x) print *, '6' call heapsort(x,y,z,index) print *, '7' i=0 j=0 k=0 do i = 1,n+2 index(i)=i end do print *, '8' do i =2,n+1 j = i+1 minD2 = 10.0 do while(((x(i)-x(j))**2).LT.minD2) if(((x(i)-x(j))**2).LT.minD2) then d = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 if (d .LT. minD2) then minD2 = d ir(index(i)) = index(j) r (index(i)) = minD2 end if k = k + 1 j = j + 1 end if end do j=i-1 do while(((x(i)-x(j))**2).LT.minD2) if(((x(i)-x(j))**2).LT.minD2) then d = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 if (d .LT. minD2) then minD2 = d ir(index(i)) = index(j) r (index(i)) = minD2 end if k = k + 1 j = j - 1 end if end do end do print *, '9' end if print *, '10' do i =2,n+1 ! write(1,*) index(i),jj,x(i), y(i), z(i),x(ir(i)), y(ir(i)), z(ir(i)), ir(i), r(i) end do print *, '11' call cpu_time(t2) write(2,*) n, t2-t1 print *, '12' deallocate(index) deallocate(x) deallocate(y) deallocate(z) deallocate(ir) deallocate(r) print *, '13' end do close(1) close(2) contains ! sorting along double precision with integer is in addition subroutine heapsort(a,b,c,ia) double precision, intent(in out) :: a(0 double precision, intent(in out) :: b(0 double precision, intent(in out) :: c(0 integer, intent(in out) :: ia(0 integer :: start, n, bottom double precision :: atemp,btemp,ctemp integer :: itemp n = size(a) do start = (n - 2) / 2, 0, -1 call siftdown(a,b,c, ia, start, n); end do do bottom = n - 1, 1, -1 atemp = a(0) btemp = b(0) ctemp = c(0) itemp = ia(0) a(0) = a(bottom) b(0) = b(bottom) c(0) = c(bottom) ia(0) = ia(bottom) a(bottom) = atemp; b(bottom) = btemp; c(bottom) = ctemp; ia(bottom) = itemp; call siftdown(a,b,c, ia, 0, bottom) end do end subroutine heapsort subroutine siftdown(a,b,c, ia, start, bottom) double precision, intent(in out) :: a(0 double precision, intent(in out) :: b(0 double precision, intent(in out) :: c(0 integer, intent(in out) :: ia(0 integer, intent(in) :: start, bottom integer :: child, root double precision :: atemp,btemp,ctemp integer :: itemp root = start do while(root*2 + 1 < bottom) child = root * 2 + 1 if (child + 1 < bottom) then if (a(child) < a(child+1)) child = child + 1 end if if (a(root) < a(child)) then atemp = a(child) btemp = b(child) ctemp = c(child) itemp = ia(child) a(child) = a (root) b(child) = b (root) c(child) = c (root) ia(child) = ia (root) a(root) = atemp b(root) = btemp c(root) = ctemp ia(root) = itemp root = child else return end if end do end subroutine siftdown end program main
Wow, I am a little impressed. I wrote the original code and I still have a hard time reading that Fortan. I was a little worried that the part where I extend the size of the vector would give you some trouble, but you apparently had none. Result is not quite O(NlogN), is it? But at least it is a lot better than O(N^2), see picture. I can't judge, are you somewhat pleased or disappointed.
If the algorithm is not the fastest available I excuse myself with the fact that I developed it in an hour. I noticed you changed some things in the code to suit your ensemble, that really tells me you understand what you are doing. This will find the m nearest neighbours of all n particles and their relative distances. Code (csharp): clear n = 10000; m = 10; % m < n x = rand(n,3); clear r k = 0; tic x(1:n,4) = 1:n; x = sortrows(x); x = [-4 0 0 0; x]; x(n+2,1) = 4; for i = 2:n+1 j = i + 1; rtemp = ones(m,2) * 3; b = 3; while (x(i,1)-x(j,1))^2 < b d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < b rtemp = [rtemp(1:m,1:2); d x(j,4)]; rtemp = sortrows(rtemp); b = max(rtemp(1:m,1)); end k = k + 1; j = j + 1; end r(x(i,4),1:m,1:2) = rtemp(1:m,1:2); j = i - 1; while (x(i,1)-x(j,1))^2 < b d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; if d < b rtemp = [rtemp(1:m,1:2); d x(j,4)]; rtemp = sortrows(rtemp); b = max(rtemp(1:m,1)); end k = k + 1; j = j - 1; end r(x(i,4),1:m,1:2) = rtemp(1:m,1:2); end toc % r k % k is not essential to the script. It is number of particles checked. Checking for 5 nearest neighbours of 6 particles yields: Code (csharp): >> NearestNeighbours Elapsed time is 0.007401 seconds. r(:,:,1) = 0.2305 0.4006 0.4595 0.7971 1.2752 0.2305 0.4034 0.4647 0.9743 1.0744 0.2411 0.3653 0.4034 0.7971 1.0028 0.3653 0.5359 0.7394 1.0744 1.2752 0.2411 0.4006 0.4127 0.4647 0.5359 0.4127 0.4595 0.7394 0.9743 1.0028 r(:,:,2) = 2 5 6 3 4 1 3 5 6 4 5 4 2 1 6 3 5 6 2 1 3 1 6 2 4 5 1 4 2 3 k = 30 Patrik
I am pretty happy with results and here I have a plot which explains everything. Unimproved algorithm seems to be working on O(N^2), while improved algorithm could be working on O(NlogN), while sorting algorithm could be working on O(N). If your law would be correct they could be working on O(N^2), O(N^1.67) and O(NlogN), so which version is correct? I am also curious what law we will get if m is getting closer to n? Do we come back to O(N^2) or not?
well, it seems like you are right and this algorithm works with O(N^1.67), as I plotted const*x*log10(x) together with sort law on this plot.
Intreresting work. What will you use this for? NlogN would have been nice but this works, and is fairly easy to implement. Maybe you could improve it furter or find the NlogN algorithm. I have assumed that m << n but it works also for m = n-1 but I don't know at what O. I guess it is not good for large m.
I found some stunning results - performance in this algorithm depends on how many space dimensions we use (I just commented z component when calculating d in 2D space (x,y plane) and y component in 1D space). It seems like laws are the following: 3D O(N^1.67) 2D O(N^1.33(?)) 1D O(N^1.00(?)) so in 1D we really have O(NlogN) (not O(N)) if we account sorting. So if we repeat 1D case for y and z coordinates separately, maybe we can get something even more close to O(NlogN). But I am not sure if it would be possible to combine these results for final r. Do you have any ideas about that? Here is that updated plot again for 2D and 1D. Oh, I am creating RTS Battle system, which looks for nearest neighbour (NN) to attack: http://forum.unity3d.com/threads/223334-3-9k-units-battle-system it currently involves the 'inefficient' i-j looping to get NN. I also included skipping j loops frequency based i-NN from previous r calculations. But everything got dependent on NN problem. So if we would manage to get O(NlogN) - it would be endless number of units used in my system compared to what it is now.
Cool. I was just thinking that if my calculations are correct it would be O(N^1.5) in 2D. I calculated (again not by book but by thinking about it) that it should be O(N^(2-1/d)) where d is number of dimensions. I cannot think of an easy way to combine three 1D runs. But then again nothing I say here have been published in Nature so consider everuthing I say as potentially wrong. RTS is fairly flat? Then 2D could be good enough? I am also feeling something about a flat 3D problem maybe could be better than O(N^1.67). Post if you find more interesting results and/or make improvements. P.S. In 1D you just need sort plus checking two scalars, right? Did you try it for 1D? Is that actual times?
Saw your system, cool. Now that I have seen the application I may have something. Since you have more information than I assumed, how about this: Get the m closest once. Next time look for closest neighbours among your closest neihbours and their closest neighbours. Statistically, with a well chosen m, it will be correct most of the time. This, I beleive could be about O(NM^2). Which for fixed m is O(N). That would be something! Am I right, is it doable?
well, I still have a hope to try here so called k-d tree which is binary search tree and is searching in k-dimensional space: http://en.wikipedia.org/wiki/Kd-tree quite good implementation seems to be implemented in C++ here: https://code.google.com/p/kdtree/downloads/list I still didn't find out how to run this code and what it gives. Another problem is that it could be not so easy to translate it from C++ to C#, as Unity uses C#. Yes, RTS is pretty flat and I can use 2D in case if there would be no other options. However as I am building 3D game, I would start experiencing problems in mountain areas. I am also thinking about using these algorithms in collision detecting problems, which is now working I guess as O(N^2). I don't like to apply current collision detection also because it applies force along r if r is close enough for NN instead of keeping NN at the constant r if critical value of r is reached. So there are several cases where I would apply NN problem. well, I assuming that in 1D the only change in your code would be that instead of d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2; you put d = (x(i,1)-x(j,1))^2; everything else remaining the same.