/* Program to perform K-group EM/BP community detection using the
 * degree-corrected SBM on an arbitrary network read from a GML file, with
 * discrete (enumerative) metadata stored in the "label" field
 *
 * Written by Mark Newman  28 NOV 2014
 * Modified to work with any number of groups  22 DEC 2014
 * Modified to initialize the parameters properly so that they add up
 *   to one  2 MAR 2015
 * Modified to calculate the final log-likelihood more accurately, using
 *   the Bethe free energy  14 MAY 2015
 */

/* Program control */

#define DEBUG
#define VERBOSE
#undef METASTART       // Set to start at the metadata solution
#undef DIAGSTART       // Set to start with gamma diagonally dominant
#define NOCONVERGE     // Set to abort if the solution doesn't converge
#define BPFIRST        // Set to start with an E-step
#define PRINTLL        // Set to print the log-likelihood at the end

/* Inclusions */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>

#include "EMPlugin.h"
#include "PluginManager.h"

//#include "readgml.h"

/* Constants */


/* Globals */
#include "network.h"

NETWORK G;            // Struct storing the network
int twom;             // Twice the number of edges

int *x;               // Metadata
int *nx;              // Number of nodes with each distinct metadata value
double **nrx;         // Expected number in each group with value
char **mlabel;        // Metadata strings
int nmlabels;         // Number of distinct metadata strings

double **gmma;        // Prior parameters (spelled "gmma" because "gamma"
                      //   is a reserved word in C math.h)
double omega[K][K];   // Mixing parameters

double ***eta;        // Messages
double **q;           // One-point marginals

gsl_rng *rng;         // Random number generator


/* Get metadata from the labels */

#define MAXMETA 300


// Functions to read a network stored in a GML file into a NETWORK struct
//
// Written by Mark Newman  11 AUG 06
// Changed to allow node labels containing the word "node", which previously
//   confused the (rather simple) code for counting network nodes  3 DEC 14
//
// To use this package, #include "readgml.h" at the head of your program
// and then call the following:
//
// Function calls:
//   int read_network(NETWORK *network, FILE *stream)
//     -- Reads a network from the FILE pointed to by "stream" into the
//        structure "network".  For the format of NETWORK structs see file
//        "network.h".  Returns 0 if read was successful.
//   void free_network(NETWORK *network)
//     -- Destroys a NETWORK struct again, freeing up the memory

#include "EMPlugin.h"

// Inclusions

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

//#include "network.h"

// Constants

#define LINELENGTH 1000

// Types

typedef struct line {
  char *str;
  struct line *ptr;
} LINE;

// Globals

LINE *first;
LINE *current;


// Function to read one line from a specified stream.  Return value is
// 1 if an EOF was encountered.  Otherwise 0.

int read_line(FILE *stream, char line[LINELENGTH])
{
  if (fgets(line,LINELENGTH,stream)==NULL) return 1;
  line[strlen(line)-1] = '\0';   // Erase the terminating NEWLINE
  return 0;
}


// Function to read in the whole file into a linked-list buffer, so that we
// can do several passes on it, as required to read the GML format
// efficiently

int fill_buffer(FILE *stream)
{
  int length;
  char line[LINELENGTH];
  LINE *previous;

  if (read_line(stream,line)!=0) {
    first = NULL;                // Indicates empty buffer
    return 1;
  }
  length = strlen(line) + 1;
  first = (LINE*) malloc(sizeof(LINE));
  first->str = (char*) malloc(length*sizeof(char));
  strcpy(first->str,line);

  previous = first;
  while (read_line(stream,line)==0) {
    length = strlen(line) + 1;
    previous->ptr = (LINE*) malloc(sizeof(LINE));
    previous = previous->ptr;
    previous->str = (char*) malloc(length*sizeof(char));
    strcpy(previous->str,line);
  }
  previous->ptr = NULL;          // Indicates last line

  return 0;
}


// Function to free up the buffer again

