Speed up nested loops

I am trying to speed these loops up but I couldn't. please help me.

this code to calculate the shortest distance between lines and the distance between the mid point of the line and the the point of minimum distance vector between them then compare it with 0, if it is less than zero it will beark the second loop and change the orientation of the new line and check it with ALL previous lines. I tried 2000 lines with 100 trails the PC takes more than 48 hrs.

n=2000;
maxtrial=100;
Rx=rand(n,1)*0.2;
Ry=rand(n,1)*0.2;
Rz=rand(n,1)*0.2;
MT=0;
k=1;
while k<=n
    
    % Create the angels for fiber oreintation
    %For Phi
Phi_min=-0.1745329;   %Angle Phi minimum value
Phi_max=0.1745329; %Angle Phi maximum value
Phi=Phi_min+rand(1,1)*(Phi_max-Phi_min);

%For theta
Z = (-1) + (1-(-1)) * rand(1,1); % value of Z to use in angle theta which will be within +- 10 degree

theta = acos (Z);

d1= sin(theta)*sin(Phi);
d2= sin(theta)*cos(Phi);
d3=cos(theta);

%first point coordinats
x2= Rx(k)+(Lf*0.5*d1);
y2= Ry(k)+(Lf*0.5*d2);
z2= Rz(k)+(Lf*0.5*d3);

%second point coordinats
x3= Rx(k)-(Lf*0.5*d1);
y3= Ry(k)-(Lf*0.5*d2);
z3= Rz(k)-(Lf*0.5*d3);

P=[x3-x2 y3-y2 z3-z2]; %orientation vector
P_all(k,:)=P;

if k>=2 && k<=n
    t=k-1;
    for t=1:k-1
    normal_vector=((cross(P_all(k-t,:),P_all(k,:))/(norm(cross(P_all(k-t,:),P_all(k,:)))))); % unit vector normal to both lines
    
    min_distance= norm(dot((R_all(k-t,:)-R_all(k,:)),normal_vector))-Df; % the minimum distance between two lines
    
    L_ij=(-(dot((R_all(k-t,:)-R_all(k,:)),P_all(k-t,:)))+(dot((R_all(k-t,:)-R_all(k,:)),P_all(k,:)))*dot(P_all(k-t,:),P_all(k,:)))/(1-(dot(P_all(k-t,:),P_all(k,:)))^2); % distance between the center of line and the point that minimum distance occure at
    
    if min_distance<0 && L_ij<=Lf/2
        k=k-1;
        MT=MT+1;
        break
    else
        MT=0;
        
    end
    
    end
end

if MT==maxtrial
        x2=0;
        x3=0;
        y2=0;
        y3=0;
        z2=0;
        z3=0;
   break
end


G1(k,:)= [z2 x2 y2]; %first points coordinates matrix
G2(k,:)= [z3 x3 y3]; %second points coordinates matrix

k=k+1;

end

[Adam]

Have you run the profiler

doc profile

to understand which lines are taking the most time? It's not perfect, given its disabling of certain run-time optimisations, but it usually gives a very good guide for where to focus speedup work.

[Stephen23]

https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

[Steven Lord]

When you profile your program, I recommend doing so for a smaller number of lines and fewer trials. Doing so will not only let you see the performance bottlenecks sooner but will also let you see if there are any places your code assumes the larger number of lines or trials. If you decrease the number of lines to say 200 and receive an error about a line of code trying to access element 201 of an array, you should probably examine that line of code (and the code that creates or manipulates the variables used by that line) more closely.

[Jan]

You see warnings in the editor already, which explain the disadvantages of iteratively growing arrays. Use a proper pre-allocation in any way.

A simple code will not run faster, but it is easier to read: Compare

Z = (-1) + (1-(-1)) * rand(1,1)

with

Z = 2 * rand - 1;

I cannot run your code, because Lf, Df and R_all are undefined. It looks strange, that MT is not reset and comparing MT==maxtrial can match 1 time only.

Please post a running code.

Answer 1

[Jan]

This line is expensive:

normal_vector = ((cross(P_all(k-t,:),P_all(k,:)) / ...
    (norm(cross(P_all(k-t,:),P_all(k,:))))));

Matlab's cross() is not efficient and norm() can be accelerated also. Calling cross() twice for the same data is a waste of time in addition. P_all(k, :) is available as P also. Faster and nicer:

v = cross(P_all(k-t, :), P);
normal_vector = v / norm(v);

Or use a faster M-function:

normal_vector = NormCross(P_all(k-t, :), P);

function c = NormCross(a, b)
c = [a(2) * b(3) - a(3) * b(2), ...
     a(3) * b(1) - a(1) * b(3), ...
     a(1) * b(2) - a(2) * b(1)];
c = c / sqrt(c(1) * c(1) + c(2) * c(2) + c(3) * c(3));
end

Equivalent improvements can be applied to the DOT commands also.

By the way: The 2 leading and trailing parentheses in "normal_vector = ((" decrease the readability of the code only.

A short run time test:

P = rand(100, 3);
t = 100;
tic
for r = 1:1e5
   for k = 1:100
      v = cross(P(k, : ), P(t, :)) / norm(cross(P(k, :), P(t, :)));
   end
end
toc

tic
for r = 1:1e5
   Pt = P(t, :);
   for k = 1:100
      v = NormCross(P(k, : ), Pt);
   end
end
toc

% Elapsed time is 11.122618 seconds.
% Elapsed time is 2.686646 seconds.

4 times faster just by avoiding calling cross twice and norm.

