#include <stdio.h>
#include <string.h>
#include "Graph.h"

main(int argc,char *argv[])
{
  MLINK link;			/* The MathLink connection */
  char *func;			/* Function to be performed */
  long nargs;			/* no. of arguments to func */
  int done=0;			/* Are we done? */
  
  fprintf(stderr,"Graph trying to link...\n");
  link=MLOpen(argc,argv);
  if (link == NULL) return 1;
  MLPutString(link,"-- linked to Graph --");
  MLFlush(link);

  while (!done) {
    MLGetFunction(link, &func, &nargs);
    if (0==strcmp(func,"MLCanonicalForm")) MLCanonicalForm(link);
    // else if (0==strcmp(func,"MLPhi")) MLPhi(link);
    else if (0==strcmp(func,"MLSplitGraph")) MLSplitGraph(link);
    else if (0==strcmp(func,"MLExit")) done=1;
    else {
      MLNewPacket(link);
      MLPutSymbol(link,"UnknownOperation");
      MLFlush(link);
    }
    MLDisownSymbol(link, func);
  }
}

Graph::~Graph()
{
  int i;

  if (N>0) {
    for (i=0; i<N; ++i) delete v[i];
    delete vd;
    delete vt;
    delete vp;
    delete vai;
    delete v;
  }
}

int Graph::PermuteSign(permutation *perm)
{
  int sign=1;
  static permutation invperm[MAXN];     /* The inverse of perm */
  static int oddvs[MAXN];		/* list of odd vertices */
  int l=0;				/* The length of oddvs */
  int i;

  for (i=0; i<n; ++i) invperm[perm[i]]=i;
  for (i=0; i<n; ++i) {
    if (vt[i]==Arrow || vt[i]==OV) sign*=Signature(vd[i],v[i],perm);
    if (vp[i]) oddvs[l++]=i;
  }
  sign*=Signature(l,oddvs,perm);
  return sign;
}

MLINK operator>>(MLINK link, Graph& G)
{
  char *vname;		/* name of current vertex */
  char *vrn;		/* reduced name of current vertex */
  int markfound=0;	/* Whether a Marker was found yet */
  long int longint;	/* Just for communication with MathLink */
  int i,j,k;

  G.coeff=1;
  MLGetFunction(link, &vname, &longint);
  MLDisownSymbol(link, vname);
  G.N=(int) longint;
  G.vd=new int[G.N];
  G.vt=new VertexType[G.N];
  G.vp=new int[G.N];
  G.vai=new int[G.N];
  G.v=new int*[G.N];
  G.n=0;
  for (i=0; i<G.N; ++i) {
    MLGetFunction(link, &vname, &longint); G.vd[i]=(int) longint;
    if (markfound) G.vt[i]=Marker; else {
      ++G.n;
      if (vname[0]=='O' && vname[1]=='d' && vname[2]=='d') {
        G.vp[i]=1;
        vrn=vname+3;
      } else {
        G.vp[i]=0;
        vrn=vname;
      }
      if (0==strcmp(vrn,"Marker")) {
        --G.n;
        G.vt[i]=Marker;
        markfound=1;
      } else if (0==strcmp(vrn,"UOV")) G.vt[i]=UOV;
      else if (0==strcmp(vrn,"OV")) G.vt[i]=OV;
      else if (0==strcmp(vrn,"Dashed")) G.vt[i]=Dashed;
      else if (0==strcmp(vrn,"Arrow")) G.vt[i]=Arrow;
      else if (0==strcmp(vrn,"UV")) G.vt[i]=UV;
      else if (0==strcmp(vrn,"Blue")) G.vt[i]=Blue;
      else if (0==strcmp(vrn,"Green")) G.vt[i]=Green;
      else if (0==strcmp(vrn,"Red")) G.vt[i]=Red;
      else if (0==strcmp(vrn,"X")) G.vt[i]=X;
      else if (0==strcmp(vrn,"Y")) G.vt[i]=Y;
      else if (0==strcmp(vrn,"Z")) G.vt[i]=Z;
      else if (*vrn=='C') {
        G.vt[i]=Ci;
        G.vai[i]=0;
        while (*(++vrn) != 0) G.vai[i]=10*G.vai[i]+(*vrn-'0');
      } else G.vt[i]=UOV;
    }
    G.v[i]=new int[G.vd[i]];
    for (j=0; j<G.vd[i]; ++j) {
      MLGetInteger(link, &k);
      G.v[i][j]=k-1;
    }
    MLDisownSymbol(link, vname);
  }
  return link;
}