void free_buffer()
{
  LINE *thisptr;
  LINE *nextptr;

  thisptr = first;
  while (thisptr!=NULL) {
    nextptr = thisptr->ptr;
    free(thisptr->str);
    free(thisptr);
    thisptr = nextptr;
  }
}


// Function to reset to the start of the buffer again

void reset_buffer()
{
  current = first;
}


// Function to get the next line in the buffer.  Returns 0 if there was
// a line or 1 if we've reached the end of the buffer.

int next_line(char line[LINELENGTH])
{
  if (current==NULL) return 1;
  strcpy(line,current->str);
  current = current->ptr;
  return 0;
}



// Function to establish whether the network read from a given stream is
// directed or not.  Returns 1 for a directed network, and 0 otherwise.  If
// the GML file contains no "directed" line then the graph is assumed to be
// undirected, which is the GML default behavior.

int is_directed()
{
  int result=0;
  char *ptr;
  char line[LINELENGTH];

  reset_buffer();

  while (next_line(line)==0) {
    ptr = strstr(line,"directed");
    if (ptr==NULL) continue;
    sscanf(ptr,"directed %i",&result);
    break;
  }

  return result;
}


// Function to count the vertices in a GML file.  Returns number of vertices.

int count_vertices()
{
  int result=0;
  char *nonspace;
  char line[LINELENGTH];

  reset_buffer();

  while (next_line(line)==0) {
    for (nonspace=line; *nonspace==' '; nonspace++);
    if (strncmp(nonspace,"node",4)==0) result++;
  }

  return result;
}


// Function to compare the IDs of two vertices

int cmpid(VERTEX *v1p, VERTEX *v2p)
{
  if (v1p->id>v2p->id) return 1;
  if (v1p->id<v2p->id) return -1;
  return 0;
}


// Function to allocate space for a network structure stored in a GML file
// and determine the parameters (id, label) of each of the vertices.

void create_network(NETWORK *network)
{
  int i;
  int length;
  char *ptr;
  char *start,*stop;
  char *nonspace;
  char line[LINELENGTH];
  char label[LINELENGTH];

  // Determine whether the network is directed

  network->directed = is_directed();

  // Count the vertices

  network->nvertices = count_vertices();

  // Make space for the vertices

  network->vertex = (VERTEX*) calloc(network->nvertices,sizeof(VERTEX));

  // Go through the file reading the details of each vertex one by one

  reset_buffer();
  for (i=0; i<network->nvertices; i++) {

    // Skip to next "node" entry.  Check that the start of the line is "node"
    // explicitly, to avoid other lines that might contain "node"

    do {
      next_line(line);
      for (nonspace=line; *nonspace==' '; nonspace++);
    } while (strncmp(nonspace,"node",4)!=0);

    // Read in the details of this vertex

    do {

      // Skip leading spaces

      for (nonspace=line; *nonspace==' '; nonspace++);

      // Look for ID

      if (strncmp(nonspace,"id",2)==0) {
	sscanf(nonspace,"id %i",&network->vertex[i].id);
      }

      // Look for label

      if (strncmp(nonspace,"label",5)==0) {
	start = strchr(line,'"');
	if (start==NULL) {
	  sscanf(nonspace,"label %s",&label);
	} else {
	  stop = strchr(++start,'"');
	  if (stop==NULL) length = strlen(line) - (start-line);
	  else length = stop - start;
	  strncpy(label,start,length);
	  label[length] = '\0';
	  network->vertex[i].label = (char*) malloc((length+1)*sizeof(char));
	  strcpy(network->vertex[i].label,label);
	}
      }

      // If we see a closing square bracket we are done

      if (strstr(line,"]")!=NULL) break;

    } while (next_line(line)==0);

  }

  // Sort the vertices in increasing order of their IDs so we can find them
  // quickly later

  qsort(network->vertex,network->nvertices,sizeof(VERTEX),/*(void*)*/(__compar_fn_t) cmpid);
}


