socialNet.h File Reference

#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <fstream>
#include <stack>
#include <string.h>
#include <math.h>

Functions
void	print_int_array (int **array, int array_size)
void	print_double_array (double **array, int array_size)
int **	initialise_int_array (int array_size)
int **	initialise_int_square_array (int x_size, int y_size)
double **	initialise_double_array (int array_size)
void	initialise_adjacency_array_WS (int **adjacency, int n, int k, double p)
void	initialise_adjacency_array_random (int **adjacency, int array_size)
void	initialise_adjacency_array (int **adjacency, int array_size, int k, double p)
void	initialise_weight_array_WS (double **weight, int array_size)
void	initialise_weight_array_random (double **weight, int array_size)
void	initialise_weight_array_ingroups (double **weight, int array_size, int numberOfGroups, double nRewiring, double threshold, double seed)
void	rewire_WS (int **adjacency, double **weight, int n, int k, double p)
void	rewire_WS2 (int **adjacency, double **weight, int n, int k, double p)
void	initialise_weight_array (double **weight, int array_size)
void	initialise_weight_array2 (double **weight, int array_size)
void	rewire (int **adjacency, double **weight, int array_size, int k, double p)
void	calculate_shortest_paths (double **shortest_path, double **shortest_path_ideal, int **adjacency, double **weight, int array_size)
double	calculate_energy (double **shortest_path, int array_size)
double	calculate_energy_local (int node, int **adjacency, double **weight, int array_size)
double	calculate_E_glob (double **shortest_path, double **shortest_path_ideal, int array_size)
double	calculate_E_loc (int **adjacency, double **weight, int array_size)
double	gamma (double l)
double	calculate_cost (int **adjacency, double **weight, int array_size)
int	load_double_array (double **array, int array_size, const char *filename)
int	load_double_array (double **array, const char *filename)
void	generate_adjacency (double **weightMat, int **adjacencyMat, double threshold, int array_size)
void	print_array_int (int **array, int array_size)
void	print_array_double (double **array, int array_size)
void	assign_distance_not_working (int current, int previous, double *distance, int d, int target, int *assigned, int **adjacency, int **pred, int array_size)
void	assign_distance (double *distance, int d, int *leaf, int **adjacency, int **pred, int *predNum, int array_size, std::stack< int > &stoiva, int *sigma)
void	calculate_edge_b (double **edge_bw, double *distance, int **pred, int **adjacency, int *predNum, int array_size, std::stack< int > &stoiva, int *sigma, double *delta, int current, int *leaf)
void	calculate_node_b (double *betw, double *distance, int **pred, int *predNum, int array_size, std::stack< int > &stoiva, int *sigma, double *delta, int current)
void	calculate_edge_betweenness (double **result_betw, int **adjacency, int array_size)
void	calculate_node_betweenness (double *result_betw, int **adjacency, int array_size)
int	assignToAGroup (int current, int currentGroup, int **adjacency, int **groups, int *numberOfMembers, int array_size, bool *assigned)
int	getGroups (int **adjacency, int **groups, int *numberOfMembers, int array_size)
bool	isInGroup (int node, int *group, int numberOfMembers)
void	splitNetwork (int **adjacency, double *betw, int array_size)
double	splitNetwork_Threshold (int **adjacency, double *betw, int array_size, double modThreshold)
double	splitNetwork_edge_Threshold (int **adjacency, int **tempadjacency, double **betw, int array_size, double modThreshold)
bool	areInTheSameGroup (int node1, int node2, int **groups, int numberOfGroups, int *numberOfMembers)
void	printGroups (int numberOfGroups, int **groups, int *numberOfMembers, int array_size)
bool	readAdjacency (string inputFile, double **weightMat, int **adjacencyMat, int numHosts, double threshold)

---

Function Documentation

bool areInTheSameGroup	(	int	node1,
		int	node2,
		int **	groups,
		int	numberOfGroups,
		int *	numberOfMembers
	)

                                                                                                      {
        
        bool result=false;
        for (int k=0;k<numberOfGroups;k++) {
                
                if (isInGroup(node1,groups[k],numberOfMembers[k]))
                        if (isInGroup(node2,groups[k],numberOfMembers[k])) {
                                result=true;
                                //break;
                        }
                }
        
        return result;
}

void assign_distance	(	double *	distance,
		int	d,
		int *	leaf,
		int **	adjacency,
		int **	pred,
		int *	predNum,
		int	array_size,
		std::stack< int > &	stoiva,
		int *	sigma
	)

void assign_distance_not_working	(	int	current,
		int	previous,
		double *	distance,
		int	d,
		int	target,
		int *	assigned,
		int **	adjacency,
		int **	pred,
		int	array_size
	)

                                                                                                                                                           {
                        if (current!=target) {
                                for (int k=0;k<array_size; k++)
                                if ((k!=current)&&(k!=previous)) {                                                                      
                                        if (adjacency[current][k]==1) {                                 
                                                if (assigned[k]==1) {
                                                        distance[k]=d+1;
                                                        assigned[k]=0;
                                                }       
                                                else if (distance[k]==d) {
                                                        //cout<<"Assign distance "<<d<<"\n";
                                                        //assign_distance(k,current,distance,d,target,assigned,adjacency,pred,array_size); 
                                                }
                                        }
                                }//if ((k!=current))
                }
}       

