Fix -V

[clav.git] / mutation-finder.c

/*
 * Copyright (c) 2018, S. Gilles <sgilles@math.umd.edu>
 *
 * Permission to use, copy, modify, and/or distribute this software
 * for any purpose with or without fee is hereby granted, provided
 * that the above copyright notice and this permission notice appear
 * in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS *ALL
 * WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING AL*L IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EV*ENT SHALL THE
 * AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
 * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <errno.h>
#include <limits.h>
#include <locale.h>
#include <signal.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include "macros.h"
#include "quiver.h"

/* Flag for whether to print info */
static volatile sig_atomic_t status_requested = 0;

/* Signal handler */
static void
handle(int sig);

/* Signal handler */
static void handle(int sig)
{
        /* XXX: I'd like to check for SIGINFO here */
        if (sig == SIGUSR1) {
                status_requested = 1;
                signal(sig, handle);
        }
}

/* Print usage */
static void usage(char *argv0)
{
        printf("Usage: %s -s <startfile> -e <endfile>\n", argv0);
        printf("       [ -v ] # Be verbose\n");
        printf("       [ -R ] # Require all non-frozen in sequence\n");
        printf("       [ -V ] # Test via valency heuristic, not exactness\n");
        printf("       [ -c <edgecap> ]\n");
        printf("       [ -f <frozen> [ -f <frozen> [...] ] ]\n");
        printf("       [ -l <startlen> ] [ -m <maxlen> ]\n");
        printf("\n");
        printf("       startfile: a quiver file representing the start "
               "state\n");
        printf("       endfile:   a quiver file representing the desired end "
               "state\n");
        printf("       edgecap:   ignore mutations if the cause edges with"
               "weights above this\n");
        printf("       startlen:  start with sequences of this length\n");
        printf("       maxlen:    halt after exhausting sequences of this "
               "length\n");
        printf("       frozen:    a vertex with this name will not be "
               "mutated\n");
}

static int cmp(const void *a, const void *b)
{
        struct vertex **va = (struct vertex **) a;
        struct vertex **vb = (struct vertex **) b;

        return strcmp((*va)->name, (*vb)->name);
}

/* Dumb valency-based checksum for graphs. */
static uintmax_t valency_checksum(struct quiver const * const q)
{
        uintmax_t c_a = 0;
        uintmax_t c_b = 0;
        uintmax_t c_c = 0;
        uintmax_t c_d = 0;
        const size_t N = sizeof (uintmax_t) * sizeof (char) / 4;
        struct rational r = { 0 };

        /*
         * This looks really, really stupid, but there is a bit of
         * thought put into it. Most graphs on which this will run
         * have, say, <30 vertices, so each sum will probably fit
         * in quarter of a uintmax_t easily. That makes the or-ing
         * slightly relevant. Also, the whole point of this heuristic
         * is that we don't know how the graph could be shuffled,
         * so dependence on indices is frowned upon, which leaves
         * most serialization-based checksums out.
         */
        for (size_t j = 0; j < q->v_num; ++j) {
                for (size_t k = 0; k < q->v_num; ++k) {
                        quiver_get_sigma(q, j, k, &r);
                        c_a = c_a + r.p;
                        c_b = c_b + r.p * r.p;
                        c_c = c_c + r.p * r.q;
                        c_d = c_d + r.q * r.q;
                }
        }

        return (c_a << (3 * N)) | (c_b << (2 * N)) | (c_c << N) | c_d;
}

