Re: super slow loop

[smateria]

Hi, thank you very much for your reply. 
I will try with the profiler as you suggested, but my feeling is that you
are right when you say I have no chance with the gamfit in the inner loop.
Though, why does gamfit work fine in the first loop?

Somehow it seems you do not see the edits I had made on my first post,
anyhow this is the complete code, working until the last loop becomes
incredibly slow.
 

clear all

mkoctfile ../statistics-1.3.0/__gammainc_lentz.cc
more off
 
lat = ncread ('precip_multimodel_ens01_SouEU.1989-2005.nc','lat');
lon = ncread ('precip_multimodel_ens01_SouEU.1989-2005.nc','lon');
time = ncread ('precip_multimodel_ens01_SouEU.1989-2005.nc','time');
 
% read the monthly time series of each ensemble precipitation

A=dir("precip_mul*"); 
prate_multimodel=NaN(74,25,14,1224);
for ens=1:74
     prate_multimodel(ens,:,:,:) = ncread(A(ens).name,'prate');
end
 
prate_multimodel=permute(prate_multimodel,[1 4 3 2]);
pre_mod=  reshape(pre_mod(:,2:6:1224,:,:),74,204,14*25);

% read the monthly time series of precipitation

prets = ncread ('cru_SouEUprec.1989-2005.nc','pre'); 
prets = permute(prets,[3 2 1]);
obsmean = reshape(prets,204,14*25); obsmean(obsmean > 1e10) = NaN;

spi_obs = NaN(192,350); 
spi_mod_res = NaN(192,25900);   
pre_mod_res = reshape(pre_mod,204,25900);

for i = 1:350 
    if isfinite(obsmean(:,i))==1; 
        spi_obs(:,i) = SPI(obsmean(:,i),3,12);     
    end 
end 

for ens = 1:25900 
    spi_mod_res(:,ens) = SPI(pre_mod_res(:,ens),3,12); 
end 

The first loop takes around 50 seconds to complete, or about 30 if I
suppress a few printouts from the function nmsmax.m. 
The second loop slows down at the 2nd iteration, and then starts proceeding
superslowly (if I ^C after five minutes, I see we are at the third iteration
:O). 

Here below, the SPI Matlab function by T.Lee (2009), perfectly working on
Matlab 

/function [Z]=SPI(Data,scale,nseas) 

erase_yr = ceil(scale/12); 

% Data setting to scaled dataset 

A1 = []; 
for is = 1:scale, A1=[A1,Data(is:length(Data)-scale+is)]; end 
XS = sum(A1,2); 

if (scale>1), XS(1:nseas*erase_yr-scale+1)=[]; end 

for is = 1:nseas 
    tind = is:nseas:length(XS); 
    Xn = XS(tind); 
    [zeroa] = find(Xn==0); 
    Xn_nozero = Xn; Xn_nozero(zeroa)=[]; 
    q = length(zeroa)/length(Xn); 
    parm = gamfit(Xn_nozero); 
    Gam_xs = q+(1-q)*gamcdf(Xn,parm(1),parm(2)); 
    Z(tind) = norminv(Gam_xs); 
end/


Thanks again everybody,
Stefano




--
View this message in context: 
http://octave.1599824.n4.nabble.com/super-slow-loop-using-SPI-function-tp4683040p4683055.html
Sent from the Octave - General mailing list archive at Nabble.com.