int assignToAGroup	(	int	current,
		int	currentGroup,
		int **	adjacency,
		int **	groups,
		int *	numberOfMembers,
		int	array_size,
		bool *	assigned
	)

                                                                                                                                      {
        
        for (int k=0;k<array_size;k++) {                
                if (adjacency[current][k]==1) {
                        if (current!=k) {
                                if (assigned[k]==false) {
                                        numberOfMembers[currentGroup-1]=numberOfMembers[currentGroup-1]+1;
                                        groups[currentGroup-1][numberOfMembers[currentGroup-1]-1]=k+1;
                                        assigned[k]=true;
                                        assignToAGroup(k,currentGroup,adjacency, groups, numberOfMembers, array_size, assigned);        
                                }
                        }
                }
        }
}

double calculate_cost	(	int **	adjacency,
		double **	weight,
		int	array_size
	)

{       double numerator = 0.0;
        double denominator = 0.0;

        for (int i = 0; i < array_size; i++)
          for (int j = 0; j < array_size; j++)
          {       if (i != j)
                        { denominator += gamma(weight[i][j]);
                          if (adjacency[i][j] != 0)
                        numerator += gamma(weight[i][j]);
                        }
          }

          return (numerator/denominator);
}

double calculate_E_glob	(	double **	shortest_path,
		double **	shortest_path_ideal,
		int	array_size
	)

{
        return (calculate_energy(shortest_path, array_size)/calculate_energy(shortest_path_ideal, array_size));
}       

double calculate_E_loc	(	int **	adjacency,
		double **	weight,
		int	array_size
	)

{       double sum = 0.0;

    for (int i = 0; i < array_size; i++)
                sum += calculate_energy_local(i, adjacency, weight, array_size);

        return (sum/array_size);
}       

void calculate_edge_b	(	double **	edge_bw,
		double *	distance,
		int **	pred,
		int **	adjacency,
		int *	predNum,
		int	array_size,
		std::stack< int > &	stoiva,
		int *	sigma,
		double *	delta,
		int	current,
		int *	leaf
	)

void calculate_edge_betweenness	(	double **	result_betw,
		int **	adjacency,
		int	array_size
	)

{       
        stack<int> stoiva;
        int* b=new int[array_size];
        int* leaf= new int[array_size];
        int* predNum= new int[array_size];
        int* sigma= new int[array_size];
        double* distance= new double[array_size];
        double* delta=new double[array_size];
        double **edge_bw=initialise_double_array(array_size);
        
        //Initialize edge_bw
        for(int k=0;k<array_size; k++) {
                for (int r=0;r<array_size;r++) {                
                        edge_bw[k][r]=0; 
                        result_betw[k][r]=0;
                }
        }
        
        int d;
        for (int i=0; i<array_size;i++) { 
                leaf[i]=1;
                distance[i]=15;
                //cout<<distance[i]<<"\n";
        
        }
                
        double **dist= initialise_double_array(array_size);
        int **pred= initialise_int_array(array_size);
                                
        for (int i=0; i<array_size; i++) {
        
                //Initialize predecessor array
                for(int l=0; l<array_size; l++)
                        for(int j=0; j<array_size; j++)
                                pred[l][j]=-1;
                
                // Initialize other arrays      
                for (int s=0; s<array_size;s++) { 
                        leaf[s]=1;
                        distance[s]=16;
                        predNum[s]=0;
                        sigma[s]=0;
                        delta[s]=0;
                }
                //for (int j=0; j<array_size; j++) {
                leaf[i]=0;
                distance[i]=0;
                sigma[i]=1;
                for (int d=0;d<15;d++) {
                        // For distances from 0 to 15 assign distances
                        assign_distance(distance,d,leaf,adjacency,pred,predNum,array_size,stoiva,sigma);        
                }                               
        
                calculate_edge_b (edge_bw, distance, pred, adjacency, predNum,array_size,stoiva,sigma,delta,i,leaf);
                
                //Add results to betweenness and reset temporal counter edge_bw
                for (int z=0;z<array_size;z++) {
                        for (int p=0;p<array_size;p++) {
                                result_betw[z][p]+=(edge_bw[z][p]/2);
                                edge_bw[z][p]=0;
                        }
                }       
                                
        }//end for (int i<0; i<array_size; i++) 

        //delete data structures
        for (int i = 0; i < array_size; i++) {
                delete edge_bw[i];
                delete pred[i];
                delete dist[i];
        }
        delete[] delta;
        delete[] edge_bw;
        delete[] pred;
        delete[] dist;
        delete[] b;
        delete[] distance;
        delete[] leaf;
        delete[] predNum;
}

double calculate_energy	(	double **	shortest_path,
		int	array_size
	)

{
  double sum = 0.0;

  for (int i = 0; i < array_size; i++)
          for (int j = 0; j < array_size; j++)
                  if (i != j && shortest_path[i][j] != infinity)
                    sum += 1.0/shortest_path[i][j];

  // Whilst, strictly, we need this to calculate the energy
  // because it is normalised, we can remove it.
  //
  // sum /= array_size * (array_size - 1);
  return(sum);
}

double calculate_energy_local	(	int	node,
		int **	adjacency,
		double **	weight,
		int	array_size
	)

