assign_lisa.w

模拟器提供了一个简单易用的平台

|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...@>=
long* col_mate; /* the column matching a given row, or $-1$ */
long* row_mate; /* the row matching a given column, or $-1$ */
long* parent_row; /* ancestor of a given column's mate, or $-1$ */
long* unchosen_row; /* node in the forest */
long t; /* total number of nodes in the forest */
long 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 */
long* slack_row; /* where the |slack| in a given column can be found */
long unmatched; /* this many rows have yet to be matched */

@ @<Allocate the intermediate data structures@>=
col_mate=gb_typed_alloc(m,long,working_storage);
row_mate=gb_typed_alloc(n,long,working_storage);
parent_row=gb_typed_alloc(n,long,working_storage);
unchosen_row=gb_typed_alloc(m,long,working_storage);
row_dec=gb_typed_alloc(m,long,working_storage);
col_inc=gb_typed_alloc(n,long,working_storage);
slack=gb_typed_alloc(n,long,working_storage);
slack_row=gb_typed_alloc(n,long,working_storage);
if (gb_trouble_code) {
  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,row_mate[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,row_mate[l]<0)) {
      o,col_mate[k]=l;
      o,row_mate[l]=k;
      if (verbose>1) printf(" matching col %ld==row %ld\n",l,k);
      goto row_done;
    }
  o,col_mate[k]=-1;
  if (verbose>1) printf("  node %ld: unmatched row %ld\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,col_mate[k]<0) {
    if (verbose>1) printf("  node %ld: unmatched row %ld\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 %ld mems I've matched %ld 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,row_mate[l]<0) goto breakthru;
          o,slack[l]=0; /* this column will now be chosen */
          o,parent_row[l]=k;
          if (verbose>1) printf("  node %ld: row %ld==col %ld--row %ld\n",
                                    t,row_mate[l],l,k);
          oo,unchosen_row[t++]=row_mate[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 %ld of %ld!\n",q,t);
while (1) {
  o,j=col_mate[k];
  o,col_mate[k]=l;
  o,row_mate[l]=k;
  if (verbose>1) printf(" rematching col %ld==row %ld\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 might be several columns tied for smallest slack. If any of them
leads to a breakthrough, 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 might 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. Either we 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 %ld produces zero at [%ld,%ld]\n",
      s,k,l);
  if (o,row_mate[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 %ld: row %ld==col %ld--row %ld\n",
                                  t,row_mate[l],l,k);
    oo,unchosen_row[t++]=row_mate[l];
  }
}

@ The code in the present section is redundant, unless cosmic
radiation has caused 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.)
@^discussion of \\{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=col_mate[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_typed_alloc(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| */
long 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(" [%ld,%ld]\n",@|
     transposed? col_mate[k]:k,@|
     transposed? k:col_mate[k]);
}

@* Encapsulated PostScript.
A special output file called \.{lisa.eps} is written if the user has
selected the \.{-P} option. This file contains 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 with the
@.dvips@>
associated \.{epsf.tex} macros, one can say
$$\.{\\epsfxsize=10cm \\epsfbox\{lisa.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("lisa.eps","w");
  if (!eps_file) {
    fprintf(stderr,"Sorry, I can't open the file `lisa.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 %ld %ld\n",n+1,m+1);
    fprintf(eps_file,"/buffer %ld string def\n",n);
    fprintf(eps_file,"%ld %ld 8 [%ld 0 0 -%ld 0 %ld]\n",n,m,n,m,m);
    fprintf(eps_file,"{currentfile buffer readhexstring pop} bind\n");
    fprintf(eps_file,"gsave %ld %ld 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 can
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 long x;
  for (l=0; l<n; l++) {
    x=(long)(conv*(float)(compl?d-aa(k,l):aa(k,l)));
    fprintf(eps_file,"%02lx",x>255?255L: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," %ld %ld bx\n",@|
     transposed? k:col_mate[k],@|
     transposed? n-1-col_mate[k]:m-1-k);
  fclose(eps_file);
}

@* Index. As usual, we close with a list of identifier definitions and uses.