ref: a409cfe16978d8e98a3754c6053cd00b09d32b26
dir: /matching.c/
/*
* Implementation of matching.h.
*/
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include "puzzles.h"
#include "matching.h"
struct scratch {
/*
* Current contents of the in-progress matching. LtoR is an array
* of nl integers, each of which holds a value in {0,1,...,nr-1},
* or -1 for no current assignment. RtoL is exactly the reverse.
*
* Invariant: LtoR[i] is non-empty and equal to j if and only if
* RtoL[j] is non-empty and equal to i.
*/
int *LtoR, *RtoL;
/*
* Arrays of nl and nr integer respectively, giving the layer
* assigned to each integer in the breadth-first search step of
* the algorithm.
*/
int *Llayer, *Rlayer;
/*
* Arrays of nl and nr integers respectively, used to hold the
* to-do queues in the breadth-first search.
*/
int *Lqueue, *Rqueue;
/*
* An augmenting path of vertices, alternating between L vertices
* (in the even-numbered positions, starting at 0) and R (in the
* odd positions). Must be long enough to hold any such path that
* never repeats a vertex, i.e. must be at least 2*min(nl,nr) in
* size.
*/
int *augpath;
/*
* Track the progress of the depth-first search at each
* even-numbered layer. Has one element for each even-numbered
* position in augpath.
*/
int *dfsstate;
/*
* Store a random permutation of the L vertex indices, if we're
* randomising the dfs phase.
*/
int *Lorder;
};
size_t matching_scratch_size(int nl, int nr)
{
size_t n;
int nmin = (nl < nr ? nl : nr);
n = (sizeof(struct scratch) + sizeof(int)-1)/sizeof(int);
n += nl; /* LtoR */
n += nr; /* RtoL */
n += nl; /* Llayer */
n += nr; /* Rlayer */
n += nl; /* Lqueue */
n += nr; /* Rqueue */
n += 2*nmin; /* augpath */
n += nmin; /* dfsstate */
n += nl; /* Lorder */
return n * sizeof(int);
}
int matching_with_scratch(void *scratchv,
int nl, int nr, int **adjlists, int *adjsizes,
random_state *rs, int *outl, int *outr)
{
struct scratch *s = (struct scratch *)scratchv;
int L, R, i, j;
/*
* Set up the various array pointers in the scratch space.
*/
{
int *p = scratchv;
int nmin = (nl < nr ? nl : nr);
p += (sizeof(struct scratch) + sizeof(int)-1)/sizeof(int);
s->LtoR = p; p += nl;
s->RtoL = p; p += nr;
s->Llayer = p; p += nl;
s->Rlayer = p; p += nr;
s->Lqueue = p; p += nl;
s->Rqueue = p; p += nr;
s->augpath = p; p += 2*nmin;
s->dfsstate = p; p += nmin;
s->Lorder = p; p += nl;
}
/*
* Set up the initial matching, which is empty.
*/
for (L = 0; L < nl; L++)
s->LtoR[L] = -1;
for (R = 0; R < nr; R++)
s->RtoL[R] = -1;
while (1) {
/*
* Breadth-first search starting from the unassigned left
* vertices, traversing edges from left to right only if they
* are _not_ part of the matching, and from right to left only
* if they _are_. We assign a 'layer number' to all vertices
* visited by this search, with the starting vertices being
* layer 0 and every successor of a layer-n node being layer
* n+1.
*/
int Lqs, Rqs, layer, target_layer;
for (L = 0; L < nl; L++)
s->Llayer[L] = -1;
for (R = 0; R < nr; R++)
s->Rlayer[R] = -1;
Lqs = 0;
for (L = 0; L < nl; L++) {
if (s->LtoR[L] == -1) {
s->Llayer[L] = 0;
s->Lqueue[Lqs++] = L;
}
}
layer = 0;
while (1) {
bool found_free_R_vertex = false;
Rqs = 0;
for (i = 0; i < Lqs; i++) {
L = s->Lqueue[i];
assert(s->Llayer[L] == layer);
for (j = 0; j < adjsizes[L]; j++) {
R = adjlists[L][j];
if (R != s->LtoR[L] && s->Rlayer[R] == -1) {
s->Rlayer[R] = layer+1;
s->Rqueue[Rqs++] = R;
if (s->RtoL[R] == -1)
found_free_R_vertex = true;
}
}
}
layer++;
if (found_free_R_vertex)
break;
if (Rqs == 0)
goto done;
Lqs = 0;
for (j = 0; j < Rqs; j++) {
R = s->Rqueue[j];
assert(s->Rlayer[R] == layer);
if ((L = s->RtoL[R]) != -1 && s->Llayer[L] == -1) {
s->Llayer[L] = layer+1;
s->Lqueue[Lqs++] = L;
}
}
layer++;
if (Lqs == 0)
goto done;
}
target_layer = layer;
/*
* Vertices in the target layer are only interesting if
* they're actually unassigned. Blanking out the others here
* will save us a special case in the dfs loop below.
*/
for (R = 0; R < nr; R++)
if (s->Rlayer[R] == target_layer && s->RtoL[R] != -1)
s->Rlayer[R] = -1;
/*
* Choose an ordering in which to try the L vertices at the
* start of the next pass.
*/
for (L = 0; L < nl; L++)
s->Lorder[L] = L;
if (rs)
shuffle(s->Lorder, nl, sizeof(*s->Lorder), rs);
/*
* Now depth-first search through that layered set of vertices
* to find as many (vertex-)disjoint augmenting paths as we
* can, and for each one we find, augment the matching.
*/
s->dfsstate[0] = 0;
i = 0;
while (1) {
/*
* Find the next vertex to go on the end of augpath.
*/
if (i == 0) {
/* In this special case, we're just looking for L
* vertices that are not yet assigned. */
if (s->dfsstate[i] == nl)
break; /* entire DFS has finished */
L = s->Lorder[s->dfsstate[i]++];
if (s->Llayer[L] != 2*i)
continue; /* skip this vertex */
} else {
/* In the more usual case, we're going through the
* adjacency list for the previous L vertex. */
L = s->augpath[2*i-2];
j = s->dfsstate[i]++;
if (j == adjsizes[L]) {
/* Run out of neighbours of the previous vertex. */
i--;
continue;
}
if (rs && adjsizes[L] - j > 1) {
int which = j + random_upto(rs, adjsizes[L] - j);
int tmp = adjlists[L][which];
adjlists[L][which] = adjlists[L][j];
adjlists[L][j] = tmp;
}
R = adjlists[L][j];
if (s->Rlayer[R] != 2*i-1)
continue; /* skip this vertex */
s->augpath[2*i-1] = R;
s->Rlayer[R] = -1; /* mark vertex as visited */
if (2*i-1 == target_layer) {
/*
* We've found an augmenting path, in the form of
* an even-sized list of vertices alternating
* L,R,...,L,R, with the initial L and final R
* vertex free and otherwise each R currently
* connected to the next L. Adjust so that each L
* connects to the next R, increasing the edge
* count in the matching by 1.
*/
for (j = 0; j < 2*i; j += 2) {
s->LtoR[s->augpath[j]] = s->augpath[j+1];
s->RtoL[s->augpath[j+1]] = s->augpath[j];
}
/*
* Having dealt with that path, and already marked
* all its vertices as visited, rewind right to
* the start and resume our DFS from a new
* starting L-vertex.
*/
i = 0;
continue;
}
L = s->RtoL[R];
if (s->Llayer[L] != 2*i)
continue; /* skip this vertex */
}
s->augpath[2*i] = L;
s->Llayer[L] = -1; /* mark vertex as visited */
i++;
s->dfsstate[i] = 0;
}
}
done:
/*
* Fill in the output arrays.
*/
if (outl) {
for (i = 0; i < nl; i++)
outl[i] = s->LtoR[i];
}
if (outr) {
for (j = 0; j < nr; j++)
outr[j] = s->RtoL[j];
}
/*
* Return the number of matching edges.
*/
for (i = j = 0; i < nl; i++)
if (s->LtoR[i] != -1)
j++;
return j;
}
int matching(int nl, int nr, int **adjlists, int *adjsizes,
random_state *rs, int *outl, int *outr)
{
void *scratch;
int size;
int ret;
size = matching_scratch_size(nl, nr);
scratch = malloc(size);
if (!scratch)
return -1;
ret = matching_with_scratch(scratch, nl, nr, adjlists, adjsizes,
rs, outl, outr);
free(scratch);
return ret;
}
void matching_witness(void *scratchv, int nl, int nr, int *witness)
{
struct scratch *s = (struct scratch *)scratchv;
int i, j;
for (i = 0; i < nl; i++)
witness[i] = s->Llayer[i] == -1;
for (j = 0; j < nr; j++)
witness[nl + j] = s->Rlayer[j] == -1;
}