{
        double res = 0.0;
        int        idx_array_size = 0;

        static int    *idx_array                        = new int[array_size];
        static int    **temp_adjacency                  = initialise_int_array(array_size);
        static double **temp_weight                     = initialise_double_array(array_size);
        static double **temp_shortest_path              = initialise_double_array(array_size);
        static double **temp_shortest_path_ideal        = initialise_double_array(array_size);

        // Calculate local energy for the graph formed by this node's immediate neighbours. Note that
        // this graph specifically excludes the node itself so as to measure fault tolerance -- how
        // efficient communication amongst the node's neighbours is if the node fails.
        //

        // work out the number of neighbours and their positions
        // the idx array can then be used to indirect through for further calcaulations
        //
        for (int i = 0; i < array_size; i++)
                if (adjacency[node][i] != 0 && i != node)               // i is an immediate neighbour of node
                        idx_array[idx_array_size++] = i;

        // Fill in temporary arrays
        // Only need to fill in the entries we're interested in, which are in idx array.
        //
        for (int i = 0; i < idx_array_size; i++)
                for (int j = 0; j < idx_array_size; j++)
                {       int idx1 = idx_array[i];
                        int idx2 = idx_array[j];

                        temp_weight[i][j]        = weight[idx1][idx2];
                        temp_adjacency[i][j] = adjacency[idx1][idx2];
                }

        // Using only these nodes calculate shortest paths
        //
        calculate_shortest_paths(temp_shortest_path, temp_shortest_path_ideal, temp_adjacency, temp_weight, idx_array_size);

        // And work out the normalised global energy for this restricted set
        //
        double energy           = calculate_energy(temp_shortest_path, idx_array_size) ;
        double energy_ideal = calculate_energy(temp_shortest_path_ideal, idx_array_size);

        if (energy_ideal != 0)
                return (energy/energy_ideal);
        else
                return (0.0);
}

void calculate_node_b	(	double *	betw,
		double *	distance,
		int **	pred,
		int *	predNum,
		int	array_size,
		std::stack< int > &	stoiva,
		int *	sigma,
		double *	delta,
		int	current
	)

void calculate_node_betweenness	(	double *	result_betw,
		int **	adjacency,
		int	array_size
	)

                                                                                      {

        int* b=new int[array_size];
        int* leaf= new int[array_size];
        int* predNum= new int[array_size];
        double* distance= new double[array_size];
        double* betw= new double[array_size];
        double* current_betw=new double[array_size];

        stack<int> stoiva;
        int* sigma= new int[array_size];
        double* delta=new double[array_size];
        
        int d;
        for (int i=0; i<array_size;i++) { 
                leaf[i]=1;
                distance[i]=15;
                betw[i]=0;
                current_betw[i]=0;
        }
                
        double **dist= initialise_double_array(array_size);
        int **pred= initialise_int_array(array_size);
                                
        for (int i=0; i<array_size; i++) {
        
                        //Initialize predecessor array
                        for(int l=0; l<array_size; l++)
                                for(int j=0; j<array_size; j++)
                                        pred[l][j]=-1;
                        
                        // Initialize other arrays      
                        for (int s=0; s<array_size;s++) { 
                                                leaf[s]=1;
                                                distance[s]=16;
                                                predNum[s]=0;
                                                sigma[s]=0;
                                                delta[s]=0;
                                                //betw[s]=0;
                        }

                        leaf[i]=0;
                        distance[i]=0;
                        sigma[i]=1;
                        for (int d=0;d<15;d++) {
                                // For distances from 0 to 15 assign distances
                                assign_distance(distance,d,leaf,adjacency,pred,predNum,array_size,stoiva,sigma);        
                        }                       
                
                        calculate_node_b (betw,distance, pred, predNum,array_size,stoiva,sigma,delta,i);                                
         }

        for (int z=0;z<array_size;z++)          
                result_betw[z]=betw[z]/2;

        //delete data structures
        for (int i = 0; i < array_size; i++) {
                delete pred[i];
                delete dist[i];
        }
        delete[] pred;
        delete[] dist;
        delete[] delta;
        delete[] b;
        delete[] current_betw;
        delete[] distance;
        delete[] leaf;
        delete[] predNum;
        delete[] betw;
}

void calculate_shortest_paths	(	double **	shortest_path,
		double **	shortest_path_ideal,
		int **	adjacency,
		double **	weight,
		int	array_size
	)

{       int i, j, k;

        // ??? MIGHT BE BETTER TO USE DIJKSTRA HERE SINCE IT IS PROBABLE THAT THE NUMBER OF
        // EDGES << N^2
        
        // Use Floyd's algorithm
        //
        for (i = 0; i < array_size; i++)
                for (j = 0; j < array_size; j++)
                {       shortest_path[i][j]       = (adjacency[i][j] != 0 ? weight[i][j] : infinity);
                        shortest_path_ideal[i][j] = weight[i][j];
                }       
                
        for (k = 0; k < array_size; k++)
                for (i = 0; i < array_size;i++)
                {       if (i == k)
                                continue;

                        for (j = 0; j < array_size; j++)
                        {   if (j == i || j == k)
                                  continue;


                          if (shortest_path[i][k] != infinity && shortest_path[k][j] != infinity)
                          { if (shortest_path[i][j] == infinity || (shortest_path[i][k] + shortest_path[k][j] < shortest_path[i][j]))
                                        shortest_path[i][j] = shortest_path[i][k] + shortest_path[k][j];

                                if (shortest_path_ideal[i][j] == infinity || (shortest_path_ideal[i][k] + shortest_path_ideal[k][j] < shortest_path_ideal[i][j]))
                                        shortest_path_ideal[i][j] = shortest_path_ideal[i][k] + shortest_path_ideal[k][j];
                          }

                        }

                }
}

