Can someone help us with choosing the best syntax and nomenclature for this code?

I have this code we are trying to run to create vectors from n closest points. However, it seems to be getting too complicated. Can someone help us to have a cleaner, more efficient code given the attched file (fpep.mat)? Below is the current code. Thanks in advance for your help.

close all;
clearvars;
% load('F_points.mat');
load('fpep.mat');
n = 10;
for pt = 1 : 953
        dist = sqrt((fpep(:,7)-fpep(pt,7)).^2 + (fpep(:,8)-fpep(pt,8)).^2); 
        [~, ascendIdx] = sort(dist); % Not really sure what this does or why there is a ~ used here. Can someone explain this?
        xyNearest(pt,:,:) = fpep(ascendIdx(1:n),7:8);
        adirector_x(pt,:,:) = ([fpep(ascendIdx(1:n),1) - xyNearest(pt,:,1)]);
        adirector_y(pt,:,:) = ([fpep(ascendIdx(1:n),2) - xyNearest(pt,:,2)]);
        bdirector_x(pt,:,:) = ([fpep(ascendIdx(1:n),3) - xyNearest(pt,:,1)]);
        bdirector_y(pt,:,:) = ([fpep(ascendIdx(1:n),4) - xyNearest(pt,:,2)]);
        cdirector_x(pt,:,:) = ([fpep(ascendIdx(1:n),5) - xyNearest(pt,:,1)]);
        cdirector_y(pt,:,:) = ([fpep(ascendIdx(1:n),6) - xyNearest(pt,:,2)]);
        ascendIdx(ascendIdx==pt) = [];  %remove the source point. Not sure why the open brackets are used here. Can someone explain this?
end

chordx = ([xyNearest(:,:,1) - fpep(:,7)]);
chordy = ([xyNearest(:,:,2) - fpep(:,8)]);
chord = sqrt((chordx).^2 + (chordy).^2);

[Rik]

This time I edited your question for you. Next time, please use the tools explained on this page to make your question more readable.

As to your question: I don't think your code will become more clear if you rewrite this.

For the questions in your comments: the tilde supresses an output you're not going to use. So in this case your code only uses the second output of the sort function. The variable name gives a hint what that is, but if you want to know for sure you need to read the documentation.

Assigning [] to something is setting an empty array. That will remove the elements from an array. That is basic Matlab syntax that any crash course should teach you. However, since the variable is not used in the next loop iteration (where it is overwritten), this doesn't actually do something usefull.

Also, you should load to a struct. That way it is clear where your variables come from and would prevent the call to clearvars. (and because you're not opening any figures here, the close all isn't required either)

[Steve]

Thanks for your feedback Rik. The problem is when I expand on the code in its current configuration, I run into issues--for example, I'm having an problem with A & B being different sizes on the dot product; this is one reason why I think a simpler code would be (at least) much easier to work with and debug (at least for a MATLAB novice like me). *Please see below:

close all;
clearvars;
% load('F_points.mat');
load('fpep.mat');
n = 10;
for pt = 1 : 953
        dist = sqrt((fpep(:,7)-fpep(pt,7)).^2 + (fpep(:,8)-fpep(pt,8)).^2); 
        [~, ascendIdx] = sort(dist);
        ascendIdx(ascendIdx==pt) = [];  %remove the source point
        xyNearest(pt,:,:) = fpep(ascendIdx(1:n),7:8);
        adirector_x(pt,:,:) = ([fpep(ascendIdx(1:n),1) - xyNearest(pt,:,1)]);
        adirector_y(pt,:,:) = ([fpep(ascendIdx(1:n),2) - xyNearest(pt,:,2)]);
        bdirector_x(pt,:,:) = ([fpep(ascendIdx(1:n),3) - xyNearest(pt,:,1)]);
        bdirector_y(pt,:,:) = ([fpep(ascendIdx(1:n),4) - xyNearest(pt,:,2)]);
        cdirector_x(pt,:,:) = ([fpep(ascendIdx(1:n),5) - xyNearest(pt,:,1)]);
        cdirector_y(pt,:,:) = ([fpep(ascendIdx(1:n),6) - xyNearest(pt,:,2)]);   
end


    chordx = ([xyNearest(:,:,1) - fpep(:,7)]);
    chordy = ([xyNearest(:,:,2) - fpep(:,8)]);
    chord = sqrt((chordx).^2 + (chordy).^2);
    adirector = sqrt((adirector_x).^2 + (adirector_y).^2);
    bdirector = sqrt((bdirector_x).^2 + (bdirector_y).^2);
    cdirector = sqrt((cdirector_x).^2 + (cdirector_y).^2);
    CosTheta_a = dot(adirector,chord)/(norm(adirector)*norm(chord));
    CosTheta_b = dot(bdirector,chord)/(norm(bdirector)*norm(chord));
    CosTheta_c = dot(cdirector,chord)/(norm(cdirector)*norm(chord));

*If anyone can help me to properly find (and store) all the angles between the vectors: "adirector", "bdirector", "cdirector" and the vector "chord" (for all the sets of points) it would be very much appreciated.

[Rik]

I have no clue what you are trying to do. Your explanation is a bit too short and general for me. The only thing I notice when running this code is that your first input to dot is 953x10x10, while the second is 953x10(x1). I'm not sure if you want to discard the third dimension, or if you wanted to have the same implicit expansion as you have in your for-loop (where you substract 1xn from nx1).

