void cal_centroid(double centroid[], node_type *over_node, node_type *extra_node){ int i, j; int total_size; if (extra_node->attribute == LEAF) { 	for (j=0; j < dim; j++) {		centroid[j] = extra_node->a[j];		total_size = 1;			for (i=0; i < M; i++) {			centroid[j] += over_node->ptr[i]->a[j]; 			total_size++;		}	}	centroid[j] = centroid[j] * 1.0 / total_size; } else { 	for (j=0; j < dim; j++) {		centroid[j] = extra_node->centroid[j];		total_size = 1;		for (i=0; i < M; i++) {			centroid[j] += over_node->ptr[i]->centroid[j] * over_node->ptr[i]->total_size;			total_size += over_node->ptr[i]->total_size;		}		centroid[j] = centroid[j] * 1.0 / total_size;	} } return;} /* cal_centroid */void reinsert(node_type *over_node, int over_level, node_type *extra_node, node_type *root) {  node_type *node_found;  node_type **overnode;  double *value;   int *sorted_index;  int i, start, stop;  double centroid[dim];  //printf("h5\n");  value = (double *)malloc(sizeof(double) * M+1);  sorted_index = (int *)malloc(sizeof(int) * M+1);  overnode = (node_type **)malloc((M+1) * sizeof(node_type *));  for (i=0; i<M; i++) {        overnode[i] = over_node->ptr[i];  }  overnode[M] = extra_node;  overnode[M]->parent = over_node;  cal_centroid(centroid, over_node, extra_node);  if (extra_node->attribute == LEAF) {  	for (i=0; i<M+1; i++) {		value[i] = dist_centroid(overnode[i]->a, centroid);		sorted_index[i] = i;	}  } else {  	for (i=0; i<M+1; i++) {		//value[i] = Dist2(overnode[i], over_node);		value[i] = dist_centroid(overnode[i]->centroid, centroid);		sorted_index[i] = i;	}  }  quicksort(sorted_index, value, 0, M);  for (i=0; i<dim; i++) {	over_node->a[i] = overnode[sorted_index[0]]->a[i];	over_node->b[i] = overnode[sorted_index[0]]->b[i];  }  over_node->ptr[0] = overnode[sorted_index[0]];  stop = M+1-reinsert_p;  for (i=1; i<stop; i++) {	over_node->ptr[i] = overnode[sorted_index[i]];	cal_MBR_node_node(over_node->a, over_node->b, over_node, overnode[sorted_index[i]]);  }  for (i=stop; i<M; i++)	over_node->ptr[i] = NULL;  over_node->vacancy = reinsert_p - 1;  adjust_MBR_delete(over_node);  update_centroid_node(over_node);  update_radius_node(over_node);  update_centroid(over_node);  start = M + 1 - reinsert_p;  for (i=start; i<M+1; i++) {	node_found = root;		choose_leaf_level(&node_found, root, 0, overnode[sorted_index[i]], over_level+extra_level); 	/* Test whether the node has room or not */	if(node_found->vacancy!=0) { 		/* have room to insert the entry */		overnode[sorted_index[i]]->parent = node_found;	        node_found->ptr[M - node_found->vacancy] = overnode[sorted_index[i]];        	node_found->vacancy--;	        adjust_MBR(overnode[sorted_index[i]]);  		update_centroid(overnode[sorted_index[i]]);	}	else {		overflow(node_found, over_level, over_level, overnode[sorted_index[i]], root);	}  }          //printf("h6\n");  return;  } /* reinsert */void overflow(node_type *over_node, int over_level, int old_level, node_type*extra_node, node_type *root){  node_type *node1, *node2;  //printf("h3\n");  if (over_level < old_level && over_level != 0) {	reinsert(over_node, over_level, extra_node, root);  }  else {	tree_node_allocate(&node1);	tree_node_allocate(&node2);        split(over_node, extra_node, node1, node2);        adjust_tree(over_node, over_level, old_level, node1, node2, root);	//printf("overflow, go out adjust, extra_node->id %d\n", extra_node->id);  }          //printf("h4\n");    return;  }void insert_node(node_type *root, double *data, int seq_id){  node_type *node_found, *new_node, *data_node;  int level_found;  int i;  /******/  /* I1 */  /******/  //printf("//////////////////////////////////////\n");  node_found = root;  extra_level=0;  tree_node_allocate(&data_node);  for (i=0; i<dim; i++) {	data_node->a[i] = data[i];	data_node->b[i] = data[i];  }  data_node->attribute = LEAF;  level_found = choose_leaf(&node_found, root, 0, data_node);  /* Now, node_found is the pointer to the leaf chosen */  /******/  /* I2 */  /******/    /********************************************************/  /* Add record to leaaf node:                            */  /* If L has room for another entry, install the entry   */  /* Otherwise invoke split() to split the node */  /********************************************************/      /* Make a leaf node */  tree_leaf_node_allocate(&new_node);  for(i=0; i < dim; i++) {        	new_node->a[i]=data[i];       	new_node->b[i]=data[i];  }  tree_node_deallocate(data_node);  new_node->id = seq_id;		//data[dim] is the seq_id     new_node->attribute=LEAF;  new_node->vacancy=M;  new_node->parent=node_found;  new_node->radius = 0.0;		// Ling  new_node->total_size = 0;		// Ling  new_node->ptr = NULL;  new_node->centroid = NULL;  /*  for(i=0; i < M; i++) {     new_node->ptr[i]=NULL;     //new_node->child_radius[i] = 0.0;     //new_node->child_size[i] = 0;  }  */  /* Test whether the node has room or not */  if(node_found->vacancy!=0) { 	//printf("insert new->id %d\n", seq_id);	//printf("h11\n");	/* have room to insert the entry */	node_found->ptr[M - node_found->vacancy] = new_node;      	node_found->vacancy--;     	adjust_MBR(new_node);	update_centroid(new_node);	//update_centroid(new_node, M-node_found->vacancy);	//printf("h12\n");  }  else {	//printf("h13\n");	overflow(node_found, level_found, level_found+1, new_node, root);	//printf("h14\n");  }  //printf("root->radius = %f, vacancy = %d, size = %d\n", root->radius, root->vacancy, root->total_size);  return;} /* insert_node */  /********************//* build_tree():    *//* build the R-tree *//********************/      void build_tree(node_type **root, double **data, int num_data){   int i, j;  /* make tree root node first */  tree_node_allocate(root);  for(i=0; i < dim; i++) { 	(*root)->b[i] = (double)(-1 * INT_MAX);       	(*root)->a[i] = (double)(INT_MAX);  }  (*root)->id = NO_ID;		//NO_ID  (*root)->attribute=ROOT;  (*root)->vacancy=M;     (*root)->parent=NULL;   (*root)->radius=0.0;  (*root)->total_size=0;  for(j=0; j < M; j++) {     	(*root)->ptr[j]=NULL;	//(*root)->child_radius[j]=0.0;	//(*root)->child_size[j]=0; }  /* add data to the tree */  for(i=0; i<num_data; i++) {        #ifdef DEBUG  	if ((i+1)%10==0)		printf("insert data %d\n", i+1);       #endif      	insert_node(*root, data[i], i);  // i is the seq id.       #ifndef DEBUG	free(data[i]);       #endif  } #ifndef DEBUG  free(data); #endif} /*build_tree */int make_data(char *datafile, double ***data){   int num_data;  int i,j;  FILE *fp_position;  fp_position=fopen(datafile,"r");    num_data = no_histogram;  (*data) = (double **)malloc(sizeof(double*) * num_data);  for(i=0; i<num_data; i++)        (*data)[i] = (double *)malloc(sizeof(double) * dim);  for(i=0; i<num_data; i++)	for (j=0; j<dim; j++)		fscanf(fp_position,"%lf",&((*data)[i][j]));  fclose(fp_position);  return(num_data);} /* make_data */void write_leaf_node(node_type *node, FILE *fp){  int i;                  for (i = 0; i<dim; i++)        fprintf(fp, "%f\n", (node->a)[i]);   for (i = 0; i<dim; i++)        fprintf(fp, "%f\n", (node->b)[i]);   fprintf(fp, "%d\n", node->attribute);  fprintf(fp, "%d\n", node->id);  fprintf(fp, "%d\n", node->vacancy);  return;}void write_inter_node(node_type *node, FILE *fp){  int i, count;  for (i = 0; i<dim; i++)	fprintf(fp, "%f\n", (node->a)[i]);  for (i = 0; i<dim; i++)          fprintf(fp, "%f\n", (node->b)[i]);  fprintf(fp, "%d\n", node->attribute);  for (i = 0; i<dim; i++)          fprintf(fp, "%f\n", (node->centroid)[i]);  fprintf(fp, "%d\n", node->vacancy);  fprintf(fp, "%f\n", node->radius);  fprintf(fp, "%d\n", node->total_size);  count = M - node->vacancy;  for (i=0; i<count; i++) {	if (node->ptr[i]->attribute != LEAF)		write_inter_node(node->ptr[i], fp);	else		write_leaf_node(node->ptr[i], fp);  }  return;}  void check_leaf_node(node_type *node, FILE *fp, int level, char point[]){  int i, j;  for (j=0; j < level; j++)                fprintf(fp, "  ");  for (i = 0; i<dim; i++) {        fprintf(fp, "[%.3f] ", (node->a)[i]);  }  fprintf(fp,"\n");  for (j=0; j < level; j++)                fprintf(fp, "  ");  fprintf(fp, "attribute: %d  ", node->attribute);  fprintf(fp, "ID: %d  ", node->id);  fprintf(fp, "vacancy: %d\n", node->vacancy);  point[node->id] = '1';  return;}int within_rectangle(double parent_a[], double parent_b[], double child_a[], double child_b[]){ int i;  for (i=0; i < dim; i++) {	if (child_a[i] < parent_a[i] || child_a[i] > parent_b[i] ||		child_b[i] > parent_b[i]) 		return FALSE;  }  return TRUE;}/*int within_sphere(double parent_centroid[], double parent_radius, double child_centroid[], double child_radius){ int i; double centroid_centroid = 0.0; for (i=0; i < dim; i++) {	centroid_centroid += pow(parent_centroid[i] - child_centroid[i], 2.0); } if (centroid_centroid + child_radius > parent_radius)	return FALSE; else	return TRUE;}*/int within_sphere(double parent_centroid[], double parent_radius, double point[]){ int i=dim; double centroid_centroid = 0.0; while (i--)        centroid_centroid += pow(parent_centroid[i] - point[i], 2.0);  centroid_centroid = sqrt(centroid_centroid); if (centroid_centroid - _EPSILON > parent_radius) {        printf("Error: centroid_centroid = %f, parent_radius = %f\n", centroid_centroid, parent_radius);        fprintf(stderr, "Error: centroid_centroid = %f, parent_radius = %f\n", centroid_centroid, parent_radius);        return FALSE; } else        return TRUE;}void check_inter_node(node_type *node, FILE *fp, int level, char point[], double **data, int total_level){  char local_point[no_histogram];  int i, count;  int j;  memset(local_point, '0', sizeof(local_point));  for (j=0; j < level; j++)                fprintf(fp, "  ");  fprintf(fp, " [a] ");  for (i = 0; i<dim; i++) {        fprintf(fp, "[%.3f,", (node->a)[i]);        fprintf(fp, "%.3f]  ", (node->b)[i]);  }  fprintf(fp, "\n");  for (j=0; j < level; j++)                fprintf(fp, "  ");  fprintf(fp, " [centroid] ");  for (i = 0; i<dim; i++)        fprintf(fp, "%.3f ", (node->centroid)[i]);  fprintf(fp, "\n");  for (j=0; j < level; j++)                fprintf(fp, "  ");  fprintf(fp, "attribute: %d  ", node->attribute);  fprintf(fp, "vacancy: %d  ", node->vacancy);  fprintf(fp, "radius: %.3f  ", node->radius);  fprintf(fp, "size: %d\n", node->total_size);  count = M - node->vacancy;  for (i=0; i<count; i++) {        fprintf(fp, "\n");        if (!within_rectangle(node->a, node->b, node->ptr[i]->a, node->ptr[i]->b)) {                fprintf(fp, "\nError!  Child %d [%p] of [%p] is out of rectangle bound\n", i, node->ptr[i], node);                fprintf(stderr, "\nError!  Child %d [%p] of [%p] is out of rectangle bound\n", i, node->ptr[i], node);        }        for (j=0; j < level+1; j++)                fprintf(fp, "--");        if (node->ptr[i]->attribute != LEAF) {                fprintf(fp, "[%p] NODE level %d\n", node->ptr[i], level+1);                check_inter_node(node->ptr[i], fp, level+1, local_point, data, total_level);        } else {                fprintf(fp, "[%p] LEAF level %d\n", node->ptr[i], level+1);                check_leaf_node(node->ptr[i], fp, level+1, local_point);                if (level+1 != total_level-1) {                        fprintf(fp, "Error!  leaf level = %d, total level = %d, node = %p\n", level+2, total_level, node);                        fprintf(stderr, "Error!  leaf level = %d, total level = %d, node = %p\n", level+2, total_level, node);                }        }  }  for (i=0; i < no_histogram; i++) {        if (local_point[i] == '1') {                within_sphere(node->centroid, node->radius, data[i]);                point[i] = '1';        }  }  return;}void print_node_address(node_type *node, FILE *fp){  int count = M - node->vacancy;  int i;  fprintf(fp, "[%p] : ", node);  for (i=0; i<count; i++)        fprintf(fp, "[%p][%d]%d ", node->ptr[i], node->ptr[i]->id, node->ptr[i]->attribute);  fprintf(fp, "\n");  if (node->ptr[0]->attribute != LEAF)        for (i=0; i < count; i++)                print_node_address(node->ptr[i], fp);  return;}void check_srtree(node_type *root, double **data){  int i;  char point[no_histogram];  int total_level=1;  FILE *fp;  node_type *node;  memset(point, '0', sizeof(point));  node = root;  while (node->attribute != LEAF) {        total_level++;        node = node->ptr[0];  }  if (!(fp=fopen(CHECK_TREE_LOG, "w"))) {        printf("Can't open file, %s\n", CHECK_TREE_LOG);        exit(EXIT_FAILURE);  }  fprintf(fp, "\nTree height: %d\n", total_level);  fprintf(fp, "\nROOT [%p]\n", root);  check_inter_node(root, fp, 0, point, data, total_level);  for (i=0; i < no_histogram; i++) {        if (point[i] == '1') {                within_sphere(root->centroid, root->radius, data[i]);        }  }  fprintf(fp, "\n");  print_node_address(root, fp);  fclose(fp);  return;}void save_srtree(node_type *root, char save_tree_file[]){  FILE *fp;  //fp = fopen(SAVE_SRTREE_FILE, "w");  fp = fopen(save_tree_file, "w");  if (!fp) {	printf("Can't write to %s\n", save_tree_file);	exit(EXIT_FAILURE);  }    write_inter_node(root, fp);  fclose(fp);   return;}void free_tree(node_type *node){  int i;  if (!node) return;  free(node->a);  free(node->b);  if (node->attribute != LEAF)  	free(node->centroid);  for (i=0; i < M-node->vacancy; i++)	free_tree(node->ptr[i]);  free(node);  return;}void destroy(double **data){ int i; for (i=0; i < no_histogram; i++)	free(data[i]); free(data); return;}int main(int argc, char *argv[]){   int num_data;  double **data;  config_type config;  node_type *root;  float  userTime, sysTime;  struct rusage myTime1, myTime2;  ////////////////////////////  if (argc != 2) {	printf("Build SR-tree\n\n");        printf("Usage: %s <config>\n", argv[0]);        exit(EXIT_FAILURE);  }  ////////////////////////////  getrusage(RUSAGE_SELF,&myTime1);  printf("initialize\n");  initialize(&config, argv[1]);     printf("make_data\n");  num_data=make_data(config.datafile, &data);  printf("build_tree\n");  build_tree(&root, data, num_data);  getrusage(RUSAGE_SELF,&myTime2);  ////////////////////////////  #ifdef DEBUG  printf("check_srtree\n");  check_srtree(root, data);  destroy(data);  #endif  printf("save_srtree\n");  save_srtree(root, config.save_tree_file);  printf("free_tree\n");  free_tree(root);  ////////////////////////////  userTime =                ((float) (myTime2.ru_utime.tv_sec  - myTime1.ru_utime.tv_sec)) +                ((float) (myTime2.ru_utime.tv_usec - myTime1.ru_utime.tv_usec)) * 1e-6;   sysTime =                ((float) (myTime2.ru_stime.tv_sec  - myTime1.ru_stime.tv_sec)) +                ((float) (myTime2.ru_stime.tv_usec - myTime1.ru_stime.tv_usec)) * 1e-6;  fprintf(stdout, "User time : %f seconds\n",userTime);  fprintf(stdout, "System time : %f seconds\n",sysTime);  fprintf(stdout, "Total time : %f seconds\n",userTime+sysTime);  return(0);  } /* main */