puzzles/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;
}

#ifdef STANDALONE_MATCHING_TEST

/*
 * Diagnostic routine used in testing this algorithm. It is passed a
 * pointer to a piece of scratch space that's just been used by
 * matching_with_scratch, and extracts from it a labelling of the
 * input graph that acts as a 'witness' to the maximality of the
 * returned matching.
 *
 * The output parameter 'witness' should be an array of (nl+nr)
 * integers, indexed such that witness[L] corresponds to an L-vertex (for
 * L=0,1,...,nl-1) and witness[nl+R] corresponds to an R-vertex (for
 * R=0,1,...,nr-1). On return, this array will assign each vertex a
 * label which is either 0 or 1, and the following properties will
 * hold:
 *
 *  + all vertices not paired up by the matching are type L0 or R1
 *  + every L0->R1 edge is used by the matching
 *  + no L1->R0 edge is used by the matching.
 *
 * The mere existence of such a labelling is enough to prove the
 * maximality of the matching, because if there is any larger matching
 * then its symmetric difference with this one must include at least
 * one 'augmenting path', which starts at a free L-vertex and ends at
 * a free R-vertex, traversing only unused L->R edges and only used
 * R->L edges. But that would mean it starts at an L0, ends at an R1,
 * and never follows an edge that can get from an 0 to a 1.
 */
static 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;
}

/*
 * Standalone tool to run the matching algorithm.
 */

#include <string.h>
#include <ctype.h>
#include <time.h>

#include "tree234.h"

int nl, nr, count;
int **adjlists, *adjsizes;
int *adjdata, *outl, *outr, *witness;
void *scratch;
random_state *rs;

void allocate(int nl_, int nr_, int maxedges)
{
    nl = nl_;
    nr = nr_;
    adjdata = snewn(maxedges, int);
    adjlists = snewn(nl, int *);
    adjsizes = snewn(nl, int);
    outl = snewn(nl, int);
    outr = snewn(nr, int);
    witness = snewn(nl+nr, int);
    scratch = smalloc(matching_scratch_size(nl, nr));
}

void deallocate(void)
{
    sfree(adjlists);
    sfree(adjsizes);
    sfree(adjdata);
    sfree(outl);
    sfree(outr);
    sfree(witness);
    sfree(scratch);
}

void find_and_check_matching(void)
{
    int i, j, k;

    count = matching_with_scratch(scratch, nl, nr, adjlists, adjsizes,
                                  rs, outl, outr);
    matching_witness(scratch, nl, nr, witness);

    for (i = j = 0; i < nl; i++) {
        if (outl[i] != -1) {
            assert(0 <= outl[i] && outl[i] < nr);
            assert(outr[outl[i]] == i);
            j++;

            for (k = 0; k < adjsizes[i]; k++)
                if (adjlists[i][k] == outl[i])
                    break;
            assert(k < adjsizes[i]);
        }
    }
    assert(j == count);

    for (i = j = 0; i < nr; i++) {
        if (outr[i] != -1) {
            assert(0 <= outr[i] && outr[i] < nl);
            assert(outl[outr[i]] == i);
            j++;
        }
    }
    assert(j == count);

    for (i = 0; i < nl; i++) {
        if (outl[i] == -1)
            assert(witness[i] == 0);
    }
    for (i = 0; i < nr; i++) {
        if (outr[i] == -1)
            assert(witness[nl+i] == 1);
    }
    for (i = 0; i < nl; i++) {
        for (j = 0; j < adjsizes[i]; j++) {
            k = adjlists[i][j];

            if (outl[i] == k)
                assert(!(witness[i] == 1 && witness[nl+k] == 0));
            else
                assert(!(witness[i] == 0 && witness[nl+k] == 1));
        }
    }
}

struct nodename {
    const char *name;
    int index;
};

int compare_nodes(void *av, void *bv)
{
    const struct nodename *a = (const struct nodename *)av;
    const struct nodename *b = (const struct nodename *)bv;
    return strcmp(a->name, b->name);
}

