# super slow loop using SPI function

3 messages
Open this post in threaded view
|

## super slow loop using SPI function

 This post was updated on . Hi everyone, I have been trying to use a function to calculate standardized precipitation index (SPI) that works fine in Matlab, but seems having huge problems in Octave. In fact, after two integrations, the loop starts proceeding incredibly slowly. I have vectorized my matrices as much as I could but the problem persists. Here is the code: [...] size (obsmean) ans =    204   350 size (pre_mod) ans =    74   204   350 spi_obs = NaN(192,350); spi_mod_res = NaN(192,25900);   pre_mod_res = reshape(pre_mod,204,25900); % reshaped from a 74ens x 204time x 350latlon 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) % Standardized Precipitation Index % Input Data % Data : Monthly Data vector not matrix (monthly or seasonal precipitation) % scale : 1,3,12,48 % nseas : number of season (monthly=12) % Example % Z=SPI(gamrnd(1,1,1000,1),3,12); 3-monthly scale, % Notice that the rest of the months of the first year are removed. % eg. if scale =3, fist year data 3-12 SPI values are not estimated. % if row vector then make coloumn vector % if (sz==1); Data(:,1) = Data;end 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); endThanks in advance for any help Stefano