MLINK operator<<(MLINK link, Graph &G)
{
  char vname[16];	/* name of current vertex */
  int i,j;

  MLPutFunction(link, "List",2);
  MLPutInteger(link, G.coeff);
  MLPutFunction(link, "Graph", G.N);
  for (i=0; i<G.N; ++i) {
    MLPutFunction(link, VertexName(vname,G.vp[i],G.vt[i],G.vai[i]), G.vd[i]);
    for (j=0; j<G.vd[i]; ++j) MLPutInteger(link, G.v[i][j]+1);
  }
  return link;
}

char *VertexName(char *vname, int vp, VertexType vt, int vai)
{
  char *vrn;            /* reduced name */
  int i,j,k;
  char c;

  if (vp) {
    strcpy(vname,"Odd");
    vrn=vname+3;
  } else vrn=vname;
  switch(vt) {
  case Marker: strcpy(vrn,"Marker"); break;
  case UOV: strcpy(vrn,"UOV"); break;
  case OV: strcpy(vrn,"OV"); break;
  case Dashed: strcpy(vrn,"Dashed"); break;
  case Arrow: strcpy(vrn,"Arrow"); break;
  case UV: strcpy(vrn,"UV"); break;
  case Blue: strcpy(vrn,"Blue"); break;
  case Green: strcpy(vrn,"Green"); break;
  case Red: strcpy(vrn,"Red"); break;
  case X: strcpy(vrn,"X"); break;
  case Y: strcpy(vrn,"Y"); break;
  case Z: strcpy(vrn,"Z"); break;
  case Ci:
    strcpy(vrn++,"C");
    j=vai; k=0;
    while (j!=0) {
      *(vrn++) = (j % 10) + '0';
      j/=10; ++k;
    }
    for (j=1; j<=k/2; ++j) {
      c=*(vrn-j);
      *(vrn-j)=*(vrn-k+j-1);
      *(vrn-k+j-1)=c;
    }
    *vrn=0;
    break;
  }
  return vname;
}

ostream& operator<<(ostream& s, Graph &G)
{
  char vname[16];	/* name of current vertex */
  char *vrn;		/* reduced name of current vertex */
  int i,j,k;
  char c;

  if (G.coeff != 1) s << G.coeff << "*";
  s << "Graph[";
  for (i=0; i<G.N; ++i) {
    s << VertexName(vname,G.vp[i],G.vt[i],G.vai[i]);
    s << "[";
    for (j=0; j<G.vd[i]; ++j) {
      s << G.v[i][j]+1;
      if (j<G.vd[i]-1) s << ",";
    }
    s << "]";
    if (i<G.N-1) s<< ",";
  }
  s << "]";
  return s;
}

