// main.c
// Copyright (C) 2023 Mark Bowron

// Verify spaces are topological and pairwise nonhomeomorphic.

// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.

// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.


#include    <stdio.h>
#include    <stdlib.h>
#include    <inttypes.h>
#include    <stdbool.h>
#include    <string.h>
#include    "set_operators.h"
#include    "set_families.h"
#include    "set_relations.h"
#include    "collapses.h"
#include    "topologies.h"


void        swap(int i, int j, int *elts);
void        permute(FILE *input_a, FILE *input_b, int spacesize);


int main(int argc, const char * argv[])
{
    FILE       *input_a = NULL, *input_b = NULL;
    int         j, k, space = 0, spacesize = 6;
    int       **adj, **tc;
    uint32_t   *point_closures;

    adj = (int **)malloc(spacesize * sizeof(int *));
    if (adj==0)
    {
        printf("\nCould not allocate space for adj.\n");
        exit(0);
    }
    tc = (int **)malloc(spacesize * sizeof(int *));
    if (tc==0)
    {
        printf("\nCould not allocate space for tc.\n");
        exit(0);
    }
    for (j=0; j<spacesize; j++)
    {
        adj[j] = (int *)malloc(spacesize * sizeof(int));
        if (adj[j]==0)
        {
            printf("\nCould not allocate space for adj[%d].\n", j);
            exit(0);
        }
        tc[j] = (int *)malloc(spacesize * sizeof(int));
        if (tc[j]==0)
        {
            printf("\nCould not allocate space for tc[%d].\n", j);
            exit(0);
        }
        for (k=0; k<spacesize; k++)
        {
            adj[j][k] = 0;
            tc[j][k] = 0;
        }
    }

// Open the input file.

    switch (spacesize)
    {
        case 3:
            input_a = fopen("3-spaces.txt", "r");
            if (input_a==0)
            {
                printf("\nCould not open file 3-spaces.txt for input_a.\n");
                exit(0);
            }
            input_b = fopen("3-spaces.txt", "r");
            if (input_b==0)
            {
                printf("\nCould not open file 3-spaces.txt for input_b.\n");
                exit(0);
            }
            break;
        case 4:
            input_a = fopen("4-spaces.txt", "r");
            if (input_a==0)
            {
                printf("\nCould not open file 4-spaces.txt for input_a.\n");
                exit(0);
            }
            input_b = fopen("4-spaces.txt", "r");
            if (input_b==0)
            {
                printf("\nCould not open file 4-spaces.txt for input_b.\n");
                exit(0);
            }
            break;
        case 5:
            input_a = fopen("5-spaces.txt", "r");
            if (input_a==0)
            {
                printf("\nCould not open file 5-spaces.txt for input_a.\n");
                exit(0);
            }
            input_b = fopen("5-spaces.txt", "r");
            if (input_b==0)
            {
                printf("\nCould not open file 5-spaces.txt for input_b.\n");
                exit(0);
            }
            break;
        case 6:
            input_a = fopen("6-spaces.txt", "r");
            if (input_a==0)
            {
                printf("\nCould not open file 6-spaces.txt for input_a.\n");
                exit(0);
            }
            input_b = fopen("6-spaces.txt", "r");
            if (input_b==0)
            {
                printf("\nCould not open file 6-spaces.txt for input_b.\n");
                exit(0);
            }
            break;
        default:
            printf("\nError: input file not found.\n\n");
            exit(0);
    }

// Sets are represented by unsigned 32-bit binary integers with leading
// zeros.  Assuming spacesize = n we could say that X = {0,1,...,n-1} and
// that an element k is in A iff the coefficient of 2^k in the binary integer
// representing A is 1.

    point_closures = (uint32_t *)malloc(spacesize * sizeof(uint32_t));
    if (point_closures==0)
    {
        printf("\nCould not allocate space for point_closures.\n");
        exit(0);
    }

    while (space < spaces[spacesize])
    {

// This program interprets each line of input as the closures of singletons.
// We could alternatively interpret each line as the smallest open neighborhoods
// as this still yields all nonhomeomorphic 6-spaces, but these two (dual)
// interpretations do not in general yield the same space for a given line of
// input and thus call for different algorithms to analyze a given space.

// We now read the current line of input into the point_closures array.

        for (j=0; j<spacesize; j++) fscanf(input_a, "%" PRIu32, &point_closures[j]);

// To verify we have a topological space we convert the point_closures array to
// matrix form, apply Warshall's algorithm to a copy of it, then compare the two
// to verify transitivity (we have a topological space iff they are equal).

        PointClosuresToMatrix(point_closures,adj,spacesize);
        matrixCopy(adj,tc,spacesize);
        Warshall(tc,spacesize);
        if (!matrixEqual(adj,tc,spacesize))
        {
            printf("Error: adj[][] does not represent a transitive relation.\n\n");
            exit(0);
        }
        space++;
    }

// If the program is still running then we know that every line in the input file
// represents a topological space.

    printf("\nAll %d lines represent topological spaces.\n", space);

// Since the number of nonhomeomorphic topological spaces on n points is currently known
// for n < 17, to show that our input file contains all topological spaces on n points
// (up to homeomorphism), it suffices to show that all pairs in it are nonhomeomorphic.
// Two finite spaces represented by adjacency matrices A and B are homeomorphic iff each
// can be obtained from the other by relabeling points iff B = PAP^{-1} for some
// permutation matrix P of dimension spacesize. Thus we need to show that B != PAP^T
// for all such permutation matrices P and all distinct A and B from the input file.
// We do this inside the permute() function.

    permute(input_a, input_b, spacesize);

// If the program is still running then we know that all of the spaces in the input file
// are pairwise nonhomeomorphic.

    printf("\nAll %d spaces are pairwise nonhomeomorphic.\n\n", space);

    fclose(input_a);
    fclose(input_b);
    for (j=0; j<spacesize; j++)
    {
        free(adj[j]);
        free(tc[j]);
    }
    free(adj);
    free(tc);
    free(point_closures);
    return EXIT_SUCCESS;
}