Instead of comparing min_distance=norm(...), you can compare the squared distance to save the computational time for the square root:

v = rand(1, 3);
min_distance = norm(v) - Df;
if min_distance < 0
    ...
end

Df2 = Df ^ 2;
...
v = rand(1, 3);
min_dist2 = v * v' - Df2;  % Faster than: sum(v .* v)
if min_dist2 < 0
    ...
end

[Majeed]

Dear Mr. Jan,

this is the working code as your request,,

clc
close all

% Define the fiber diameter and numbers
Df=200e-6; % fiber diameter
n=100;    % number of fibers
nfiber=n;

% Define the dimension of RVE
L= 250*Df;
W= 250*Df;
H= 250*Df;

% Fiber length in term of Df
Lf= L * 1.5;    % fiber length

%RVE size
RVE_size=L*W*H;

% Define the max. no. of trials
maxtrial=3;


% Create random points inside the RVE for centers of fibers
Rx=rand(n,1)*L;
Ry=rand(n,1)*W;
Rz=rand(n,1)*H;


R_all=[Rx Ry Rz];

%Creats the centers of fibers
j=1;

while j<n
    i=j+1;
    while i<=n

         k  = norm(R_all(j,:)-R_all(i,:)); %function of distance between points
         
         if k < 1.5*Df  % check the distance between all points; should not be < 3 Df
             R_all(i,:)=[];
             n=n-1;
         end
         
         
         i=i+1;
    end
     j=j+1; 
end

% scatter3(Rx, Ry, Rz); %Fibers' centers figure

G1=zeros(n,3);  %first points coordinates matrix
G2=zeros(n,3);  %second points coordinates matrix

P_all(n,3)=0; %direction vectors

counter=0;
MT=0;
k=1;
tic
while k<=n
    
% Create the angels for fiber oreintation
%For Phi
Phi_min=-0.087266463;   %Angle Phi minimum value
Phi_max=0.087266463; %Angle Phi maximum value
Phi=Phi_min+rand(1,1)*(Phi_max-Phi_min);

%For theta
Z = (-1) + (1-(-1)) * rand(1,1); % value of Z to use in angle theta which will be within +- 10 degree

% angle =[acos(Z) -acos(Z)]; % To get +- angles
% pos = randi(length(angle)); % To choose position in previous array
% theta = angle(pos); % To choose angle value

theta = acos (Z);

d1= sin(theta)*sin(Phi);
d2= sin(theta)*cos(Phi);
d3=cos(theta);

%first point coordinats
x2= Rx(k)+(Lf*0.5*d1);
y2= Ry(k)+(Lf*0.5*d2);
z2= Rz(k)+(Lf*0.5*d3);

%second point coordinats
x3= Rx(k)-(Lf*0.5*d1);
y3= Ry(k)-(Lf*0.5*d2);
z3= Rz(k)-(Lf*0.5*d3);

P=[x3-x2 y3-y2 z3-z2]; %orientation vector
P_all(k,:)=P;
% keyboard

if k>=2 && k<=n
    t=k-1;
    for t=1:k-1
    normal_vector=((cross(P_all(k-t,:),P_all(k,:))/(norm(cross(P_all(k-t,:),P_all(k,:)))))); % unit vector normal to both fibers
    
    min_distance= norm(dot((R_all(k-t,:)-R_all(k,:)),normal_vector))-Df; % the minimum distance between two fibers
    
    L_ij=(-(dot((R_all(k-t,:)-R_all(k,:)),P_all(k-t,:)))+(dot((R_all(k-t,:)-R_all(k,:)),P_all(k,:)))*dot(P_all(k-t,:),P_all(k,:)))/(1-(dot(P_all(k-t,:),P_all(k,:)))^2); % distance between the center of fiber and the point that minimum distance occure at
    
    if min_distance<0 && L_ij<=Lf/2

        k=k-1;
        MT=MT+1;
        counter=counter+1;
        break
    else
        MT=0;
    end
    
    end
    
end

if MT==maxtrial
        x2=0;
        x3=0;
        y2=0;
        y3=0;
        z2=0;
        z3=0;
   continue
end


G1(k,:)= [z2 x2 y2]; %first points coordinates matrix
G2(k,:)= [z3 x3 y3]; %second points coordinates matrix

k=k+1;

% line ( [x2,x3] , [y2,y3] , [z2,z3]) %fibers distribution figure

end

toc
 % open text file for writing points for each fiber for spaceclaim
 
fid = fopen('uper_fibers.txt','wt');
fid2= fopen('Normal_fibers.txt','wt');
    for m=1:n
        if G1(m,:)~=0
            if G1(m,3)>0.08 && G2(m,3)>0.08
     fprintf(fid,'%i %i %i\n%d %d %d\n\n',G1(m,:),G2(m,:));
        else
     fprintf(fid2,'%i %i %i\n%d %d %d\n\n',G1(m,:),G2(m,:));       
            end
        end
    end
  fclose(fid);

[Jan]

Did you use the profiler already to find the bottleneck? I offered a function to speed up the computation of norm(cross()). Did you try it? If so with which effect?

[Jan]

@Abdulmajeed Altassan: The runtime is not the problem of the code. When it works, it runs in less than 0.1 seconds. But for some inputs the code produces an infnite loop with unlimited runtime. So my first idea to optimze the code does not hit the point, because the code does not work reliably.

[Majeed]

Dear Jan;

Thank you for your answers, it was jamming with me when I set n=2000 and maxtrial=5 it tooks 48 hr and it dose not finished.

I will follow you idea thank you very much