// Function to find a vertex with a specified ID using binary search.
// Returns the element in the vertex[] array holding the vertex in question,
// or -1 if no vertex was found.

int find_vertex(int id, NETWORK *network)
{
  int top,bottom,split;
  int idsplit;

  top = network->nvertices;
  if (top<1) return -1;
  bottom = 0;
  split = top/2;

  do {
    idsplit = network->vertex[split].id;
    if (id>idsplit) {
      bottom = split + 1;
      split = (top+bottom)/2;
    } else if (id<idsplit) {
      top = split;
      split = (top+bottom)/2;
    } else return split;
  } while (top>bottom);

  return -1;
}
    

// Function to determine the degrees of all the vertices by going through
// the edge data

void get_degrees(NETWORK *network)
{
  int s,t;
  int vs,vt;
  char *ptr;
  char line[LINELENGTH];

  reset_buffer();

  while (next_line(line)==0) {

    // Find the next edge entry

    ptr = strstr(line,"edge");
    if (ptr==NULL) continue;

    // Read the source and target of the edge

    s = t = -1;

    do {

      ptr = strstr(line,"source");
      if (ptr!=NULL) sscanf(ptr,"source %i",&s);
      ptr = strstr(line,"target");
      if (ptr!=NULL) sscanf(ptr,"target %i",&t);

      // If we see a closing square bracket we are done

      if (strstr(line,"]")!=NULL) break;

    } while (next_line(line)==0);

    // Increment the degrees of the appropriate vertex or vertices

    if ((s>=0)&&(t>=0)) {
      vs = find_vertex(s,network);
      network->vertex[vs].degree++;
      if (network->directed==0) {
	vt = find_vertex(t,network);
	network->vertex[vt].degree++;
      }
    }

  }

  return;
}


// Function to read in the edges

void read_edges(NETWORK *network)
{
  int i;
  int s,t;
  int vs,vt;
  int *count;
  double w;
  char *ptr;
  char line[LINELENGTH];

  // Malloc space for the edges and temporary space for the edge counts
  // at each vertex

  for (i=0; i<network->nvertices; i++) {
    network->vertex[i].edge = (EDGE*) malloc(network->vertex[i].degree*sizeof(EDGE));
  }
  count = (int*) calloc(network->nvertices,sizeof(int));

  // Read in the data

  reset_buffer();

  while (next_line(line)==0) {

    // Find the next edge entry

    ptr = strstr(line,"edge");
    if (ptr==NULL) continue;

    // Read the source and target of the edge and the edge weight

    s = t = -1;
    w = 1.0;

    do {

      ptr = strstr(line,"source");
      if (ptr!=NULL) sscanf(ptr,"source %i",&s);
      ptr = strstr(line,"target");
      if (ptr!=NULL) sscanf(ptr,"target %i",&t);
      ptr = strstr(line,"value");
      if (ptr!=NULL) sscanf(ptr,"value %lf",&w);

      // If we see a closing square bracket we are done

      if (strstr(line,"]")!=NULL) break;

    } while (next_line(line)==0);

    // Add these edges to the appropriate vertices

    if ((s>=0)&&(t>=0)) {
      vs = find_vertex(s,network);
      vt = find_vertex(t,network);
      network->vertex[vs].edge[count[vs]].target = vt;
      network->vertex[vs].edge[count[vs]].weight = w;
      count[vs]++;
      if (network->directed==0) {
	network->vertex[vt].edge[count[vt]].target = vs;
	network->vertex[vt].edge[count[vt]].weight = w;
	count[vt]++;
      }
    }

  }

  free(count);
  return;
}


// Function to read a complete network

int read_network(NETWORK *network, FILE *stream)
{
  fill_buffer(stream);
  create_network(network);
  get_degrees(network);
  read_edges(network);
  free_buffer();

  return 0;
}


// Function to free the memory used by a network again

