90 year data simulation in need of performance increase

Question

I developed some code for my masters project that will simulate 90 years daily data using 1000 different data sets. The code is working fine and gives the correct output that I wanted but the processing time are very high. It took about 8 hours to finish the simulation.

    tic
%% importing the csv file with selected column
files=dir('*_scen_*.csv');
for i=1:length(files);
    LHR=importcsv(files(i).name);

%% Definable variables
% Define These Value
TAW=-216;     %total available water
RAW=-129;     %readily available water
KC=1.0;       %crop coefficient
IRL=15;       %intense rain level
RC=(80/100);  %percentage of recharge 
RO=(1-RC);    %percentage of runoff

% The very first row of Soil Moisture Deficit
for j=1
    SMD(j,i)=(LHR.RAIN(j)-LHR.PET(j));
 if     SMD(j,i)>0;
        SMD(j,i)=0;
 elseif SMD(j,i)<RAW;
        SMD(j,i)=(LHR.RAIN(j)-(LHR.PET(j)*((TAW-SMD(j-1))/(TAW-RAW))));
 end
end

%for the following SMD Calculation
for k=2:(length(LHR.RAIN));
    SMD(k,i)=SMD(k-1,i)+(LHR.RAIN(k)-LHR.PET(k));
    % The SMD conditions
    if SMD(k,i)>0;
       SMD(k,i)=0;
 elseif SMD(k,i)<RAW;
     SMD(k,i)=SMD(k-1,i)+(LHR.RAIN(k)-(LHR.PET(k)*((TAW-SMD(k-1,i))/(TAW-RAW))));
    end
end
%Convert negative SMD to Positive
SMD=abs(SMD);

 %%Evapotranspiration Calculation
 for l=1:(length(SMD));
    if SMD(l,i)<abs(RAW);
        AET(l,i)=LHR.PET(l);
    elseif SMD(l,i)>abs(RAW);
        AET(l,i)=KC*LHR.PET(l)*((abs(TAW)-(SMD(l,i)))/(abs(TAW)-abs(RAW)));
    end
end
for m=2:(length(SMD));
    if SMD(m,i)<abs(RAW);
        AET(m,i)=LHR.PET(m);
    elseif SMD(m,i)>abs(RAW);
       AET(m,i)=KC*LHR.PET(m)*((abs(TAW)-(SMD(m-1,i)))/(abs(TAW)-abs(RAW)));
    end
end
%% HER calculation
for n=1:length(SMD);
    if SMD(n,i)<(LHR.RAIN(n)-AET(n,i));
        HER(n,i)=(LHR.RAIN(n)-AET(n,i)-SMD(n,i));
    elseif SMD(n,i)>(LHR.RAIN(n)-AET(n,i));
        HER(n,i)=0;
    end
end
%% Calculation of recharge anf runoff
for o=1:(length(HER));
    if (HER(o,i)+(abs(TAW)-SMD(o,i)))<abs(TAW);
        RUNOFF(o,i)=0;
    elseif (HER(o,i)+(abs(TAW)-SMD(o,i)))>abs(TAW);
        if HER(o,i)>IRL;
            RUNOFF(o,i)=RO*HER(o,i);
        elseif HER(o,i)<IRL;
            RUNOFF(o,i)=0;    
        end
    end
    if (HER(o,i)+(abs(TAW)-SMD(o,i)))<abs(TAW);
        RECHARGE(o,i)=0;
    elseif (HER(o,i)+(abs(TAW)-SMD(o,i)))>abs(TAW);
        if HER(o,i)>IRL;
            RECHARGE(o,i)=RC*HER(o,i);
        elseif HER(o,i)<IRL;
            RECHARGE(o,i)=HER(o,i);
        end
    end
end
%% rainfall
for p=1:length(LHR.RAIN);
RAINFALL(p,i)=LHR.RAIN(p);
PET(p,i)=LHR.PET(p);
end

end
clear i
clear j
clear k
clear l
clear m
clear n
clear o
clear p

toc

