% 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=}\, \.{m0=}\, 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 */ @# @@; main(argc,argv) int argc; /* the number of command-line arguments */ char *argv[]; /* an array of strings containing those arguments */ {@+@; @; 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) @; if (PostScript) @; mems=0; @; if (printing) @; if (PostScript) @; 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 @= 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: @= 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? */ @ @= 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; } } @ @= for (k=0;k0$. 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= { for (l=0; l0) for (k=0;k= 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. @= 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 */ @ @= 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) */ @= t=0; /* the forest starts out empty */ for (l=0; l1) 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. @= t=0; for (l=0; l1) 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)$. @= @; 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; q++; } @; } breakthru: @; if(--unmatched==0) goto done; @; } done: @; @ @= { o,k=unchosen_row[q]; o,s=row_dec[k]; for (l=0; l1) 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. @= 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. @= s=INF; for (l=0; l; } 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. @= { 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; j1) 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.) @= for (k=0;km) { 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. @= if (m>n) @@; else transposed=0; @; if (compl==0) for (k=0; k; @; @ @= { 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= long* tmtx; /* the transpose of |mtx| */ int transposed; /* has the data been transposed? */ @ @= { printf("The following entries produce an optimum assignment:\n"); for (k=0; k 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. @= { 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; fprintf(eps_file,"grestore\n"); } } @ @= 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. @= {@+register float conv=255.0/(float)d; register int x; for (l=0; l255?255:x); if ((l&0x1f)==0x1f) fprintf(eps_file,"\n"); } if (n&0x1f) fprintf(eps_file,"\n"); } @ @= { 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