int node_index(tree234 *n2i, tree234 *i2n, const char *name)
{
    struct nodename *nn, *nn_prev;
    char *namedup = dupstr(name);

    nn = snew(struct nodename);
    nn->name = namedup;
    nn->index = count234(n2i);

    nn_prev = add234(n2i, nn);
    if (nn_prev != nn) {
        sfree(nn);
        sfree(namedup);
    } else {
        addpos234(i2n, nn, nn->index);
    }

    return nn_prev->index;
}

struct edge {
    int L, R;
};

int compare_edges(void *av, void *bv)
{
    const struct edge *a = (const struct edge *)av;
    const struct edge *b = (const struct edge *)bv;
    if (a->L < b->L) return -1;
    if (a->L > b->L) return +1;
    if (a->R < b->R) return -1;
    if (a->R > b->R) return +1;
    return 0;
}

void matching_from_user_input(FILE *fp, const char *filename)
{
    tree234 *Ln2i, *Li2n, *Rn2i, *Ri2n, *edges;
    char *line = NULL;
    struct edge *e;
    int i, lineno = 0;
    int *adjptr;

    Ln2i = newtree234(compare_nodes);
    Rn2i = newtree234(compare_nodes);
    Li2n = newtree234(NULL);
    Ri2n = newtree234(NULL);
    edges = newtree234(compare_edges);

    while (sfree(line), lineno++, (line = fgetline(fp)) != NULL) {
        char *p, *Lname, *Rname;

        p = line;
        while (*p && isspace((unsigned char)*p)) p++;
        if (!*p)
            continue;

        Lname = p;
        while (*p && !isspace((unsigned char)*p)) p++;
        if (*p)
            *p++ = '\0';
        while (*p && isspace((unsigned char)*p)) p++;

        if (!*p) {
            fprintf(stderr, "%s:%d: expected 2 words, found 1\n",
                    filename, lineno);
            exit(1);
        }

        Rname = p;
        while (*p && !isspace((unsigned char)*p)) p++;
        if (*p)
            *p++ = '\0';
        while (*p && isspace((unsigned char)*p)) p++;

        if (*p) {
            fprintf(stderr, "%s:%d: expected 2 words, found more\n",
                    filename, lineno);
            exit(1);
        }

        e = snew(struct edge);
        e->L = node_index(Ln2i, Li2n, Lname);
        e->R = node_index(Rn2i, Ri2n, Rname);
        if (add234(edges, e) != e) {
            fprintf(stderr, "%s:%d: duplicate edge\n",
                    filename, lineno);
            exit(1);
        }
    }

    allocate(count234(Ln2i), count234(Rn2i), count234(edges));

    adjptr = adjdata;
    for (i = 0; i < nl; i++)
        adjlists[i] = NULL;
    for (i = 0; (e = index234(edges, i)) != NULL; i++) {
        if (!adjlists[e->L])
            adjlists[e->L] = adjptr;
        *adjptr++ = e->R;
        adjsizes[e->L] = adjptr - adjlists[e->L];
    }

    find_and_check_matching();

    for (i = 0; i < nl; i++) {
        if (outl[i] != -1) {
            struct nodename *Lnn = index234(Li2n, i);
            struct nodename *Rnn = index234(Ri2n, outl[i]);
            printf("%s %s\n", Lnn->name, Rnn->name);
        }
    }
}

