//
// Reads adjacency matrix and returns various measurements thereof
// (clustering coeff., mean minimum path, assortativity, degree distribution,
// mixing distribution, etc.)
//
// Input matrix should come as 3 columns: node i, node j, adjacency(i,j)
// (alternatively, appropriate changes must be made as of line 152).
// If data is, e.g., in file stuff.txt, info is put in files
// p_stuff.txt (distribution), mix_stuff.txt (mixing) and sw_stuff.txt 
// (various magnitudes and their expected values for an equivalent network
// in the Erdös-Rényi and in the Configuration ensembles). 
//
// Example compilation: gcc nat.c -o nat.exe -lm
// Example execution: ./nat.exe stuff.txt
//
// Developed by Samuel Johnson. Thanks to Joaquín J. Torres for providing
// codes for some functions. For any queries, please write to
// samuel@onsager.ugr.es.
//
// Instituto Carlos I de Física Teórica y Computacional,
// Dpto. de Electromagnetismo y Física de la Materia,
// Universidad de Granada, Granada, Spain.
// www.ugr.es/~sam


#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define theta(X,Y)    ((X>=Y)? 1.:0.)
#define MIN(X,Y)	((X>Y)? Y:X)
#define ABSO(X)     ((X>=0.0)? X :-X)
#define NR_END 1
#define FREE_ARG char*
#define PARAMETROS  1
#define ARCHIVO1 1
#define ARCHIVO2 2
#define Cmax (50) // max # of characters for name.dat
int *ivector(long, long);
double **dmatrix(long, long, long, long);
float *vector(long, long);
double *dvector(long, long);
void nrerror(char []);
float **matrix(long,long,long,long);
int **imatrix(long, long, long, long);
void free_ivector(int *v, long nl, long nh);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_dvector(double *v, long nl, long nh);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);

void short_path();
int Dim;
double L, Lm, L2m, sigL;
int count, count_disc;
double mean_disc, mean_disc_m;

void clustering_coeff();
double C, Cm, C2m, sigC;
double L_m, L2_m, sigL_m;
int tria;
int q, q1, q2, q3;
int **los;

double dispersion();
double L_conf,C_conf,L_rand,C_rand,C_reg,L_reg,disp;
double cosa, nei_conf;

void neighbours();
double *nei, nei_med, *nei_m;
double r, nei_k2, r_conf;


double nu, nuc, sigma, eta, mini;
double **n_bas;
int **adj;
int *kt;

double tauref=0;
int **epsi,*k,*n;
double *p;
int j,q,n_max,con;
double k_med, k2_med, k3_med;
int Nbig=12000;   // Input matrix can't be bigger than Nbig x Nbig
int N;
long i,ii;	
FILE *ar,*br;
FILE *deg_dis;
FILE *sw;
FILE *mix;

char p_data[Cmax];
char sw_data[Cmax];
char mix_data[Cmax];
double xvar1,xvar2;

int got_los=0;

double pear;