void free_network(NETWORK *network)
{
  int i;

  for (i=0; i<network->nvertices; i++) {
    free(network->vertex[i].edge);
    free(network->vertex[i].label);
  }
  free(network->vertex);
}


void get_metadata()
{
  int u,i;

  /* Make space for the metadata numbers and labels */

  x = (int*) malloc(G.nvertices*sizeof(int));
  mlabel = (char**) malloc(MAXMETA*sizeof(char*));
  nmlabels = 0;

  /* Go through the vertices */

  for (u=0; u<G.nvertices; u++) {

    /* Check to see if this label is already in the list of labels */

    for (i=0; i<nmlabels; i++) {
      if (strcmp(G.vertex[u].label,mlabel[i])==0) break;
    }

    /* If not, add it */

    if (i==nmlabels) {
      mlabel[nmlabels++] = G.vertex[u].label;  // Just set pointers equal
    }

    /* Record this as the metadata type for this vertex */

    x[u] = i;
  }

  /* Count how many nodes there are in each metadata group */

  nx = (int*) calloc(nmlabels,sizeof(int));
  for (u=0; u<G.nvertices; u++) nx[x[u]]++;

#ifdef VERBOSE
  fprintf(stderr,"Found %i distinct metadata values:\n",nmlabels);
  for (i=0; i<nmlabels; i++) fprintf(stderr," %i %s\n",i,mlabel[i]);
#endif
}


/* Function to generate d numbers at random that add up to unity */

void random_unity(int d, double *x)
{
  int k;
  double sum=0.0;

  for (k=0; k<d; k++) {
    x[k] = gsl_rng_uniform(rng);
    sum += x[k];
  }
  for (k=0; k<d; k++) x[k] /= sum;
}


/* Do BP */

#define BP_ACC 1e-4
#define BP_MAXSTEPS 20
#define SMALL 1.0e-100