Graph &Permute(Graph &G, permutation *perm)
{
  Graph *PG;				/* the permuted G */
  static permutation invperm[MAXN];	/* The inverse of perm */
  int oddvs[MAXN];			/* list of odd vertices */
  int l=0;				/* The length of oddvs */
  int i,j;

  PG=new Graph;
  PG->N=G.N; PG->n=G.n; PG->coeff=G.coeff;
  PG->vd=new int[G.N];
  PG->vt=new VertexType[G.N];
  PG->vp=new int[G.N];
  PG->vai=new int[G.N];
  PG->v=new int*[G.N];
  for (i=0; i<G.n; ++i) invperm[perm[i]]=i;
  for (i=0; i<G.n; ++i) {
    PG->vd[i]=G.vd[perm[i]];
    PG->vt[i]=G.vt[perm[i]];
    PG->vp[i]=G.vp[perm[i]];
    PG->vai[i]=G.vai[perm[i]];
    PG->v[i]=new int[PG->vd[i]];
    switch (PG->vt[i]) {
    case Ci:
    case UV:
    case Blue:
    case Green:
    case Red:
    case X:
    case Y:
    case Z:
      PG->v[i][0]=invperm[G.v[perm[i]][0]];
      break;
    case Dashed:
    case UOV:
      for (j=0; j<PG->vd[i]; ++j) PG->v[i][j]=invperm[G.v[perm[i]][j]];
      SignedSort(PG->vd[i], PG->v[i]);
      break;
    case Arrow:
    case OV:
      for (j=0; j<PG->vd[i]; ++j) PG->v[i][j]=invperm[G.v[perm[i]][j]];
      PG->coeff*=SignedSort(PG->vd[i], PG->v[i]);
      break;
    default:
      for (j=0; j<PG->vd[i]; ++j) PG->v[i][j]=invperm[G.v[perm[i]][j]];
    }
    if (PG->vp[i]) oddvs[l++]=i;
  }
  for (i=G.n; i<G.N; ++i) {
    PG->vd[i]=G.vd[i];
    PG->vt[i]=G.vt[i];
    PG->vp[i]=G.vp[i];
    PG->vai[i]=G.vai[i];
    PG->v[i]=new int[G.vd[i]];
    for (j=0; j<PG->vd[i]; ++j) PG->v[i][j]=invperm[G.v[i][j]];
  }
  PG->coeff*=Signature(l,oddvs,perm);
  return *PG;
}

void MLCanonicalForm(MLINK link)
{
  Graph G;                      /* Current graph considered */
  link >> G;
  link << CanonicalForm(G);
}

Graph &CanonicalForm(Graph &G)
{
  nvector orbits[MAXN];
  static optionblk options =	/* The nauty options. */
  {
    TRUE,		/* getcanon: make canong and canonlab? */
    TRUE,		/* digraph: multiple edges or loops? */
    FALSE,		/* writeautoms: write automorphisms? */
    FALSE,		/* writemarkers: write stats on pts fixed, etc.? */
    FALSE,		/* defaultptn: set lab,ptn,active for single cell? */
    FALSE,		/* cartesian: use cartesian rep for writing automs? */
    CONSOLWIDTH,	/* linelength: max chars/line (excl. '\n') for output */
    (FILE*)NULL,	/* outfile: file for output, if any */
    NILFUNCTION,	/* userrefproc: replacement for usual refine procedure*/
    (void(*)()) OnAutomorphism,
			/* userautomproc: procedure called for each */
			/* automorphism */
    NILFUNCTION,	/* userlevelproc: procedure called for each level */
    NILFUNCTION,	/* usernodeproc: procedure called for each node */
    NILFUNCTION,	/* usertcellproc: replacement for targetcell procedure*/
    NILFUNCTION,	/* invarproc: procedure to compute vertex-invariant */
    0,			/* tc_level: max level for smart target cell choosing */
    0,			/* mininvarlevel: min level for invariant computation */
    0,			/* maxinvarlevel: max level for invariant computation */
    0,			/* invararg: value passed to (*invarproc)() */
    (groupblk*)NULL	/* groupopts: placeholder for future group options */
  };
  statsblk stats;		/* Statistics returned by nauty */
  setword workspace[50*MAXM];	/* nauty's workspace */
  nautycoloring coloring;	/* nauty's coloring of G */
  int m;			/* nauty's m */
  static graph g[MAXN*MAXM];	/* The nauty form of G */
  static graph cg[MAXN*MAXM];	/* The canonical form of g */

  nautycolors(G,coloring);
  m=(coloring.n+WORDSIZE-1)/WORDSIZE;
  tonauty(G,g);
  OnAutomorphismGraph = &G;
  nauty(
    g,coloring.lab,coloring.ptn,NILSET,orbits,
    &options,&stats,workspace,50*MAXM,m,coloring.n,cg
  );
  return Permute(G,coloring.lab);
}

extern "C" void OnAutomorphism(
  int count,
  permutation *perm,
  nvector *orbits,
  int numorbits,
  int stabvertex,
  int n
)
{
  if (OnAutomorphismGraph->PermuteSign(perm) != 1)
    OnAutomorphismGraph->SetCoefficient(0);
}

