higher-ordering cluster的C语言实现

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

#define INITIAL_SIZE	100
#define INCREMENT_SIZE 	100

int vertax_size;
int edge_size;
char filename_edge[20];
typedef struct Node{
	double value;
	int index;
}Node;
Node* eigenvalue;
Node* ordering;
double** eigenvector;
double* second_small_eigenvector;
double** D_inverse_matrix;
int** edge_info;
int** adjacentMatrix;
typedef struct Motif_Node{
	int x;
	int y;
	int z;
}Motif_Node;
typedef struct Motif_Set{
	Motif_Node* data;
	int length;
	int capacity;
}Motif_Set;
Motif_Set motif_set;

void readEigenvalue();
void readEigenVector();
void generate_sigma();
void read_D_inverse_matrix();
void readEdge();

void find_motifs();
int* get_neighbors(int, int*);
int* get_common(int*, int*, int, int);
void initial_motif_set();
void identify_motif(int, int, int*);
void test();
int calcu_cut(int*, int, int*, int);
int searchfor(int, int*, int);
void calcu_volume(int*, int, int*, int*, int, int*);
void test_1();

void singleCluster();

main(int argc, char* argv[])
{
	if( argc != 4 )
	{
		printf("This algorithm require 2 parameters"
				"\n\t1, the size of vertax"
				"\n\t2, the file name that contain edge inforamtion"
				"\n\t3, the size of edge"
				"\n");
		exit(0);
	}
	vertax_size = atoi(argv[1]);
	strcat(filename_edge, argv[2]);
	edge_size = atoi(argv[3]);
	readEigenvalue();
 	readEigenVector();
 	read_D_inverse_matrix();
 	generate_sigma();

	readEdge();
	find_motifs();
	singleCluster();
	return 0;
}