int bp()
{
  int i,j;
  int u,v;
  int r,s;
  int steps;
  double deltaeta,maxdelta;
  double logeta,neweta,sum,norm,largest;
  double d[K];
  double logpre[K];
  double logqun[K];
  double ***logetaun;

  // Make space for the new log-etas, which are called "logetaun" because
  // the are initially calculated in unnormalized form

  logetaun = (double***) malloc(G.nvertices*sizeof(double**));
  for (u=0; u<G.nvertices; u++) {
    logetaun[u] = (double**) malloc(G.vertex[u].degree*sizeof(double*));
    for (i=0; i<G.vertex[u].degree; i++) {
      logetaun[u][i] = (double*) malloc(K*sizeof(double));
    }
  }

  // Main BP loop

  steps = 0;
  do {

    /* Calculate the expected group degrees */
    
    for (r=0; r<K; r++) {
      d[r] = 0.0;
      for (u=0; u<G.nvertices; u++) d[r] += q[u][r]*G.vertex[u].degree;
    }
  
    /* Calculate the log-prefactors (without the leading factor of d_i or
     * the prior) */

    for (r=0; r<K; r++) {
      logpre[r] = 0.0;
      for (s=0; s<K; s++) logpre[r] -= omega[r][s]*d[s];
    }

    /* Calculate new values for the one-vertex marginals */

#ifdef VERBOSE
    fprintf(stderr,"Calculating one-vertex marginals...    \r");
#endif
    for (u=0; u<G.nvertices; u++) {
      for (r=0; r<K; r++) {
	logqun[r] = log(gmma[r][x[u]]) + G.vertex[u].degree*logpre[r];
	for (i=0; i<G.vertex[u].degree; i++) {
	  sum = 0.0;
	  for (s=0; s<K; s++) sum += eta[u][i][s]*omega[r][s];
	  if (sum<SMALL) sum = SMALL;
	  logqun[r] += log(sum);
	}
	if (r==0) largest = logqun[r];
	else if (logqun[r]>largest) largest = logqun[r];
      }

      /* Normalize */

      norm = 0.0;
      for (r=0; r<K; r++) {
	logqun[r] -= largest;
	norm += exp(logqun[r]);
      }
      for (r=0; r<K; r++) q[u][r] = exp(logqun[r])/norm;
    }

    /* Calculate (unnormalized) new values for the (log) messages */

#ifdef VERBOSE
    fprintf(stderr,"Calculating messages...              \r");
#endif
    for (u=0; u<G.nvertices; u++) {
      for (i=0; i<G.vertex[u].degree; i++) {
	v = G.vertex[u].edge[i].target;
	for (r=0; r<K; r++) {
	  logeta = log(gmma[r][x[v]]) + G.vertex[v].degree*logpre[r];
	  for (j=0; j<G.vertex[v].degree; j++) {
	    if (G.vertex[v].edge[j].target!=u) {
	      sum = 0.0;
	      for (s=0; s<K; s++) sum += eta[v][j][s]*omega[r][s];
	      if (sum<SMALL) sum = SMALL;   // Prevent -Inf
	      logeta += log(sum);
	    }
	  }
	  logetaun[u][i][r] = logeta;
	}
      }
    }

    /* Normalize and calculate largest change */

#ifdef VERBOSE
    fprintf(stderr,"Normalizing messages...\r");
#endif
    maxdelta = 0.0;
    for (u=0; u<G.nvertices; u++) {
      for (i=0; i<G.vertex[u].degree; i++) {
	norm = 0.0;
	largest = logetaun[u][i][0];
	for (r=1; r<K; r++) {
	  if (logetaun[u][i][r]>largest) largest = logetaun[u][i][r];
	}
	for (r=0; r<K; r++) {
	  logetaun[u][i][r] -= largest;
	  norm += exp(logetaun[u][i][r]);
	}	  
	for (r=0; r<K; r++) {
	  neweta = exp(logetaun[u][i][r])/norm;
	  deltaeta = fabs(neweta-eta[u][i][r]);
	  if (deltaeta>maxdelta) maxdelta = deltaeta;
	  eta[u][i][r] = neweta;
	}
      }
    }

#ifdef VERBOSE
    fprintf(stderr,"BP steps %i, max change = %g                   \r",
	    steps,maxdelta);
#endif

  } while ((maxdelta>BP_ACC)&&(++steps<=BP_MAXSTEPS));

#ifdef VERBOSE
  fprintf(stderr,"\n");
#endif

  // Free space

  for (u=0; u<G.nvertices; u++) {
   for (i=0; i<G.vertex[u].degree; i++) free(logetaun[u][i]);
    free(logetaun[u]);
  }
  free(logetaun);

  return steps;
}


// Function to calculate new values of the parameters

