How can I speed up my Code?

Good evening,

I have a problem with the performance of my Code, as Im new to matlab.

I hope the comments help you to understand everything.

The main problem is, that most of the time I have no points at the x-coordinate. Thats why I interpolate the closest and the secound closest. But that I have to do for every AeroDataId, so every Mach/CL combination in my database. To speed it up a little bit, I only use the points in a certin mach and Cl delta.

load('new_upper_airfoil');


%% searched point
Ma=0.7;
C_l=0.65;

%% delta of cl and Mach
cl_delta=0.015;
Ma_delta=0.015;

%% delta of the points in the original pressure distribution
Delta=2/height(originalDV);

DV=zeros(height(originalDV),2);
Temp=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
Temp_update=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
i=1;
%% load only the points within a Mach Number and Lift coefficient delta
new_upper_airfoil=new_upper_airfoil(new_upper_airfoil.Ma < (Ma+Ma_delta) & new_upper_airfoil.Ma > (Ma-Ma_delta) & new_upper_airfoil.CL_res < (C_l+cl_delta) & new_upper_airfoil.CL_res > (C_l-cl_delta),:);
new_lower_airfoil=new_lower_airfoil(new_lower_airfoil.Ma < (Ma+Ma_delta) & new_lower_airfoil.Ma > (Ma-Ma_delta) & new_lower_airfoil.CL_res < (C_l+cl_delta) & new_lower_airfoil.CL_res > (C_l-cl_delta),:);

%% Interpolation of upper_airfoil

for x_c=0:Delta:1
    %% Vector with all unique Ids (muss nicht nochmal für Unterseite erstellt werden)
    DataId_unique=unique(new_upper_airfoil.AeroDataId); %könnte gespeichert werden für alle Daten
   
    for a=1:1:size(DataId_unique)
        %% temp table with all the information for one ID (AeroData ID is no longer needed)
        new_upper_airfoil_unique=new_upper_airfoil(new_upper_airfoil.AeroDataId==DataId_unique(a,:),2:5); % könnte später gespeichert werden
        
        %% try to find x/c
        FIND=new_upper_airfoil_unique(new_upper_airfoil_unique.X_C==x_c,:);    
        
        if(isempty(FIND))
            %% find the closest and secound closest x-coordinate and interpolate between them
            [d,sortIndex]=sort(abs(new_upper_airfoil_unique.X_C-x_c));
            
            closest=new_upper_airfoil_unique(sortIndex(1),:);
            secound_closest=new_upper_airfoil_unique(sortIndex(2),:);
            
            average=interp1([closest.X_C,secound_closest.X_C], [closest.CP,secound_closest.CP], x_c, 'linear','extrap');
            
            %% saved with closest Ma and Cl as they are the same as in secound closest
            Temp=[Temp;{x_c, average, closest.Ma, closest.CL_res}];
        else
            Temp=[Temp;FIND];
        end  
        
    end
    %% Interpolation
    cp_inter=IDW(Temp.Ma, Temp.CL_res, Temp.CP, Ma, C_l, -2, 'ng', 4);
    
    Temp(:,:)=[];
    %% Result
    DV(i,:)=[x_c,cp_inter];
    i=i+1; 
end

[Philippe Lebel]

could you update the question with the file that the following line loads please?

load('new_upper_airfoil');

In this code is there any funcitons that you made yourself or this is all standard functions? (will i be able to run this on my computer)

[Gerrit]

