Reputation: 307
I have a set of n
complex numbers that move through the complex plane from time step 1
to nsampl
. I want to plot those numbers and their trace over time (y-axis shows imaginary part, x-axis the real part). The numbers are stored in a n x nsampl
vector. However in each time step the order of the n
points is random. Thus in each time step I pick a point in the last time step, find its nearest neighbor in the current time step and put it at the same position as the current point. Then I repeat that for all other n-1
points and go on to the next time step. This way every point in the previous step is associated with exactly one point in the new step (1:1 relation). My current implementation and an example are given below. However my implementation is terribly slow (takes about 10s for 10 x 4000 complex numbers). As I want to increase both, the set size n
and the time frames nsampl
this is really important to me. Is there a smarter way to implement this to gain some performance?
Example with n=3 and nsampl=2:
%manually create a test vector X
X=zeros(3,2); % zeros(n,nsampl)
X(:,1)=[1+1i; 2+2i; 3+3i];
X(:,2)=[2.1+2i; 5+5i; 1.1+1.1i]; % <-- this is my vector with complex numbers
%vector sort algorithm
for k=2:nsampl
Xlast=[real(X(:,k-1)) imag(X(:,k-1))]; % create vector with x/y-coords from last time step
Xcur=[real(X(:,k)) imag(X(:,k))]; % create vector with x/y-coords from current time step
for i=1:size(X,1) % loop over all n points
idx = knnsearch(Xcur(i:end,:),Xlast(i,:)); %find nearest neighbor to Xlast(i,:), but only use the points not already associated, thus Xcur(i:end,:) points
idx = idx + i - 1;
Xcur([i idx],:) = Xcur([idx i],:); %sort nearest neighbor to the same position in the vector as it was in the last time step
end
X(:,k) = Xcur(:,1)+1i*Xcur(:,2); %revert x/y coordinates to a complex number
end
Result:
X(:,2)=[1.1+1.1i; 2.1+2i; 5+5i];
Can anyone help me to speed up this code?
Upvotes: 2
Views: 532
Reputation: 36720
The problem you are tying to solve is combinatorial optimizartion which is solved by the hungarian algorithm (aka munkres
). Luckily there is a implementation for matlab available for download. Download the file and put it either on your search path or next to your function. The code to use it is:
for k=2:size(X,2)
%build up a cost matrix, here cost is the distance between two points.
pairwise_distances=abs(bsxfun(@minus,X(:,k-1),X(:,k).'));
%let the algorithm find the optimal pairing
permutation=munkres(pairwise_distances);
%apply it
X(:,k)=X(permutation,k);
end
Upvotes: 1