double gamma

(

double

l

)

{ return l; }

void generate_adjacency	(	double **	weightMat,
		int **	adjacencyMat,
		double	threshold,
		int	array_size
	)

                                                                                                   {
        
        for (int i=0; i<array_size; i++)  
                for (int j=0; j<array_size; j++) {
                        if (weightMat[i][j]>threshold)
                                adjacencyMat[i][j]=1; 
                        else
                                adjacencyMat[i][j]=0;
                }
}

int getGroups	(	int **	adjacency,
		int **	groups,
		int *	numberOfMembers,
		int	array_size
	)

                                                                                    {
        
        int numberOfGroups=0;
        bool assigned[array_size];
        
        for (int i=0; i<array_size; i++) {
                assigned[i]=false;      
        } 

        for (int i=0; i<array_size; i++) {
                
                if (assigned[i]==false) {
                        numberOfGroups++;
                        numberOfMembers[numberOfGroups-1]=numberOfMembers[numberOfGroups-1]+1;
                        groups[numberOfGroups-1][numberOfMembers[numberOfGroups-1]-1]=i+1;
                        assigned[i]=true;
                        assignToAGroup(i,numberOfGroups,adjacency, groups, numberOfMembers, array_size, assigned);      
                }                       
        }
        return numberOfGroups;
}

void initialise_adjacency_array	(	int **	adjacency,
		int	array_size,
		int	k,
		double	p
	)

{       
        // initialise_adjacency_array_random(adjacency, array_size);
        initialise_adjacency_array_WS(adjacency, array_size, k, p);

}

void initialise_adjacency_array_random	(	int **	adjacency,
		int	array_size
	)

{       for (int i = 0; i < array_size; i++)
                for (int j = 0; j < array_size; j++)
                        if (i == j)
                                adjacency[i][j] = 1;
                        else if (i < j)
                        { adjacency[i][j] = (((double) rand())/RAND_MAX < 0.5);
                          adjacency[j][i] = adjacency[i][j];
                        }
}

void initialise_adjacency_array_WS	(	int **	adjacency,
		int	n,
		int	k,
		double	p
	)

{       int i, j;

        for (i = 0; i < n; i++)
                for (j = 0; j < n; j++)
                        if (i == j)
                                adjacency[i][j] = 1;
                        else    
                                adjacency[i][j] = 0;

        // Create the regular lattice -- note that this is a circular lattice
        // If we iterate over the whole of the arrayand put in the k/2 forward links, then,
        // because this is symmetric, we will end up with the nearest k links in.
        //
        for (i = 0; i < n; i++)
                for (j = 1; j <= k/2; j++)
                {   int idx = (i+j) % n;                // Make it circular

                        adjacency[i][idx] = 1;
                        adjacency[idx][i] = 1;
                }

}

double** initialise_double_array

(

int

array_size

)

{       double **result = new double * [array_size];

        for (int i = 0; i < array_size; i++) {
                result[i] = new double [array_size];
                for(int j = 0; j < array_size; j++)
                        result[i][j]=0;
        }
        
        return result;
}

int** initialise_int_array

(

int

array_size

)

                                           {
        
        int **result = new int * [array_size];
        
        for (int i = 0; i < array_size; i++)
                result[i] = new int [array_size];

        return result;
}

int** initialise_int_square_array	(	int	x_size,
		int	y_size
	)

                                                          {

        int **result = new int * [y_size];

        for (int i = 0; i < y_size; i++)
                result[i] = new int [x_size];

        return result;
}

void initialise_weight_array	(	double **	weight,
		int	array_size
	)

{ 
        // initialise_weight_array_random(weight, array_size);
        initialise_weight_array_WS(weight, array_size);
}

void initialise_weight_array2	(	double **	weight,
		int	array_size
	)

                                                                {
        int seed = 5;
        srand(seed);
                for (int i=0;i<array_size;i++)
                        for (int j=0;j<array_size;j++) {
                                if (i==j) weight[i][j]=1;
                                if (i<j){
                                        weight[i][j]=rand()/(RAND_MAX+1.0);
                                        weight[j][i]=weight[i][j];
                                }
                        }
}