double params()
{
  int u,v;
  int i,j;
  int r,s;
  double norm,esum,quvrs;
  double d[K];
  double term[K][K];
  double sum[K][K];
  double L;

  // Calculate some basics

  for (r=0; r<K; r++) {
    d[r] = 0.0;
    for (i=0; i<nmlabels; i++) nrx[r][i] = 0.0;
  }

  for (u=0; u<G.nvertices; u++) {
    for (r=0; r<K; r++) {
      d[r] += q[u][r]*G.vertex[u].degree;
      nrx[r][x[u]] += q[u][r];
    }
  }

  // Calculate new values of the gammas

  for (r=0; r<K; r++) {
    for (i=0; i<nmlabels; i++) gmma[r][i] = nrx[r][i]/nx[i];
  }

  // Calculate the new values of the omegas

  // Zero out the sum variables

  for (r=0; r<K; r++) {
    for (s=0; s<K; s++) sum[r][s] = 0.0;
  }
  esum = 0.0;

  // Perform the sums

  for (u=0; u<G.nvertices; u++) {
    for (i=0; i<G.vertex[u].degree; i++) {
      v = G.vertex[u].edge[i].target;

      // Find which edge leads back from v to u

      for (j=0; j<G.vertex[v].degree; j++) {
	if (G.vertex[v].edge[j].target==u) break;
      }
      if (j==G.vertex[v].degree) {
	fprintf(stderr,"Error!\n");
	exit(23);
      }

      // Calculate the terms and the normalization factor

      norm = 0.0;
      for (r=0; r<K; r++) {
	for (s=0; s<K; s++) {
	  term[r][s] = omega[r][s]*eta[u][i][r]*eta[v][j][s];
	  norm += term[r][s];
	}
      }

      // Add to the running sums

      for (r=0; r<K; r++) {
	for (s=0; s<K; s++) {
	  quvrs = term[r][s]/norm;
	  sum[r][s] += quvrs;
	  esum += quvrs*log(quvrs);
	}
      }
    }
  }

  // Calculate the new values of the omega variables (after calculating
  // the likelihood using the old omegas)

  for (r=0; r<K; r++) {
    for (s=0; s<K; s++) omega[r][s] = sum[r][s]/(d[r]*d[s]);
  }

  // Calculate the expected log-likelihood

  // Internal energy first

  L = 0.0;
  for (r=0; r<K; r++) {
    for (s=0; s<K; s++) L += 0.5*sum[r][s]*log(omega[r][s]);
    for (i=0; i<nmlabels; i++) {
      if (gmma[r][i]>0.0) L += nx[i]*gmma[r][i]*log(gmma[r][i]);
    }
  }

  // Now the entropy

  L -= 0.5*esum;
  for (u=0; u<G.nvertices; u++) {
    for (r=0; r<K; r++) {
      if ((q[u][r]>0.0)&&(G.vertex[u].degree>0)) {
	L += (G.vertex[u].degree-1)*q[u][r]*log(q[u][r]);
      }
    }
  }

  return L;
}


#define EM_MAXSTEP 100
#define EM_ACC 1e-4

//main(int argc, char *argv[])
void EMPlugin::input(std::string file)
{
  int /*u,*/v/*,i*/;/*r,s;*/
  //int step;
  double norm;// ,deltac, maxdelta;
  double ru[K];

  // Initialize random number generator

  rng = gsl_rng_alloc(gsl_rng_mt19937);
  gsl_rng_set(rng,time(NULL));

  // Read the network and the metadata from stdin

#ifdef VERBOSE
  fprintf(stderr,"Reading network...\n");
#endif
  FILE* inputfile = fopen(file.c_str(), "r");
  if (!inputfile) printf("OPEN OF %s FAILED.\n", file.c_str());
  read_network(&G,inputfile);
  //read_network(&G,stdin);
  for (u=twom=0; u<G.nvertices; u++) twom += G.vertex[u].degree;
  get_metadata();

  // Make space for the marginals and initialize to random initial values,
  // unless METASTART is set, in which case initialize to the metadata
  // values

  eta = (double***) malloc(G.nvertices*sizeof(double**));
  q = (double**) malloc(G.nvertices*sizeof(double*));

  for (u=0; u<G.nvertices; u++) {
    q[u] = (double*) malloc(K*sizeof(double));
#ifdef METASTART
    for (r=0; r<K; r++) q[u][r] = (r==x[u]?1:0);
#else
    random_unity(K,q[u]);
#endif
  }

  // Make space for the messages and initialize to the same values as the
  // marginals

  for (u=0; u<G.nvertices; u++) {
    eta[u] = (double**) malloc(G.vertex[u].degree*sizeof(double*));
    for (i=0; i<G.vertex[u].degree; i++) {
      eta[u][i] = (double*) malloc(K*sizeof(double));
      v = G.vertex[u].edge[i].target;
      for (r=0; r<K; r++) eta[u][i][r] = q[v][r];
    }
  }

  nrx = (double**) malloc(K*sizeof(double*));
  for (r=0; r<K; r++) nrx[r] = (double*) malloc(nmlabels*sizeof(double));

  // Malloc space for parameters gmma and choose random initial values

  gmma = (double**) malloc(K*sizeof(double*));
  for (r=0; r<K; r++) gmma[r] = (double*) malloc(nmlabels*sizeof(double));
  for (i=0; i<nmlabels; i++) {
#ifdef DIAGSTART
    norm = 0.0;
    for (r=0; r<K; r++) {
      if (r==i) gmma[r][i] = gsl_rng_uniform(rng);
      else gmma[r][i] = 0.1*gsl_rng_uniform(rng);
      norm += gmma[r][i];
    }
    for (r=0; r<K; r++) gmma[r][i] /= norm;
#else
    random_unity(K,ru);
    for (r=0; r<K; r++) gmma[r][i] = ru[r];
#endif
  }

  // Choose random values for the omegas, but with a bias toward
  // assortative choices (change if necessary for other networks)

  for (r=0; r<K; r++) {
    for (s=0; s<K; s++) {
      if (r==s) c[r][s] = 1 + gsl_rng_uniform(rng);
      else if (r<s) c[r][s] = gsl_rng_uniform(rng);
      else c[r][s] = c[s][r];
      omega[r][s] = c[r][s]/twom;
    }
  }
}