int main (int argc,char **argv){

    epsi=imatrix(0,Nbig-1,0,Nbig-1);
    k=ivector(0,Nbig-1);
    n=ivector(0,Nbig-1);
    p=dvector(0,Nbig-1);
    los=imatrix(0,Nbig-1,0,Nbig-1);
    nei=dvector(0,Nbig-1);
    nei_m=dvector(0,Nbig-1);

    adj=imatrix(0,Nbig-1,0,Nbig-1);
    kt=ivector(0,Nbig-1);

    n_bas=dmatrix(0,Nbig-1,0,Nbig-1);

    for (i=0;i<Nbig;i++)
    {
	k[i]=0;
	for (j=0;j<Nbig;j++) 
	{
	    epsi[i][j]=0;
	}
    }  


    if (argc<=PARAMETROS) {
	puts("Command parameters");
	puts("-----------------------");
	puts("File to read the data");
	exit(0);
    } 
    
 
    if (!(ar=fopen(argv[ARCHIVO1],"r"))) {
	printf("%s can't be read\n",argv[ARCHIVO1]);
	exit(0);
    }


    sprintf(p_data,"p_simp_%s",argv[ARCHIVO1]);
    sprintf(sw_data,"sw_simp_%s",argv[ARCHIVO1]);
    sprintf(mix_data,"mix_simp_%s",argv[ARCHIVO1]);
    deg_dis = fopen(p_data,"w");
    sw = fopen(sw_data,"w");
    mix=fopen(mix_data,"w");

    printf("# %s -> %s, %s, %s\n",argv[ARCHIVO1],sw_data,p_data,mix_data);
    
    n_max=0;


//  Read data:
    for(i=1;;i++) {

	fscanf(ar,"%lf\n", &xvar1);  // first node
	j=xvar1;
	fscanf(ar,"%lf\n", &xvar1);  // second node
	q=xvar1;
	fscanf(ar,"%lf\n", &xvar1);  // adjacency of both
	con=xvar1;	

// ATTENTION: Ignores multiple edges:
	if (xvar1>0) con=1;
	else con=0;

	epsi[j][q]=con;  // adjacency matrix
	epsi[q][j]=con;

	if (j>n_max) n_max=j;
	if (q>n_max) n_max=q;

	if ((feof(ar)!=0)) break;
	
    }	
   
    N=n_max;

    for (i=0;i<n_max;i++)
    {
	for (j=0;j<n_max;j++)
	{
	    if (epsi[i][j]==1) k[i]++;
	}
    }

    k_med=0.0;
    k2_med=0.0;
    k3_med=0.0;
    for (i=0;i<n_max;i++)
    {
	q=k[i];
	n[q]++;
	k_med+=k[i];
	cosa=k[i]*k[i];
//	k2_med+=k[i]*k[i]; 
	k2_med+=cosa;
//	k3_med+=k[i]*k[i]*k[i];
 	k3_med+=cosa*k[i];
    }
    k_med=k_med/(1.0*n_max);
    k2_med=k2_med/(1.0*n_max);
    k3_med=k3_med/(1.0*n_max);
    for (i=0;i<n_max;i++) p[i]=n[i]/(1.0*n_max);

    fprintf(deg_dis,"# k\tp(k)\n");
    for (i=0;i<n_max;i++)
    {
	fprintf(deg_dis,"%i\t%f\t%i\n",i,p[i],n[i]);  
    }	

    disp=dispersion();
    printf("# <k>= %f\tN= %i\t\tstd.dev[p(k)]= %f\n",k_med,n_max,disp);
    printf("C_rand= %f\tl_rand= %f\n",C_rand,L_rand);
    printf("C_conf= %f\tl_conf= %f\n",C_conf,L_conf);

    neighbours();
    printf("k_nei_med= %f\t k_nei_med_conf= %f\n",nei_med,nei_conf);
    for (i=0; i<N; i++)
    {
	fprintf(mix,"%i\t%i\t%f\t%f\t # i, k(i), k_nn(i), <k_nn(k)>\n",i,k[i],nei[i],nei_m[i]);
    }

    clustering_coeff();
    printf("C= %f\tsig(C)= %f\n",Cm,sigC);

    short_path();
    printf("l= %f\tsig(l)= %f\n",Lm,sigL);

    fprintf(sw,"# File: %s\n",argv[ARCHIVO1]);

    for (i=0; i<N; i++)
    {
	kt[i]=k[i];
	for (j=0; j<N; j++) adj[i][j]=epsi[i][j];
    }

	
    fprintf(sw,"<k>= %f\tN= %i\t\tstd.dev.[p(k)]= %f\n",k_med,n_max,disp);
    fprintf(sw,"C= %f\tstd.dev.= %f\tC_conf= %f\tC_rand= %f\n",Cm,sigC,C_conf,C_rand);
    fprintf(sw,"l= %f\tstd.dev.= %f\tl_conf= %f\tl_rand= %f\n",Lm,sigL,L_conf,L_rand);
    fprintf(sw,"<knn>= %f\t<knn>_conf= %f\n",nei_med,nei_conf);
    fprintf(sw,"r= %f\n",r);

    fprintf(mix,"# %f\t%f # -> <g>, g_conf\n",nei_med,nei_conf);

    free_imatrix(epsi,0,Nbig-1,0,Nbig-1);    
    free_ivector(k,0,Nbig-1);
    free_ivector(n,0,Nbig-1);
    free_dvector(p,0,Nbig-1);
    free_dvector(nei,0,Nbig-1);
    free_dvector(nei_m,0,Nbig-1);
    free_imatrix(adj,0,Nbig-1,0,Nbig-1);    
    free_ivector(kt,0,Nbig-1);
    free_dmatrix(n_bas,0,Nbig-1,0,Nbig-1);

    fclose(ar);
    fclose(deg_dis);
    fclose(sw);

   return 1;
}