void test_subsets(void)
{
    int b = 8;
    int n = 1 << b;
    int i, j, nruns, expected_size;
    int *adjptr;
    int *edgecounts;
    struct stats {
        int min, max;
        double n, sx, sxx;
    } *stats;
    static const char seed[] = "fixed random seed for repeatability";

    /*
     * Generate a graph in which every subset of [b] = {1,...,b}
     * (represented as a b-bit integer 0 <= i < n) has an edge going
     * to every subset obtained by removing exactly one element.
     *
     * This graph is the disjoint union of the corresponding graph for
     * each layer (collection of same-sized subset) of the power set
     * of [b]. Each of those graphs has a matching of size equal to
     * the smaller of its vertex sets. So we expect the overall size
     * of the output matching to be less than n by the size of the
     * largest layer, that is, to be n - binomial(n, floor(n/2)).
     *
     * We run the generation repeatedly, randomising it every time,
     * and we expect to see every possible edge appear sooner or
     * later.
     */

    rs = random_new(seed, strlen(seed));

    allocate(n, n, n*b);
    adjptr = adjdata;
    expected_size = 0;
    for (i = 0; i < n; i++) {
        adjlists[i] = adjptr;
        for (j = 0; j < b; j++) {
            if (i & (1 << j))
                *adjptr++ = i & ~(1 << j);
        }
        adjsizes[i] = adjptr - adjlists[i];
        if (adjsizes[i] != b/2)
            expected_size++;
    }

    edgecounts = snewn(n*b, int);
    for (i = 0; i < n*b; i++)
        edgecounts[i] = 0;

    stats = snewn(b, struct stats);

    nruns = 0;
    while (nruns < 10000) {
        nruns++;
        find_and_check_matching();
        assert(count == expected_size);

        for (i = 0; i < n; i++)
            for (j = 0; j < b; j++)
                if ((i ^ outl[i]) == (1 << j))
                    edgecounts[b*i+j]++;

        if (nruns % 1000 == 0) {
            for (i = 0; i < b; i++) {
                struct stats *st = &stats[i];
                st->min = st->max = -1;
                st->n = st->sx = st->sxx = 0;
            }

            for (i = 0; i < n; i++) {
                int pop = 0;
                for (j = 0; j < b; j++)
                    if (i & (1 << j))
                        pop++;
                pop--;

                for (j = 0; j < b; j++) {
                    if (i & (1 << j)) {
                        struct stats *st = &stats[pop];
                        int x = edgecounts[b*i+j];
                        if (st->max == -1 || st->max < x)
                            st->max = x;
                        if (st->min == -1 || st->min > x)
                            st->min = x;
                        st->n++;
                        st->sx += x;
                        st->sxx += (double)x*x;
                    } else {
                        assert(edgecounts[b*i+j] == 0);
                    }
                }
            }
        }
    }

    printf("after %d runs:\n", nruns);
    for (j = 0; j < b; j++) {
        struct stats *st = &stats[j];
        printf("edges between layers %d,%d:"
               " min=%d max=%d mean=%f variance=%f\n",
               j, j+1, st->min, st->max, st->sx/st->n,
               (st->sxx - st->sx*st->sx/st->n) / st->n);
    }
}

int main(int argc, char **argv)
{
    static const char stdin_identifier[] = "<standard input>";
    const char *infile = NULL;
    bool doing_opts = true;
    enum { USER_INPUT, AUTOTEST } mode = USER_INPUT;

    while (--argc > 0) {
        const char *arg = *++argv;

        if (doing_opts && arg[0] == '-' && arg[1]) {
            if (!strcmp(arg, "--")) {
                doing_opts = false;
            } else if (!strcmp(arg, "--random")) {
                char buf[64];
                int len = sprintf(buf, "%lu", (unsigned long)time(NULL));
                rs = random_new(buf, len);
            } else if (!strcmp(arg, "--autotest")) {
                mode = AUTOTEST;
            } else {
                fprintf(stderr, "matching: unrecognised option '%s'\n", arg);
                return 1;
            }
        } else {
            if (!infile) {
                infile = (!strcmp(arg, "-") ? stdin_identifier : arg);
            } else {
                fprintf(stderr, "matching: too many arguments\n");
                return 1;
            }
        }
    }

    if (mode == USER_INPUT) {
        FILE *fp;

        if (!infile)
            infile = stdin_identifier;

        if (infile != stdin_identifier) {
            fp = fopen(infile, "r");
            if (!fp) {
                fprintf(stderr, "matching: could not open input file '%s'\n",
                        infile);
                return 1;
            }
        } else {
            fp = stdin;
        }

        matching_from_user_input(fp, infile);

        if (infile != stdin_identifier)
            fclose(fp);
    }

    if (mode == AUTOTEST) {
        if (infile) {
            fprintf(stderr, "matching: expected no filename argument "
                    "with --autotest\n");
            return 1;
        }

        test_subsets();
    }

    return 0;
}

#endif /* STANDALONE_MATCHING_TEST */