#include "misc.h"
#include "sets.h"
#include "graph.h"
#include "queue.h"
#include <assert.h>
#include <string.h>

/*************************************************************************
**
**                            The Basics
**
*************************************************************************/

GRAPH *GraphAlloc(unsigned int n)
{
    static Boolean startup = 1;
    int i;
    GRAPH *G = Calloc(1, sizeof(GRAPH));
    if(startup)
    {
	startup = 0;
	SetStartup();
    }
    G->n = n;
    G->A = Calloc(n, sizeof(G->A[0]));
    for(i=0; i<n; i++)
	G->A[i] = SetAlloc(n);
    G->degree = Calloc(n, sizeof(G->degree[0]));
    G->neighbors = Calloc(n, sizeof(G->neighbors[0]));
    return G;
}

void GraphFree(GRAPH *G)
{
    int i;
    for(i=0; i<G->n; i++)
    {
	Free(G->neighbors[i]);
	SetFree(G->A[i]);
    }
    Free(G->neighbors);
    Free(G->degree);
    Free(G->A);
    Free(G);
}

/* If Gc == NULL, create duplicate.  Otherwise just copy G's info into Gc. */
GRAPH *GraphCopy(GRAPH *Gc, GRAPH *G)
{
    int i;
    if(Gc == G) return G;
    if(!Gc)
	Gc = GraphAlloc(G->n);

    if(G->n > Gc->n)
    {
	Gc->A = Realloc(Gc->A, G->n * sizeof(Gc->A[0]));
	Gc->degree = Realloc(Gc->degree, G->n * sizeof(Gc->degree[0]));
	Gc->neighbors = Realloc(Gc->neighbors, G->n * sizeof(Gc->neighbors[0]));
	/* reallaoc doesn't zero out new entries, damn it */
	for(i=Gc->n; i < G->n; i++)
	{
	    Gc->A[i] = NULL;
	    Gc->neighbors[i] = NULL;
	}
    }

    for(i=0; i < G->n; i++)
    {
	if(Gc->A[i] && Gc->n != G->n)
	{
	    SetFree(Gc->A[i]);
	    Gc->A[i] = NULL;
	}
	Gc->A[i] = SetCopy(Gc->A[i], G->A[i]);
	if(G->degree[i] > Gc->degree[i])
	    Gc->neighbors[i] = Realloc(Gc->neighbors[i], G->degree[i] * sizeof(Gc->neighbors[i][0]));
	memmove(Gc->neighbors[i], G->neighbors[i], G->degree[i] * sizeof(G->neighbors[i][0]));
	Gc->degree[i] = G->degree[i];
    }

    Gc->n = G->n;
    return Gc;
}

GRAPH *GraphConnect(GRAPH *G, int i, int j)
{   
    assert(0 <= i && i < G->n && 0 <= j && j < G->n);
    if(GraphAreConnected(G, i, j))
	return G;
    SetAdd(G->A[i], j);
    SetAdd(G->A[j], i);
    ++G->degree[i];
    ++G->degree[j];
    G->neighbors[i] = Realloc(G->neighbors[i], G->degree[i]*sizeof(G->neighbors[i][0]));
    G->neighbors[j] = Realloc(G->neighbors[j], G->degree[j]*sizeof(G->neighbors[j][0]));
    G->neighbors[i][G->degree[i]-1] = j;
    G->neighbors[j][G->degree[j]-1] = i;
    return G;
}

GRAPH *GraphEdgesAllDelete(GRAPH *G)
{
    int i;
    for(i=0; i < G->n; i++)
    {
	SetEmpty(G->A[i]);
	G->degree[i] = 0;
	/* Don't need to realloc/free neighbors, it'll happen automatically once we start re-adding edges */
    }
    return G;
}