void initialise_weight_array_ingroups	(	double **	weight,
		int	array_size,
		int	numberOfGroups,
		double	nRewiring,
		double	threshold,
		double	seed
	)

                                                                                                                                                {
        
        srand((int)seed);
        int groupId=0;

        int **groups=initialise_int_array(array_size);
        int numberOfMembers[numberOfGroups];
        
        for (int i=0;i<numberOfGroups;i++) {
                numberOfMembers[i]=0;
        }
        
        //Initial assignment of nodes to groups
        for (int i=0;i<array_size;i++) {
                do {
                        groupId=rand()%numberOfGroups;
                } while(numberOfMembers[groupId]!=floor(i/numberOfGroups));
                //int groupId=i%numberOfGroups;
                //cout<<"\n";
                //cout<<"Group id of host "<<i<<" is "<<groupId<<".\n";
                groups[groupId][numberOfMembers[groupId]]=i+1;
                numberOfMembers[groupId]+=1;    
        }
        
        for (int i=0;i<array_size;i++)
                for (int j=0;j<=i;j++) {
                         if (areInTheSameGroup (i+1,j+1,groups,numberOfGroups,numberOfMembers)==true) {
                                        if(i==j) {// if it's a diagonal element
                                                weight[i][j]=1;
                                        }
                                        else if ((rand()/(RAND_MAX+1.0))<probRewiring) {        //rewiring
                                                bool found=false;
                                                for (int z=0;z<array_size;z++) {
                                                        if ((areInTheSameGroup (i+1,z+1,groups,numberOfGroups,numberOfMembers)==false)&&(weight[i][z]==0)&&(found==false)) {
                                                                weight[i][z]=1.0-(rand()*threshold/(RAND_MAX+1.0));
                                                                weight[z][i]=weight[i][z];
                                                                found=true;
                                                        }
                                                }
                                                weight[i][j]=(rand()*threshold/(RAND_MAX+1.0));
                                                weight[j][i]=weight[i][j];
                                        }
                                        else {//no rewiring 
                                                weight[i][j]=1.0-(rand()*threshold/(RAND_MAX+1.0));
                                                weight[j][i]=weight[i][j];
                                        }
                         } else if(weight[i][j]==0){//the hosts are not in the same cluster
                                weight[i][j]=(rand()*threshold/(RAND_MAX+1.0));
                                weight[j][i]=weight[i][j];
                         }                      
                }
        
        //deletion of the variables
        for (int i=0;i<array_size;i++) 
                delete[] groups[i];
        delete[] groups;
}

void initialise_weight_array_random	(	double **	weight,
		int	array_size
	)

{       for (int i = 0; i < array_size; i++)
                for (int j = 0; j < array_size; j++)
                        if (i == j)
                                weight[i][j] = 0;
                        else if (i < j)
                        { weight[i][j] = ((double) rand())/RAND_MAX;

                          // Impose arbitrary lower limit to ensure that reciprocals don't get
                          // unmanageably large
                          //
                          while (weight[i][j] < 0.01)
                                weight[i][j] = ((double) rand())/RAND_MAX;

                          weight[j][i] = weight[i][j];
                        }
}

void initialise_weight_array_WS	(	double **	weight,
		int	array_size
	)

{
        for (int i = 0; i < array_size; i++)
                for (int j = 0; j < array_size; j++)
                        if (i == j)
                                weight[i][j] = 0;
                        else
                                weight[i][j] = 1;
}

bool isInGroup	(	int	node,
		int *	group,
		int	numberOfMembers
	)

                                                          {
        
        bool result=false;
        
        for (int k=0;k<numberOfMembers;k++)     
                if (group[k]==node) {
                        result=true;
                        //break;
                }
        return result;
}

int load_double_array	(	double **	array,
		const char *	filename
	)

                                                            {
        
        //no checks
        int i,j=0;
        int r;
        double value;
        //int * valueint;
        //FILE *fopen (const char *filename, const char *mode)
        FILE* file1=fopen (filename,"r");
        
        int array_size;
        //Segmentation fault
        r=fscanf(file1,"%d",&array_size);
        
        //<<"Loading matrix from file "<<filename<<"...\n";
        
        for (int i=0; i<array_size; i++) {
                for (int j=0; j<array_size; j++) {
                        r=fscanf (file1, "%lf", &value);
                        array[i][j]=value;
                }
        }
        
        fclose (file1);
        return array_size;
}

int load_double_array	(	double **	array,
		int	array_size,
		const char *	filename
	)

                                                                           {
        
        //no checks
        int i,j=0;
        int r;
        double value;
        //int * valueint;
        //FILE *fopen (const char *filename, const char *mode)
        FILE* file1=fopen (filename,"r");
        
        //cout<<"Loading matrix from file "<<filename<<"...\n";
        
        for (int i=0; i<array_size; i++) {
                for (int j=0; j<array_size; j++) {
                        r=fscanf (file1, "%lf", &value);
                        array[i][j]=value;
                }
        }
        
        fclose (file1);
}

void print_array_double	(	double **	array,
		int	array_size
	)

                                                         {
        for (int i=0; i<array_size; i++) {
                for (int j=0; j<array_size; j++) { 
                        printf ("%.2f  ",array[i][j]);
                }
                cout<<"\n";
        }
}

void print_array_int	(	int **	array,
		int	array_size
	)

                                                   {
        
        for (int i=0; i<array_size; i++) {
                for (int j=0; j<array_size; j++) { 
                        printf ("%i  ",array[i][j]);
                }
                cout<<"\n";
        }
}

void print_double_array	(	double **	array,
		int	array_size
	)

{       for (int i = 0; i < array_size; i++)
        { cout << array[i][0];
          for (int j = 1; j < array_size; j++)
                cout << ", " << array[i][j];
          cout << "\n";
        }
        cout << "\n"; fout.flush();
}

void print_int_array	(	int **	array,
		int	array_size
	)

{       for (int i = 0; i < array_size; i++)
        { cout << array[i][0];
          for (int j = 1; j < array_size; j++)
                cout << ", " << array[i][j];
          cout << "\n";
        }
        cout << "\n"; fout.flush();
}