All standard function except IDW. Its an interpolation method.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% INVERSE DISTANCE WEIGHT %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Vint]=IDW(xc,yc,vc,x,y,e,r1,r2)
%%%%%%%%%%%%%%%%%
%%% INPUTS
%xc = stations x coordinates (columns) [vector]
%yc = stations y coordinates (rows) [vector]
%vc = variable values on the point [xc yc]
%x = interpolation points  x coordinates [vector]
%y = interpolation points y coordinates [vector]
%e = distance weight
%r1 --- 'fr' = fixed radius ;  'ng' = neighbours
%r2 --- radius lenght if r1 == 'fr' / number of neighbours if  r1 =='ng'
%%% OUTPUTS
%Vint --- Matrix [length(y),length(x)] with interpolated  variable values
%%% EXAMPLES
%%% --> V_spa=IDW(x1,y1,v1,x,y,-2,'ng',length(x1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Simone Fatichi -- [e-mail]
%   Copyright 2009
%   $Date: 2009/06/19 $ 
%   $Updated 2012/02/24 $ 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vint=zeros(length(y),length(x));
xc=reshape(xc,1,length(xc));
yc=reshape(yc,1,length(yc));
vc=reshape(vc,1,length(vc));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if  strcmp(r1,'fr')                                                         % compareStrings, fixed radius oder nearest neighbour
    if  (r2<=0)                                                             % Kontrolle, ob Radius positiv 
        disp('Error: Radius must be positive')
        return
    end
    for i=1:length(x)                                                       
        for j=1:length(y)
            D=[]; V=[]; wV =[]; vcc=[];
            D= sqrt((x(i)-xc).^2 +(y(j)-yc).^2);                            % Euklidische Norm des Abstands zwischen bekannten und gesuchten Punkten
            if min(D)==0                                                    % Kontrolle, ob Abstand bei einem bekannten Punkt 0 wird
                disp('Error: One or more stations have the coordinates of an interpolation point')
                return
            end
            vcc=vc(D<r2);                                                    % wählt vc-Werte der Wertepaare aus, bei denen D<r2 gilt 
            D=D(D<r2);                                                      % Euklidische Norm, die kleiner als der gesetzte Radius ist; alle dazugehörigen Z-Werte  
            V = vcc.*(D.^e);                                                % Zähler 
            wV = D.^e;                                                      % Nenner 
            if isempty(D)                                                   % Wenn Abstandsvektor D leer ist V=NaN   
                V=NaN;                                                      
            else
                V=sum(V)/sum(wV);                                           % Enformel, Summe über alle Einträge
            end
            Vint(j,i)=V;                                                    % Ergebnis
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else                                                                        % r1 ist gleich 'ng'
    if (r2 > length(vc)) || (r2<1)                                          % Kontrolle, ob so viele Nachbarn wie gesetzt verfügbar sind
        disp('Error: Number of neighbours not congruent with data')
        return
    end
    for i=1:length(x)           
        for j=1:length(y)
            D=[]; V=[]; wV =[];vcc=[];                                      % erstellen der Vektoren
            D= sqrt((x(i)-xc).^2 +(y(j)-yc).^2);                            % Berechnen der euklidischen Norm
            if min(D)==0                                                    % Kontrolle, ob gesuchter Punkt in der Menge enthalten 
                disp('Error: One or more stations have the coordinates of an interpolation point') 
                return
            end
            [D,I]=sort(D);                                                  % sortiert D nach aufsteigender Reihenfolge, I ist Reihenfolge von vorher                                                
            vcc=vc(I);                                                      % vcc ist der Vektor mit den Z- Werten der Nachbarn in der Reihenfolge wie D
            V = vcc(1:r2).*(D(1:r2).^e);                                    % vcc(1:r2): Vektor mit z-Werten der nächsten Nachbarn; V = Zähler                                
            wV = D(1:r2).^e;                                                % Nenner der Endformel
            V=sum(V)/sum(wV);                                               % Endformel, Summe über alle Einträge
            Vint(j,i)=V;                                                    % Ergebnis
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return

[Gerrit]

I use another function now, so its save that I use the left and right neighours and not extrapolate.

t0=tic();


load('new_upper_airfoil');

load('new_lower_airfoil');



%% Eingabe der Parameter

Ma=0.7;
C_l=0.65;

% x_delta=0.008;
cl_delta=0.02;
Ma_delta=0.02;

%% Interpolation of upper_airfoil
% Delta=2/height(originalDV);
DV=zeros(202,2);
Temp=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
Temp_update=array2table(zeros(0,4),'VariableNames',{'X_C','CP','Ma','CL_res'});
i=1;

%% load only the points within a Mach Number and Lift coefficient delta
new_upper_airfoil=new_upper_airfoil(new_upper_airfoil.Ma < (Ma+Ma_delta) & new_upper_airfoil.Ma > (Ma-Ma_delta) & new_upper_airfoil.CL_res < (C_l+cl_delta) & new_upper_airfoil.CL_res > (C_l-cl_delta),:);
new_lower_airfoil=new_lower_airfoil(new_lower_airfoil.Ma < (Ma+Ma_delta) & new_lower_airfoil.Ma > (Ma-Ma_delta) & new_lower_airfoil.CL_res < (C_l+cl_delta) & new_lower_airfoil.CL_res > (C_l-cl_delta),:);


for x_c=0:0.01:1
    %% Vector with all unique Ids (muss nicht nochmal für Unterseite erstellt werden)
    DataId_unique=unique(new_upper_airfoil.AeroDataId); %könnte gespeichert werden für alle Daten
    t_start=tic();
    for a=1:1:size(DataId_unique)
        %% temp table with all the information for one ID (AeroData ID is no longer needed)
        new_upper_airfoil_unique=new_upper_airfoil(new_upper_airfoil.AeroDataId==DataId_unique(a,:),2:5); % könnte später gespeichert werden
        
        %% try to find x/c
        FIND=new_upper_airfoil_unique(abs(new_upper_airfoil_unique.X_C-x_c)<0.0005,:); 
        
        if(isempty(FIND))
            %% find left and right neighbour
            leftNeighbourIndex=nearestpoint(x_c,new_upper_airfoil_unique.X_C,'previous');
            rightNeighbourIndex=nearestpoint(x_c,new_upper_airfoil_unique.X_C,'next');
            
            if isnan(leftNeighbourIndex)
                leftNeighbourIndex=rightNeighbourIndex+1;
            end
            if isnan(rightNeighbourIndex)
                rightNeighbourIndex=leftNeighbourIndex-1;
            end
            %% interpolation of the coordinate
            average=interp1([new_upper_airfoil_unique.X_C(leftNeighbourIndex,:),new_upper_airfoil_unique.X_C(rightNeighbourIndex,:)], [new_upper_airfoil_unique.CP(leftNeighbourIndex,:),new_upper_airfoil_unique.CP(rightNeighbourIndex,:)], x_c, 'linear','extrap');
            new_upper_airfoil_unique.X_C(1,:)=x_c;
            new_upper_airfoil_unique.CP(1,:)=average; 
            
            Temp=[Temp;new_upper_airfoil_unique(1,:)];
        else
            FIND.X_C=x_c;
            Temp=[Temp;FIND];
        end  
        
    end
    t_stuetz=toc(t_start);
    fprintf('Dauer für Stützstellen: %6.5f s\n',t_stuetz);
     
    cp_inter=IDW(Temp.Ma, Temp.CL_res, Temp.CP, Ma, C_l, -2, 'ng', 4);

    Temp(:,:)=[];
    DV(i,:)=[x_c,cp_inter];
    i=i+1;
    
end

nearestpoint function: https://de.mathworks.com/matlabcentral/fileexchange/8939-nearestpoint-x-y-m

[Adam]

Have you used the profiler

doc profile

to understand where it is slow? Should always be the first step before doing speedup work. No point speeding up something that takes only 1% of the total run time anyway!

[Gerrit]

No I didnt. The function are not slow. I thnik the only way to get some speed, would be to somehow vectorize the loop. It just takes so long, because there is a loop, which adds up all the small time amounts of the function.

[darova]

First of all: pay attention to preallocation

You didn't attached new_lower_airfoil.mat

[Gerrit]

Hi,

you do not need this file to run this part of the code. Preallocation is slower in this case.

[Stephen23]

"Preallocation is slower in this case."

Please upload the code that you tried.

[Gerrit]

The problem here is, that Temp is getting filled up again with every loop step.