void short_path()
{
    double xr, ia, a1, a2;
    int i,j,ka;    
    double **I;
    double **A;

    int DIM,DIM1,DIM2,DIM3,NNEU;

    Dim=n_max;
    DIM1=Dim;
    
    A=dmatrix(0,DIM1,0,DIM1);
    I=dmatrix(0,DIM1,0,DIM1);

    for (i=0;i<DIM1;i++){
	for (j=0;j<DIM1;j++){
	    
	    if (epsi[i][j]==1) A[i][j]=1.0;
	    else A[i][j]=0.0;
	}
    }
 

   
    /* sequential floyd's algorithm */
    
    
    for (i=0;i<DIM1;i++){
	for (j=0;j<DIM1;j++){
	    
	    I[i][j]=10000000000000.;
	    
	    ia=A[i][j];
	    
	    if (j==i)I[i][j]=0.;
	    else if (ia>0.5) I[i][j]=1.0;
	    
	}
    }
    
    
    for (ka=0;ka<DIM1;ka++){
	for (i=0;i<DIM1;i++){
	    for (j=0;j<DIM1;j++){
		
		a1=I[i][ka]+I[ka][j];
		a2=I[i][j];
		I[i][j]=MIN(a2,a1);
	    }
	}	 
    }


// Get mean short-path and variance thereof:
   

    Lm=0.0;
    L2m=0.0;
    sigL=0.0;
    count=0;
    count_disc=0;
    for (i=0; i<DIM1; i++) for (j=0;j<DIM1;j++)
    {

	if (I[i][j]<=2*N && i!=j)  // not to count the "infinite" paths of disconnected bits
	{
	    Lm=Lm+I[i][j];
	    L2m=L2m+I[i][j]*I[i][j];
	    count++;
	}    

	if (I[i][j]>2*N) count_disc++;

    }

    Lm=Lm/(1.0*count);
    L2m=L2m/(1.0*count);  
    sigL=sqrt(ABSO(L2m-Lm*Lm ));

    mean_disc+=count_disc;

    free_dmatrix(A,0,DIM1,0,DIM1);
    free_dmatrix(I,0,DIM1,0,DIM1);
}




void clustering_coeff()
{    
// Get los[][]:

    int i,j,l,neu;

    for (i=0; i<n_max; i++) k[i]=k[i];

    if (got_los==0)
    {
	for (i=0; i<n_max; i++) for (q=0; q<n_max; q++) los[i][q]=-1;     
	for (i=0; i<n_max; i++) 
	{
	    q=0;
	    for (j=0; j<n_max; j++)
	    {
		if (epsi[i][j]==1)
		{		
		    los[i][q]=j;
		    q++;
		}    
	    }
	}
    }
    Dim=n_max;
   
    Cm=0.0;
    C2m=0.0;
    sigC=0.0;
    for (i=0; i<Dim; i++)
    {
	C=0.0;
	tria=0;

	
	for (q1=0; q1<k[i]; q1++)  // i's neighbours = j
	{
	    j=los[i][q1];	    
	    if (j>=0)
	    {
		for (q2=0; q2<k[j]; q2++)   // j's neighbs = l
		{
		    l=los[j][q2];
		    if (l>-1)  // l's neighbs = neu
		    {
			for (q3=0; q3<Dim; q3++)
			{
			    neu=los[l][q3]; 
			    if (neu==i) tria++;// if neu==i, we have a triangle
			}
		    }
		}
	    }
	}
	
	
	if (k[i]>1) C=1.0*tria/(1.0*k[i]*(k[i]-1.0));  // C(i)
	else C=0.0;
	Cm=Cm+C;
	C2m=C2m+C*C;
    }
    Cm=Cm/(1.0*Dim);
    C2m=C2m/(1.0*Dim);
    sigC=sqrt(ABSO(C2m-Cm*Cm));
}

