Monday, April 13, 2009

Speeding Up Gams

The following calculations took more than 4 hours:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =
    sum(phzone$CountryLGSxTRCrops(faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr), 1);

scalar max;
max = smax((faoqcountry,faocrop,grid,irrg_key,lgsxtr), count(faoqcountry,faocrop,grid,irrg_key, lgsxtr));
display max;

The parameter CountryLGSxTRCrops is very large (but sparse): 2.5 million elements. The reason for this fragment to be so slow is two-fold:

  • In the assignment to parameter count we are doing things out of order. GAMS prefers to loop over sets at the end of the index list, but here we loop over an index in the middle (phzone). Sometimes GAMS is able to reorder things automatically to do things faster, but in this case apparently not.
  • smax is executed dense, as zeroes can be significant. In GAMS zero and “does not exist” is the same. This can be exploited when sparse processing is used. Here GAMS is conservative and reverts to dense processing as the implicit zeros may be important in calculating the correct SMAX value.

This a direct rewrite into an explicit loop that executes in 30 seconds:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
scalar mx /0/;
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr)=0;
loop((faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr)$CountryLGSxTRCrops(faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr),
    count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =  count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) + 1;
    mx = max(mx,count(faoqcountry,faocrop,grid,irrg_key,lgsxtr));
);

Explicit loops are expensive in GAMS but still much cheaper compared to dense processing. Even faster is:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
parameter ReorderedCrops(faoqcountry,faocrop,grid,irrg_key,lgsxtr,phzone);
option ReorderedCrops < CountryLGSxTRCrops;

count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =
     sum(phzone$reorderedCrops(faoqcountry,faocrop,grid,irrg_key,lgsxtr,phzone), 1);

scalar max;
max = smax((faoqcountry,faocrop,grid,irrg_key,lgsxtr)$count(faoqcountry,faocrop,grid,irrg_key,lgsxtr),
             count(faoqcountry,faocrop,grid,irrg_key,lgsxtr));
display max;

Here we use the (undocumented)  < option to reorder the big parameter such that the phzone index is last. This is the representation that GAMS likes most when looping over phzone. Now the calculation of count is very fast. Next we force the calculation of max to be performed over the nonzero elements in parameter count only. This formulation takes less than 5 seconds. We use some more memory by duplicating CountryLGSxTRCrops.