void readEigenvalue()
{
	FILE* fread;
	if( NULL == (fread = fopen("eigenvalue.data", "r")))
	{
		printf("function: readEigenvalue: open file(eigenvalue.data) error");
		exit(0);
	}
	int i;
	if( !(eigenvalue = (Node*)malloc(sizeof(Node) * (vertax_size + 1))))
	{
		printf("function: readEigenvalue: eigenvalue* malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
	{
		if( 1 != fscanf(fread, "%lf ", &eigenvalue[i].value))
		{
			printf("function: readEigenvalue: fscanf error: %d", i);
			exit(0);
		}
	}
}
void readEigenVector()
{
	FILE* fread;
	if( NULL == (fread = fopen("eigenvector.data", "r")))
	{
		printf("function: readEigenVector: open file(eigenvector.data) error");
		exit(0);
	}
	if( !(eigenvector = (double**)malloc(sizeof(double*) * (vertax_size + 1))) )
	{
		printf("function: readEigenvector: eigenvector** malloc error");
		exit(0);
	}
	int i;
	for( i = 1; i <= vertax_size; i++ )
	{
		if( !( eigenvector[i] = (double*)malloc(sizeof(double) * (vertax_size + 1)) ) )
		{
			printf("function: readEigenvector: eigenvector[%d]* malloc error: ", i);
			exit(0);
		}
	}
	int j;
	for( i = 1; i <= vertax_size; i++ )
	{
		for( j = 1; j <= vertax_size; j++ )
		{
			if( 1 != fscanf(fread, "%lf ", &eigenvector[i][j]))
			{
				printf("function: readEigenvector: fscanf error: (%d, %d)", i, j);
				exit(0);
			}
		}
	}
}
void read_D_inverse_matrix()
{
	FILE* fread;
	if( NULL == (fread = fopen("D_inverse_matrix.data", "r")) )
	{
		printf("function: read_D: open file error");
		exit(0);
	}
	if( !(D_inverse_matrix = (double**)malloc(sizeof(double*) * (vertax_size + 1))) )
	{
		printf("function: read_D: D_inverse_matrix** malloc error");
		exit(0);
	}
	int i;
	for( i = 1; i <= vertax_size; i++ )
		if( !(D_inverse_matrix[i] = (double*)malloc(sizeof(double) * (vertax_size + 1))) )
		{
			printf("function: read_D: D_inverse_matrix[%d]* malloc error");
			exit(0);
		}
	int j;
	for( i = 1; i <= vertax_size; i++ )
	{
		for( j = 1; j <= vertax_size; j++ )
		{
			if( 1 != fscanf(fread, "%lf ", &D_inverse_matrix[i][j]))
			{
				printf("function: read_D_inverse_matrix: fscanf error (%d, %d)", i, j);
				exit(0);
			}
		}
	}
	fclose(fread);
}
void generate_sigma()
{
	int i, j;
	double smaller_value, temp;
	int smaller_index, temp_index;
	for( i = 1; i <= vertax_size; i++ )
		eigenvalue[i].index = i;
	for( i = 1; i < vertax_size; i++ )
	{
		smaller_index = i;
		smaller_value = eigenvalue[i].value;
		for( j = i + 1; j <= vertax_size; j++ )
			if( eigenvalue[j].value < eigenvalue[smaller_index].value )
				smaller_index = j;
		if( smaller_index != i )
		{
			temp = eigenvalue[smaller_index].value;
			temp_index = eigenvalue[smaller_index].index;
			eigenvalue[smaller_index].value = eigenvalue[i].value;
			eigenvalue[smaller_index].index = eigenvalue[i].index;
			eigenvalue[i].value = temp;
			eigenvalue[i].index = temp_index;
		}
	}

	if( !(second_small_eigenvector = (double*)malloc(sizeof(double) * (vertax_size + 1))) )
	{
		printf("function: generate_sigma: second_small_eigenvector* malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
		second_small_eigenvector[i] = eigenvector[i][eigenvalue[2].index];

	if( !(ordering = (Node*)malloc(sizeof(Node) * (vertax_size + 1))) )
	{
		printf("function: generate_sigma: ordering* malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
		ordering[i].index = i;
	for( i = 1;  i <= vertax_size; i++ )
	{
		for( j = 1; j <= vertax_size; j++ )
		{
			ordering[i].value += D_inverse_matrix[i][j] * second_small_eigenvector[j];
		}
	}

	for( i = 1; i < vertax_size; i++ )
	{
		smaller_index = i;
		smaller_value = ordering[i].value;
		for( j = i + 1; j <= vertax_size; j++ )
			if( ordering[j].value < ordering[smaller_index].value )
				smaller_index = j;
		if( smaller_index != i )
		{
			temp = ordering[i].value;
			ordering[i].value = ordering[smaller_index].value;
			ordering[smaller_index].value = temp;
			temp_index = ordering[i].index;
			ordering[i].index = ordering[smaller_index].index;
			ordering[smaller_index].index = temp_index;
		}
	}

	printf("show the sigma value: \n");
	for( i = 1; i <= vertax_size; i++ )
		printf("%d ", ordering[i].index);
	printf("\n");
}

void readEdge()
{
	FILE* fread;
	if( NULL == (fread = fopen(filename_edge, "r")))
	{
		printf("function: readEdge: open file(%s) error", filename_edge);
		exit(0);
	}
	if( !(edge_info = (int**)malloc(sizeof(int*) * (edge_size + 1))) )
	{
		printf("function: readEdge: edge_info** malloc error");
		exit(0);
	}
	int i, j;
	for( i = 1; i <= edge_size; i++ )
	{
		if( !(edge_info[i] = (int*)malloc(sizeof(int) * (2 + 1))) )
		{
			printf("function: readEdge: edge_info[%d]* malloc error", i);
			exit(0);
		}
	}
	for( i = 1; i <= edge_size; i++ )
	{
		if( 2 != fscanf(fread, "%d\t%d", &edge_info[i][1], &edge_info[i][2]))
		{
			printf("function: readEdge: fscanf error: %d", i);
			exit(0);
		}
	}
	fclose(fread);
	
	/*	create the adjacent matrix	*/
	if( !(adjacentMatrix = (int**)malloc(sizeof(int*) * (vertax_size + 1))) )
	{
		printf("function: readEdge: adjacentMatrix** malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
		if( !(adjacentMatrix[i] = (int*)malloc(sizeof(int) * (vertax_size + 1))) )
		{
			printf("f: readEdge: adjacentMatrix[%d]* malloc error");
			exit(0);
		}
	for( i = 1; i <= vertax_size; i++ )
		for( j = 1; j <= vertax_size; j++ )
			adjacentMatrix[i][j] = 0;
	for( i = 1; i <= edge_size; i++ )
	{
		adjacentMatrix[edge_info[i][1]][edge_info[i][2]] = 1;
		//adjacentMatrix[edge_info[i][2]][edge_info[i][1]] = 1;
	}
}


void find_motifs()
{
	int i, j;
	int neighbor_length_0 = 0;
	int neighbor_length_1 = 0;
	int* neighbor_0;
	int* neighbor_1;
	int* neighbor_2;
	int* common_neighbors;
	initial_motif_set();
	for( i = 1; i <= vertax_size; i++ )
	{
		neighbor_0 = get_neighbors(i, &neighbor_length_0);			//得到节点i的所有的邻居和邻居的个数
		printf("Node: %3d ---- the length : %d\n", i, neighbor_length_0);
		for( j = 1; j <= neighbor_length_0; j++ )		
		{
			neighbor_1 = get_neighbors(neighbor_0[j], &neighbor_length_1);	//得到i的第j个邻居的所有邻居和邻居的个数
			common_neighbors = get_common(neighbor_0, neighbor_1, neighbor_length_0, neighbor_length_1);//求共同邻居
			identify_motif(i, neighbor_0[j], common_neighbors);		//求motif的数量
		}
	}
}
void initial_motif_set()
{
	if( !(motif_set.data = (Motif_Node*)malloc(sizeof(Motif_Node) * (INITIAL_SIZE + 1))) )
	{
		printf("F: initial_motif_set: motif_set.data* malloc error");
		exit(0);
	}
	motif_set.length = 0;
	motif_set.capacity = INITIAL_SIZE;
}
int* get_neighbors(int vertaxID, int* length)
{
	int* neighborSet;
	int i;
	*length = 0;
	for( i = vertaxID + 1; i <= vertax_size; i++ )
		if( 1 == adjacentMatrix[vertaxID][i] || 1 == adjacentMatrix[i][vertaxID] )
			(*length)++;
	//printf("length: %d\n", *length);
	if( !(neighborSet = (int*)malloc(sizeof(int) * (*length + 1))) )
	{
		printf("fun: get_neighbors: neighborSet* malloc error");
		exit(0);
	}
	int counter = 1;
	for( i = vertaxID + 1; i <= vertax_size; i++ )
		if( 1 == adjacentMatrix[vertaxID][i] || 1 == adjacentMatrix[i][vertaxID] )
			neighborSet[counter++] = i;

	return neighborSet;
}

int* get_common(int* set_1, int* set_2, int length_1, int length_2)
{
	int* lib;
	if( !(lib = (int*)malloc(sizeof(int) * (vertax_size + 1))) )
	{
		printf("function: get_common: lib* malloc error");
		exit(0);
	}
	int i, counter, common_counter;
	counter = common_counter = 0;
	for( i = 1; i <= vertax_size; i++ )
		lib[i] = 0;
	for( i = 1; i <= length_1; i++ )
		lib[set_1[i]]++;
	for( i = 1; i <= length_2; i++ )
		lib[set_2[i]]++;
	for( i = 1; i <= vertax_size; i++ )
		if( 2 == lib[i] )
			common_counter++;
	int* common_neighbor;
	if( !(common_neighbor = (int*)malloc(sizeof(int) * (common_counter + 1))) )
	{
		printf("f: get_common: common_neighbor* malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
		if( 2 == lib[i] )
			common_neighbor[++counter] = i;
	common_neighbor[0] = counter;
	return common_neighbor;
}
void identify_motif(int vertaxID_1, int vertaxID_2, int* commonSet)
{
	int i, addFlag;
	for( i = 1; i <= commonSet[0]; i++ )
	{
		addFlag = 0;
		if( adjacentMatrix[vertaxID_1][vertaxID_2] * adjacentMatrix[vertaxID_2][vertaxID_1] )
		{
			if( adjacentMatrix[vertaxID_1][commonSet[i]] * adjacentMatrix[vertaxID_2][commonSet[i]] && !(adjacentMatrix[commonSet[i]][vertaxID_1] + adjacentMatrix[commonSet[i]][vertaxID_2]) )
				addFlag = 1;
		}
		else if( adjacentMatrix[vertaxID_1][commonSet[i]] * adjacentMatrix[commonSet[i]][vertaxID_1] )
		{
			if( adjacentMatrix[vertaxID_1][vertaxID_2] * adjacentMatrix[commonSet[i]][vertaxID_2] && !(adjacentMatrix[vertaxID_2][vertaxID_1] + adjacentMatrix[vertaxID_2][commonSet[i]]) )
				addFlag = 1; 
		}
		else if( adjacentMatrix[vertaxID_2][commonSet[i]] * adjacentMatrix[commonSet[i]][vertaxID_2] )
		{
			if( adjacentMatrix[vertaxID_2][vertaxID_1] * adjacentMatrix[commonSet[i]][vertaxID_1] && !(adjacentMatrix[vertaxID_1][vertaxID_2] + adjacentMatrix[vertaxID_1][commonSet[i]]) )
				addFlag = 1; 
		}
		if( 1 == addFlag )
		{
			if( motif_set.length >= motif_set.capacity )
			{
				if( !(motif_set.data = (Motif_Node*)realloc(motif_set.data, sizeof(Motif_Node) * (motif_set.capacity + INCREMENT_SIZE))) )
				{
					printf("F: indentify_motif: motif_set.data* realloc error");
					exit(0);
				}
				motif_set.capacity += INCREMENT_SIZE;
			}
			motif_set.length++;
			motif_set.data[motif_set.length].x = vertaxID_1;
			motif_set.data[motif_set.length].y = vertaxID_2;
			motif_set.data[motif_set.length].z = commonSet[i];
		}
	}
}


calcu_cut(int* set_1, int length_1, int* set_2, int length_2)
{
	int cut_size = 0;
	int i, j;
	int small_index;
	int small_value;
	int temp;

	length_1 -= 1;
	length_2 -= 1;

	for( i = 1; i < length_1; i++ )		//all the operator of -1 is because of useless index 0
	{
		small_index = i;
		small_value = set_1[i];
		for( j = i + 1; j <= length_1; j++ )
			if( set_1[j] < set_1[small_index] )
				small_index = j;
		if( i != small_index )
		{
			temp = set_1[i];
			set_1[i] = set_1[small_index];
			set_1[small_index] = temp;
		}
	}

	for( i = 1; i <= motif_set.length; i++ )
	{
		if( 3 != searchfor(motif_set.data[i].x, set_1, length_1) + searchfor(motif_set.data[i].y, set_1, length_1) + searchfor(motif_set.data[i].z, set_1, length_1) && 0 != searchfor(motif_set.data[i].x, set_1, length_1) + searchfor(motif_set.data[i].y, set_1, length_1) + searchfor(motif_set.data[i].z, set_1, length_1) )
		{
			cut_size++;
		}
	}
	return cut_size;
}
/*
 * search the neighbors int the @set with length @length
 * */
searchfor(int vertaxID, int* set, int length)
{
	int low, high, mid;
	low = 1;
	high = length;		//because of my hibit of coding, (do not use index 0 of array), so have to -1
	while( low <= high )
	{
		mid = (low + high) / 2;
		if( set[mid] == vertaxID )
			return 1;
		else if( set[mid] > vertaxID )
			high = mid - 1;
		else
			low = mid + 1;
	}
	if( low > high )
		return 0;
}

void calcu_volume(int* set_1, int length_1, int* result_1, int* set_2, int length_2, int* result_2)
{
	length_1 -= 1;
	length_2 -= 1;

	int i;
	int* base;
	if( !(base = (int*)malloc(sizeof(int) * (vertax_size + 1))) )
	{
		printf("F: calcu_volume: base* malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
		base[i] = 0;
	*result_1 = *result_2 = 0;
	for( i = 1; i <= motif_set.length; i++ )
	{
		if( searchfor(motif_set.data[i].x, set_1, length_1) && !base[motif_set.data[i].x] )
		{
			(*result_1)++;
			base[motif_set.data[i].x] = 1;
		}
		else if( !base[motif_set.data[i].x] )
		{
			(*result_2)++;
			base[motif_set.data[i].x] = 1;
		}
		if( searchfor(motif_set.data[i].y, set_1, length_1) && !base[motif_set.data[i].y] )
		{
			(*result_1)++;
			base[motif_set.data[i].y] = 1;
		}
		else if( !base[motif_set.data[i].y] )
		{
			(*result_2)++;
			base[motif_set.data[i].y] = 1;
		}
		if( searchfor(motif_set.data[i].z, set_1, length_1) && !base[motif_set.data[i].z] )
		{
			(*result_1)++;
			base[motif_set.data[i].z] = 1;
		}
		else if( !base[motif_set.data[i].z] )
		{
			(*result_2)++;
			base[motif_set.data[i].z] = 1;
		}
	}
}

void singleCluster()
{
	int* set_1;
	int* set_2;
	int i, j;
	int cut, vol_1, vol_2;
	cut = vol_1 = vol_2 = 0;
	double* conductance;
	if( !(conductance = (double*)malloc(sizeof(double) * (vertax_size + 1))) )
	{
		printf("F: singleCluster: conductance* malloc error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
	{
		if( !(set_1 = (int*)malloc(sizeof(int) * (i + 1))) )
		{
			printf("F: singleCluster: set_1* malloc error");
			exit(0);
		}
		if( !(set_2 = (int*)malloc(sizeof(int) * (vertax_size - i + 1))) )
		{
			printf("F: singleCluster: set_2* malloc error");
			exit(0);
		}
		for( j = 1; j <= i; j++ )
			set_1[j] = ordering[j].index;
		for( j = i + 1; j <= vertax_size; j++ )
			set_2[j - i] = ordering[j].index;
		cut = calcu_cut(set_1, i + 1, set_2, vertax_size - i + 1);
		calcu_volume(set_1, i + 1, &vol_1, set_2, vertax_size - i + 1, &vol_2);
		conductance[i] = (double)cut/(double)(vol_1 < vol_2 ? vol_1 : vol_2);
		free(set_1);
		set_1 = NULL;
		free(set_2);
		set_2 = NULL;
	}
	
	FILE* fwrite;
	if( NULL == (fwrite = fopen("motif_conductance.data", "w")) )
	{
		printf("F: single_cluster: open file(motif_conductance.data) error");
		exit(0);
	}
	for( i = 1; i <= vertax_size; i++ )
	{
		printf("%f ", conductance[i]);
		fprintf(fwrite, "%d\t%f\n", i, conductance[i]);
	}
	fclose(fwrite);
}