super slow loop using SPI function

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

super slow loop using SPI function

smateria
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
Reply | Threaded
Open this post in threaded view
|

Re: super slow loop

Kire Pudsje


On Wed, Apr 26, 2017 at 4:35 PM, smateria <[hidden email]> wrote:

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:

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
Reply | Threaded
Open this post in threaded view
|

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