## Monday, August 11, 2008

### 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.462big 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.