void EMPlugin::run() {
  // EM loop

#ifdef VERBOSE
  fprintf(stderr,"Starting EM algorithm...\n");
#endif
  step = 0;
  do {

    // Run BP to calculate the messages and one-vertex marginals

#ifdef BPFIRST
    bpsteps = bp();
#endif

    // Calculate the new values of the parameters

    L = params();

    // Calculate the new values of the c variables

    for (r=0; r<K; r++) {
      for (s=0; s<K; s++) {
	oldc[r][s] = c[r][s];
	c[r][s] = omega[r][s]*twom;
      }
    }

#ifndef BPFIRST
    bpsteps = bp();
#endif

    // Find the largest change in any of the c's

    maxdelta = 0.0;
    for (r=0; r<K; r++) {
      for (s=0; s<K; s++) {
	deltac = fabs(c[r][s]-oldc[r][s]);
        if (deltac>maxdelta) maxdelta = deltac;
      }
    }

    // Print out new values of the parameters

#ifdef VERBOSE
    fprintf(stderr,"EM step %i, max change = %g\n",step,maxdelta);
    fprintf(stderr,"gamma =\n");
    for (r=0; r<K; r++) {
      for (i=0; i<nmlabels; i++) fprintf(stderr," %.6f",gmma[r][i]);
      fprintf(stderr,"\n");
    }

    fprintf(stderr,"c =\n");
    for (r=0; r<K; r++) {
      for (s=0; s<K; s++) fprintf(stderr," %.6f",c[r][s]);
      fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");
#endif

    if (++step>EM_MAXSTEP) {
#ifdef NOCONVERGE
      fprintf(stderr,"Solution failed to converge in %i EM steps\n",
	      EM_MAXSTEP);
#endif
      break;
    }

  } while (maxdelta>EM_ACC);

#ifdef NOCONVERGE
  if (bpsteps>BP_MAXSTEPS) {
    fprintf(stderr,"BP failed converge on final EM step\n");
  }
#endif

#ifdef VERBOSE
  fprintf(stderr,"Log-likelihood = %g\n",L);
#endif
}

void EMPlugin::output(std::string file) {
  // Output the results
  FILE* outputfile = fopen(file.c_str(), "w");

#ifdef PRINTLL
  fprintf(outputfile, "%g\n",L);
#endif
  for (u=0; u<G.nvertices; u++) {
    fprintf(outputfile, "%i %s",u,mlabel[x[u]]);
    for (r=0; r<K; r++) fprintf(outputfile, " %.6f",q[u][r]);
    fprintf(outputfile,"\n");
  }
}



PluginProxy<EMPlugin> EMPluginProxy = PluginProxy<EMPlugin>("EM", PluginManager::getInstance());
