How to enhance the performance of for-loops and cell-arrays (related to statistical calculations)?

The following code calculates some performance measures out of different periodes. No error messages occured, but the processing time is very long.

F = 'runoff.txt'; % name of the file
D = 'C:\Users\heute\model\results\model_standalone\'; % absolute or relative path of base directory
S = dir(fullfile(D,'results*'));
X = [S.isdir] & ~ismember({S.name},{'.','..'});
N = {S(X).name};
L = cell(size(N));
C = cell(size(N));
for k = 1:numel(N)
    T = fullfile(D,N{k},F);
    fid = fopen(T,'rt');
    fmt = ['%s',repmat('%f',1,6)];
    opt = {'HeaderLines',1,'CollectOutput',true};
    Z = textscan(fid,fmt,opt{:});
    fclose(fid);
    L{k} = Z{1}; % timestamp
    C{k} = Z{2}; % data
    %        
    Qs = C{k}(:,6); % define the simulated runoff, as column 6 in each cell array
    %
    % define the periodes for computing performance measures
    sdatelim_neu = [datenum(2013,10,01,00,00,00) datenum(2016,10,01,00,00,00)];
    dt = 1/24;
    date = sdatelim_neu(1):dt:sdatelim_neu(2);
    date_runoff = transpose(date);
    %
    sdatelim1 = [datenum(2014,05,01,00,00,00) datenum(2014,10,01,00,00,00)];
    dt = 1/24;
    sdate_sdatelim1 = sdatelim1(1):dt:sdatelim1(2);
    %
    sdatelim2 = [datenum(2015,05,01,00,00,00) datenum(2015,10,01,00,00,00)];
    sdate_sdatelim2 = sdatelim2(1):dt:sdatelim2(2);
    %
    sdatelim3 = [datenum(2016,05,01,00,00,00) datenum(2016,10,01,00,00,00)];
    sdate_sdatelim3 = sdatelim3(1):dt:sdatelim3(2);
    %
    % loop over the different periodes
    for s = 1:length(sdate_sdatelim1);
        for a = 1:length(sdate_sdatelim2);
            for b = 1:length(sdate_sdatelim3);
                j = find(date_runoff >= sdate_sdatelim1(s) & date_runoff < sdate_sdatelim1(k)+dt) & find(date_runoff >= sdate_sdatelim2(a) & date_runoff < sdate_sdatelim2(a)+dt) & find(date_runoff >= sdate_sdatelim3(b) & date_runoff < sdate_sdatelim3(b)+dt);
                f_1k = 1-cov(Qs(j) - Qo)/var(Qo); %NSE
                f_2k = sqrt(mean((Qs(j) - Qo).^2)); %RMSE
                f_3k = abs(mean(Qs(j)- Qo)); %BIAS
                %Qo is the observed runoff -> imported from file
                %
                % write into matrix YA -> for use in further analysis
                YA = [f_1k, f_2k, f_3k];
            end
        end
    end
end

As a test case, I ran this code for two inputfiles (each of them has 26280 rows in column 6). In the end however several 1000 input-files should be processed.

How can I reduce the computing time?

or is there an error within the for-loop over the different periods? or is this:

Qs = C{k}(:,6); % define the simulated runoff, as column 6 in each cell array

an inefficient command?

(I use Matlab R2012a)

[dpb]

How about telling us what it is you're trying to compute instead, so we can think about how to best solve the problem rather than try to figure out how to speed up existing code?

A small sample of the data sufficient to illustrate the computations never hurts...

[Glazio]

@dpb:

I have some 1000 txt files that contain simulated runoff values. This files have all the same name and each file is saved within a subfolder called results and a consecutive number (that means results1, results2, results3,....).

The simulated runoff values are read from each txt-file and written to a cell-array (called C{k} = Z{2}).

Now, the task is to compute some performance measures (NSE, RMSE, BIAS) from the comparison of observed and simulated values. For each of the runoff.txt files one NSE, RMSE and BIAS should be calculated at the end.

Within a single file (runoff.txt), only certain periods should be considered to calculate the performance measures.

These periods are:

• 1. May 2014 to 30. September 2014
• 1. May 2015 to 30. September 2015
• 1. May 2016 to 30. September 2016

In the end the Matrix YA should contain 3 columns (NSE, RMSE, BIAS) and the number of rows should correspond to the number of processed runoff.txt-files (e.g. for 1000 files -> 1000 lines).

How can this be solved efficiently?

There are attached two example files:

1.) runoff.txt -> simulation results for one parameter combination 2.) runoff_observed.txt -> measurements of runoff

[Walter Roberson]

What is the code trying to calculate? What is the goal?

[Adam]