void printGroups	(	int	numberOfGroups,
		int **	groups,
		int *	numberOfMembers,
		int	array_size
	)

                                                                                       {
        
        for (int i=0;i<numberOfGroups;i++) {
                cout<<"\n";
                cout<<"The members of group "<<i+1<<" are: ";
                for (int j=0;j<numberOfMembers[i];j++)
                        cout<<groups[i][j]<<" ";
                cout<<"\n";
        }       
}

bool readAdjacency	(	string	inputFile,
		double **	weightMat,
		int **	adjacencyMat,
		int	numHosts,
		double	threshold
	)

                                                                                                             {

        int adj;
        int linecnt=0;
        int cnt=0;
        string line;
        
        //Input file
        ifstream infile(inputFile.c_str());

        //Opening file
        if (!infile) {
                cout << "There was a problem opening file " << inputFile << " for reading." << endl;
                return false;
        }
        
        //Count line number
        while (getline(infile,line)) {
                linecnt++;
        }
        
        //Compare line number with number of hosts
        if(linecnt!=numHosts) {
                cout << "Number of hosts in input file does not comply with given input" << endl;
                return false;
        }
        
        //Move pointer to the beggining of file
        infile.clear();
        infile.seekg(0,ios::beg);
        
        //Read values
        for (int i=0; i<numHosts; i++) {
                for (int j=0; j<numHosts; j++) {
                        if (infile >> adj) {
                                if(adj >= 1) {
                                        adjacencyMat[i][j]=1;
                                        if(weightMat[i][j]==0) {
                                                weightMat[i][j]=1.0-rand()*threshold/(RAND_MAX+1.0);
                                                weightMat[j][i]=weightMat[i][j];
                                        }
                                }
                                else {
                                                adjacencyMat[i][j]=0;
                                                weightMat[i][j]=0;
                                }
                        }
                        else {
                                cout << "Problem when reading adjacency matrix" << endl;
                                return false;
                        }
                }
        }

        return true;
}

void rewire	(	int **	adjacency,
		double **	weight,
		int	array_size,
		int	k,
		double	p
	)

{
        rewire_WS2(adjacency, weight, array_size, k, p);
}

void rewire_WS	(	int **	adjacency,
		double **	weight,
		int	n,
		int	k,
		double	p
	)

{
        //
        // OK, now rewire.
        // Start with the edges closest to each node (in a clockwise direction) and rewire with
        // probability p to give a new endpoint, provided that (i) it is not wired to us and
        // (ii) there is not already a link between us and that node. Then move out to the edges
        // connecting us to the next furthest away nodes (again in a clockwise direction) and repeat
        //
        for (int i = 1; i <= k/2; i++)
                for (int j = 0; j < n ; j++)
                        if (((double) rand())/RAND_MAX <= p)
                        {       int idx = (i+j) % n;    // Make it circular

                                // Choose new endpoint
                                //
                                int new_end = rand() % n;

                                while (new_end == idx || adjacency[j][new_end] != 0)
                                        new_end = rand() % n;

                                adjacency[j][idx]         = 0;
                                adjacency[idx][j]         = 0;
                                adjacency[j][new_end] = 1;
                                adjacency[new_end][j] = 1;
                        }
}

void rewire_WS2	(	int **	adjacency,
		double **	weight,
		int	n,
		int	k,
		double	p
	)

{
        //
        // OK, now rewire.
        // Start with the edges closest to each node (in a clockwise direction) and rewire with
        // probability p to give a new endpoint, provided that (i) it is not wired to us and
        // (ii) there is not already a link between us and that node. Then move out to the edges
        // connecting us to the next furthest away nodes (again in a clockwise direction) and repeat
        //
        for (int i = 1; i <= k/2; i++)
                for (int j = 0; j < n ; j++)
                        if (((double) rand())/RAND_MAX <= p)
                        {       int idx = (i+j) % n;    // Make it circular

                                // Choose new endpoint
                                //
                                int new_end = rand() % n;

                                while (new_end == idx || adjacency[j][new_end] != 0)
                                        new_end = rand() % n;

                                adjacency[j][idx]         = 0;
                                adjacency[idx][j]         = 0;
                                adjacency[j][new_end] = 1;
                                adjacency[new_end][j] = 1;

                                weight[j][new_end]        = 3;
                                weight[new_end][j]    = 3;
                        }
}

void splitNetwork	(	int **	adjacency,
		double *	betw,
		int	array_size
	)

                                                                 { 

        int best=0;
        double bestBetw=0;
        
        double betw2[array_size];
        
        //calculate highest betweenness
        for (int i=0;i<array_size;i++) {
                if (betw[i]>bestBetw) { 
                        bestBetw=betw[i];
                        best=i; 
                }       
        }
        
        bool oneFound=false;
        double minBetw=bestBetw;
        int linkToDelete=0;
        for (int i=0;i<array_size;i++) {
                if ((adjacency[best][i]==1)&&(best!=i)){
                                        adjacency[best][i]=0;
                                        adjacency[i][best]=0;   
                                        calculate_node_betweenness(betw2, adjacency, array_size);
                                        if (betw2[best]<minBetw) {
                                                minBetw=betw2[best];
                                                linkToDelete=i;
                                        }       
                                        adjacency[best][i]=1;
                                        adjacency[i][best]=1;
                                        
                }
        }

        adjacency[best][linkToDelete]=0;
        adjacency[linkToDelete][best]=0;
        
        /*bool oneFound=false;
        for (int i=0;i<array_size;i++) {
                if ((adjacency[best][i]==1)&&(best!=i)&&(oneFound==false)){
                                        adjacency[best][i]=0;
                                        adjacency[i][best]=0;   
                                        oneFound=true;
                }
        }*/

        
        cout<<"The highest value of betweenness is "<<betw[best]<< " (node "<<best+1<<")\n";
        cout<<"Link to delete is between link "<<(best+1)<<" and "<<linkToDelete+1<<endl;
}

