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

#define NORMAL 0
#define CHOSEN 1

#define MAX_NINDEX 45 /* Max number of subgraph size boxes to be considered */
#define MAX_CANDLINKS 4000000 /* Max number of candidate links for the next node */

#define NR_END 1
#define FREE_ARG char*

struct node
{
        int data;
        struct node *next;
} **neigh_list;

int *neighbors;
int labels[MAX_NINDEX];

void nrerror(char error_text[]);
int *ivector(long nl, long nh);
int **imatrix(long nrl, long nrh, long ncl, long nch);
void free_ivector(int *v, long nl, long nh);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void add_node(int i, struct node** headRef);
void clear_lists(int cluster_size);
void calc_dens_hist(char *, int, int, int **);
int read_network(char *);

int int_cmp(const void *a, const void *b)
{
        int x=*((int *)a);
        int y=*((int *)b);
        if(x==y) return 0;
        else if(x<y) return -1;
        else return 1;
}


int main(argc,argv)
int argc;
char *argv[];
{
	double ks,d1,d2,dt,fn1,fn2;
	FILE *fp;
	int RUNS,cluster_size,cluster_size2,i,maxnple,**dens_hist,**dens_hist2,MAX_n,j1,j2;

	if(argc==5)
	{
		RUNS=atoi(argv[3]);
		MAX_n=atoi(argv[4]);
	}
	else
        {
                printf("Usage: %s filename1 filename2 runs max_n\n",argv[0]);
		return(-1);
	}
	
	/* Read the network size of the two clusters */
	fp=fopen(argv[1],"r");
	if(fp==NULL)
	{
		printf("File %s does not exist\n",argv[1]);
		exit(-1);
	}
	fscanf(fp,"%d\n",&cluster_size);
	fclose(fp);
	fp=fopen(argv[2],"r");
	if(fp==NULL)
	{
		printf("File %s does not exist\n",argv[2]);
		exit(-1);
	}
	fscanf(fp,"%d\n",&cluster_size2);
	fclose(fp);
	
	if(cluster_size<cluster_size2) i=cluster_size;
	else i=cluster_size2;
	if(i<10) maxnple=i;
	else if(i<100) maxnple=10*(i/10);
	else if(i<1000) maxnple=100*(i/100);
	else if(i<10000) maxnple=1000*(i/1000);
	else if(i<100000) maxnple=10000*(i/10000);
	else maxnple=100000;
	if(MAX_n>0 && MAX_n<maxnple) maxnple=MAX_n;

	dens_hist=imatrix(0,MAX_NINDEX-1,0,RUNS-1);
	dens_hist2=imatrix(0,MAX_NINDEX-1,0,RUNS-1);
	for(i=0;i<MAX_NINDEX;i++)
	{
		labels[i]=0;
	}
	
	calc_dens_hist(argv[1],RUNS,maxnple,dens_hist);
	calc_dens_hist(argv[2],RUNS,maxnple,dens_hist2);
	
	i=0;
	while(labels[i]>0)
	{
		qsort(dens_hist[i],RUNS,sizeof(int),int_cmp);
		qsort(dens_hist2[i],RUNS,sizeof(int),int_cmp);

		fn1=0.0;
		fn2=0.0;
		j1=0;
		j2=0;
		ks=0.0;
		while(j1<RUNS && j2<RUNS)
		{
			if((d1=dens_hist[i][j1]) <= (d2=dens_hist2[i][j2])) fn1=j1++/(double)RUNS;
			if(d2<=d1) fn2=j2++/(double)RUNS;
			if((dt=fabs(fn2-fn1))> ks) ks=dt;
		}
		printf("%d\t%g\n",labels[i],ks);
		i++;
	}	

	free_imatrix(dens_hist,0,MAX_NINDEX-1,0,RUNS-1);
	free_imatrix(dens_hist2,0,MAX_NINDEX-1,0,RUNS-1);
	
	return(0);
}