void tonauty(const Graph& G, graph *g)
{
  int n;	/* The order of G and g */
  int m;	/* nauty's m */
  set *gv;	/* A single vertex in g */
  int vd;	/* The length of gv */
  int i,j;

  n=G.Order();
  m=(n+WORDSIZE-1)/WORDSIZE;
  for (i=0; i<n; ++i) {
    gv =  GRAPHROW(g,i,m);
    EMPTYSET(gv,m);
    vd=G.VertexDegree(i);
    for (j=0; j<vd; ++j) ADDELEMENT(gv, G.Neighbor(i,j) % n);
  }
}

void nautycolors(const Graph& G, nautycoloring& c)
{
  int p,np;		/* Priority and next priority */
  int vp[MAXN];		/* List of vertex priorities */
  int i,j;

  c.n=G.Order();
  p=0;
  for (i=0; i<c.n; ++i) {
    vp[i]=G.VertexPriority(i);
    p=(p>vp[i] ? p : vp[i]);
  }
  j=0;
  while (p>0) {
    np=0;
    for (i=0; i<c.n; ++i) if (vp[i]<=p) {
      if (vp[i]==p) {
        c.lab[j]=i;
        c.ptn[j++]=1;
      } else np=(np>vp[i] ? np : vp[i]);
    }
    c.ptn[j-1]=0;
    p=np;
  }
}

int Signature(int l, int *list, permutation *perm)
{
  int sign=1;
  int sl[MAXN];	/* sorted `list' */
  int i;

  for(i=0; i<l; ++i) sl[i]=list[i];
  SignedSort(l,sl);
  for (i=0; i<l; ++i) sl[i]=perm[sl[i]];
  sign*=SignedSort(l,sl);
  return sign;
}

int SignedSort(int l, int *list)
{
  int sign=1;
  int lastflip,newlastflip;
  int i,t;

  lastflip=l-1;
  while(lastflip>0) {
    newlastflip=0;
    for(i=0; i<lastflip; ++i) if (list[i]==list[i+1]) sign=0;
    else if (list[i]>list[i+1]) {
      sign*=-1;
      newlastflip=i;
      t=list[i]; list[i]=list[i+1]; list[i+1]=t;
    }
    lastflip=newlastflip;
  }
  return sign;
}

int SplitGraph(Graph &G, Graph *&glist)
{
  char vname[16];		/* name of current vertex */
  int i,j,k,m;			/* Temporary variables */
  Graph *g;			/* Current component */

  int *a=new int[G.n];
  int *b=new int[G.n];
  int *c=new int[G.n];
  for (i=0; i<G.n; ++i) a[i]=i;
  int done;
  do {
    done=1;
    for (i=0; i<G.n; ++i) for (j=0; j<G.vd[i]; ++j)
      if (a[G.v[i][j]]<a[i]) {done=0; a[i]=a[G.v[i][j]];}
  } while (!done);
  int l=0;
  for (i=0; i<G.n; ++i) if (a[i]==i) {
    ++l; b[i]=0; c[i]=1;
  }
  else b[i]=c[a[i]]++;
  g=glist=new Graph[l];
  for (i=0; i<G.n; ++i) if (a[i]==i) {
    g->N=g->n=c[i]; g->coeff=1;
    g->vd=new int[c[i]]; g->vp=new int[c[i]];
    g->vai=new int[c[i]]; g->v=new int*[c[i]];
    g->vt=new VertexType[c[i]];
    m=0;
    for (j=i; j<G.n; ++j) if (a[j]==i) {
      g->vp[m]=G.vp[j]; g->vt[m]=G.vt[j]; g->vai[m]=G.vai[j];
      g->vd[m]=G.vd[j];
      g->v[m]=new int[G.vd[j]];
      for (k=0; k<G.vd[j]; ++k) g->v[m][k]=b[G.v[j][k]];
      ++m;
    }
    (*g)=CanonicalForm(*g);
    ++g;
  }
  delete a; delete b; delete c;
  return l;
}

void MLSplitGraph(MLINK link)
{
  Graph G;                      /* Current graph considered */
  Graph *glist;			/* List returned by SplitGraph */
  link >> G;
  int l=SplitGraph(G,glist);
  MLPutFunction(link, "Times", l);
  for (int i=0; i<l; ++i) link << glist[i];
  delete glist;
}