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.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.