GRAPH *GraphDisconnect(GRAPH *G, int i, int j)
{   
    int k;
    assert(0 <= i && i < G->n && 0 <= j && j < G->n);
    if(!GraphAreConnected(G, i, j))
	return G;
    SetDelete(G->A[i], j);
    SetDelete(G->A[j], i);
    --G->degree[i];
    --G->degree[j];

    /* now find and delete each other's neighbors */
    k=0;
    while(G->neighbors[i][k] != j)
	k++;
    assert(k <= G->degree[i] && G->neighbors[i][k] == j); /* this is the new degree, so using "<=" is correct */
    memmove(&(G->neighbors[i][k]), &(G->neighbors[i][k+1]), (G->degree[i]-k)*sizeof(G->neighbors[i][0]));
    
    k=0;
    while(G->neighbors[j][k] != i)
	k++;
    assert(k <= G->degree[j] && G->neighbors[j][k] == i);
    memmove(&(G->neighbors[j][k]), &(G->neighbors[j][k+1]), (G->degree[j]-k)*sizeof(G->neighbors[j][0]));

    return G;
}

#ifndef GraphAreConnected
Boolean GraphAreConnected(GRAPH *G, int i, int j)
{
    assert(0 <= i && i < G->n && 0 <= j && j < G->n);
    if(SetIn(G->A[i],j))
    {
	assert(SetIn(G->A[j],i));
	return true;
    }
    else
	return false;
}
#endif

void GraphPrintAdjMatrix(FILE *fp, GRAPH *G)
{
    int i, j;
    fprintf(fp, "%d\n", G->n);
    for(i=0; i<G->n; i++)
    {
	for(j=0; j<G->n - 1; j++)
	    fprintf(fp, "%d ", !!GraphAreConnected(G,i,j));
	fprintf(fp, "%d\n", !!GraphAreConnected(G,i,j));
    }
}


GRAPH *GraphReadAdjMatrix(FILE *fp)
{
    int i,j,n;
    GRAPH *G;
    if(fscanf(fp, "%d", &n) != 1)
	Fatal("GraphReadAdjMatrix: reading 'n' failed");
    assert(n >= 0);
    G = GraphAlloc(n);
    for(i=0; i<n; i++) for(j=0; j<n; j++)
    {
	int connected;
	if(fscanf(fp, "%d", &connected) != 1)
	    Fatal("GraphReadAdjMatrix: reading entry(%d,%d) failed", i, j);
	if(connected)
	    GraphConnect(G,i,j);
    }
    return G;
}


void GraphPrintAdjList(FILE *fp, GRAPH *G)
{
    int i, j;
    fprintf(fp, "%d\n", G->n);
    for(i=0; i<G->n; i++)
    {
	fprintf(fp, "%d ", G->degree[i]);
	for(j=0; j<G->degree[i] - 1; j++)
	    fprintf(fp, "%d ", G->neighbors[i][j]);
	fprintf(fp, "%d\n", G->neighbors[i][j]);
    }
}


GRAPH *GraphReadAdjList(FILE *fp)
{
    GRAPH *G;
    int n, i, j, d;
    if(fscanf(fp, "%d", &n) != 1)
	Fatal("GraphReadAdjList: failed to read 'n'");
    assert(n >= 0);
    G = GraphAlloc(n);
    for(i=0; i<n; i++)
    {
	if(fscanf(fp, "%d", &d) != 1)
	    Fatal("node %d: expecting degree, but couldn't find an integer", i);
	for(j=0; j<d; j++)
	{
	    int neigh;
	    if(fscanf(fp, "%d", &neigh) != 1)
		Fatal("node %d, degree %d, ran out of integers on neighbor #%d", i, d, j);
	    else
	    {
		assert(0 <= neigh && neigh < n);
		assert(neigh != i);
		GraphConnect(G, i, neigh);
	    }
	}
    }
    return G;
}

void GraphPrintConnections(FILE *fp, GRAPH *G)
{
    int i, j;
    fprintf(fp, "%d\n", G->n);
    for(i=0; i<G->n; i++) for(j=0; j<G->degree[i]; j++)
	fprintf(fp, "%d %d\n", i, G->neighbors[i][j]);
}


GRAPH *GraphReadConnections(FILE *fp)
{
    GRAPH *G;
    int n, i, j, d;
    if(fscanf(fp, "%d", &n) != 1)
	Fatal("GraphReadAdjList: failed to read 'n'");
    assert(n >= 0);
    G = GraphAlloc(n);
    
    while((d=fscanf(fp, "%d %d", &i, &j)) == 2)
    {
	if(i==-1 && j==-1)
	{
	    d=0;
	    break;
	}
	assert(0 <= i && i < G->n);
	assert(0 <= j && j < G->n);
	GraphConnect(G, i, j);
    }
    if(d > 0)
	Fatal("expecting no more integers, but got %d integers", d);
    return G;
}


