// topologies.h
// Copyright (C) 2023 Mark Bowron

// Various functions dealing with topologies.

// 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.

// Sets are represented by unsigned 32-bit integers with leading zeros.

void    BaseToClosures(int **base, int **closures, int spacesize);
// Computes and stores the closures of singletons given a base for the topology.

void    MatrixToPointClosures(int **closures, uint32_t *point_closures, int spacesize);
// Converts each row of a 0-1 matrix of int variables (these rows usually represent the
// closures of singletons) into a uint32_t variable.

void    BaseToPointClosures(int **base, int **closures, uint32_t *point_closures, int spacesize);
// Converts a 0-1 base matrix into a uint32_t array of point closures.

void    PointClosuresToMatrix(uint32_t *point_closures, int **closures, int spacesize);
// This function is the inverse of MatrixToPointClosures().

int     KFGOrdering(KZ *kz, KFG *fam, int spacesize, uint32_t *point_closures);
// Determines which of the 12 KFG orderings a space satisfies.

// #    index to KFGOrdering()
// -    ----------------------
// 1    D
// 2    P
// 3    EO 1 (fb <= id not true)
// 4    ED 1 (fb <= id not true)
// 5    OU 1 (fbg <= id not true)
// 6    GE 1 (bg <= aibi not true)
// 7    KD 1 (fbg <= id not true)
// 8    EO 2
// 9    ED 2
// 10   OU 2
// 11   GE 2
// 12   KD 2

int     ClosuresToKFGOrdering(KZ *kz, KFG *kfg, int **closures, uint32_t *point_closures, int spacesize);
// Returns the KFG ordering of a space represented by a 0-1 matrix of point closures.

int     BaseToKFGOrdering(KZ *kz, KFG *kfg, int **base, int **closures, uint32_t *point_closures, int spacesize);
// Returns the KFG ordering of a space represented by a 0-1 base matrix.

void    BaseToPsiNumbers(KFZ *kfz, int **base, int **closures, uint32_t *point_closures, int spacesize);
// Computes the psi-number of every subset of a space represented by a 0-1 base matrix
// then prints every psi-number that occurred.

int     KFZtoKNum(KFZ *kfz, int spacesize);
// Computes the k-number of a space represented by a 0-1 base matrix.

int     KFZtoKFNum(KFZ *kfz, int spacesize);
// Computes the k_f-number of a space represented by a 0-1 base matrix.

int     AllPhiNumbers(KZ *kz, int **base, int **closures, uint32_t *point_closures, int spacesize);
// Computes the phi-number of every subset of a space represented by a 0-1 base matrix
// then verifies that every phi-number occurred.

void    Warshall(int **adj, int dim);
// Computes the transitive closure of a preorder represented by and adjacency matrix.

int     matrixEqual(int **a, int **b, int dim);
int     matrixSubset(int **a, int **b, int dim);
void    matrixIntersect(int **a, int **b, int **intersection, int dim);
void    matrixCopy(int **from, int **to, int n);
void    matrixErase(int **matrix, int dim);
void    matrixMultiply(int **a, int **b, int **result, int dim);
void    getTransitiveReduction(int **adj, int **copy, int **tc, int **ab, int **tr, int dim);
// Computes the transitive reduction of a partial order represented by an adjacency
// matrix.

// Note: swapping the terms "open" and "closed" results in dual topologies on the
// underlying set when it is finite.  The BaseToClosures() function therefore works
// equally well as a "ClosuresToBase()" function that takes the closures of singletons
// as input (in the base[][] matrix) and outputs (in the closures[][] matrix) their
// smallest open neighborhoods (under the same topology).

void BaseToClosures(int **base, int **closures, int spacesize)
{
    int j, k, m;

// We initially set the closure of each singleton to the whole space.

    for (j=0; j<spacesize; j++)
        for (k=0; k<spacesize; k++)
            closures[j][k] = 1;

// Now we remove points in base sets not containing the singleton.

    for (j=0; j<spacesize; j++)
        for (k=0; k<spacesize; k++)
        {

// base[k] is the smallest open neighborhood of k.  If it contains j we continue.

            if (base[k][j] == 1) continue;

// Since base[k] is open and does not contain j, it is disjoint from the closure of j.
// Thus every point in base[k] needs to be removed from closures[j].

            for (m=0; m<spacesize; m++)
            {
                if (base[k][m]) closures[j][m] = 0;
            }
        }
}