Is there any improvement scope for this code that might reduce the processing time? Sorry if the code looks unprofessional; I am in the beginner stage for MATLAB programming.

Answer 1

This cleanup code at the end is indicative of the problem with this program:

clear i
clear j
clear k
clear l
clear m
clear n
clear o
clear p

Nobody is going to reverse-engineer your minified / obfuscated code to understand it. You will probably have a hard time understanding it yourself if you come back to it after a few weeks. The program is therefore unmaintainable.

I wouldn't have much confidence in the correctness of the results, either. The code should be broken down into functions, each of which has a single purpose, specific inputs and outputs, and is separately testable.

Answer 2

You can define functions in MATLAB by creating files named like functionname.m with content like this:

% the 'g' is the return value, 'sigmoid' is the function name and
% 'z' is the only parameter
function g = sigmoid(z)

g = zeros(size(z));
g = 1./ (1 + exp(-z));

end

One other general idea that can speed up MATLAB quite a lot is vectorization: When you have two arrays and you multiply their contents and add them up, that can also be done with a vector multiplication (and similar matrix multiplication). But although the result is the same, the execution time with the vector variant is MUCH better.

(I agree with 200_success - you should name your variables with something that makes sense so that we can actually understand what is going on)

Another problem is indentation: At the first glance I did not see the outermost for loop.

Answer 3

First off:

1. Indentation! If you have written the code and figure out afterwards that your indentation is wrong / missing / messy, Matlab makes it very easy for you! Ctrl+A followed by Ctrl+I. Simple as that, auto-indentation.

2. Variable names. Choose descriptive, unambigous variable names. Remember that other people might need to read and understand your code!

3. Don't use i and j as variable names in Matlab.

4. Spaces! Use spaces. Expressions such as x=y*z>3+1 are hard to read.

5. tic/toc. timeit is the preferred choice when benchmarking Matlab processes. If you're using an old version of Matlab, you can get the function from the File Exchange.

6. numel is better than length. It's more robust, but also much faster

7. All the constants can be declared before the loop.

---

You want it faster?

Pre-allocate memory for SMD! SMD = zeros(tot_rows, numel_files).

Vectorization is the way to go!

Let's have a look at the first loop:

for j=1
    SMD(j,i)=(LHR.RAIN(j)-LHR.PET(j));
 if     SMD(j,i)>0;
        SMD(j,i)=0;
 elseif SMD(j,i)<RAW;
        SMD(j,i)=(LHR.RAIN(j)-(LHR.PET(j)*((TAW-SMD(j-1))/(TAW-RAW))));
 end
end

Wait, what?! This isn't even a loop, for j = 1. Get rid of the for-part. And, it doesn't work if the elseif part get's evaluated: SMD(j-1) for j=1 will cause an error as Matlab is 1-indexed. Therefore, I'm assuming that will never happen. Rewrite it:

SMD(1,:) = LHR.RAIN(1)-LHR.PET(1) 
SMD(1, SMD(1,:) > 0) = 0;

In fact, let's get rid of the entire outer loop, you're never using i to anything useful anyway. We'll keep the inner loops for now:

num_rain = numel(LHR.RAIN);
for k=2:num_rain
    SMD(k,:)=SMD(k-1,i)+(LHR.RAIN(k)-LHR.PET(k));
    % The SMD conditions

    SMD(k, (SMD(k,:) > 0)) = 0;
    SMD(k, (SMD(k,:) < RAW) = SMD(k-1,:) + (LHR.RAIN(k)-(LHR.PET(k)).*(TAW-SMD(k-1,:))./(TAW-RAW)
end

Voila, you don't need the outer loop anymore! A lot of time saved! Also, all the if's and elses are gone with the wind!

Can we get rid of the inner loops too? All rows are dependent on the previous rows, and unfortunately that's hard to vectorize. It could be possible to use cumsum and subtract the appropriate terms, but I think that'll be a hazzle.

Anyway, you'll have a much cleaner and much faster function if you follow these tips!

Good luck! =)