GRAPH *GraphComplement(GRAPH *Gbar, GRAPH *G)
{
    int i, j;
    assert(Gbar != G);
    if(!Gbar)
	Gbar = GraphAlloc(G->n);

    assert(Gbar->n == G->n);

    for(i=0; i < G->n; i++) for(j=i+1; j < G->n; j++)
	if(!GraphAreConnected(G, i, j))
	    GraphConnect(Gbar, i, j);
    return Gbar;
}


GRAPH *GraphUnion(GRAPH *dest, GRAPH *G1, GRAPH *G2)
{
    int i, j, n = G1->n;

    if(G1->n != G2->n)
	return NULL;
    
    if(dest)
	dest->n = n;
    else
	dest = GraphAlloc(n);

    for(i=0; i < n; i++) for(j=i+1; j < n; j++)
	if(GraphAreConnected(G1, i, j) || GraphAreConnected(G2, i, j))
	    GraphConnect(dest, i ,j);

    return dest;
}


int GraphBFS(GRAPH *G, int root, int distance, int *nodeArray, int *distArray)
{
    QUEUE *BFSQ;
    int i, count = 0;

    assert(0 <= root && root < G->n);
    assert(distance >= 0);
    assert(nodeArray != NULL);
    assert(distArray != NULL);

    if(distance == 0) /* We could let the rest of the routine run, but why bother? */
    {
	nodeArray[0] = root;
	distArray[root] = 0;
	return 1;
    }

    for(i=0; i<G->n; i++)
	nodeArray[i] = distArray[i] = -1;

    distArray[root] = 0;
    BFSQ = QueueAlloc(G->n);
    QueuePut(BFSQ, (foint)root);
    while(QueueSize(BFSQ) > 0)
    {
	int v = QueueGet(BFSQ).i;

	/* At this point, distArray[v] should be assigned (when v was appended
	 * onto the queue), but v hasn't been "visited" or "counted" yet.
	 */

	assert(0 <= v && v < G->n);
	assert(0 <= distArray[v] && distArray[v] < G->n);

	assert(nodeArray[count] == -1);
	nodeArray[count] = v;
	count++;

	if(distArray[v] < distance) /* v's neighbors will be within BFS distance */
	{
	    int j;
	    for(j=0; j < G->degree[v]; j++)
		if(distArray[G->neighbors[v][j]] == -1) /* some of the neighbors might have already been visited */
		{
		    distArray[G->neighbors[v][j]] = distArray[v] + 1;
		    QueuePut(BFSQ, (foint)G->neighbors[v][j]);
		}
	}
    }
    QueueFree(BFSQ);
    return count;
}

/* doesn't allow Gv == G */
GRAPH *GraphInduced(GRAPH *Gv, GRAPH *G, SET *V)
{
    unsigned array[G->n], nV = SetToArray(array, V), i, j;
    assert(Gv != G);
    if(Gv)
    {
	assert(Gv->n == nV);
	GraphEdgesAllDelete(Gv);
    }
    else
	Gv = GraphAlloc(nV);
    for(i=0; i < nV; i++) for(j=i+1; j < nV; j++)
	if(GraphAreConnected(G, array[i], array[j]))
	    GraphConnect(Gv, i, j);
    return Gv;
}

/* allows Gv == G */
GRAPH *GraphInduced_NoVertexDelete(GRAPH *Gv, GRAPH *G, SET *V)
{
    unsigned array[G->n], nV = SetToArray(array, V), i, j;
    GRAPH *GGv = GraphAlloc(G->n);

    for(i=0; i < nV; i++) for(j=i+1; j < nV; j++)
	if(GraphAreConnected(G, array[i], array[j]))
	    GraphConnect(GGv, array[i], array[j]);
    GraphCopy(Gv, GGv);
    GraphFree(GGv);
    return Gv;
}

/*
** A reasonably fast search for triangles.
*/

/*
** This is a helper function for Contains Kn, but you can use it.  It
** tells you if adding this edge will cause a triangle.  Note that
** this function ONLY works if self-loops do not exist!
*/
Boolean GraphConnectingCausesK3(GRAPH *G, int i, int j)
{
    int numIntersect;
    SET *C = SetAlloc(G->n);
    SetIntersect(C, G->A[i], G->A[j]);
    numIntersect = SetCardinality(C);
    SetFree(C);
    return numIntersect != 0;
}