Anyway, below is my suggestion to slightly clean up your code. It still lacks a lot of comments, but that is up to you to write.

%load input matrix
S=load('fpep.mat');fpep=S.fpep;
n = 10;

%pre-allocate loop outputs
xyNearest=  zeros( size(fpep,1) , n , 2);
adirector_x=zeros( size(fpep,1) , n , n);
adirector_y=zeros( size(fpep,1) , n , n);
bdirector_x=zeros( size(fpep,1) , n , n);
bdirector_y=zeros( size(fpep,1) , n , n);
cdirector_x=zeros( size(fpep,1) , n , n);
cdirector_y=zeros( size(fpep,1) , n , n);

for pt = 1 : size(fpep,1)
    dist = sqrt((fpep(:,7)-fpep(pt,7)).^2 + (fpep(:,8)-fpep(pt,8)).^2);
    [~, ascendIdx] = sort(dist);
    ascendIdx(ascendIdx==pt) = [];  %remove the source point
    xyNearest(pt,:,:) = fpep(ascendIdx(1:n),7:8);
    adirector_x(pt,:,:) = fpep(ascendIdx(1:n),1) - xyNearest(pt,:,1);
    adirector_y(pt,:,:) = fpep(ascendIdx(1:n),2) - xyNearest(pt,:,2);
    bdirector_x(pt,:,:) = fpep(ascendIdx(1:n),3) - xyNearest(pt,:,1);
    bdirector_y(pt,:,:) = fpep(ascendIdx(1:n),4) - xyNearest(pt,:,2);
    cdirector_x(pt,:,:) = fpep(ascendIdx(1:n),5) - xyNearest(pt,:,1);
    cdirector_y(pt,:,:) = fpep(ascendIdx(1:n),6) - xyNearest(pt,:,2);
end

chordx = xyNearest(:,:,1) - fpep(:,7);
chordy = xyNearest(:,:,2) - fpep(:,8);
chord = sqrt((chordx).^2 + (chordy).^2);
adirector = sqrt((adirector_x).^2 + (adirector_y).^2);
bdirector = sqrt((bdirector_x).^2 + (bdirector_y).^2);
cdirector = sqrt((cdirector_x).^2 + (cdirector_y).^2);
CosTheta_a = dot(adirector,chord)/(norm(adirector)*norm(chord));%probably needs ./ and .* instead
CosTheta_b = dot(bdirector,chord)/(norm(bdirector)*norm(chord));%probably needs ./ and .* instead
CosTheta_c = dot(cdirector,chord)/(norm(cdirector)*norm(chord));%probably needs ./ and .* instead

[Steve]

Thanks again so much for your help Rik. I am not really sure why the two vectors are different sizes, or (more importantly), how to get them to be the same size. Below are pictures of what I'm trying to achieve. Hope this helps.

[Rik]

How did you generate this plot? That might help in getting what you mean.

[Steve]

The "zoomed in" plot was manufactured in Paint. However, the other plot was obtained with the following code:

close all;
clearvars;
load('F_points.mat');
load('fpep.mat');
%xy = cell2mat(F_points)';  %[n x 2] matrix of (x,y) coordinates
xy = F_points';
pt = 499; % choose a point in xy and we'll find the 10 nearest neighbors. 
% Euclidean distance between xy(p,:) and all other points
dist = sqrt((xy(:,1)-xy(pt,1)).^2 + (xy(:,2)-xy(pt,2)).^2); 
% Find the n closest values (excluding the point selected)
n = 10;
dist_ep = sqrt((fpep(:,1)-xy(pt,1)).^2 + (fpep(:,2)-xy(pt,2)).^2);
[Z,I] = sort(dist);
[~, ascendIdx_ep] = sort(dist_ep);
ascendIdx_ep(ascendIdx_ep==pt) = [];  %remove the pt point
epNearest = fpep(ascendIdx_ep(1:n),:);
ptNearest = xy(I(1:n),:);
[~, ascendIdx] = sort(dist); 
ascendIdx(ascendIdx==pt) = [];  %remove the pt point
xyNearest = xy(ascendIdx(1:n),:);
IZ = ([I Z xy(I,1) xy(I,2)]); %% This shows the index number associated with each of the closest Fermat points
IZNearest = IZ(I(1:n),:);
% fpepNearest2 = fpep(I(1:n),:);
chordx = ([xyNearest(:,1) - xy(pt,1)]);
chordy = ([xyNearest(:,2) - xy(pt,2)]);
chord = sqrt((chordx).^2 + (chordy).^2);
director_x = ([fpep(:,1) - xy(pt,1)]);
director_y = ([fpep(:,2) - xy(pt,2)]);
director = sqrt((director_x).^2 + (director_y).^2);
figure()
for ii = 1 : length(chordx)
    xvals = [xy(pt,2),xy(pt,2)+chordy(ii)];
    yvals = [xy(pt,1),xy(pt,1)+chordx(ii)];
plot(xvals,yvals,'r-')
hold on
end
plot(xy(:,2),xy(:,1),'b.')
hold on
% Show the selected point
plot(xy(pt,2),xy(pt,1),'b.','MarkerFaceColor', 'y')
% Show nearest 'n' dots
plot(xyNearest(:,2),xyNearest(:,1),'r+')
plot(fpep(:,2),fpep(:,1),'k.')
plot(fpep(:,4),fpep(:,3),'k.')
plot(fpep(:,6),fpep(:,5),'k.')
axis equal