double splitNetwork_edge_Threshold	(	int **	adjacency,
		int **	tempadjacency,
		double **	betw,
		int	array_size,
		double	modThreshold
	)

                                                                                                                               {

        int numberOfGroups_before=0;
        int numberOfGroups_after=0;
        double numberOfLinks_before=0;
        double bestbw=0;
        int bestx=-1;
        int besty=-1;

        int **groups_before=initialise_int_array(array_size);
        int **groups_after=initialise_int_array(array_size);
        
        double **totalResults=initialise_double_array(array_size);
        int **adjacency_before=initialise_int_array(array_size);        
        
        //adjacency_before is the connectivity matrix before the splitting procedure
        for (int i=0;i<array_size;i++) {
                for (int j=0;j<array_size;j++) {
                        adjacency_before[i][j]=adjacency[i][j];
                        totalResults[i][j]=0;
                }
        }
        
        int numberOfMembers_before[array_size];
        int numberOfMembers_after[array_size];

        for (int i=0;i<array_size; i++) {
                numberOfMembers_before[i]=0;
                numberOfMembers_after[i]=0;
        }
         
        numberOfGroups_before=getGroups(adjacency, groups_before,numberOfMembers_before,array_size);
        
        //If there is a link to remove (the number of nodes is not the same as the groups)
        if(numberOfGroups_before<array_size) {
                //Split network by removing the edge with the higher betweenness
                //First finding edge with highest betweenness
                for (int i=0;i<array_size;i++) {
                        for (int j=0;j<array_size;j++) {
                                if(betw[i][j]>bestbw) {
                                        bestbw=betw[i][j];
                                        bestx=i;besty=j;
                                }
                        }
                }
        
                //Removing edge
                tempadjacency[bestx][besty]=0;
                tempadjacency[besty][bestx]=0;
                
        }
        //adjacency is not renamed but it must be considered adjacency_after
        
        numberOfGroups_after=getGroups(tempadjacency, groups_after ,numberOfMembers_after,array_size);
        
        double **modularity=initialise_double_array(array_size);
        double **modularity_num=initialise_double_array(array_size);
        
        for (int i=0;i<numberOfGroups_after;i++) 
                for (int j=0;j<numberOfGroups_after;j++) {
                        modularity_num[i][j]=0;         
                }
        
        //Measuring total nymber of links       
        for (int i=0;i<array_size;i++)
                for (int j=0;j<array_size;j++) {
                         if(adjacency_before[i][j]==1 && i!=j)
                                numberOfLinks_before=numberOfLinks_before+1;    
        }

        for (int i=0;i<array_size;i++) 
                for (int j=0;j<array_size;j++) {
                        if (adjacency_before[i][j]==1 && i!=j)  {
                                for (int g1=0;g1<numberOfGroups_after;g1++) {
                                        if (isInGroup(i+1,groups_after[g1],numberOfMembers_after[g1])) {
                                                for (int g2=0;g2<numberOfGroups_after;g2++)
                                                        if (isInGroup(j+1,groups_after[g2],numberOfMembers_after[g2]))
                                                        {                                               
                                                                modularity_num[g1][g2]=modularity_num[g1][g2]+1;
                                                        }
                                        }
                                }
                        }
        }               
        
        for (int i=0;i<numberOfGroups_after;i++) 
                for (int j=0;j<numberOfGroups_after;j++) {
                        modularity[i][j]=((double)modularity_num[i][j]/(double)(numberOfLinks_before));         
                }       
                
        //cout<<"The total number of links is:"<<numberOfLinks_before<<"\n";
        //cout<<"The e matrix is the following:\n";
        //print_array_double(modularity_num,numberOfGroups_after);

        //calculation of the modularity Q
        //Q= Tr e - ||e^2||
        //where ||x|| indicates the sum of all elements of x  

        //calculation of the trace of the matrix e
        double total1=0;   
        for (int i=0;i<numberOfGroups_after;i++)
                total1=total1+modularity[i][i];
  
        for (int i=0;i<numberOfGroups_after;i++) {
                for (int j=0;j<numberOfGroups_after;j++) {
                        double partial=0;
                        for (int t=0;t<numberOfGroups_after;t++)
                                totalResults[i][j]+=partial+modularity[i][t]*modularity[t][j];
                }
        }         
  
        //calculation of ||e^2|| 
        double total2=0;
        for (int i=0;i<numberOfGroups_after;i++)
                for (int j=0;j<numberOfGroups_after;j++)
                        total2=total2+totalResults[i][j];
                
        //sum of the two parts Tr e and ||e^2|| 

        double q=total1-total2;

        //cout<<"The scalar value of q is "<<q<<"\n";

        for (int i=0;i<array_size;i++) {
                delete[] groups_before[i];
                delete[] groups_after[i];
                delete[] totalResults[i];
                delete[] adjacency_before[i];
                delete[] modularity[i];
                delete[] modularity_num[i];
        }

        delete[] groups_before;
        delete[] groups_after;
        delete[] totalResults;
        delete[] adjacency_before;
        delete[] modularity;
        delete[] modularity_num;
        
        if (numberOfGroups_after==array_size) {
        //cout<<"It is not possible to split the networks, since the number of groups is equal to the number of hosts.\n";
        return -1;
    }

    if (2*q>=modThreshold || q==0) {

        //cout<<"Successful splitting!\n";
        return q;

    }
    else {
        for (int i=0;i<array_size;i++)
            for (int j=0;j<array_size;j++) {
                tempadjacency[i][j]=adjacency_before[i][j];
        }
        return -1;
    }
}

