Thursday, July 2, 2009

GAMS: Fill parameters with random sparse data

For testing I need sometimes very large, but very sparse data to play with. To generate random data, the easiest is something like:

a(i,j,k)$(uniform(0,1)<fracsparse) = uniform(-100,100);

However this will loop over every (i,j,k). If very sparse this can be become relatively expensive. It is faster in that case to loop over the expected number of nonzero elements, and place them randomly in the parameter. For this we use the lag operator for indices.

The two methods are not exactly the same: the second method can never generate more than nn nonzeroes, while the first method can. However for testing algorithms in GAMS the generated data is often “random enough” to give an indication about performance etc.

$ontext

  
fill a(i,j,k) and b(i,j,k) with random sparse data.
  
Both sparsity pattern as the data should be random.

  
We demonstrate two methods. The second method is efficient if
  
the parameters are *very* sparse.

$offtext

$set n 500
set i /i1*i%n%/;
scalar n /%n%/;
alias (i,j,k);


scalar fracsparse 'fraction of nonzero elements' / 0.0001 /;


*
* method 1
* this is slow because we loop over (i,j,k)
*

parameter a(i,j,k);
a(i,j,k)$(uniform(0,1)<fracsparse) = uniform(-100,100);


*
* method 2
* loop over number of nonzero elements only
*

scalar nn; nn = fracsparse*n*n*n;
display nn;

parameter b(i,j,k);
scalar kk;
loop(i$sameas(i,'i1'),
 
for(kk=1 to nn,
     b(i+uniformint(0,n-1),i+uniformint(0,n-1),i+uniformint(0,n-1)) = uniform(-100,100);
  );
);


parameter statistics(*,*);
statistics(
'a','sum')   = sum((i,j,k),a(i,j,k));
statistics(
'a','count') = sum((i,j,k)$a(i,j,k),1);
statistics(
'b','sum')   = sum((i,j,k),b(i,j,k));
statistics(
'b','count') = sum((i,j,k)$b(i,j,k),1);
display statistics;


The timings are as follows:

  seconds
assignment a 18.237
assignment b (including loops) 0.031

Note1: the outer loop for b is not doing anything. It is looping over one element. The reason is we cannot use a lead or lag with a single element b(‘i1’+…).

Note2: This is not completely original: I think I have seen similar code in some model from Tom Rutherford. Sorry: I don’t have a reference for that.