void neighbours()
{
// With los:

    if (got_los==0)
    {
	for (i=0; i<n_max; i++) for (q=0; q<n_max; q++) los[i][q]=-1;     
	for (i=0; i<n_max; i++) 
	{
	    q=0;
	    for (j=0; j<n_max; j++)
	    {
		if (epsi[i][j]==1)
		{		
		    los[i][q]=j;
		    q++;
		}    
	    }
	}
    }

    nei_med=0.0;
    for (i=0; i<n_max; i++)
    {
	nei[i]=0.0;	
	for (q=0; q<n_max; q++)	// for (q=0; q<k[i]; q++) // should be the same
	{
	    j=los[i][q];
	    if (j>=0) nei[i]+=k[j];
	}

	if (k[i]>0) 
	{
	    nei[i]=nei[i]/(1.*k[i]);
	    nei_med=nei_med+nei[i];
	}
    }

    for (j=0; j<N; j++) nei_m[i]=0.0;
    for (i=0; i<N; i++)
    {
	j=k[i];
	nei_m[j]+=nei[i];
    }

    for (j=0; j<n_max; j++)
    {
	if (n[j]>0) nei_m[j]=nei_m[j]/(1.*n[j]);  // <g(i|k)> = g(k)
    }


    nei_med=nei_med/(1.*n_max);

    nei_k2=0.0;
    for (i=0; i<N; i++)
    {
	nei_k2+=k[i]*k[i]*nei[i];
    }
    nei_k2=nei_k2/(1.*N);

    r_conf=(k_med*k2_med*nei_conf-k2_med*k2_med)/(k_med*k3_med-k2_med*k2_med);


//  Pearson's coefficient (ki*kj, averaged over links):
    pear=0.0;
    for (i=0; i<N; i++) for (j=0; j<N; j++)
    {
	if (epsi[i][j]==1) pear+=k[i]*k[j];
    }
    pear=pear/(k_med*N);
    r=(pear-k2_med/k_med*k2_med/k_med)/(k3_med/k_med-k2_med/k_med*k2_med/k_med);

    printf("r= %f\n",r);

}


double dispersion()
{
    double k2,var,vari,di,km;
    km=0.0;
    k2=0.0;
    var=0.0;

    for (i=0; i<N; i++)
    {
	km=km+k[i];
	k2=k2+k[i]*k[i];
    }
    km=km/(1.0*N);
    k2=k2/(1.0*N);
    var=k2-km*km;
    di=sqrt(var);
    C_reg=km/(1.0*N);

// Configuration ensemble:
    vari=var;
    cosa=vari/km/km+1.-1./km;
    C_conf=cosa*cosa*km/(1.*N);
    L_conf=1.+2.*(log(N)-log(km))/(log(N)+log(km)+log(C_conf));
    nei_conf=k2/km;

// Erdös-Rényi ensemble (Poisson distrib => var=<k>):
    vari=km;
    cosa=vari/km/km+1.-1./km;
    C_rand=cosa*cosa*km/(1.*N);
    L_rand=1.+2.*(log(N)-log(km))/(log(N)+log(km)+log(C_rand));

    return di;
}  // end short_path()




float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        float **m;
 
        /* allocate pointers to rows */
        m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m += NR_END;
        m -= nrl;
 
        /* allocate rows and set pointers to them */
        m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
        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 nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
        fprintf(stderr,"Numerical Recipes run-time error...\n");
        fprintf(stderr,"%s\n",error_text);
        fprintf(stderr,"...now exiting to system...\n");
        exit(1);
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
        double *v;
 
        v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
        if (!v) nrerror("allocation failure in dvector()");
        return v-nl+NR_END;
}
float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
        float *v;
 
        v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
        if (!v) nrerror("allocation failure in vector()");
        return v-nl+NR_END;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        double **m;
 
        /* allocate pointers to rows */
        m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m += NR_END;
        m -= nrl;
 
        /* allocate rows and set pointers to them */
        m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
        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;
}
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_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

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_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-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));
}