Remove typos in license header

[clav.git] / check-graph-iso.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 ALL IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT 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 <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

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

#define zmax(x, y) (((x) > (y)) ? (x) : (y))

static int
cmp(const void *a, const void *b)
{
        const size_t *as = (const size_t *) a;
        const size_t *bs = (const size_t *) b;

        if (*as < *bs) {
                return -1;
        } else if (*as > *bs) {
                return 1;
        }

        return 0;
}

int
check_iso(struct quiver *a, struct quiver *b, const char **fix, size_t fix_len)
{
        int ret = -1;
        size_t n = a->v_num;
        size_t *a_i = 0;
        size_t *b_i = 0;
        size_t i_idx = 0;
        size_t ia = 0;
        size_t ib = 0;
        size_t a_next = 0;
        size_t b_next = 0;
        size_t next_highest = 0;
        size_t next_highest_idx = 0;
        size_t t = 0;

        /*
           a_i is { f₁, f₂, …, fₘ, u₁, u₂, …, uₙ₋ₘ }, where f's are
           fixed indices and u₁ < u₂ < … < uₙ₋ₘ .

           b_i will be the same, but as the sequence progresses the
           u's will be permuted, and we require always that a and
           b, when restricted to indices less than i_idx, are
           isomorphic.
         */
        if (!(a_i = calloc(n, sizeof *a_i))) {
                perror(L("calloc"));
                goto done;
        }

        if (!(b_i = calloc(n, sizeof *b_i))) {
                perror(L("calloc"));
                goto done;
        }

        for (size_t j = 0; j < fix_len; ++j) {
                const char *name = fix[j];

                for (size_t k = 0; k < n; ++k) {
                        if (!(strcmp(a->v[k].name, name))) {
                                a_i[j] = k;
                        }

                        if (!(strcmp(b->v[k].name, name))) {
                                b_i[j] = k;
                        }
                }
        }

        /* Fill out the permutation to make u's increasing */
        for (i_idx = fix_len; i_idx < n; ++i_idx) {
                /* This is incredibly inefficient, but it's readable */
                int seen_a_next = 0;
                int seen_b_next = 0;

                /* Find smallest a_next value that's not already used */
                for (a_next = 0; a_next < n; ++a_next) {
                        seen_a_next = 0;

                        for (size_t k = 0; k < i_idx; ++k) {
                                if (a_i[k] == a_next) {
                                        seen_a_next = 1;
                                }
                        }

                        if (!seen_a_next) {
                                a_i[i_idx] = a_next;
                                break;
                        }
                }

                /* The same; for b */
                for (b_next = 0; b_next < n; ++b_next) {
                        seen_b_next = 0;

                        for (size_t k = 0; k < i_idx; ++k) {
                                if (b_i[k] == b_next) {
                                        seen_b_next = 1;
                                }
                        }

                        if (!seen_b_next) {
                                b_i[i_idx] = b_next;
                                break;
                        }
                }
        }

        i_idx = fix_len;

        /*
         * Require that, restricted to the fixed parts, we have a
         * graph isomorphism.
         */
        for (size_t j = 0; j < fix_len; ++j) {
                size_t ja = a_i[j];
                size_t jb = b_i[j];

                if (a->v[ja].fatness != b->v[jb].fatness) {
                        fprintf(stderr,
                                "Fixed vertex \u00ab%s\u00bb has inconsistent fatness %d, %d\n",
                                a->v[ja].name, a->v[ja].fatness,
                                b->v[ja].fatness);
                        goto done;
                }

                for (size_t k = 0; k < fix_len; ++k) {
                        size_t ka = a_i[k];
                        size_t kb = b_i[k];
                        struct rational *ea = &a->e[ja * a->v_len + ka];
                        struct rational *eb = &b->e[jb * b->v_len + kb];

                        if (ea->p != eb->p ||
                            ea->q != eb->q) {
                                fprintf(stderr,
                                        "S: \u00ab%s\u00bb \u2192 \u00ab%s\u00bb, %d/%d\n",
                                        a->v[ja].name, a->v[ka].name, ea->p,
                                        ea->q);
                                fprintf(stderr,
                                        "T: \u00ab%s\u00bb \u2192 \u00ab%s\u00bb, %d/%d\n",
                                        b->v[jb].name, b->v[kb].name, eb->p,
                                        eb->q);
                                fprintf(stderr,
                                        "These are fixed vertices; so ");
                                fprintf(stderr, "no isomorphism is possible\n");
                                goto done;
                        }
                }
        }

        /* The long haul */
        while (1) {
check_partial_isomorphism:

                /*
                 * At this point, a_i and b_i are valid permutations
                 * (completely, up to the last position). Furthermore, the
                 * map Gₐ[a_i[k] <-> Gᵦ[b_i[k]], restricted to k < i_idx,
                 * gives a graph isomorphism. We seek to check if incrementing
                 * i_idx by one works.
                 */
                ia = a_i[i_idx];
                ib = b_i[i_idx];

                if (a->v[ia].fatness != b->v[ib].fatness) {
                        goto increment_i_idxth_place;
                }

                for (size_t j = 0; j < i_idx; ++j) {
                        size_t ja = a_i[j];
                        size_t jb = b_i[j];
                        struct rational *ea = &a->e[ia * a->v_len + ja];
                        struct rational *eb = &b->e[ib * b->v_len + jb];
                        struct rational *fa = &a->e[ja * a->v_len + ia];
                        struct rational *fb = &b->e[jb * b->v_len + ib];

                        if (ea->p != eb->p ||
                            ea->q != eb->q ||
                            fa->p != fb->p ||
                            fa->q != fb->q) {
                                goto increment_i_idxth_place;
                        }
                }

                /*
                 * Expanding the isomorphism works. Good, let's try
                 * to expand a bit more
                 */
                i_idx++;

                if (i_idx < n) {
                        goto check_partial_isomorphism;
                }

                /* The whole thing is an isomorphism. Good, we're done. */
                printf("Isomorphism found\n");
                printf("     S      |      T     \n");
                printf("------------|------------\n");

                for (size_t j = 0; j < n; ++j) {
                        size_t ja = a_i[j];
                        size_t jb = b_i[j];

                        printf("%11s | %s\n", a->v[ja].name, b->v[jb].name);
                }

                ret = 0;
                goto done;
increment_i_idxth_place:

                /*
                 * Try to get the next in the orderings of b_i.
                 * It's a permutation, so we're in the easy case
                 * of the incrementation algorithm.
                 */
                next_highest = (size_t) -1;
                next_highest_idx = i_idx;

                for (size_t j = i_idx + 1; j < n; ++j) {
                        if (b_i[j] < next_highest &&
                            b_i[j] > b_i[i_idx]) {
                                next_highest = b_i[j];
                                next_highest_idx = j;
                        }
                }

                if (next_highest == (size_t) -1) {
                        if (i_idx <= fix_len) {
                                /* We've checked all we can */
                                printf("No isomorphism possible\n");
                                goto done;
                        }

                        i_idx--;
                        goto increment_i_idxth_place;
                }

                t = b_i[i_idx];
                b_i[i_idx] = b_i[next_highest_idx];
                b_i[next_highest_idx] = t;
                qsort(b_i + i_idx + 1, n - i_idx - 1, sizeof *b_i, cmp);
        }

done:
        free(a_i);
        free(b_i);

        return ret;
}