Best way to speed up numeric calculation code that uses cell arrays is to not use cell arrays, in general! Numeric arrays are vastly faster than cell arrays so if you can rework your data structures to use them instead it should help for starters.

[Stephen23]

• Those nested loops could likely be vectorized.
• It does not seem to be necessary store the data from all files, as you only seem to process the data from the current file.
• https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

[Glazio]

@Walter Roberson: The code is trying to calculate the Root Mean Square Error, BIAS and NSE for each inputfile (runoff.txt) and should consider only certain periods for calculation.

The goal is a matrix YA which contains the performance measure combination for each runoff-file.

[Glazio]

@Stephen Cobeldick, thanks for your help.

What exactly do you mean with:

"It does not seem to be necessary store the data from all files, as you only seem to process the data from the current file." ?

In the end, all results should be in YA.

Answer 1

[dpb]

S = 'runoff.txt';
O = 'runoff_observed.txt';
D = 'C:\Users\heute\model\results\model_standalone\';
d = dir(fullfile(D,'results*'));                          % list of directories
fmtS = ['%{dd.MM.yyyy-HH:mm}D' repmat('%*f',1,5) '%f'];   % simulated format string
fmtO = ['%{dd.MM.yyyy HH:mm}D' %f'];                      % observed format string
L=length(d);                                              % number sudirs
YA=zeros(L,3);                                            % preallocate
for k = 1:L                                               % iterate over subdirs
  fid = fopen(fullfile(D,d{k}.name,S),'rt');              % open simulated
  Z=textscan(fid,fmtS,'headerlines',2,'collectoutput',1); % read simulated
  fclose(fid);
  dtS=Z{:,1};                                             % timestamp simulated (datetime)
  Qs=Z{:,2};                                              % simulated data
  fid = fopen(fullfile(D,d.name{k},O),'rt');              % open observed
  Z=textscan(fid,fmtO,'headerlines',1,'collectoutput',1); % read observed
  fclose(fid);
  dtO=Z{:,1};                                             % timestamp observed (datetime)
  Qo=Z{:,2};                                              % observed data
  % define the periods for computing performance measures
  yr1=2014; yr2=2016;                                     % years to compute over
output
  ix=isbetween(dtO,datenum(yr(1),05,01),datenum(yr(1),10,01)); % first year
  for yr=yr1+1:yr2  % subsequent years
    ix=ix | isbetween(dtO,datenum(yr,05,01),datenum(yr,10,01));
  end
  YA(k,:) = [f_1k, f_2k, f_3k];
end

ADDENDUM

Cleaned up to incorporate changes from conversation below excepting for opening the reference file--treat that as need to. Above should then return the L records in the output array.

ERRATUM

NB: Remove the (1) index from year reference yr in the loop to get the subsequent years after first...inadvertently left it in there when copied line.

[Glazio]

@dpb: Thanks for your help!

The following code:

fid = fopen(fullfile(D,d.name{k},S),'rt');

returns this error message:

"Field reference for multiple structure elements that is followed by more reference blocks is an error."

What could be the reason?

[Stephen23]

It is a mistake in dpb's code. Replace this:

d.name{k}

with the correct way to access a non-scalar structure:

d(k).name

[Jan]

@Galzio: Are you able to fix this typo by your own? did you try it?

"d.name{k}" ==> "d(k).name"

[dpb]

I did note it was "air code'" :)

I seem to make a lot of typos and transpositions any more; that's one I don't recall having done before, though. Invention reigns...

[Glazio]

@Stephen Cobeldick, Jan Simon, dpb: Thanks for your help, yes I've tried to fix this and other typos.

Which alternative is suitable for the function isbetween, because this function isn't defined in R2012a?

ix=isbetween(dtO,date1,date2)

[dpb]

It's just "syntactic sugar" for the comparison...my utility iswithin would be the same syntax/result..

function flg=iswithin(x,lo,hi)
% returns T for values within range of input
% SYNTAX:
%  [log] = iswithin(x,lo,hi)
%      returns T for x between lo and hi values, inclusive

flg= (x>=lo) & (x<=hi);
>>

[Glazio]

@dpb: Thanks, now the code runs without any error messages.

For me, however, it seems that the program does not provide exactly what is desired. In my test case I use two inputfiles with simulated results and a file with the reference data.

The resulting matrix should now have two rows and three columns:

• number of rows = number of used input files with simulated data (here 2 files)
• number of columns = number of used performance measures (NSE, RMSE, BIAS).

But in this case, YA has 3 rows and 3 columns (so one line to much).

Could it be that the program outputs the performance measures separately for each individual period?

How could the problem be fixed?

[Glazio]

Faith to have found a solution to my last comment.

Have the following code part:

