Sorting in GAMS

Although there is an external tool that can help, sorting in GAMS remains somewhat problematic. In this case I needed to sort each column in a matrix p(i,j). The length of each column is 10-15 while the number of columns can be hundreds. To confirm I had a reasonable implementation I tried for implementations:

• Writing sorting directly with loops
• Use smin with other loops
• Use card(j) times a call to gdxrank
• Make a very large one dimensional vector and sort that

   1  set
   2    i /i1*i15/
   3    j /j1*j1000/
   4  *  i /i1*i5/
   5  *  j /j1*j2/
   6  ;
   7  
   8  
   9  alias(i,ii);
  10  
  11  parameter p(i,j);
  12  p(i,j) = uniform(0,1);
  13  
  14  parameter timing(*);
  15  scalar t;
  16  
  17  parameter order1(i,j),p1(i,j);
  18  p1(i,j) = p(i,j);
  19  *display p;
  20  
  21  t = jnow;
  22  
  23  scalar pmin;
  24  loop(j,
  25  
  26     loop(ii,
  27        pmin = 9999999.99;
  28        loop(i$(p1(i,j)<pmin),
  29           order1(ii,j) = ord(i);
  30           pmin = p1(i,j);
  31        );
  32        p1(i,j)$(order1(ii,j)=ord(i)) =  9999999.99;
  33     );
  34  
  35  );
  36  
  37  timing('loop algorithm') = (jnow-t)*24*60*60;
  38  
  39  t = jnow;
  40  
  41  parameter order2(i,j),p2(i,j);
  42  p2(i,j) = p(i,j);
  43  *display p;
  44  
  45  scalar pmin,first;
  46  loop(j,
  47  
  48     loop(ii,
  49        pmin = smin(i, p2(i,j));
  50        first = 1;
  51        loop(i$(p2(i,j)=pmin and first),
  52           order2(ii,j) = ord(i);
  53           first = 0;
  54        );
  55        p2(i,j)$(order2(ii,j)=ord(i)) =  9999.9999;
  56     );
  57  
  58  );
  59  
  60  timing('loop/smin algorithm') = (jnow-t)*24*60*60;
  61  
  62  t = jnow;
  63  
  64  parameter order3(i,j),psort(i),pindex(i);
  65  
  66  loop(j,
  67     psort(i) = p(i,j);
  68     execute_unload "rank_in.gdx", psort;
  69     execute "=gdxrank rank_in.gdx rank_out.gdx";
  70     execute_load  "rank_out.gdx", pindex=psort;
  71     loop((i,ii)$(pindex(i)=ord(ii)),
  72         order3(ii,j) = ord(i);
  73     );
  74  );
  75  
  76  timing('gdxrank loop algorithm') = (jnow-t)*24*60*60;
  77  
  78  t = jnow;
  79  
  80  set k /1*15000/;
  81  parameter psort4(k),pindex4(k),order4(i,j);
  82  scalar pmax,pmin;
  83  pmax = smax((i,j),p(i,j));
  84  pmin = smin((i,j),p(i,j));
  85  set pmap(i,j,k);
  86  pmap(i,j,k)$(ord(k)=card(i)*(ord(j)-1)+ord(i)) = yes;
  87  *display pmap;
  88  
  89  psort4(k) = sum(pmap(i,j,k), [p(i,j)-pmin]/(pmax-pmin)+ord(j));
  90  *display p,psort5;
  91  execute_unload "rank_in.gdx", psort4;
  92  execute "=gdxrank rank_in.gdx rank_out.gdx";
  93  execute_load  "rank_out.gdx", pindex4=psort4;
  94  
  95  loop(j,
  96     pindex(i) = sum(pmap(i,j,k),pindex4(k)) - (ord(j)-1)*card(i);
  97  *   display pindex;
  98     loop((i,ii)$(pindex(i)=ord(ii)),
  99         order4(ii,j) = ord(i);
 100     );
 101  );
 102  
 103  timing('big gdxrank algorithm') = (jnow-t)*24*60*60;
 104  
 105  
 106  *display p,order1,order2,order3,order4;
 107  
 108  display timing;

The timing results are:

----    108 PARAMETER timing 
 
loop algorithm          0.641,    loop/smin algorithm     0.867,    gdxrank loop algorithm 19.462
big gdxrank algorithm  54.200
 

Note that for the last example we spend all time in the calculation of the pmap set.

Conclusion: GAMS would benefit from a built-in sorting facility.