econ_order.w
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:13k
源码类别:

通讯编程

开发平台:

Visual C++

  1. % This file is part of the Stanford GraphBase (c) Stanford University 1993
  2. @i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
  3. @i gb_types.w
  4. deftitle{ECON_,ORDER}
  5. def<#1>{$langle${rm#1}$rangle$}
  6. prerequisite{GB_,ECON}
  7. @* Near-triangular ordering.
  8. This demonstration program takes a matrix of data
  9. constructed by the {sc GB_,ECON} module and permutes the economic sectors
  10. so that the first sectors of the ordering tend to be producers of
  11. primary materials for other industries, while the last sectors
  12. tend to be final-product
  13. industries that deliver their output mostly to end users.
  14. More precisely, suppose the rows of the matrix represent the outputs
  15. of a sector and the columns represent the inputs. This program attempts
  16. to find a permutation of rows and columns that minimizes the sum of
  17. the elements below the main diagonal. (If this sum were zero, the
  18. matrix would be upper triangular; each supplier of a sector would precede
  19. it in the ordering, while each customer of that sector would follow it.)
  20. The general problem of finding a minimizing permutation is NP-complete;
  21. it includes, as a very special case, the {sc FEEDBACK ARC SET} problem
  22. discussed in Karp's classic paper [{sl Complexity of Computer
  23. @^Karp, Richard Manning@>
  24. Computations} (Plenum Press, 1972), 85--103].
  25. But sophisticated ``branch and cut'' methods have been developed that work
  26. well in practice on problems of reasonable size.
  27. Here we use a simple heuristic downhill method
  28. to find a permutation that is locally optimum, in the sense that
  29. the below-diagonal sum does not decrease if any individual
  30. sector is moved to another position while preserving the relative order
  31. of the other sectors. We start with a random permutation and repeatedly
  32. improve it, choosing the improvement that gives the least positive
  33. gain at each step. A primary motive for the present implementation
  34. was to get further experience with this method of cautious descent, which
  35. was proposed by A. M. Gleason in {sl AMS Proceedings of Symposia in Applied
  36. @^Gleason, Andrew Mattei@>
  37. Mathematics/ bf10} (1958), 175--178. (See the comments following
  38. the program below.)
  39. @ As explained in {sc GB_,ECON}, the subroutine call |econ(n,2,0,s)|
  40. constructs a graph whose |n<=79| vertices represent sectors of the
  41. U.S. economy and whose arcs $uto v$ are assigned numbers corresponding to the
  42. flow of products from sector~|u| to sector~|v|. When |n<79|, the
  43. |n| sectors are obtained from a basic set of 79 sectors by
  44. combining related commodities. If |s=0|, the combination is done in
  45. a way that tends to equalize the row sums, while if |s>0|, the combination
  46. is done by choosing a random subtree of a given 79-leaf tree;
  47. the ``randomness'' is fully determined by the value of~|s|.
  48. This program uses two random number seeds, one for |econ| and one
  49. for choosing the random initial permutation. The former is called~|s|
  50. and the latter is called~|t|. A further parameter, |r|, governs the
  51. number of repetitions to be made; the machine will try |r|~different
  52.  starting permutations
  53. on the same matrix. When |r>1|, new solutions are displayed only when
  54. they improve on the previous best.
  55. By default, |n=79|, |r=1|, and |s=t=0|. The user can change these
  56. default parameters by specifying options
  57. on the command line, at least in a UNIX/ implementation, thereby
  58. obtaining a variety of special effects. The relevant
  59. command-line options are .{-n}<number>, .{-r}<number>,
  60. .{-s}<number>, and/or .{-t}<number>. Additional options
  61. .{-v} (verbose), .{-V} (extreme verbosity), and .{-g}
  62. (greedy or steepest descent instead of cautious descent) are also provided.
  63. @^UNIX dependencies@>
  64. Here is the overall layout of this CEE/ program:
  65. @p
  66. #include "gb_graph.h" /* the GraphBase data structures */
  67. #include "gb_flip.h" /* the random number generator */
  68. #include "gb_econ.h" /* the |econ| routine */
  69. @h@#
  70. @<Global variables@>@;
  71. main(argc,argv)
  72.   int argc; /* the number of command-line arguments */
  73.   char *argv[]; /* an array of strings containing those arguments */
  74. {@+unsigned long n=79; /* the desired number of sectors */
  75.   long s=0; /* random \{seed} for |econ| */
  76.   long t=0; /* random \{seed} for initial permutation */
  77.   unsigned long r=1; /* the number of repetitions */
  78.   long greedy=0; /* should we use steepest descent? */
  79.   register long j,k; /* all-purpose indices */
  80.   @<Scan the command-line options@>;
  81.   g=econ(n,2L,0L,s);
  82.   if (g==NULL) {
  83.     fprintf(stderr,"Sorry, can't create the matrix! (error code %ld)n",
  84.              panic_code);
  85.     return -1;
  86.   }
  87.   printf("Ordering the sectors of %s, using seed %ld:n",g->id,t);
  88.   printf(" (%s descent method)n",greedy?"Steepest":"Cautious");
  89.   @<Put the graph data into matrix form@>;
  90.   @<Print an obvious lower bound@>;
  91.   gb_init_rand(t);
  92.   while (r--)
  93.     @<Find a locally optimum permutation and report the below-diagonal sum@>;
  94.   return 0; /* normal exit */
  95. }
  96. @ Besides the matrix $M$ of input/output coefficients, we will find it
  97. convenient to use the matrix $Delta$, where $Delta_{jk}=M_{jk}-M_{kj}$.
  98. @d INF 0x7fffffff /* infinity (or darn near) */
  99. @<Global...@>=
  100. Graph *g; /* the graph we will work on */
  101. long mat[79][79]; /* the corresponding matrix */
  102. long del[79][79]; /* skew-symmetric differences */
  103. long best_score=INF; /* the smallest below-diagonal sum we've seen so far */
  104. @ @<Scan the command-line options@>=
  105. while (--argc) {
  106. @^UNIX dependencies@>
  107.   if (sscanf(argv[argc],"-n%lu",&n)==1) ;
  108.   else if (sscanf(argv[argc],"-r%lu",&r)==1) ;
  109.   else if (sscanf(argv[argc],"-s%ld",&s)==1) ;
  110.   else if (sscanf(argv[argc],"-t%ld",&t)==1) ;
  111.   else if (strcmp(argv[argc],"-v")==0) verbose=1;
  112.   else if (strcmp(argv[argc],"-V")==0) verbose=2;
  113.   else if (strcmp(argv[argc],"-g")==0) greedy=1;
  114.   else {
  115.     fprintf(stderr,"Usage: %s [-nN][-rN][-sN][-tN][-g][-v][-V]n",argv[0]);
  116.     return -2;
  117.   }
  118. }
  119. @ The optimum permutation is a function only of the $Delta$ matrix, because
  120. we can subtract any constant from both $M_{jk}$ and $M_{kj}$ without changing
  121. the basic problem.
  122. @<Put the graph data into matrix form@>=
  123. {@+register Vertex *v;
  124.   register Arc *a;
  125.   n=g->n;
  126.   for (v=g->vertices;v<g->vertices+n;v++)
  127.     for (a=v->arcs;a;a=a->next)
  128.       mat[v-g->vertices][a->tip-g->vertices]=a->flow;
  129.   for (j=0;j<n;j++)
  130.     for (k=0;k<n;k++)
  131.       del[j][k]=mat[j][k]-mat[k][j];
  132. }
  133. @ Nontrivial lower bounds that can be made strong enough to find provably
  134. optimum solutions to the ordering problem can be based on linear programming,
  135. as shown for example by Gr"otschel, J"unger, and Reinelt
  136. @^Gr"otschel, Martin@>
  137. @^J"unger, Michael@>
  138. @^Reinelt, Gerhard@>
  139. [{sl Operations Research bf32} (1984), 1195--1220].
  140. The basic idea is to formulate the problem as
  141. the task of minimizing $sum M_{jk}x_{jk}$ for integer variables $x_{jk}ge0$,
  142. subject to the conditions $x_{jk}+x_{kj}=1$ and $x_{ik}le x_{ij}+x_{jk}$
  143. for all triples $(i,j,k)$ of distinct subscripts; these conditions are
  144. necessary and sufficient. Relaxing the integrality constraints gives a
  145. lower bound, and we can also add additional inequalities such as
  146. $x_{14}+x_{25}+x_{36}+x_{42}+x_{43}+x_{51}+x_{53}+x_{61}+x_{62}le7$.
  147. The interesting story of inequalities like this has been surveyed by
  148. P. C. Fishburn [{sl Mathematical Social Sciences/ bf23} (1992), 67--80].
  149. @^Fishburn, Peter Clingerman@>
  150. However, our goal is more modest---we just want
  151. to study two of the simplest heuristics. So we will be happy with a trivial
  152. bound based only on the constraints $x_{jk}+x_{kj}=1$.
  153. @<Print an obvious lower bound@>=
  154. {@+register long sum=0;
  155.   for (j=1;j<n;j++)
  156.     for (k=0;k<j;k++)
  157.       if (mat[j][k]<=mat[k][j]) sum+=mat[j][k];
  158.       else sum+=mat[k][j];
  159.   printf("(The amount of feed-forward must be at least %ld.)n",sum);
  160. }
  161. @* Descent.
  162. At each stage in our search, |mapping| will be the current permutation;
  163. in other words, the sector in row and column~|k| will be
  164. |g->vertices+mapping[k]|. The current below-diagonal sum will be
  165. the value of |score|. We will not actually have to permute anything
  166. inside of |mat|.
  167. @d sec_name(k) (g->vertices+mapping[k])->name
  168. @<Glob...@>=
  169. long mapping[79]; /* current permutation */
  170. long score; /* current sum of elements above main diagonal */
  171. long steps; /* the number of iterations so far */
  172. @ @<Find a locally optimum perm...@>=
  173. {
  174.   @<Initialize |mapping| to a random permutation@>;
  175.   while(1) {
  176.     @<Figure out the next move to make; |break| if at local optimum@>;
  177.     if (verbose) printf("%8ld after step %ldn",score,steps);
  178.     else if (steps%1000==0 && steps>0) {
  179.       putchar('.');
  180.       fflush(stdout); /* progress report */
  181.     }
  182.     @<Take the next step@>;
  183.   }
  184.   printf("n%s is %ld, found after %ld step%s.n",@|
  185.          best_score==INF?"Local minimum feed-forward":
  186.                "Another local minimum",@|
  187.          score,steps,steps==1?"":"s");
  188.   if (verbose || score<best_score) {
  189.     printf("The corresponding economic order is:n");
  190.     for (k=0;k<n;k++) printf(" %sn",sec_name(k));
  191.   if (score<best_score) best_score=score;
  192.   }
  193. }
  194. @ @<Initialize |mapping| to a random permutation@>=
  195. steps=score=0;
  196. for (k=0; k<n; k++) {
  197.   j=gb_unif_rand(k+1);
  198.   mapping[k]=mapping[j];
  199.   mapping[j]=k;
  200. }
  201. for (j=1; j<n; j++) for (k=0;k<j;k++) score+=mat[mapping[j]][mapping[k]];
  202. if (verbose>1) {
  203.   printf("nInitial permutation:n");
  204.   for (k=0;k<n;k++) printf(" %sn",sec_name(k));
  205. }
  206. @ If we move, say, |mapping[5]| to |mapping[3]| and shift the previous
  207. entries |mapping[3]| and |mapping[4]| right one, the score decreases by
  208. $$hbox{|del[mapping[5]][mapping[3]]+del[mapping[5]][mapping[4]]|},.$$
  209. Similarly, if we move |mapping[5]| to |mapping[7]| and shift the previous
  210. entries |mapping[6]| and |mapping[7]| left one, the score decreases by
  211. $$hbox{|del[mapping[6]][mapping[5]]+del[mapping[7]][mapping[5]]|},.$$
  212. The number of possible moves is $(n-1)^2$. Our job is to find the
  213. one that makes the score decrease, but by as little as possible (or, if
  214. |greedy!=0|, to make the score decrease as much as possible).
  215. @<Figure out the next move to make; |break| if at local optimum@>=
  216. best_d=greedy? 0: INF;
  217. best_k=-1;
  218. for (k=0;k<n;k++) {@+register long d=0;
  219.   for (j=k-1;j>=0;j--) {
  220.     d+=del[mapping[k]][mapping[j]];
  221.     @<Record the move from |k| to |j|, if |d| is better than |best_d|@>;
  222.   }
  223.   d=0;
  224.   for (j=k+1;j<n;j++) {
  225.     d+=del[mapping[j]][mapping[k]];
  226.     @<Record the move...@>;
  227.     }
  228.   }
  229. if (best_k<0) break;
  230. @ @<Record the move...@>=
  231. if (d>0 && (greedy? d>best_d: d<best_d)) {
  232.   best_k=k;
  233.   best_j=j;
  234.   best_d=d;
  235. }
  236. @ @<Glob...@>=
  237. long best_d; /* best improvement seen so far on this step */
  238. long best_k,best_j; /* moving |best_k| to |best_j| improves by |best_d| */
  239. @ @<Take the next step@>=
  240. if (verbose>1)
  241.   printf("Now move %s to the %s, pastn",sec_name(best_k),
  242.           best_j<best_k? "left": "right");
  243. j=best_k;
  244. k=mapping[j];
  245. do@+{
  246.   if (best_j<best_k) mapping[j]=mapping[j-1],j--;
  247.   else mapping[j]=mapping[j+1],j++;
  248.   if (verbose>1) printf("    %s (%ld)n",sec_name(j),@|
  249.            best_j<best_k?del[mapping[j+1]][k]:
  250.                          del[k][mapping[j-1]]);
  251. }@+while(j!=best_j);
  252. mapping[j]=k;
  253. score-=best_d;
  254. steps++;
  255. @ How well does cautious descent work? In this application, it
  256. is definitely too cautious. For example, after lots of computation with the
  257. default settings, it comes up
  258. with a pretty good value (457342), but only after taking 39,418 steps!
  259. Then (if |r>1|) it tries again and stops with 461584 after 47,634 steps.
  260. The greedy algorithm with the same starting permutations obtains the
  261. local minimum 457408 after only 93 steps, then 460411 after 83 steps.
  262. The greedy algorithm tends to find solutions that are a bit inferior,
  263. but it is so much faster that it allows us to run many
  264. more experiments. After 20 trials with the default settings, it finds
  265. a permutation with only 456315 below the diagonal,
  266. and after about 250 more it reduces this upper bound to 456295.
  267. (Gerhard Reinelt has proved, via branch-and-cut,
  268. that 456295 is in fact optimum.)
  269. @^Reinelt, Gerhard@>
  270. The method of stratified greed, which is illustrated in the {sc
  271. FOOTBALL} module, should do better than the ordinary greedy algorithm;
  272. and interesting results can be expected when stratified greed is compared
  273. also to other methods like simulated annealing and genetic breeding.
  274. Comparisons should be made by seeing which method can come up with the
  275. best upper bound after calculating for a given number of mems (see
  276. {sc MILES_,SPAN}). The upper bound obtained in any run is a random
  277. variable, so several independent trials of each method should be~made.
  278. Question: Suppose we divide the vertices into two subsets and prescribe
  279. a fixed permutation on each subset. Is it NP-complete to find the
  280. optimum way to merge these two permutations---i.e., to find a
  281. permutation, extending the given ones, that has the smallest
  282. below-diagonal sum?
  283. @* Index. We close with a list that shows where the identifiers of this
  284. program are defined and used.