yr1=2014; yr2=2016;                                     % years to compute over
  YA=zeros(yr2-yr1+1,3);                                  % preallocate for output
  idx=0;                                                  % output index
  for yr=yr1:yr2
    date1 = datenum(yr,05,01);                            % dates to consider
    date2 = datenum(yr,10,01);
    ix=isbetween(dtO,date1,date2);                        % inclusive dates logical vector
    f_1k = 1-cov(Qs(ix)-Qo(ix))/var(Qo(ix));              %NSE
    f_2k = rms(Qs(ix)-Qo(ix));                            %RMSE
    f_3k = abs(mean(Qs(ix)-Qo(ix));                       %BIAS
    % write into matrix YA -> for use in further analysis
    idx=idx+1;
    YA(idx,:) = [f_1k, f_2k, f_3k];
  end

slightly modified to:

Y1 = zeros(k,3);
    idx=0;                                                  % output index
    for yr=yr1:yr2
      date1 = datenum(yr,05,01,00,00,00);                   % dates to consider
      date2 = datenum(yr,10,01,00,00,00);
      %ix=iswithin(dto,date1,date2);                         % inclusive dates logical vector
      ix = iswithin(dto,date1,date2);
      f_1k = 1-cov(Qs(ix)-Qo(ix))/var(Qo(ix));              %NSE
      f_2k = rms(Qs(ix)-Qo(ix));                            %RMSE
      f_3k = abs(mean(Qs(ix)-Qo(ix)));                      %BIAS
      % write into matrix YA -> for use in further analysis
      idx=idx+1;
      Y1(idx,:) = [f_1k, f_2k, f_3k];
      Y2 = sum(Y1);

      end
      YA(k,:) = Y2;

The program now works much faster :-)

A general question:

Is it more meaningful to preallocate the output YA to a fixed known file size or to change it over k?

[dpb]

Well, it wasn't easy to decipher what was really intended vis a vis what was coded...I took a guess. :)

The above assumed there was a reference file for each year as well; if that's not so then move the code to read reference data out of the loop and use whatever logic it is that generates the proper file name for it.

And yes, the logic test is for each year in the loop over years; if it is to compare all in "one swell foop" then you'll have to build the composite logical vector to select all at once--that would just be or ing those together

...
ix=isbetween(dtO,datenum(yr(1),05,01),datenum(yr(1),10,01)); % first year
for yr=yr1+1:yr2  % subsequent years
  ix=ix | isbetween(dtO,datenum(yr,05,01),datenum(yr,10,01));
end
f_1k = 1-cov(Qs(ix)-Qo(ix))/var(Qo(ix));              %NSE
f_2k = rms(Qs(ix)-Qo(ix));                            %RMSE
f_3k = abs(mean(Qs(ix)-Qo(ix));                       %BIAS
% write into matrix YA -> for use in further analysis
idx=idx+1;
YA(idx,:) = [f_1k, f_2k, f_3k];
...

[dpb]

...
Y1(idx,:) = [f_1k, f_2k, f_3k];
Y2 = sum(Y1);

It may be faster, but it's not the same answer by any stretch.

sum([rms(x(yr1)) rms(x(yr2)) rms(x(yr3))])

is NOT same as

rms([x(yr1) x(yr2) x(yr3)])

[Glazio]

@dpb: Yes, I just saw the mistake and of course it is not the same ;-)

If I use the following code:

...
ix=isbetween(dtO,datenum(yr(1),05,01),datenum(yr(1),10,01)); % first year
for yr=yr1+1:yr2  % subsequent years
  ix=ix | isbetween(dtO,datenum(yr(1),05,01),datenum(yr(1),10,01));
end
f_1k = 1-cov(Qs(ix)-Qo(ix))/var(Qo(ix));              %NSE
f_2k = rms(Qs(ix)-Qo(ix));                            %RMSE
f_3k = abs(mean(Qs(ix)-Qo(ix));                       %BIAS
% write into matrix YA -> for use in further analysis
idx=idx+1;
YA(idx,:) = [f_1k, f_2k, f_3k];
...

I get only one row for all two files, but it would have to be two lines.

What could be the problem?

[dpb]

Because don't need idx any longer; it was just to increment for each year as thought was the intent. Use k instead and can remove idx entirely.

[Glazio]

@dpb: Thanks, now the first row has only 0's and relevant results are only visible in the last row.

What could be the reason?

[Glazio]

@dpb: Thanks for your help and patience, now it seems as if the right results are delivered :-)

[dpb]

Well, not if you didn't make the fixup I noted above--it'll run but will process the first year three times as was.

MORAL: ALWAYS debug thoroughly; running w/o error doesn't guarantee correctness!