// This function converts each row of a 0-1 matrix of int variables (these rows usually
// represent the closures of singletons) into a uint32_t variable.

void MatrixToPointClosures(int **closures, uint32_t *point_closures, int spacesize)
{
    int j, k;

// Convert each row of the closures matrix to uint32_t.

    for (j=0; j<spacesize; j++)
    {
        point_closures[j] = 0;
        for (k=0; k<spacesize; k++)
        {
            if (closures[j][k]) point_closures[j] = point_closures[j] + (one<<(uint32_t)k);
        }
    }
}


// This function converts a 0-1 base matrix into a uint32_t array of point closures.

void BaseToPointClosures(int **base, int **closures, uint32_t *point_closures, int spacesize)
{
    BaseToClosures(base, closures, spacesize);
    MatrixToPointClosures(closures, point_closures, spacesize);
}


// This function is the inverse of MatrixToPointClosures().

void PointClosuresToMatrix(uint32_t *point_closures, int **closures, int spacesize)
{
    int j, k;

// Convert point_closures array to adjacency matrix.

    for (j=0; j<spacesize; j++)
        for (k=0; k<spacesize; k++)
        {
            if (point_closures[j] & (one<<(uint32_t)k)) closures[j][k] = 1;
            else closures[j][k] = 0;
        }
}


// This function determines which of the 12 KFG orderings a space satisfies.
// Never feed it a blank fam variable; it only works if fam has been computed.

int KFGOrdering(KZ *kz, KFG *fam, int spacesize, uint32_t *point_closures)
{
    int  ordering;

// To obtain the GE monoid we extract kz from fam.

    CopyKZfromKFG(kz,fam,spacesize);
    ordering = GE(kz,spacesize,point_closures);

// Five of the seven GE monoids require further testing to determine whether or not a
// certain inequality holds.

    if (ordering == 3 || ordering == 4) // X is EO or ED
    {
        if (OpSubset(fam->_fb,fam->_id,spacesize)) ordering = ordering + 5;
    }
    if (ordering == 5 || ordering == 7) // X is OU or KD
    {
        if (OpSubset(fam->_fbg,fam->_id,spacesize)) ordering = ordering + 5;
    }
    if (ordering == 6)          // X is GE
    {
        if (OpSubset(fam->_bg,fam->_aibi,spacesize)) ordering = ordering + 5;
    }
    return ordering;
}


// This function returns the KFG ordering of a space represented by a 0-1 matrix of point closures.

int ClosuresToKFGOrdering(KZ *kz, KFG *kfg, int **closures, uint32_t *point_closures, int spacesize)
{
    MatrixToPointClosures(closures,point_closures,spacesize);
    ComputeKFG(kfg,spacesize,point_closures);
    return KFGOrdering(kz,kfg,spacesize,point_closures);
}


// This function returns the KFG ordering of a space represented by a 0-1 base matrix.

int BaseToKFGOrdering(KZ *kz, KFG *kfg, int **base, int **closures, uint32_t *point_closures, int spacesize)
{
    BaseToClosures(base,closures,spacesize);
    return ClosuresToKFGOrdering(kz,kfg,closures,point_closures,spacesize);
}


// This function computes the psi-number of every subset of a space represented by a 0-1
// base matrix then prints every psi-number that occurred.

void BaseToPsiNumbers(KFZ *kfz, int **base, int **closures, uint32_t *point_closures, int spacesize)
{
    uint32_t    set;
    int         psi, psi_number[71] = {0};

    BaseToClosures(base,closures,spacesize);
    MatrixToPointClosures(closures,point_closures,spacesize);
    ComputeKFZ(kfz,spacesize,point_closures);
    for (set=1; set<(one<<spacesize); set++)
    {
        psi_number[Psi(set,spacesize,kfz)] = 1;
    }
    for (psi=1; psi<69; psi++)
    {
        if (psi_number[psi]) printf("%d ",psi);
    }
}


// This function computes the psi-number of every subset of a space represented by a 0-1
// base matrix then returns the largest k-number that occurred.

int KFZtoKNum(KFZ *kfz, int spacesize)
{
    uint32_t    set;
    int         knum = 0, psi_number[71] = {0}, psi;

    for (set=1; set<(one<<spacesize); set++)
    {
        psi_number[Psi(set,spacesize,kfz)] = 1;
    }
    for (psi=68; psi>0; psi--)
    {
        if (psi_number[psi]) knum = k_num[psi];
    }
    if (knum) return knum;
    return 2;
}