/* Load start/end quivers, figure out frozen vertices, order things, etc. */
static int setup(char *start_file, char *end_file, char **frozen, size_t
                 num_frozen, struct quiver *out_qstart, struct quiver *out_qend,
                 size_t *out_num_nonfrozen)
{
        struct quiver a = { 0 };
        struct quiver b = { 0 };
        FILE *fa = 0;
        FILE *fb = 0;
        const char *errstr = 0;
        int ret = 0;
        struct vertex **averts = 0;
        struct vertex **bverts = 0;
        size_t *a_to_start = 0;
        size_t *b_to_end = 0;

        if (!out_qstart ||
            !out_qend) {
                fprintf(stderr, L("invalid quivers"));
                ret = EINVAL;
                goto done;
        }

        if (out_qstart->v_num ||
            out_qend->v_num) {
                fprintf(stderr, L("output quivers are not empty\n"));
                ret = EINVAL;
                goto done;
        }

        if (!(fa = fopen(start_file, "r"))) {
                ret = errno;
                perror(L("fopen"));
                fprintf(stderr, "cannot open file %s\n", start_file);
                goto done;
        }

        if (!(fb = fopen(end_file, "r"))) {
                ret = errno;
                perror(L("fopen"));
                fprintf(stderr, "cannot open file %s\n", end_file);
                goto done;
        }

        if ((ret = quiver_load_from_file(&a, fa, &errstr))) {
                fprintf(stderr, "%s\n", errstr);
                goto done;
        }

        if ((ret = quiver_load_from_file(&b, fb, &errstr))) {
                fprintf(stderr, "%s\n", errstr);
                goto done;
        }

        if (a.v_num != b.v_num) {
                fprintf(stderr, L(
                                "Quivers do not have same number of vertices\n"));
                ret = EINVAL;
                goto done;
        }

        /*
         * Make sure there are no duplicate vertices, but don't sort in
         * place - that would change ordering, and ordering of vertices
         * must match user input - tuning input so that known prefixes
         * arise early is sometimes crucial.
         */
        if (!(averts = calloc(a.v_num, sizeof *averts))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        if (!(bverts = calloc(b.v_num, sizeof *bverts))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        for (size_t j = 0; j < a.v_num; ++j) {
                averts[j] = &(a.v[j]);
                bverts[j] = &(b.v[j]);
        }

        qsort(averts, a.v_num, sizeof(averts[0]), cmp);
        qsort(bverts, b.v_num, sizeof(bverts[0]), cmp);

        for (size_t j = 0; j < a.v_num; ++j) {
                if (j &&
                    !strcmp(averts[j]->name, averts[j - 1]->name)) {
                        fprintf(stderr, L("Duplicate vertex \"%s\"\n"),
                                averts[j]->name);
                        ret = EINVAL;
                        goto done;
                }

                int d = strcmp(averts[j]->name, bverts[j]->name);

                if (d > 0) {
                        fprintf(stderr, L("Vertex \"%s\" not in both files\n"),
                                bverts[j]->name);
                        ret = EINVAL;
                        goto done;
                }

                if (d < 0) {
                        fprintf(stderr, L("Vertex \"%s\" not in both files\n"),
                                averts[j]->name);
                        ret = EINVAL;
                        goto done;
                }

                if (averts[j]->fatness != bverts[j]->fatness) {
                        fprintf(stderr, L(
                                        "Vertex \"%s\" has inconsistent fatness\n"),
                                averts[j]->name);
                        ret = EINVAL;
                        goto done;
                }
        }

        /* a_to_s[i] = j <=> "vertex i in a is vertex j in out_qstart" */
        if (!(a_to_start = calloc(a.v_num, sizeof(*a_to_start)))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        if (!(b_to_end = calloc(b.v_num, sizeof(*b_to_end)))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        /*
         * Order so that frozen vertices are all at the end - required
         * for incrementation algorithm
         */
        for (uint_fast8_t want_frozen = 0; want_frozen <= 1; want_frozen++) {
                for (size_t j = 0; j < a.v_num; ++j) {
                        uint_fast8_t is_frozen = 0;

                        for (size_t k = 0; k < num_frozen; ++k) {
                                is_frozen = is_frozen ||
                                            !strcmp(a.v[j].name, frozen[k]);

                                if (is_frozen) {
                                        k = num_frozen;
                                }
                        }

                        if (is_frozen != want_frozen) {
                                continue;
                        }

                        if ((ret = quiver_add_vertex(out_qstart,
                                                     &(a_to_start[j]),
                                                     a.v[j].name,
                                                     a.v[j].fatness, 0, 0, 0,
                                                     &errstr))) {
                                fprintf(stderr, "%s\n", errstr);
                                goto done;
                        }

                        if (!want_frozen &&
                            out_num_nonfrozen) {
                                *out_num_nonfrozen = a_to_start[j] + 1;
                        }

                        size_t tj = 0;

                        if ((ret = quiver_add_vertex(out_qend, &tj, a.v[j].name,
                                                     a.v[j].fatness, 0, 0, 0,
                                                     &errstr))) {
                                fprintf(stderr, "%s\n", errstr);
                                goto done;
                        }

                        for (size_t k = 0; k < b.v_num; ++k) {
                                if (!(strcmp(b.v[k].name,
                                             out_qend->v[tj].name))) {
                                        b_to_end[k] = tj;
                                }
                        }
                }
        }

        /* Add the edges in correct order */
        for (size_t j = 0; j < a.v_num; ++j) {
                size_t tj = a_to_start[j];

                for (size_t k = 0; k < a.v_num; ++k) {
                        size_t tk = a_to_start[k];
                        size_t idx_o = tj * out_qstart->v_len + tk;
                        size_t idx_a = j * a.v_len + k;

                        out_qstart->e[idx_o] = a.e[idx_a];
                }
        }

        for (size_t j = 0; j < b.v_num; ++j) {
                size_t tj = b_to_end[j];

                for (size_t k = 0; k < b.v_num; ++k) {
                        size_t tk = b_to_end[k];
                        size_t idx_o = tj * out_qend->v_len + tk;
                        size_t idx_b = j * b.v_len + k;

                        out_qend->e[idx_o] = b.e[idx_b];
                }
        }

done:
        quiver_clean(&a);
        quiver_clean(&b);

        if (fa) {
                fclose(fa);
        }

        if (fb) {
                fclose(fb);
        }

        free(averts);
        free(bverts);
        free(a_to_start);
        free(b_to_end);

        return ret;
}

/* Main loop */
int main(int argc, char **argv)
{
        int ret = 0;
        struct quiver q = { 0 };
        struct quiver qend = { 0 };
        char *frozen[argc];
        size_t fidx = 0;
        char *start_file = 0;
        char *end_file = 0;
        long start_len_l = 0;
        long max_len_l = 0;
        size_t start_len = 0;
        size_t max_len = (size_t) -1;
        size_t warn_level = (size_t) -1;
        struct rational *target = 0;
        char *p = 0;
        const char *errstr = 0;
        uint_fast8_t verbose = 0;
        uint_fast8_t require_all_nonfrozen = 0;
        uint_fast8_t check_valency_only = 0;
        int *seen_in_seq = 0;
        uintmax_t desired_valency_checksum = 0;
        uintmax_t this_valency_checksum = 0;

        setlocale(LC_ALL, "");
        int opt = 0;

        for (int j = 0; j < argc; ++j) {
                frozen[j] = 0;
        }

        while ((opt = getopt(argc, argv, "vhRVs:e:c:m:l:f:")) != -1) {
                switch (opt) {
                case 'v':
                        verbose++;
                        break;
                case 'h':
                        usage(argv[0]);
                        ret = 0;
                        goto done;
                case 'R':
                        require_all_nonfrozen = 1;
                        break;
                case 'V':
                        check_valency_only = 1;
                        break;
                case 's':
                        start_file = optarg;
                        break;
                case 'e':
                        end_file = optarg;
                        break;
                case 'c':
                        warn_level = strtoll(optarg, &p, 0);

                        if (p &&
                            *p) {
                                fprintf(stderr, "Invalid edge cap: %s\n",
                                        optarg);
                                ret = EINVAL;
                                goto done;
                        }

                        break;
                case 'l':
                        start_len_l = strtoll(optarg, &p, 0);

                        if ((start_len_l <= 0) ||
                            (p &&
                             *p)) {
                                fprintf(stderr, "Invalid start length: %s\n",
                                        optarg);
                                ret = EINVAL;
                                goto done;
                        }

                        start_len = start_len_l;
                        break;
                case 'm':
                        max_len_l = strtoll(optarg, &p, 0);

                        if ((max_len_l <= 0) ||
                            (p &&
                             *p)) {
                                fprintf(stderr, "Invalid max length: %s\n",
                                        optarg);
                                ret = EINVAL;
                                goto done;
                        }

                        max_len = max_len_l;
                        break;
                case 'f':
                        frozen[fidx] = optarg;
                        fidx++;
                        break;
                default:
                        usage(argv[0]);
                        ret = EINVAL;
                        goto done;
                }
        }

        if (!start_file ||
            !end_file) {
                usage(argv[0]);
                ret = EINVAL;
                goto done;
        }

        size_t num_nonfrozen = 0;

        if (setup(start_file, end_file, frozen, fidx, &q, &qend,
                  &num_nonfrozen)) {
                goto done;
        }

        if (warn_level == (size_t) -1) {
                warn_level = 0;

                for (size_t k = 0; k < q.v_num; ++k) {
                        if (warn_level < q.v[k].fatness) {
                                warn_level = q.v[k].fatness;
                        }
                }

                warn_level = warn_level + 1;
        }

        if (warn_level >= UINT_MAX) {
                warn_level = UINT_MAX - 1;
        }

        if ((ret = quiver_set_warning_threshold(&q, warn_level, &errstr))) {
                fprintf(stderr, "%s\n", errstr);
                goto done;
        }

        /* Save off qend's vertices.  XXX: Why not just USE qend? */
        if ((qend.v_num * qend.v_num) / qend.v_num != qend.v_num) {
                /* Should never happen - isn't qend.v_num capped? */
                ret = errno = EOVERFLOW;
                perror(L(""));
                goto done;
        }

        if (!(target = calloc(qend.v_num * qend.v_num, sizeof *target))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        if (!(seen_in_seq = calloc(qend.v_num, sizeof *seen_in_seq))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        for (size_t i = 0; i < qend.v_num * qend.v_num; ++i) {
                target[i] = (struct rational) { .p = 0, .q = 1 };
        }

        for (size_t i = 0; i < qend.v_num; ++i) {
                for (size_t j = 0; j < qend.v_num; ++j) {
                        target[i * qend.v_num + j] = qend.e[i * qend.v_len + j];
                }
        }

        if (check_valency_only) {
                desired_valency_checksum = valency_checksum(&qend);
        }

        /* Now, the long haul*/
        size_t len = 0;

        if (start_len) {
                len = start_len - 1;
        }

        signal(SIGUSR1, handle);
        size_t s = 0;
        size_t c_idx;
        int warn = 0;
        uint_fast8_t correct = 0;

        while (len < max_len) {
                len++;
                size_t sequence[len];

                for (size_t k = 0; k < num_nonfrozen; ++k) {
                        seen_in_seq[k] = 0;
                }

                for (size_t k = 0; k < len; ++k) {
                        /* This will get pruned, it's okay */
                        sequence[k] = 0;
                        seen_in_seq[sequence[k]]++;
                }

                c_idx = 0;
have_full_sequence:

                if (status_requested) {
                        status_requested = 0;
                        printf("Current sequence %s", q.v[sequence[0]].name);

                        for (size_t k = 1; k < len; ++k) {
                                printf(", %s", q.v[sequence[k]].name);
                        }

                        printf("\n");
                        fflush(stdout);
                }

                /*
                 * Here, we have computed a full sequence that we want
                 * to check.  c_idx is in [0, len - 1], and the quiver
                 * has been mutated following the sequence up to (but
                 * not including) c_idx.
                 *
                 * This code is duplicated from increment_sequence
                 * because it is reasonable to start the whole process
                 * with a sequence that may, in fact, produce warnings.
                 */
                warn = 0;

                /*
                 * If we are looking for a flip, we're probably
                 * going to require that all non-frozen vertices
                 * participate in the mutation.
                 */
                if (require_all_nonfrozen &&
                    c_idx == len - 1) {
                        for (size_t j = 0; j < num_nonfrozen; ++j) {
                                if (!seen_in_seq[j]) {
                                        goto increment_sequence;
                                }
                        }
                }

                while (c_idx < len) {
                        /* Prune duplicates (mutations are involutions) */
                        if (c_idx &&
                            sequence[c_idx] == sequence[c_idx - 1]) {
                                goto increment_sequence;
                        }

                        /* Proceed one more mutation */
                        quiver_mutate(&q, sequence[c_idx], &warn, 0);

                        if (warn) {
                                /* That mutation caused a warning: skip it */
                                quiver_mutate(&q, sequence[c_idx], 0, 0);

                                /* And try the sequence, starting at c_idx */
                                goto increment_sequence;
                        }

                        /*
                         * Skip ... X, Y ... when X > Y and the edge
                         * from X to Y is trivial
                         */
                        if (c_idx &&
                            sequence[c_idx] < sequence[c_idx - 1]) {
                                size_t x = sequence[c_idx - 1];
                                size_t y = sequence[c_idx];
                                struct rational *xy = &(q.e[x * q.v_len + y]);
                                struct rational *yx = &(q.e[y * q.v_len + x]);

                                if (!xy->p &&
                                    !yx->p) {
                                        /* Undo mutation */
                                        quiver_mutate(&q, sequence[c_idx], 0,
                                                      0);

                                        /* And try the sequence, starting at c_idx */
                                        goto increment_sequence;
                                }
                        }

                        c_idx++;
                }

                /*
                 * Here, there's a full sequence, and the quiver has
                 * been mutated according to it (all of it).  We want
                 * to check if the result matches the target.
                 */
                if (check_valency_only) {
                        /* Take a short-cut */
                        this_valency_checksum = valency_checksum(&q);

                        if (this_valency_checksum != desired_valency_checksum) {
                                goto finished_check;
                        }
                }

                correct = 1;

                for (size_t i = 0; i < q.v_num; ++i) {
                        for (size_t j = 0; j < q.v_num; ++j) {
                                struct rational *e = &(q.e[i * q.v_len + j]);
                                struct rational *f = &(target[i * q.v_num + j]);

                                if (e->p * f->q != f->p * e->q) {
                                        correct = 0;
                                        i = q.v_num + 1;
                                        j = i;
                                }
                        }
                }

                if (correct ||
                    check_valency_only) {
                        printf("SEQUENCE:     %s", q.v[sequence[0]].name);

                        for (size_t k = 1; k < len; ++k) {
                                printf(", %s", q.v[sequence[k]].name);
                        }

                        printf("\n");
                        fflush(stdout);
                }

finished_check:
                c_idx = len - 1;
                quiver_mutate(&q, sequence[c_idx], 0, 0);
increment_sequence:

                /*
                 * Here, we have some kind of sequence. The quiver agrees
                 * with it up to (but not including) c_idx, and we want
                 * to increment it at c_idx.
                 */
                s++;
                seen_in_seq[sequence[c_idx]]--;
                sequence[c_idx]++;
                seen_in_seq[sequence[c_idx]]++;

                if (verbose &&
                    c_idx == 2) {
                        double completed = sequence[0] * num_nonfrozen *
                                           num_nonfrozen + sequence[1] *
                                           num_nonfrozen +
                                           sequence[2] - 1;
                        double total = num_nonfrozen * num_nonfrozen *
                                       num_nonfrozen;

                        printf("Length %zu: completed %.3g%%\n", len, (100.0 *
                                                                       completed)
                               / total);
                        fflush(stdout);
                }

                /* Skip over repeats: each mutation is an involution */
                if (c_idx > 0 &&
                    sequence[c_idx] == sequence[c_idx - 1]) {
                        seen_in_seq[sequence[c_idx]]--;
                        sequence[c_idx]++;
                        seen_in_seq[sequence[c_idx]]++;
                }

                /* Wrap back to 0 */
                if (sequence[c_idx] >= num_nonfrozen) {
                        if (!c_idx) {
                                goto exhausted_len;
                        }

                        c_idx--;
                        quiver_mutate(&q, sequence[c_idx], 0, 0);
                        goto increment_sequence;
                }

                if (c_idx != len - 1) {
                        /*
                         * Here, we now have a start of a sequence.
                         * The quiver agrees with our start up to (but
                         * not including) c_idx, and we want to possibly
                         * complete this to a full sequence. First,
                         * however, we need to make sure that we don't
                         * introduce any warnings at c_idx.
                         *
                         * This single-mutation testing block could be
                         * ignored, but I have a vague suspicion it will
                         * be a hot path.
                         */
                        warn = 0;
                        quiver_mutate(&q, sequence[c_idx], &warn, 0);

                        if (warn) {
                                /* That mutation caused a warning: skip it */
                                quiver_mutate(&q, sequence[c_idx], 0, 0);

                                /* And try the sequence, starting at c_idx */
                                goto increment_sequence;
                        }

                        /*
                         * Quivers almost commute: in particular if
                         * the edge v_i -> v_j is zero, then mutation
                         * at v_i and v_j commute. So if we're at X, Y
                         * in the list, and X > Y, then we've actually
                         * already computed that branch of the tree,
                         * back when it was Y, X. Therefore, increment X.
                         */
                        if (c_idx &&
                            sequence[c_idx] < sequence[c_idx - 1]) {
                                size_t x = sequence[c_idx - 1];
                                size_t y = sequence[c_idx];
                                struct rational *xy = &(q.e[x * q.v_len + y]);
                                struct rational *yx = &(q.e[y * q.v_len + x]);

                                if (!xy->p &&
                                    !yx->p) {
                                        /* Undo mutation */
                                        quiver_mutate(&q, sequence[c_idx], 0,
                                                      0);

                                        /* And try the sequence, starting at c_idx */
                                        goto increment_sequence;
                                }
                        }

                        /*
                         * Here the quiver agrees with our start up
                         * to and including c_idx, and we need to fill
                         * it out into a full sequence so that it can
                         * be tested.  We don't need to perform all the
                         * mutations of that full sequence, but we do
                         * need to generate the minimal one.
                         */
                        for (size_t i = c_idx + 1; i < len; ++i) {
                                seen_in_seq[sequence[i]]--;
                                sequence[i] = 1 - !!((i - c_idx) % 2);
                                seen_in_seq[sequence[i]]++;
                        }

                        c_idx++;
                }

                goto have_full_sequence;
exhausted_len:

                /*
                 * Here, we have c_idx = 0, and the quiver has not been
                 * mutated by anything in sequence
                 */
                printf("Exhausted length %zu\n", len);
        }

        fflush(stdout);
done:
        free(target);
        free(seen_in_seq);

        return ret;
}