static void
usage(char *argv0)
{
        printf("Usage: %s -s <file> -t <file> [ -f <v> [ -f <v> [ ... ] ] ]\n",
               argv0);
        printf("       <file>: a quiver file output by clav\n");
        printf("       <v>: a name of a vertex in both files\n");
        printf("\n");
        printf("Find (and possibly output) some graph isomorphism between s\n");
        printf("and t, fixing every vertex named by ‘-f’.\n");
        printf("\n");
        printf("Returns 0 if such an isomorphism was found, < 0 if not.\n");
}

/* Run the thing */
int
main(int argc, char **argv)
{
        int ret = -1;
        int opt = 0;
        const char *sarg = 0;
        const char *targ = 0;
        struct quiver qs = { 0 };
        struct quiver qt = { 0 };
        FILE *fs = 0;
        FILE *ft = 0;
        const char **fixnames = 0;
        const char *errstr = 0;
        size_t fidx = 0;

        if (!(fixnames = calloc(1 + (argc / 2), sizeof *fixnames))) {
                perror(L("calloc"));
                goto done;
        }

        while ((opt = getopt(argc, argv, "hs:t:f:")) != -1) {
                switch (opt) {
                case 's':
                        sarg = optarg;
                        break;
                case 't':
                        targ = optarg;
                        break;
                case 'f':
                        fixnames[fidx] = optarg;
                        fidx++;
                        break;
                case 'h':
                        ret = 0;

                /* fall through */
                default:
                        usage(argv[0]);
                        goto done;
                }
        }

        if (!sarg) {
                fprintf(stderr, "\u2018-s\u2019 is required\n");
                usage(argv[0]);
                goto done;
        }

        if (!targ) {
                fprintf(stderr, "\u2018-t\u2019 is required\n");
                usage(argv[0]);
                goto done;
        }

        if (!(fs = fopen(sarg, "r"))) {
                perror(L("fopen"));
                goto done;
        }

        if ((ret = quiver_load_from_file(&qs, fs, &errstr))) {
                fprintf(stderr, "cannot load \u00ab%s\u00bb: %s\n", sarg,
                        errstr);
                goto done;
        }

        fclose(fs);
        fs = 0;

        if (!(ft = fopen(targ, "r"))) {
                perror(L("fopen"));
                goto done;
        }

        if ((ret = quiver_load_from_file(&qt, ft, &errstr))) {
                fprintf(stderr, "cannot load \u00ab%s\u00bb: %s\n", sarg,
                        errstr);
                goto done;
        }

        fclose(ft);
        ft = 0;

        if (qs.v_num != qt.v_num) {
                fprintf(stderr,
                        "\u00ab%s\u00bb and \u00ab%s\u00bb have different |V|\n",
                        sarg,
                        targ);
                goto done;
        }

        /* Make sure all the fixed vertices are present exactly once */
        for (size_t j = 0; j < fidx; ++j) {
                const char *f = fixnames[j];
                int seen_in_s = 0;
                int seen_in_t = 0;

                for (size_t k = 0; k < qs.v_num; ++k) {
                        seen_in_s += !strcmp(qs.v[k].name, f);
                        seen_in_t += !strcmp(qt.v[k].name, f);
                }

                if (seen_in_s != 1 ||
                    seen_in_t != 1) {
                        fprintf(stderr, "\u00ab%s\u00bb appears\n", f);
                        fprintf(stderr, "  %d time%s in \u00ab%s\u00bb\n",
                                seen_in_s, (seen_in_s == 1) ? " " : "s", sarg);
                        fprintf(stderr, "  %d time%s in \u00ab%s\u00bb\n",
                                seen_in_t, (seen_in_t == 1) ? " " : "s", targ);
                        goto done;
                }
        }

        /* Make sure all vertices in all files are distinct */
        for (size_t j = 0; j < qs.v_num; ++j) {
                for (size_t k = 0; k < j; ++k) {
                        if (!(strcmp(qs.v[j].name, qs.v[k].name))) {
                                fprintf(stderr,
                                        "\u00ab%s\u00bb is not distinct in \u00ab%s\u00bb\n",
                                        qs.v[j].name, sarg);
                        }

                        if (!(strcmp(qt.v[j].name, qt.v[k].name))) {
                                fprintf(stderr,
                                        "\u00ab%s\u00bb is not distinct in \u00ab%s\u00bb\n",
                                        qs.v[j].name, sarg);
                        }
                }
        }

        /*
         * No immediate reason why there shouldn't be an isomorphism.
         * Let's try it.
         */
        ret = check_iso(&qs, &qt, fixnames, fidx);
done:
        quiver_clean(&qs);
        quiver_clean(&qt);

        if (fs) {
                fclose(fs);
        }

        if (ft) {
                fclose(ft);
        }

        free(fixnames);

        return ret;
}