How can I make this code run faster?

This code takes about 30s to run, and I will need to run it several times

Can anyone help me make it faster?

here is the code:

sum_Load = rand(1,1000);             % For demonstration purposes
Pseudo_Load_Factor = rand(500,1000); % For demonstration purposes
Disp_Factor =  rand(500,1000);       % For demonstration purposes
M = 1000;
N = 500;
Disp_Inter = zeros(M,N);
for j = 1:M
    Load_Inter = sum_Load(j);
    for i = 1:N
        idx = find(cummin(Pseudo_Load_Factor(i,:))<Load_Inter & cummax(Pseudo_Load_Factor(i,:))>Load_Inter,1,'first');
        if isempty(idx)
            Disp_Inter(j,i) = 0.10*7000;
        else
            Disp_Inter(j,i) = interp1(Pseudo_Load_Factor(i,(idx-1):idx),Disp_Factor(i,(idx-1):idx),Load_Inter); 
        end
    end
end

Thanks a lot!!

[Jan]

Please post some running code. What is the contents of sum_Load, Pseudo_Load_Factor, cumin, Load_Inter etc?

interp1 was very slow in older versions. Does this concern modern Matlab versions also?

[Adam Danz]

You could also use the profiler to look into which sections are relatively slower.

[Conrado Neto]

In my codes sum_load, Pseudo_Load_Factor and Disp_Factor are obtained/imported from .dat files.

I've updated the question with a running code, using rand for those arrays

Answer 1

[Bruno Luong]

M = 1000;
N = 500;
P = 1000;
sum_Load = rand(1,M);             % For demonstration purposes
Pseudo_Load_Factor = rand(N,P); % For demonstration purposes
Disp_Factor =  rand(N,P);       % For demonstration purposes

Disp_Inter = zeros(M,N);
for j = 1:M
    Load_Inter = sum_Load(j);
    A = Pseudo_Load_Factor-Load_Inter;
    B = A(:,1:end-1).*A(:,2:end) < 0;
    [Bm,idx] = max(B,[],2);
    idxl = (1:N)'+(idx-1)*N;
    idxr = idxl+N;
    PLFR = Pseudo_Load_Factor(idxr);
    w = (PLFR-Load_Inter)./(PLFR-Pseudo_Load_Factor(idxl));
    DPI = w.*Disp_Factor(idxl) + (1-w).*Disp_Factor(idxr);
    DPI(Bm==0) = 700;
    Disp_Inter(j,:) = DPI;
end

[Conrado Neto]

This soluton is a lot faster, but for some reason, I used isequal to compare the vector I obtain with my code versus using this code, and they are different

although, by looking at many numbers in the matrices, they are the same

is there any approximation in your solution (or mine?)?

thanks!

[Bruno Luong]

No, simply because the arithmethic order is different. The error is small and can be considered as finite precision numerical errors, they are normal.

They are similar to this

>> 1/3 -1/2 + 1/6

ans =

  -2.7756e-17

>> 1/3 + 1/6 -1/2

ans =

     0

[Conrado Neto]

I see

another thing I am not sure is included in your code, and that is very important on mine is ensure that the first point be returned from the interpolation, check the following figure

that is why I use the 'find' command with the option 'first'.

[Bruno Luong]

Yes it's the firstpoint (the idx returned by MAX command is the first one). Otherwise mine result won't match yours.

[Conrado Neto]

perfect, thanks

Answer 2

[Steven Lord]

I have some off the top of my head suggestions that come from the structure of the code, but if you explain in words what Disp_Inter represents and/or the underlying problem you're trying to solve we may be able to suggest functions to help improve your algorithm and/or its implementation at a higher level.

%snip
for j = 1:M
    Load_Inter = sum_Load(j);
    for i = 1:N
        idx = find(cummin(Pseudo_Load_Factor(i,:))<Load_Inter & cummax(Pseudo_Load_Factor(i,:))>Load_Inter,1,'first');

An obvious low-hanging fruit: since Pseudo_Load_Factor doesn't change inside the nested loops, compute the cumulative minimum and cumulative maximum of them once before entering the loop and index into those precomputed matrices inside the loop.

        if isempty(idx)
            Disp_Inter(j,i) = 0.10*7000;

Preallocate Disp_Inter to this constant value before the loop (though I'd use 700 rather than multiplying the integer value 7000 by the double precision approximation to one tenth.) That way you only need to change the value if idx is not empty.

%snip the rest of the code

[Conrado Neto]

Sure

Disp_inter is a vector with values of displacements

each Disp_inter has a corresponding Pseudo_Load_Factor value

(It is a load x displacement curve)

Then, separately, I have another vector, sum_Load, with values of load that I need to find the corresponding displacement in this load displacement curve

since sum_Load may not correspond exactly to Pseudo_load_Factor, I interpolate to find the corresponding displacement.

look at the following figure:

I hope I was able to explain clearly, if not, let me know

thanks