double splitNetwork_Threshold	(	int **	adjacency,
		double *	betw,
		int	array_size,
		double	modThreshold
	)

                                                                                                   {

        int numberOfGroups_before=0;
        int numberOfGroups_after=0;
        double numberOfLinks_before=0;

        int **groups_before=initialise_int_array(array_size);
        int **groups_after=initialise_int_array(array_size);
        
        double **totalResults=initialise_double_array(array_size);
        int **adjacency_before=initialise_int_array(array_size);        
        
        //adjacency_before is the connectivity matrix before the splitting procedure
        for (int i=0;i<array_size;i++)
                for (int j=0;j<array_size;j++) {
                        adjacency_before[i][j]=adjacency[i][j];
                }
        int numberOfMembers_before[array_size];
        int numberOfMembers_after[array_size];

        for (int i=0;i<array_size; i++) {
                numberOfMembers_before[i]=0;
                numberOfMembers_after[i]=0;
        }
         
        numberOfGroups_before=getGroups(adjacency, groups_before,numberOfMembers_before,array_size);
        
        splitNetwork(adjacency, betw, array_size);
        
        //adjacency is not renamed but it must be considered adjacency_after
        
        numberOfGroups_after=getGroups(adjacency, groups_after ,numberOfMembers_after,array_size);
        
        double **modularity=initialise_double_array(array_size);
        double **modularity_num=initialise_double_array(array_size);
        
        for (int i=0;i<numberOfGroups_after;i++) 
                for (int j=0;j<numberOfGroups_after;j++) {
                        modularity_num[i][j]=0;         
                }
                
        for (int i=0;i<array_size;i++)
                for (int j=0;j<array_size;j++) {
                        if (adjacency_before[i][j]==1)
                                numberOfLinks_before=numberOfLinks_before+1;    
        }


        for (int i=0;i<array_size;i++) 
                for (int j=0;j<array_size;j++) {
                        if (adjacency_before[i][j]==1)  {
                                for (int g1=0;g1<numberOfGroups_after;g1++) {
                                        if (isInGroup(i+1,groups_after[g1],numberOfMembers_after[g1])) {
                                                for (int g2=0;g2<numberOfGroups_after;g2++)
                                                        if (isInGroup(j+1,groups_after[g2],numberOfMembers_after[g2])) {
                                                        
                                                                modularity_num[g1][g2]=modularity_num[g1][g2]+1;
                                                        }
                                        }
                                }
                        }
        }               
        
        for (int i=0;i<numberOfGroups_after;i++) 
                for (int j=0;j<numberOfGroups_after;j++) {
                        modularity[i][j]=modularity_num[i][j]/numberOfLinks_before;             
                }
   
        //cout<<"The total number of links is:"<<numberOfLinks_before<<"\n";
        //cout<<"The e matrix is the following:\n";
        //print_array_double(modularity_num,numberOfGroups_after);

        //calculation of the modularity Q
        //Q= Tr e - ||e^2||
        //where ||x|| indicates the sum of all elements of x  

        //calculation of the trace of the matrix e
        double total1=0;   
        for (int i=0;i<numberOfGroups_after;i++)
        total1=total1+modularity[i][i];

        for (int i=0;i<numberOfGroups_after;i++) {
                for (int j=0;j<numberOfGroups_after;j++) {
                                double partial=0;
                                for (int t=0;t<numberOfGroups_after;t++)
                                        totalResults[i][j]=partial+modularity[i][t]*modularity[t][j];
                }
        }         


        //calculation of ||e^2|| 
        double total2=0;
        for (int i=0;i<numberOfGroups_after;i++)
        for (int j=0;j<numberOfGroups_after;j++)
                total2=total2+totalResults[i][j];
        
        //sum of the two parts Tr e and ||e^2|| 

        double q=total1-total2;

        //cout<<"The scalar value of q is "<<q<<"\n";


        if (numberOfGroups_after==array_size) {
                //cout<<"It is not possible to split the networks, since the number of groups is equal to the number of hosts.\n";
                return -1;      
        }       

        if (q>modThreshold) {
        
                //cout<<"Successful splitting!\n";
                return q;

        }

        else {
                for (int i=0;i<array_size;i++)
                        for (int j=0;j<array_size;j++) {
                                adjacency[i][j]=adjacency_before[i][j];
                }
                return -1;
        }

        for (int i=0;i<array_size;i++) {
                delete[] groups_before[i];
                delete[] groups_after[i];
                delete[] totalResults[i];
                delete[] adjacency_before[i];
                delete[] modularity[i];
                delete[] modularity_num[i];
        }

        delete[] groups_before;
        delete[] groups_after;
        delete[] totalResults;
        delete[] adjacency_before;
        delete[] modularity;
        delete[] modularity_num;
        
}