// This function computes the psi-number of every subset of a space represented by a 0-1
// base matrix then returns the largest k_f-number that occurred.

int KFZtoKFNum(KFZ *kfz, int spacesize)
{
    uint32_t    set;
    int         kfnum = 0, psi_number[71] = {0}, psi;

    for (set=1; set<(one<<spacesize); set++)
    {
        psi_number[Psi(set,spacesize,kfz)] = 1;
    }
    for (psi=0; psi<70; psi++)
    {
        if (psi_number[psi])
        {
            if (kfnum < kf_num[psi]) kfnum = kf_num[psi];
        }
    }
    if (kfnum) return kfnum;
    return 2;
}


// This function computes the phi-number of every subset of a space represented by a 0-1
// base matrix then verifies that every phi-number occurred.

int AllPhiNumbers(KZ *kz, int **base, int **closures, uint32_t *point_closures, int spacesize)
{
    uint32_t    set;
    int         phi, phi_number[31] = {0};

    BaseToClosures(base,closures,spacesize);
    MatrixToPointClosures(closures,point_closures,spacesize);
    ComputeKZ(kz,spacesize,point_closures);
    for (set=1; set<(one<<spacesize); set++)
    {
        phi_number[Phi(set,kz)] = 1;
    }
    for (phi=1; phi<30; phi++)
    {
        if (phi_number[phi]) continue;
        else return 0;
    }
    return 1;
}


// This function computes the transitive closure of a preorder represented by an
// adjacency matrix.

void Warshall(int **adj, int dim)
{
    int k, r, c;

    /*
     *  Apply Warshall's algorithm to adj[][]
     */

    for (k=0; k<dim; k++) adj[k][k] = 1;    //   force reflexivity

    for (k=0; k<dim; k++)
        for (r=0; r<dim; r++)
            for (c=0; c<dim; c++)
                adj[r][c] = adj[r][c] || (adj[r][k] && adj[k][c]);
}


int matrixEqual(int **a, int **b, int dim)
{
    int i, j;

    for (i=0; i<dim; i++)
        for (j=0; j<dim; j++)
            if (a[i][j]!=b[i][j]) return 0;
    return 1;
}


int matrixSubset(int **a, int **b, int dim)
{
    int i, j;

    for (i=0; i<dim; i++)
        for (j=0; j<dim; j++)
            if (a[i][j]>b[i][j]) return 0;
    return 1;
}


void matrixIntersect(int **a, int **b, int **intersection, int dim)
{
    int i, j;

    for (i=0; i<dim; i++)
        for (j=0; j<dim; j++)
        {
            if (a[i][j]*b[i][j]) intersection[i][j] = 1;
            else intersection[i][j] = 0;
        }
}


void matrixCopy(int **from, int **to, int dim)
{
    int r, c;

    for (r=0; r<dim; r++)
        for (c=0; c<dim; c++)
            to[r][c] = from[r][c];
}


void matrixErase(int **matrix, int dim)
{
    int i, j;
    for (i=0; i<dim; i++)
        for (j=0; j<dim; j++)
            matrix[i][j] = 0;
}


void matrixMultiply(int **a, int **b, int **result, int dim)
{
    int sum, i, j, k;

    for (i=0; i<dim; i++)
        for (j=0; j<dim; j++)
        {
            sum = 0;
            for (k=0; k<dim; k++)
                sum = sum + a[i][k] * b[k][j];
            result[i][j] = sum;
        }
}


// This function computes the transitive reduction of a partial order represented by an
// adjacency matrix.

void getTransitiveReduction(int **adj, int **copy, int **tc, int **ab, int **tr, int dim)
{
    int i, j;

    matrixCopy(adj, copy, dim);
    matrixCopy(adj, tc, dim);

// Apply Warshall's algorithm to get the transitive closure tc of adj

    Warshall(tc, dim);

// Transitive reduction algorithm requires that the relation be anti-reflexive

    for (i=0; i<dim; i++)
    {
        copy[i][i] = 0;
        tc[i][i] = 0;
    }

// Multiply copy by tc to get matrix product ab

    matrixMultiply(copy, tc, ab, dim);

// We can now calculate the transitive reduction of adj:
// tr[u][v] = 1 iff copy[u][v] = 1 and ab[u][v] = 0

    for (i=0; i<dim; i++)
        for (j=0; j<dim; j++)
        {
            if (copy[i][j] == 1 && ab[i][j] == 0) tr[i][j] = 1;
            else tr[i][j] = 0;
        }
}