void swap(int i, int j, int *elts)
{
    int tmp = elts[i];
    elts[i] = elts[j];
    elts[j] = tmp;
}


// The permute() function uses Heap's algorithm to check that B != PAP^T
// for every permutation matrix P of dimension spacesize and every pair of
// adjacency matrices A, B obtained from separate lines in the input file.

void permute(FILE *input_a, FILE *input_b, int spacesize)
{
    int         j, k, upper, lower;
    int        *wgt, *elts;
    int       **a, **b, **p, **pt, **papt, space_a, space_b;
    uint32_t   *closures;

    a = (int **)malloc(spacesize * sizeof(int *));
    if (a==0)
    {
        printf("\nCould not allocate space for a.\n");
        exit(0);
    }
    b = (int **)malloc(spacesize * sizeof(int *));
    if (b==0)
    {
        printf("\nCould not allocate space for b.\n");
        exit(0);
    }
    p = (int **)malloc(spacesize * sizeof(int *));
    if (p==0)
    {
        printf("\nCould not allocate space for p.\n");
        exit(0);
    }
    pt = (int **)malloc(spacesize * sizeof(int *));
    if (pt==0)
    {
        printf("\nCould not allocate space for pt.\n");
        exit(0);
    }
    papt = (int **)malloc(spacesize * sizeof(int *));
    if (papt==0)
    {
        printf("\nCould not allocate space for papt.\n");
        exit(0);
    }
    for (j=0; j<spacesize; j++)
    {
        a[j] = (int *)malloc(spacesize * sizeof(int));
        if (a[j]==0)
        {
            printf("\nCould not allocate space for a[%d].\n", j);
            exit(0);
        }
        b[j] = (int *)malloc(spacesize * sizeof(int));
        if (b[j]==0)
        {
            printf("\nCould not allocate space for b[%d].\n", j);
            exit(0);
        }
        p[j] = (int *)malloc(spacesize * sizeof(int));
        if (p[j]==0)
        {
            printf("\nCould not allocate space for p[%d].\n", j);
            exit(0);
        }
        pt[j] = (int *)malloc(spacesize * sizeof(int));
        if (pt[j]==0)
        {
            printf("\nCould not allocate space for pt[%d].\n", j);
            exit(0);
        }
        papt[j] = (int *)malloc(spacesize * sizeof(int));
        if (papt[j]==0)
        {
            printf("\nCould not allocate space for papt[%d].\n", j);
            exit(0);
        }
        for (k=0; k<spacesize; k++)
        {
            a[j][k] = 0;
            b[j][k] = 0;
            p[j][k] = 0;
            pt[j][k] = 0;
            papt[j][k] = 0;
        }
    }

    closures = (uint32_t *)malloc(spacesize * sizeof(uint32_t));
    if (closures==0)
    {
        printf("\nCould not allocate space for closures.\n");
        exit(0);
    }

    elts = (int *)malloc(spacesize * sizeof(int));
    if (elts==0)
    {
        printf("\nCould not allocate space for elts.\n");
        exit(0);
    }

    wgt = (int *)malloc(spacesize * sizeof(int));
    if (wgt==0)
    {
        printf("\nCould not allocate space for wgt.\n");
        exit(0);
    }
    
    for (upper=0; upper<spacesize; upper++)
    {
        elts[upper] = upper;
        wgt[upper] = 0;
    }

    upper = 1;                  // first swap points are 1 and 0
    while(upper < spacesize)
    {
        if (wgt[upper] < upper)
        {

// if upper odd then lower = wgt[upper] otherwise lower = 0

            lower = upper % 2 * wgt[upper];
            swap(upper, lower, elts);

    /*
     *  The permutation stored in elts[] is ready to use to create P and P^T.
     */

            matrixErase(p,spacesize);
            matrixErase(pt,spacesize);
            for (j=0; j<spacesize; j++)
            {
                p[j][elts[j]] = 1;
                pt[elts[j]][j] = 1;
            }

// Move pointer input_a to top of input file.

            rewind(input_a);

// The following for loop advances input_a to the second line.

            for (j=0; j<spacesize; j++) fscanf(input_a, "%" PRIu32, &closures[j]);
            space_a = 1;

// We check that every B above A in the input file satisfies B != PAP^T.

            while (space_a < spaces[spacesize])
            {
                for (j=0; j<spacesize; j++) fscanf(input_a, "%" PRIu32, &closures[j]);
                PointClosuresToMatrix(closures,a,spacesize);

// B comes in handy to compute PAP^T.

                matrixMultiply(p,a,b,spacesize);
                matrixMultiply(b,pt,papt,spacesize);

// Move pointer input_b to top of input file.

                rewind(input_b);
                space_b = 0;

// We now use B for its main purpose, namely, to check that B != PAP^T for all
// spaces B above A in the input file.

                while (space_b < space_a)
                {
                    for (j=0; j<spacesize; j++) fscanf(input_b, "%" PRIu32, &closures[j]);
                    PointClosuresToMatrix(closures,b,spacesize);

                    if (matrixEqual(b,papt,spacesize))
                    {
                        printf("Error: homeomorphic A and B found.\n\n");
                        exit(0);
                    }
                    space_b++;
                }
                space_a++;
            }

    /*
     *  We are now finished using the permutation stored in elts[].
     */

            wgt[upper]++;       // increase index "weight" for upper by one
            upper = 1;
        }
        else                    // wgt[upper] == upper
        {
            wgt[upper] = 0;
            upper++;
        }
    }

    for (j=0; j<spacesize; j++)
    {
        free(a[j]);
        free(b[j]);
        free(p[j]);
        free(pt[j]);
        free(papt[j]);
    }
    free(a);
    free(b);
    free(p);
    free(pt);
    free(papt);
    free(closures);
    free(elts);
    free(wgt);
}