Boolean GraphContainsK3(GRAPH *G)
{
    int i,j;
    for(i=0; i < G->n; i++)
	for(j=i+1; j < G->n; j++)
	    if(GraphAreConnected(G,i,j) && GraphConnectingCausesK3(G,i,j))
		return true;
    return false;
}



/**************************************************************************
**
** These are the Clique and Indep set (exponential time) algorithms.
** They work for general graphs.
**
**************************************************************************/

/*
** We check to see if these particular k vertices are a clique by
** bitwise-anding together their adjacency lists (with self-loop added
** to each).  Return false of it's not a clique, otherwise return the set
** of nodes.
** Note that if the graph has changed since the previous call, we may miss some cliques.
*/
SET *Graph_IsCombClique(CLIQUE *c)
{
    SET *intersect;
    int i;

    SetEmpty(c->set);
    for(i=0; i < c->cliqueSize; i++)
	SetAdd(c->set, c->inducedArray[c->combArray[i]]);

    intersect = SetCopy(NULL, c->set);

    for(i=0; i < c->cliqueSize; i++)
    {
	int node = c->inducedArray[c->combArray[i]];
	/* Temporarily add in self-loop for intersection purposes */
	SetAdd(c->G->A[node], node);
	SetIntersect(intersect, intersect, c->G->A[node]);
	SetDelete(c->G->A[node], node);
	if(!SetEq(intersect, c->set))
	{
	    SetFree(intersect);
	    return NULL;
	}
    }
    assert(SetEq(intersect, c->set));
    SetFree(intersect);
    return c->set;
}

SET *GraphKnNext(CLIQUE *c)
{
    SET *s;
    while(CombinNext(c->combin))
	if((s = Graph_IsCombClique(c)))
	    return s;
    return NULL;
}

CLIQUE *GraphKnFirst(GRAPH *G, int k)
{
    int i, nDegk = G->n;
    SET *setDegk;
    CLIQUE *c;

    assert(k <= G->n);
    if(k == 0)
	return NULL;
    
    setDegk = SetAlloc(G->n);
    c = (CLIQUE*)Calloc(1,sizeof(CLIQUE));
    c->G = GraphCopy(NULL, G);
    c->cliqueSize = k;
    c->inducedArray = Calloc(G->n, sizeof(c->inducedArray[0]));

#if 1
    /*
    ** First reduce the potential number of vertices needed to check. A
    ** necessary condition for a vertex to be a member of a Kk is to have
    ** degree >= k-1, and all edges go to other vertices of degree >= k-1.
    ** Iterate inducing subgraphs 'til we cain't induce no more...
    */
    while(nDegk >= k)
    {
	int prevNDegk = nDegk;  /* just to check for changes */
	nDegk = 0;
	SetEmpty(setDegk);
	for(i=0; i < c->G->n; i++)
	    if(c->G->degree[i] >= k-1)
	    {
		++nDegk;
		SetAdd(setDegk, i);
	    }
	if(nDegk == prevNDegk)  /* nobody eliminated */
	    break;
	GraphInduced_NoVertexDelete(c->G, c->G, setDegk);
    }

    if(nDegk < k)
    {
	Free(c->inducedArray);
	GraphFree(c->G);
	Free(c);
	SetFree(setDegk);
	return NULL;
    }
    i = SetToArray(c->inducedArray, setDegk);
    assert(i == nDegk);
#else
    nDegk = G->n;
    for(i=0; i < nDegk; i++)
	c->inducedArray[i] = i;
#endif
    SetEmpty(setDegk);
    c->set = setDegk; /* re-use setDegK */
    c->combArray = Calloc(k, sizeof(c->combArray[0]));
    c->combin = CombinZeroth(nDegk, k, c->combArray);
    if(Graph_IsCombClique(c) || GraphKnNext(c))
	return c;
    /* else */
    GraphCliqueFree(c);
    return NULL;
}

void GraphCliqueFree(CLIQUE *c)
{
    CombinFree(c->combin);
    Free(c->combArray);
    Free(c->inducedArray);
    GraphFree(c->G);
    SetFree(c->set);
    Free(c);
}

