diff options
author | Norbert Preining <norbert@preining.info> | 2019-09-02 13:46:59 +0900 |
---|---|---|
committer | Norbert Preining <norbert@preining.info> | 2019-09-02 13:46:59 +0900 |
commit | e0c6872cf40896c7be36b11dcc744620f10adf1d (patch) | |
tree | 60335e10d2f4354b0674ec22d7b53f0f8abee672 /support/graphbase/assign_mona.w |
Initial commit
Diffstat (limited to 'support/graphbase/assign_mona.w')
-rw-r--r-- | support/graphbase/assign_mona.w | 692 |
1 files changed, 692 insertions, 0 deletions
diff --git a/support/graphbase/assign_mona.w b/support/graphbase/assign_mona.w new file mode 100644 index 0000000000..3a1ab9db8a --- /dev/null +++ b/support/graphbase/assign_mona.w @@ -0,0 +1,692 @@ +% This file is part of the Stanford GraphBase (c) Stanford University 1992 +\def\title{ASSIGN\_\thinspace MONA} +@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES! +\def\<#1>{$\langle${\rm#1}$\rangle$} +\def\dash{\mathrel-\joinrel\joinrel\mathrel-} % adjacent vertices +\def\ddash{=\joinrel\joinrel=} % matched vertices + +\prerequisite{GB\_\thinspace MONA} +@* The assignment problem. +This demonstration program takes a matrix +constructed by the |gb_mona| module and chooses at most one number from +each row and column in such a way as to maximize the sum of the numbers +chosen. It also reports the number of ``mems'' (memory references) +expended during its computations, so that the algorithm it uses +can be compared with alternative procedures. + +The matrix has $m$ rows and $n$ columns. If $m\le n$, one number will +be chosen in each row; if $m\ge n$, one number will be chosen in each column. +The numbers in the matrix are brightness levels (i.e., pixel values) in +a digitized version of the Mona Lisa. + +Of course the author does not pretend that the location of ``highlights'' in +da Vinci's painting, one per row and one per column, has any application +to art appreciation. However, this program does seem to have pedagogic value, +because the relation between pixel values and shades of gray allows us +to visualize the data underlying this special case of the +assignment problem; ordinary matrices of numeric data are much harder +to perceive. The non-random nature of pixels +in a work of art may also have similarities to the ``organic'' properties +of data in real-world applications. + +This program is optionally able to produce an encapsulated PostScript file +from which the solution can be displayed graphically, with halftone shading. + +@ As explained in |gb_mona|, the subroutine call |mona(m,n,d,m0,m1,n0,n1,d0,d1, +area)| constructs an $m\times n$ matrix of integers between $0$ and~$d$, +inclusive, based on the brightness levels in a rectangular region of +a digitized Mona Lisa, where |m0|, |m1|, |n0|, and |n1| define that +region. The raw data is obtained as a sum of |(m1-m0)(n1-n0)| pixel +values between $0$ and~$255$, then scaled in such a way that sums |<=d0| +are mapped to zero, sums |>=d1| are mapped to~$d$, and intermediate sums are +mapped linearly to intermediate values. Default values |m1=360|, |n1=250|, +|m=m1-m0|, |n=n1-n0|, |d=255|, and |d1=255(m1-m0)(n1-n0)| are substituted if +any of the parameters |m|, |n|, |d|, |m1|, |n1|, or |d1| are zero. + +The user can specify the nine parameters |(m,n,d,m0,m1,n0,n1,d0,d1)| +on the command line, at least in a \UNIX\ implementation, thereby +obtaining a variety of special effects; the relevant +command-line options are \.{m=}\<number>, \.{m0=}\<number>, and so on, +with no spaces before or after the \.= signs that separate parameter +names from parameter values. Additional options are also provided: +\.{-s} (use only Mona's $16\times32$ ``smile''); +\.{-c} (complement black/white); \.{-p} (print the matrix and solution); +\.{-P} (produce a PostScript file \.{mona.eps} for graphic output); +\.{-h} (use a heuristic that applies only when $m=n$); and +\.{-v} or \.{-V} (print verbose or Very verbose commentary about the + algorithm's performance). +@^UNIX dependencies@> + +Here is the overall layout of this \Cee\ program: + +@p +#include "gb_graph.h" /* the GraphBase data structures */ +#include "gb_mona.h" /* the |mona| routine */ +@# +@<Global variables@>@; +main(argc,argv) + int argc; /* the number of command-line arguments */ + char *argv[]; /* an array of strings containing those arguments */ +{@+@<Local variables@>; + @<Scan the command line options@>; + mtx=mona(m,n,d,m0,m1,n0,n1,d0,d1,working_storage); + if (mtx==NULL) { + fprintf(stderr,"Sorry, can't create the matrix! (error code %d)\n", + panic_code); + return -1; + } + printf("Assignment problem for %s%s\n",mona_id,(compl?", complemented":"")); + sscanf(mona_id,"mona(%u,%u,%lu",&m,&n,&d); /* adjust for defaults */ + if (m!=n) heur=0; + if (printing) @<Display the input matrix@>; + if (PostScript) @<Output the input matrix in PostScript format@>; + mems=0; + @<Solve the assignment problem@>; + if (printing) @<Display the solution@>; + if (PostScript) @<Output the solution in PostScript format@>; + printf("Solved in %d mems%s.\n",mems, + (heur?" with square-matrix heuristic":"")); +} + +@ @f Vertex int /* |gb_graph| defines these data types */ +@f Arc int +@f Graph int +@f Area int + +@<Glob...@>= +Area working_storage; /* where to put the input data and auxiliary arrays */ +long *mtx; /* input data for the assignment problem */ +long mems; /* the number of memory references counted + while solving the problem */ + +@ The following local variables are related to the command-line options: + +@<Local v...@>= +unsigned m=0,n=0; /* number of rows and columns desired */ +unsigned long d=0; /* number of pixel values desired, minus~1 */ +unsigned m0=0,m1=0; /* input will be from rows $[|m0|\,.\,.\,|m1|)$ */ +unsigned n0=0,n1=0; /* and from columns $[|n0|\,.\,.\,|n1|)$ */ +unsigned long d0=0,d1=0; /* lower and upper threshold of raw pixel scores */ +int compl=0; /* should the input values be complemented? */ +int heur=0; /* should the square-matrix heuristic be used? */ +int printing=0; /* should the input matrix and solution be printed? */ +int PostScript=0; /* should an encapsulated PostScript file be produced? */ + +@ @<Scan the command line options@>= +while (--argc) { +@^UNIX dependencies@> + if (sscanf(argv[argc],"m=%u",&m)==1) ; + else if (sscanf(argv[argc],"n=%u",&n)==1) ; + else if (sscanf(argv[argc],"d=%lu",&d)==1) ; + else if (sscanf(argv[argc],"m0=%u",&m0)==1) ; + else if (sscanf(argv[argc],"m1=%u",&m1)==1) ; + else if (sscanf(argv[argc],"n0=%u",&n0)==1) ; + else if (sscanf(argv[argc],"n1=%u",&n1)==1) ; + else if (sscanf(argv[argc],"d0=%u",&d0)==1) ; + else if (sscanf(argv[argc],"d1=%u",&d1)==1) ; + else if (strcmp(argv[argc],"-s")==0) smile; /* sets |m0|, |m1|, |n0|, |n1| */ + else if (strcmp(argv[argc],"-c")==0) compl=1; + else if (strcmp(argv[argc],"-h")==0) heur=1; + else if (strcmp(argv[argc],"-v")==0) verbose=1; + else if (strcmp(argv[argc],"-V")==0) verbose=2; /* terrifically verbose */ + else if (strcmp(argv[argc],"-p")==0) printing=1; + else if (strcmp(argv[argc],"-P")==0) PostScript=1; + else { + fprintf(stderr, + "Usage: %s [param=value] [-s] [-c] [-h] [-v] [-p] [-P]\n",argv[0]); + return -2; + } +} + +@ @<Display the input matrix@>= +for (k=0;k<m;k++) { + for (l=0;l<n;l++) printf("% 4d",compl?d-*(mtx+k*n+l):*(mtx+k*n+l)); + printf("\n"); +} + +@ We obtain a crude but useful estimate of the computation time +by counting mem units, as explained in the |miles_span| program. + +@d o mems++ +@d oo mems+=2 +@d ooo mems+=3 + +@* Algorithmic overview. The {\it assignment problem\/} is the classical +problem of weighted bipartite matching, the problem of choosing +a maximum-weight set of disjoint edges in a bipartite graph. We will consider +only the case of complete bipartite graphs, when the weights are +specified by an $m\times n$ matrix. + +An algorithm is most easily developed if we begin with the assumption +that the matrix is square (i.e., that $m=n$), and if we change from +maximization to minimization. Then the assignment problem is the task +of finding a permutation $\pi[0]\ldots\pi[n-1]$ of $\{0,\ldots,n-1\}$ +such that $\sum_{k=0}^{n-1} a_{k\pi[k]}$ is minimized, where +$A=(a_{kl})$ is a given matrix of numbers $a_{kl}$ for $0\le k,l<n$. +The algorithm below works for arbitrary real numbers $a_{kl}$, but we +will assume in our implementation that the matrix entries are integers. + +One way to approach the assignment problem is to make three simple +observations: (a)~Adding a constant to any row of the matrix does not +change the solution $\pi[0]\ldots\pi[n-1]$. (b)~Adding a constant to +any column of the matrix does not change the solution. (c)~If $a_{kl}\ge0$ +for all $k$ and~$l$, and if $\pi[0]\ldots\pi[n-1]$ is a permutation +with the property that $a_{k\pi[k]}=0$ for all~$k$, then $\pi[0]\ldots\pi[n-1]$ +solves the assignment problem. + +The remarkable fact is that these three observations actually suffice. In +other words, there is always a sequence of constants $(\sigma_0,\ldots,\sigma_ +{n-1})$ and $(\tau_0,\ldots,\tau_{n-1})$ and a permutation $\pi[0]\ldots +\pi[n-1]$ such that +$$\vbox{\halign{$#$,\hfil&\quad for #\hfil\cr +a_{kl}-\sigma_k+\tau_{\,l}\ge0& $0\le k<n$ and $0\le l<n$;\cr +a_{k\pi[k]}-\sigma_k+\tau_{\pi[k]}=0& $0\le k<n$.\cr}}$$ + +@ To prove the remarkable fact just stated, we start by reviewing the +theory of {\it unweighted\/} bipartite matching. Any $m\times n$ matrix +$A=(a_{kl})$ befines a bipartite graph on the vertices $(r_0,\ldots,r_{m-1})$ +and $(c_0,\ldots,c_{n-1})$ if we say that $r_k\dash c_l$ whenever +$a_{kl}=0$; in other words, the edges of the bipartite graph are the zeroes +of the matrix. Two zeroes of~$A$ are called {\it independent\/} if they appear +in different rows and columns; this means that the corresponding edges have +no vertices in common. A set of mutually independent zeroes of the matrix +therefore corresponds to a set of mutually disjoint edges, also called a +{\it matching\/} between rows and columns. + +The Hungarian mathematicians Egerv\'ary and K\"onig proved +[{\sl Matematikai \'es Fizikai Lapok\/ \bf38} (1931), 16--28, 116--119] +that the maximum number of independent zeroes in a matrix is equal to +the minimum number of rows and/or columns that are needed to ``cover'' +every zero. In other words, if we can find $p$ independent zeroes but +not~$p+1$, then there is a way to choose $p$ lines in such a way that +every zero of the matrix is included in at least one of the chosen lines, +where a ``line'' is either a row or a column. + +Their proof was constructive, and it leads to a useful computer algorithm. +Given a set of $p$ independent zeroes of a matrix, let us write +$r_k\ddash c_l$ or $c_l\ddash r_k$ and say that $r_k$ is matched with $c_l$ +if $a_{kl}$ is one of these $p$ special +zeroes, while we continue to write $r_k\dash c_l$ or $c_l\dash r_k$ +if $a_{kl}$ is one of the nonspecial zeroes. A given set of $p$ +special zeroes defines a choice of $p$ lines in the following way: Column~$c$ +is chosen if and only if it is reachable by a path of the form +$$r_0\dash c_1\ddash r_1\dash c_2\ddash\cdots\dash c_q\ddash r_q\,,\eqno(*)$$ +where $r_0$ is unmatched, $q\ge1$, and $c=c_q$. Row~$r$ is chosen if +and only if it is matched with a column that is not chosen. Thus exactly +$p$ lines are chosen. We can now prove that the chosen lines cover +all the zeroes, unless there is a way to find $p+1$ independent zeroes. + +For if $c\ddash r$, either $c$ or $r$ has been chosen. And +if $c\dash r$, one of the following cases must arise. (1)~If $r$ and~$c$ +are both unmatched, we can increase~$p$ by matching them to each other. +(2)~If $r$ is unmatched and $c\ddash r'$, then $c$ has been chosen, so +the zero has been covered. (3)~If $r$ is matched to $c'\ne c$, then +either $r$ has been chosen or $c'$ has been chosen. In the latter case +there is a path of the form +$$r_0\dash c_1\ddash r_1\dash c_2\ddash\cdots\ddash + r_{q-1}\dash c'\ddash r\dash c\,,$$ +where $r_0$ is unmatched and $q\ge1$. +If $c$ is matched, it has therefore been chosen; otherwise we can increase $p$ +by redefining the matching to include +$$r_0\ddash c_1\dash r_1\ddash c_2\dash\cdots\dash + r_{q-1}\ddash c'\dash r\ddash c\,.$$ + +@ Now suppose $A$ is a {\it nonnegative\/} matrix. +Cover the zeroes of~$A$ with a minimum number of lines, $p$, using the +algorithm of Egerv\'ary and K\"onig. If $p<n$, some elements are still +uncovered, so those elements are positive; suppose the minimum uncovered +value is $\delta>0$. We can subtract $\delta$ from each unchosen row +and add $\delta$ to each chosen column; the net effect is to subtract~$\delta$ +from all uncovered elements and to add~$\delta$ to all doubly-covered +elements, while leaving all singly-covered elements unchanged. This +transformation causes a new zero to appear, while preserving +$p$ independent zeroes of the previous matrix (since they were each +covered only once). If we repeat the Egerv\'ary-K\"onig construction +with the same $p$ independent zeroes, we find that either $p$~is no +longer maximum or at least one more column has been chosen. +(The new zero $r\dash c$ occurs in a row~$r$ that was either unmatched +or matched to a previously chosen column, because row~$r$ was not +chosen.) Therefore if we repeat the process, we must eventually +be able to increase $p$ until finally $p=n$. This will solve the +assignment problem, proving the remarkable claim made earlier. + +@ If the given matrix $A$ has $m$ rows and $n>m$ columns, +we can extend it artificially +until it is square, by setting $a_{kl}=0$ for all $m\le k<n$ and +$0\le l<n$. The construction above will then apply. But we need not +waste time making such an extension, because it suffices to run the +algorithm on the original $m\times n$ matrix until $m$ independent zeroes +have been found. The reason is that the set of matched vertices always +grows monotonically in the Egerv\'ary-K\"onig construction: If a +column is matched at some stage, it will remain matched from that time on, +although it may well change partners. The $n-m$ dummy rows at the bottom +of~$A$ are always chosen to be part of the covering; so the dummy entries +become nonzero only in the columns that are part of some covering. +Such columns are part of some matching, so they are part of the +final matching. Therefore at most $m$ columns of the dummy entries +become nonzero during the procedure. We can always find $n-m$ independent +zeroes in the $n-m$ dummy rows of the matrix, so we need not deal with the +dummy elements explicitly. + +@ It has been convenient to describe the algorithm by saying that +we add and subtract constants to and from the colums and rows of~$A$. +But all those additions and subtractions can take a lot of time. So we will +merely pretend to make the adjustments that the method calls for; we will +represent them implicitly by two vectors $(\sigma_0,\ldots,\sigma_{m-1})$ +and $(\tau_0,\ldots,\tau_{n-1})$. Then the current value of each matrix +entry will be $a_{kl}-\sigma_k+\tau_{\,l}$, instead of $a_{kl}$. The +``zeroes'' will be positions such that $a_{kl}=\sigma_k-\tau_{\,l}$. + +Initially we will set $\tau_{\,l}=0$ for $0\le l<n$ and $\sigma_k= +\min\{a_{k0},\ldots,a_{k(n-1)}\}$ for $0\le k<m$. If $m=n$ we can also +make sure that there's a zero in every column by subtracting +$\min\{a_{0l},\ldots,a_{(n-1)l}\}$ from $a_{kl}$ for all $k$ and~$l$. +(This initial adjustment can conveniently be made to the original +matrix entries, instead of indirectly via the $\tau$'s.) Users can +discover if such a transformation is worthwhile by trying the program +both with and without the \.{-h} option. + +We have been saying a lot of things and proving a bunch of theorems, +without writing any code. Let's get back into programming mode +by writing the routine that is called into +action when the \.{-h} option has been specified: + +@d aa(k,l) *(mtx+k*n+l) /* a macro to access the matrix elements */ + +@<Subtract column minima in order to start with lots of zeroes@>= +{ + for (l=0; l<n; l++) { + o,s=aa(0,l); /* the |o| macro counts one mem */ + for (k=1;k<n;k++) + if (o,aa(k,l)<s) s=aa(k,l); + if (s>0) + for (k=0;k<n;k++) + oo,aa(k,l)-=s; /* |oo| counts two mems */ + } + if (verbose) printf(" The heuristic has cost %d mems.\n",mems); +} + +@ @<Local var...@>= +register int k; /* the current row of interest */ +register int l; /* the current column of interest */ +register int j; /* another interesting column */ +register long s; /* the current matrix element of interest */ + +@* Algorithmic details. +The algorithm sketched above is quite simple, except that we did not +discuss how to determine the chosen columns~$c_q$ that +are reachable by paths of the stated form $(*)$. It is easy to find +all such columns by constructing an unordered forest whose nodes are rows, +beginning with all unmatched rows~$r_0$ and adding a row~$r$ +for which $c\ddash r$ when $c$ is adjacent to a row already in the forest. + +Our data structure, which is based on suggestions of Papadimitriou and +Steiglitz [{\sl Combinatorial Optimization\/} (Prentice-Hall, 1982), +$\mathchar"278$11.1], will use several arrays. If row~$r$ is matched +with column~$c$ we will have |matching_col[r]=c| and |matching_row[c]=r|; +if row~$r$ is unmatched, |matching_col[r]| will be |-1|, and +if column~$c$ is unmatched, |matching_row[c]| will be |-1|. +If column~$c$ has a mate and is also reachable in a path of the form $(*)$, +we will have $|parent_row|[c]=r'$ for some $r'$ in the forest. Otherwise +column~$c$ is not chosen, and we will have |parent_row[c]=-1|. The rows +in the current forest will be called |unchosen_row[0]| through +|unchosen_row[t-1]|, where |t| is the current total number of nodes. + +The amount $\sigma_k$ subtracted from row $k$ is called |row_dec[k]|; the +amount $\tau_{\,l}$ added to row~$l$ is called |col_inc[l]|. In order to +compute the minimum uncovered element efficiently, we maintain a +quantity called |slack[l]| representing the minimum uncovered element +in each column. More precisely, if column~$l$ is not chosen, +|slack[l]| is the minimum of $a_{kl} +-\sigma_k+\tau_{\,l}$ for $k\in\{|unchosen_row|[0],\ldots, +|unchosen_row|[q-1]\}$, where $q\le t$ is the number of rows in the +forest that we have explored so far. We also remember |slack_row[l]|, +the number of a row where the stated minimum occurs. + +Column $l$ is chosen if and only if |parent_row[l]>=0|. We will arrange +things so that we also have |slack[l]=0| in every chosen column. + +@<Local var...@>= +int* matching_col; /* the column matching a given row, or $-1$ */ +int* matching_row; /* the row matching a given column, or $-1$ */ +int* parent_row; /* ancestor of a given column's mate, or $-1$ */ +int* unchosen_row; /* node in the forest */ +int t; /* total number of nodes in the forest */ +int q; /* total number of explored nodes in the forest */ +long* row_dec; /* $\sigma_k$, the amount subtracted from a given row */ +long* col_inc; /* $\tau_{\,l}$, the amount added to a given column */ +long* slack; /* minimum uncovered entry seen in a given column */ +int* slack_row; /* where the |slack| in a given column can be found */ +int unmatched; /* this many rows have yet to be matched */ + +@ @<Allocate the intermediate data structures@>= +matching_col=gb_alloc_type(m,@[int@],working_storage); +matching_row=gb_alloc_type(n,@[int@],working_storage); +parent_row=gb_alloc_type(n,@[int@],working_storage); +unchosen_row=gb_alloc_type(m,@[int@],working_storage); +row_dec=gb_alloc_type(m,@[long@],working_storage); +col_inc=gb_alloc_type(n,@[long@],working_storage); +slack=gb_alloc_type(n,@[long@],working_storage); +slack_row=gb_alloc_type(n,@[int@],working_storage); +if (gb_alloc_trouble) { + fprintf(stderr,"Sorry, out of memory!\n"); return -3; +} + +@ The algorithm operates in stages, where each stage terminates +when we are able to increase the number of matched elements. + +The first stage is different from the others; it simply goes through +the matrix and looks for zeroes, matching as many rows and columns +as it can. This stage also initializes table entries that will be +useful in later stages. + +@d INF 0x7fffffff /* infinity (or darn near) */ + +@<Do the initial stage@>= +t=0; /* the forest starts out empty */ +for (l=0; l<n; l++) { + o,matching_row[l]=-1; + o,parent_row[l]=-1; + o,col_inc[l]=0; + o,slack[l]=INF; +} +for (k=0; k<m; k++) { + o,s=aa(k,0); /* get ready to calculate the minimum entry of row $k$ */ + for (l=1; l<n; l++) if (o,aa(k,l)<s) s=aa(k,l); + o,row_dec[k]=s; + for (l=0; l<n; l++) + if ((o,s==aa(k,l)) && (o,matching_row[l]<0)) { + o,matching_col[k]=l; + o,matching_row[l]=k; + if (verbose>1) printf(" matching col %d==row %d\n",l,k); + goto row_done; + } + o,matching_col[k]=-1; + if (verbose>1) printf(" node %d: unmatched row %d\n",t,k); + o,unchosen_row[t++]=k; +row_done:; +} + +@ If a subsequent stage has not succeeded in matching every row, +we prepare for a new stage by reinitializing the forest as follows. + +@<Get ready for another stage@>= +t=0; +for (l=0; l<n; l++) { + o,parent_row[l]=-1; + o,slack[l]=INF; +} +for (k=0; k<m; k++) + if (o,matching_col[k]<0) { + if (verbose>1) printf(" node %d: unmatched row %d\n",t,k); + o,unchosen_row[t++]=k; + } + +@ Here, then, is the algorithm's overall control structure. +There are at most $m$ stages, and each stage does $O(mn)$ operations, +so the total running time is $O(m^2n)$. + +@<Do the Hungarian algorithm@>= +@<Do the initial stage@>; +if (t==0) goto done; +unmatched=t; +while(1) { + if (verbose) printf(" After %d mems I've matched %d rows.\n",mems,m-t); + q=0; + while(1) { + while (q<t) { + @<Explore node |q| of the forest; + if the matching can be increased, |goto breakthru|@>; + q++; + } + @<Introduce a new zero into the matrix by modifying |row_dec| and |col_inc|; + if the matching can be increased, |goto breakthru|@>; + } +breakthru: @<Update the matching by pairing row $k$ with column $l$@>; + if(--unmatched==0) goto done; + @<Get ready for another stage@>; +} +done: @<Doublecheck the solution@>; + +@ @<Explore node |q| of the forest; + if the matching can be increased, |goto breakthru|@>= +{ + o,k=unchosen_row[q]; + o,s=row_dec[k]; + for (l=0; l<n; l++) + if (o,slack[l]) {@+register long del; + oo,del=aa(k,l)-s+col_inc[l]; + if (del<slack[l]) { + if (del==0) { /* we found a new zero */ + if (o,matching_row[l]<0) goto breakthru; + o,slack[l]=0; /* this column will now be chosen */ + o,parent_row[l]=k; + if (verbose>1) printf(" node %d: row %d==col %d--row %d\n", + t,matching_row[l],l,k); + oo,unchosen_row[t++]=matching_row[l]; + } else { + o,slack[l]=del; + o,slack_row[l]=k; + } + } + } +} + +@ At this point, column $l$ is unmatched, and row $k$ is in +the forest. By following parent links in the forest, +we can rematch rows and columns so that a previously unmatched row~$r_0$ +gets a mate. + +@<Update the matching by pairing row $k$ with column $l$@>= +if (verbose) printf(" Breakthrough at node %d of %d!\n",q,t); +while (1) { + o,j=matching_col[k]; + o,matching_col[k]=l; + o,matching_row[l]=k; + if (verbose>1) printf(" rematching col %d==row %d\n",l,k); + if (j<0) break; + o,k=parent_row[j]; + l=j; +} + +@ If we get to this point, we have explored the entire forest; none of +the unchosen rows has led to a breakthrough. An unchosen column with +smallest |slack| will allow us to make further progress. + +@<Introduce a new zero into the matrix by modifying |row_dec| and |col_inc|; + if the matching can be increased, |goto breakthru|@>= +s=INF; +for (l=0; l<n; l++) + if (o,slack[l] && slack[l]<s) + s=slack[l]; +for (q=0; q<t; q++) + ooo,row_dec[unchosen_row[q]]+=s; +for (l=0; l<n; l++) + if (o,slack[l]) { /* column $l$ is not chosen */ + o,slack[l]-=s; + if (slack[l]==0) @<Look at a new zero, and |goto breakthru| with + |col_inc| up to date if there's a breakthrough@>; + } else oo,col_inc[l]+=s; + +@ There may be several columns tied for smallest slack. If any of them +leads to a breakthough, we are very happy; but we must finish the loop on~|l| +before going to |breakthru|, because the |col_inc| variables +need to be maintained for the next stage. + +Within column |l|, there may be several rows that produce the same slack; +we have remembered only one of them, |slack_row[l]|. Fortunately, one is +sufficient for our purposes. We either have a breakthrough, or we choose +column~|l|, regardless of which row or rows led us to consider that column. + +@<Look at a new zero, and |goto breakthru| with + |col_inc| up to date if there's a breakthrough@>= +{ + o,k=slack_row[l]; + if (verbose>1) + printf(" Decreasing uncovered elements by %d produces zero at [%d,%d]\n", + s,k,l); + if (o,matching_row[l]<0) { + for (j=l+1; j<n; j++) + if (o,slack[j]==0) oo,col_inc[j]+=s; + goto breakthru; + } else { /* not a breakthrough, but the forest continues to grow */ + o,parent_row[l]=k; + if (verbose>1) printf(" node %d: row %d==col %d--row %d\n", + t,matching_row[l],l,k); + oo,unchosen_row[t++]=matching_row[l]; + } +} + +@ The code in the present section is redundant, unless cosmic +radiation has cause the hardware to malfunction. But there is some +reassurance whenever we find that mathematics still appears to be +consistent, so the author could not resist writing these few unnecessary lines, +which verify that the assignment problem has indeed been solved optimally. +(We don't count the mems.) + +@<Doublecheck...@>= +for (k=0;k<m;k++) + for (l=0;l<n;l++) + if (aa(k,l)<row_dec[k]-col_inc[l]) { + fprintf(stderr,"Oops, I made a mistake!\n"); + return -6; /* can't happen */ + } +for (k=0;k<m;k++) { + l=matching_col[k]; + if (l<0 || aa(k,l)!=row_dec[k]-col_inc[l]) { + fprintf(stderr,"Oops, I blew it!\n"); return-66; /* can't happen */ + } +} +k=0; +for (l=0;l<n;l++) if (col_inc[l]) k++; +if (k>m) { + fprintf(stderr,"Oops, I adjusted too many columns!\n"); + return-666; /* can't happen */ +} + +@* Interfacing. +A few nitty-gritty details still need to be handled: Our algorithm +is not symmetric between rows and columns, and it works only for $m\le n$; +so we will transpose the matrix when +$m>n$. Furthermore, our algorithm minimizes, but we actually want +it to maximize (except when |compl| is nonzero). + +Hence, we want to make the following transformations to the data before +processing it with the algorithm developed above. + +@<Solve the assignment problem@>= +if (m>n) @<Transpose the matrix@>@; +else transposed=0; +@<Allocate the intermediate data structures@>; +if (compl==0) + for (k=0; k<m; k++) for (l=0; l<n; l++) + aa(k,l)=d-aa(k,l); +if (heur) @<Subtract column minima...@>; +@<Do the Hungarian algorithm@>; + +@ @<Transpose...@>= +{ + if (verbose>1) printf("Temporarily transposing rows and columns...\n"); + tmtx=gb_alloc_type(m*n,@[long@],working_storage); + if (tmtx==NULL) { + fprintf(stderr,"Sorry, out of memory!\n"); return -4; + } + for (k=0; k<m; k++) for (l=0; l<n; l++) + *(tmtx+l*m+k)=*(mtx+k*n+l); + m=n;@+n=k; /* |k| holds the former value of |m| */ + mtx=tmtx; + transposed=1; +} + +@ @<Local v...@>= +long* tmtx; /* the transpose of |mtx| */ +int transposed; /* has the data been transposed? */ + +@ @<Display the solution@>= +{ + printf("The following entries produce an optimum assignment:\n"); + for (k=0; k<m; k++) + printf(" [%d,%d]\n",@| + transposed? matching_col[k]:k,@| + transposed? k:matching_col[k]); +} + +@* Encapsulated PostScript. +A special output file called \.{mona.eps} is written if the user has +selected the \.{-P} option. This file will contain a sequence of +PostScript commands that can be used to generate an illustration +within many kinds of documents. For example, if \TeX\ is being used +with the \.{dvips} output driver from Radical Eye Software and the +@.dvips@> +associated \.{epsf.tex} macros, one can say +$$\.{\\epsfxsize=10cm \\epsfbox\{mona.eps\}}$$ +within a \TeX\ document and the illustration will be typeset in +a box that is 10 centimeters wide. + +The conventions of PostScript allow the illustration to be scaled to +any size. Best results are probably obtained if each pixel is at +least one millimeter wide (about 1/25 inch) when printed. + +The illustration is formed by first +``painting'' the input data as a rectangle of pixels, +with up to 256 shades of gray. Then the solution pixels are +framed in black, with a white trim just inside the black edges +to help make the frame visible in already-dark places. The frames are +created by painting over the original image; the +center of each solution pixel retains its original color. + +Encapsulated PostScript files have a simple format that is recognized +by many software packages and printing devices. We use a subset of +PostScript that should be easy to convert to other languages if necessary. + +@<Output the input matrix in PostScript format@>= +{ + eps_file=fopen("mona.eps","w"); + if (!eps_file) { + fprintf("Sorry, I can't open the file `mona.eps'!\n"); + PostScript=0; + } else { + fprintf(eps_file,"%%!PS-Adobe-3.0 EPSF-3.0\n"); /* 1.0 and 2.0 also OK */ + fprintf(eps_file,"%%%%BoundingBox: -1 -1 %d %d\n",n+1,m+1); + fprintf(eps_file,"/buffer %d string def\n",n); + fprintf(eps_file,"%d %d 8 [%d 0 0 -%d 0 %d]\n",n,m,n,m,m); + fprintf(eps_file,"{currentfile buffer readhexstring pop} bind\n"); + fprintf(eps_file,"gsave %d %d scale image\n",n,m); + for (k=0;k<m;k++) @<Output row |k| as a hexadecimal string@>; + fprintf(eps_file,"grestore\n"); + } +} + +@ @<Glob...@>= +FILE *eps_file; /* file for encapsulated PostScript output */ + +@ This program need not produce machine-independent output, so we may +safely use floating-point arithmetic here. At most 64 characters +(32 pixel-bytes) are output on each line. + +@<Output row |k|...@>= +{@+register float conv=255.0/(float)d; register int x; + for (l=0; l<n; l++) { + x=(int)(conv*(float)(compl?d-aa(k,l):aa(k,l))); + fprintf(eps_file,"%02x",x>255?255:x); + if ((l&0x1f)==0x1f) fprintf(eps_file,"\n"); + } + if (n&0x1f) fprintf(eps_file,"\n"); +} + +@ @<Output the solution in PostScript format@>= +{ + fprintf(eps_file, + "/bx {moveto 0 1 rlineto 1 0 rlineto 0 -1 rlineto closepath\n"); + fprintf(eps_file," gsave .3 setlinewidth 1 setgray clip stroke"); + fprintf(eps_file," grestore stroke} bind def\n"); + fprintf(eps_file," .1 setlinewidth\n"); + for (k=0; k<m; k++) + fprintf(eps_file," %d %d bx\n",@| + transposed? k:matching_col[k],@| + transposed? n-1-matching_col[k]:m-1-k); + fclose(eps_file); +} + +@* Index. As usual, we close with a list of identifier definitions and uses. + |