void calc_dens_hist(name,RUNS,maxnple,dens_hist)
char *name;
int RUNS,maxnple,**dens_hist;
{
	int i,j,runs,cluster_size,nple,nlist,**edgelist,ksum,kr,end,k;
	int *chosenlist,density,*status,nindex;
	struct node *temp;
	
	/* Read the network. cluster_size is the real network size */
	cluster_size=read_network(name);
	
	edgelist=imatrix(0,MAX_CANDLINKS-1,0,1);

	status=ivector(0,cluster_size-1);
	chosenlist=ivector(0,cluster_size-1);

	ksum=0;
	for(i=0;i<cluster_size;i++) ksum+=neighbors[i];

	for(runs=0;runs<RUNS;runs++)
	{
		srand48(runs);
		for(i=0;i<cluster_size;i++) status[i]=NORMAL;

		/* Choose a random link */
		kr=1+(int)(drand48()*ksum);
		j=0;
		k=neighbors[j];
		while(k<kr)
		{
			j++;
			k+=neighbors[j];
		}
		temp=neigh_list[j];
		for(i=0;i<k-kr;i++) temp=temp->next;
		i=temp->data;
		/* Store the initial link and nodes */
		status[i]=CHOSEN;
		status[j]=CHOSEN;
		chosenlist[0]=i;
		chosenlist[1]=j;
		nple=2;
		nindex=0;

		/* Fill in the list with candidate links */
		nlist=0; /* nlist shows the number of candidate links */
		temp=neigh_list[i];
		while(temp!=NULL)
		{
			if(status[temp->data]!=CHOSEN)
			{
				edgelist[nlist][0]=i;
				edgelist[nlist][1]=temp->data;
				nlist++;
				if(nlist>=MAX_CANDLINKS) {printf("Error: MAX POSSIBLE LINKS EXCEEDED!\n");exit(-1);}
			}
			temp=temp->next;
		}
		temp=neigh_list[j];
		while(temp!=NULL)
		{
			if(status[temp->data]!=CHOSEN)
			{
				edgelist[nlist][0]=j;
				edgelist[nlist][1]=temp->data;
				nlist++;
				if(nlist>=MAX_CANDLINKS) {printf("Error: MAX POSSIBLE LINKS EXCEEDED!\n");exit(-1);}
			}
			temp=temp->next;
		}
		
		density=1; /* density counts the number of links in the current selection */

		while(nple<maxnple)
		{
			/* Choose a random link from the list and store node in the cluster */
			kr=(int)(drand48()*nlist);
			j=edgelist[kr][1];
			status[j]=CHOSEN;
			chosenlist[nple]=j;
	
			/* Insert the new links in the candidate list */
			temp=neigh_list[j];
			while(temp!=NULL)
			{
				if(status[temp->data]!=CHOSEN)
				{
					edgelist[nlist][0]=j;
					edgelist[nlist][1]=temp->data;
					nlist++;
					if(nlist>=MAX_CANDLINKS) {printf("Error: MAX POSSIBLE LINKS EXCEEDED!\n");exit(-1);}
				}
				else density++;
				temp=temp->next;
			}
			/* Clean up list */
			end=nlist-1;
			i=0;
			while(i<=end)
			{
				while(edgelist[i][1]==j && i<=end)
				{
					edgelist[i][0]=edgelist[end][0];
					edgelist[i][1]=edgelist[end][1];
					end--;
				}
				i++;
			}
			nlist=end+1;
			nple++;
			
			/* Count the density (number of links) in the cluster */
			if(nple==cluster_size || nple<10 || (nple<100 && nple%10==0) || (nple<=1000 && nple%100==0) || (nple<=10000 && nple%1000==0) || (nple<=100000 && nple%10000==0))
			{
				if(runs==0) labels[nindex]=nple;
				dens_hist[nindex][runs]=density;
				nindex++;
			}
		}
	}

	free_imatrix(edgelist,0,MAX_CANDLINKS-1,0,1);
	free_ivector(status,0,cluster_size-1);
	free_ivector(chosenlist,0,cluster_size-1);
	clear_lists(cluster_size);

	return;
}

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
	fprintf(stderr,"Run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
	int *v;
	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl+NR_END;
}


int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	int **m;
	/* allocate pointers to rows */
	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;
	/* allocate rows and set pointers to them */
	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;
	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
	/* return pointer to array of pointers to rows */
	return m;
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

int read_network(name)
char name[100];
{
        int i,j,size;
        FILE *fp;

        fp=fopen(name,"r");
        fscanf(fp,"%d\n",&size);
        neighbors=ivector(0,size-1);
        for(i=0;i<size;i++) neighbors[i]=0;
        /* Allocate neighbor lists memory for each node */
        neigh_list=(struct node **)malloc(size*sizeof(struct node));
        for(i=0;i<size;i++) neigh_list[i]=NULL;

        while(!feof(fp))
        {
                fscanf(fp,"%d\t%d\n",&i,&j);
		add_node(j,&neigh_list[i]);
		add_node(i,&neigh_list[j]);
                neighbors[i]++;
		neighbors[j]++;
        }
        fclose(fp);

        return(size);
}

void add_node(i,headRef)
int i;
struct node** headRef;
{ 
        struct node *newnode;
 
        newnode = (struct node *)malloc(sizeof(struct node));
        newnode->data=i;
        newnode->next=*headRef;
 
 	*headRef=newnode;
	
        return;
}

void clear_lists(cluster_size)
int cluster_size;
{
	struct node *current,*cur2;
	int i;
	
	for(i=0;i<cluster_size;i++)
	{
		current=neigh_list[i];
               	while(current!=NULL)
               	{
                       	cur2=current->next;
                       	free(current);
                       	current=cur2;
	       	}
       	       	neigh_list[i]=NULL;
	}
	free(neigh_list);
	free_ivector(neighbors,0,cluster_size-1);
	return;
}