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); end Thanks in advance for any help Stefano |
On Wed, Apr 26, 2017 at 4:35 PM, smateria <[hidden email]> wrote:
With gamfit in the inner loop, you can vectorize what you want, but it will not help. I hope it is not just the second isfinite check that is killing you, and trying to fit on infinite data. The stored data variable is different from the pre-allocated data. Try to supply a complete working example that demonstrates the problem. pkg load statistics optim obsmean=rand(204, 350); pre_mod = rand(74, 204, 350); ... Finally try using the profiler (profile on/off/clear, profshow & profexplore) _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 |
Free forum by Nabble | Edit this page |