CLIQUE *GraphInFirst(GRAPH *G, int n)
{
    static GRAPH Gbar;  /* non re-entrant */
    GraphComplement(&Gbar, G);
    return GraphKnFirst(&Gbar, n);
}


Boolean GraphKnContains(GRAPH *G, int n)
{
    CLIQUE *c;
    assert(n>=0);
    if(n < 2)
	return true;
    if(n > G->n)
	return false;
    c = GraphKnFirst(G, n);
    if(c)
	GraphCliqueFree(c);
    return c != NULL;
}

Boolean GraphInContains(GRAPH *G, int n)
{
    static GRAPH Gbar;  /* non re-entrant */
    GraphComplement(&Gbar, G);
    return GraphKnContains(&Gbar, n);
}

/**************************************************************************
**
**  Graph Isomorphism
**
**************************************************************************/

static GRAPH *isoG1, *isoG2;

static Boolean _permutationIdentical(int n, int perm[n])
{
    int i, j;
    for(i=0; i<n; i++)
	if(isoG1->degree[i] != isoG2->degree[perm[i]])
	    return false;

    for(i=0; i<n; i++) for(j=i+1; j<n; j++)
	/* The !GraphAreConnected is just to turn a bitstring into a boolean */
	if(!GraphAreConnected(isoG1, i,j) !=
	    !GraphAreConnected(isoG2, perm[i], perm[j]))
	    return false;   /* non-isomorphic */
    return true;   /* isomorphic! */
}

Boolean GraphsIsomorphic(GRAPH *G1, GRAPH *G2)
{
    int i, n = G1->n, degreeCount1[n], degreeCount2[n];
    SET *degreeOnce;

    /*
    ** First some simple tests.
    */
    if(G1->n != G2->n)
	return false;
    
    if(n < 2)
	return true;

    /*
    ** Ensure each degree occurs the same number of times in each.
    */
    for(i=0; i<n; i++)
	degreeCount1[i] = degreeCount2[i] = 0;
    for(i=0; i<n; i++)
    {
	++degreeCount1[G1->degree[i]];
	++degreeCount2[G2->degree[i]];
    }
    for(i=0; i<n; i++)
	if(degreeCount1[i] != degreeCount2[i])
	    return false;

    /*
    ** Let degree d appear only once.  Then there is exactly one vertex
    ** v1 in G1 with degree d, and exactly one vertex v2 in G2 with degree d.
    ** G1 and G2 are isomorphic only if the neighborhood of v1 is isomorphic
    ** to the neighborhood of v2.
    */
    degreeOnce = SetAlloc(n);
    for(i=0; i<n; i++)
	if(degreeCount1[i] == 1)
	    SetAdd(degreeOnce, i);
    for(i=0; i<n; i++)
    {
	/* Find out if the degree of vertex i in G1 appears only once */
	if(SetIn(degreeOnce, G1->degree[i]))
	{
	    int j, degree = G1->degree[i];
	    GRAPH *neighG1i, *neighG2j;

	    /* find the (unique) vertex in G2 that has the same degree */
	    for(j=0; j < n; j++)
		if(G2->degree[j] == degree)
		    break;
	    assert(j < n);

	    neighG1i = GraphInduced(NULL, G1, G1->A[i]);
	    neighG2j = GraphInduced(NULL, G2, G2->A[j]);

	    /*
	    ** Note: this recursion works only as long as
	    ** _permutationIdentical doesn't call GraphsIsomorphic.
	    ** (if it does, isoG1 and isoG2 get messed up).
	    */

	    j = GraphsIsomorphic(neighG1i, neighG2j);
	    GraphFree(neighG1i);
	    GraphFree(neighG2j);
	    if(!j)
		return false;
	    /* Otherwise they *might* be isomorphic, so keep going */
	}
    }
    SetFree(degreeOnce);

    /*
    ** Oh well, fire up the exponential search.
    ** CombinAllPermutations will return 0 iff all permutations were
    ** tried; the function _permutationIdentical should return non-zero
    ** when it finds an identical permutation, and that non-zero value
    ** will be returned here, indicating an identical permutation was
    ** found, ie, that the graphs are isomorphic.
    */
    isoG1 = G1; isoG2 = G2;
    return !!CombinAllPermutations(n, _permutationIdentical);
}
