///////////
// Changes made:
// 1. Added another parameter: number of populations
// 2. Added a hierarchical merging step based on density change between centroids of two hyper-regions
// 3. Picked the longest dimensions between the population centroids to judge whether the two parts should be merged
// 4. Removed checking time parameter
// 5. Output error to stderr
// 6. Fixed the bug of density threshold always = 3
// 7. Added another error (select_num_bin<min_grid) || (select_num_bin>max_grid) to STDERR
// 8. Fixed a bug for 2D data by using K=e*K
// 9. Added some header files, may not be necessary
// 10. Added a lower bound (at least two) for number of populations
/***************************************************************************************************************************************
	
	FLOCK: FLOw cytometry Clustering without K (Named by: Jamie A. Lee and Richard H. Scheuermann) 
	
	Author: (Max) Yu Qian, Ph.D.
	
	Copyright: Scheuermann Lab, Dept. of Pathology, UTSW
	
	Development: November 2005 ~ Forever

	Status: July 2010: Release 2.0

	Usage: flock data_file
		    Note: the input file format must be channel values and the delimiter between two values must be a tab.


    	
****************************************************************************************************************************************/

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <unistd.h>
#include <assert.h>



#define DEBUG 0
#define LINE_LEN 1024
#define FILE_NAME_LEN 128
#define PARA_NAME_LEN 64
#define MAX_VALUE 1000000000
#define MIN_GRID 6
#define MAX_GRID 50
#define E_T 1.0

#define NORM_METHOD 2 //2 if z-score; 0 if no normalization; 1 if min-max 
#define KMEANS_TERM 10 
#define MAX_POP_NUM 128

#ifndef max
    #define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif

#ifndef min
    #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#endif

static long **Gr=0;
static long *gcID = 0;          /* grid cluster IDs */
static long *cluster_count=0;   /* count of nodes per cluster */
static long ndense=0;
static long ndim=0;
/* cid changes between depth-first searches, but is constant within a
   single search, so it goes here. */
static long cid=0;
/* Do a depth-first search for a single connected component in graph
 * G.  Start from node, tag the nodes found with cid, and record
 * the tags in grid_clusterID.  Also, record the node count in
 * cluster_count.  If we find a node that has already been assigned to
 * a cluster, that means we're merging two clusters, so zero out the
 * old cid's node count.
 *
 * Note that our graph is constructed as a DAG, so we can be
 * guaranteed to terminate without checking for cycles.  
 *
 * Note2: this function can potentially recurse to depth = ndense.
 * Plan stack size accordingly.  
 *
 * Output:
 *
 * grid_clusterID[] -- array where we tag the nodes.
 * cluster_count[]  -- count of the number of nodes per cluster.
 */
static void merge_cluster(long from, long into)
{
  int i;

  for(i=0; i<ndense; ++i)
    if(gcID[i] == from)
      gcID[i] = into;
}
   
void depth_first(long node)
{
  long i;

  if(gcID[node] == cid)         // we're guaranteed no cycles, but it is possible to reach a node 
    return;                     // through two different paths in the same cluster.  This early
                                // return saves us some unneeded work.

  /* Check to see if we are merging a cluster */
  if(gcID[node] >= 0) {
    /* We are, so zero the count for the old cluster. */
    cluster_count[ gcID[node] ] = 0;  
    merge_cluster(gcID[node], cid);
    return;
  }

  /* Update for this node */
  gcID[node] = cid;
  cluster_count[cid]++;

  /* Recursively search the child nodes */
  for(i=0; i<ndim; ++i)
    if(Gr[node][i] >= 0)      /* This is a child node */
      depth_first(Gr[node][i]);
}

void bail(const char *msg)
{
  fprintf(stderr,"%s",msg);
  exit(0);
}


static void check_clusters(long *gcID, int ndense)
{
  int i;

  for(i=0; i<ndense; ++i)
    if(gcID[i] < 0) {
      fprintf(stderr,"faulty cluster id at i= %d\n",i);
      exit(0);
    }
}



long find_connected(long **G, long n_dense_grids, long num_dm, long *grid_clusterID)
{
  long nclust=0;                  /* number of clusters found */
  long i;
  long *subfac;
  long clustid=0;
  int subval=0,nempty=0;
  
  int sz = n_dense_grids*sizeof(long); 
  cluster_count = malloc(sz);
  if(!cluster_count)
    bail("find_connected: Unable to allocate %zd bytes.\n");
  memset(cluster_count,0,sz);

  /* set up the statics that will be used in the DFS */
  Gr=G;
  gcID = grid_clusterID;
  ndense = n_dense_grids;
  ndim = num_dm;
  
  for(i=0;i<ndense;++i)
    grid_clusterID[i] = -1;

  for(i=0;i<ndense;++i) {
    if(grid_clusterID[i] < 0) {  /* grid hasn't been assigned yet */
      cid = nclust++;
      depth_first(i);
    }
  }


  
  /* At this point we probably have some clusters that are empty due to merging.
     We want to compact the cluster numbering to eliminate the empty clusters. */

  subfac = malloc(sz);
  if(!subfac)
    bail("find_connected: Unable to allocate %zd bytes.\n");
  subval=0;
  nempty=0;

  /* cluster #i needs to have its ID decremented by 1 for each empty cluster with
     ID < i.  Precaclulate the decrements in this loop: */
  for(i=0;i<nclust;++i) {
    //clustid = grid_clusterID[i];
    if(cluster_count[i] == 0) {
      subval++;
      nempty++;
    }
    subfac[i] = subval;
  }

  //printf("nempty is %d\n",nempty);
  
  /* Now apply the decrements to all of the dense grids */
  for(i=0;i<ndense;++i) {
    clustid = grid_clusterID[i];
    grid_clusterID[i] -= subfac[clustid];
  }


  
  /* correct the number of clusters found */
  nclust -= nempty;

  //printf("nclust is %d\n",nclust);

  return nclust;
}

/************************************* Read basic info of the source file **************************************/
/************************************* Read basic info of the source file **************************************/
void getfileinfo(FILE *f_src, long *file_Len, long *num_dm, char *name_string, long *time_ID)
{
  char src[LINE_LEN];
  char current_name[64];
  char prv;

  long num_rows=0;
  long num_columns=0;
  long ch='\n';
  long prev='\n';
  long time_pos=0;
  long i=0;
  long j=0;
  int sw=0;

  src[0]='\0';
  fgets(src, LINE_LEN, f_src);

  name_string[0]='\0';
  current_name[0]='\0';
  prv='\n';

  while ((src[i]==' ') || (src[i]=='\t')) //skip space and tab characters
    i++;

  while ((src[i]!='\0') && (src[i]!='\n')) //repeat until the end of the line
    {
      current_name[j]=src[i];
		
      if ((src[i]=='\t') && (prv!='\t')) //a complete word
        {
          current_name[j]='\0';

          if (0!=strcmp(current_name,"Time"))
            {
              num_columns++; //num_columns does not inlcude the column of Time
              time_pos++;
              if (sw) {
                  strcat(name_string,"\t");
              }
              strcat(name_string,current_name); 
              sw = 1;
            }
          else
            {
              *time_ID=time_pos;
            }

          current_name[0]='\0';
          j=0;			
        }		
		
      if ((src[i]=='\t') && (prv=='\t')) //a duplicate tab or space
        {
          current_name[0]='\0';
          j=0;
        }
		
      if (src[i]!='\t')
        j++;
		
      prv=src[i];
      i++;
    }
	
  //name_string[j]='\0';
  if (prv!='\t') //the last one hasn't been retrieved
    {
      current_name[j]='\0';
      if (0!=strcmp(current_name,"Time"))
        {
          num_columns++;
          strcat(name_string,"\t");
          strcat(name_string,current_name);
          time_pos++;
        }
      else
        {
          *time_ID=time_pos;
        }
    }
  if (DEBUG==1)
    {
      printf("time_ID is %ld\n",*time_ID);
      printf("name_string is %s\n",name_string);
    }

  //start computing # of rows

  while ((ch = fgetc(f_src))!= EOF )
    {
      if (ch == '\n')
        {
          ++num_rows;
        }
      prev = ch;
    }
  if (prev!='\n')
    ++num_rows;

   if (num_rows<50)
  {
    fprintf(stderr,"Number of events in the input file is too few and should not be processed!\n"); //modified on July 23, 2010
	exit(0);
  }
	
  *file_Len=num_rows;
  *num_dm=num_columns; 

  printf("original file size is %ld; number of dimensions is %ld\n", *file_Len, *num_dm);
}



/************************************* Read the source file into uncomp_data **************************************/
void readsource(FILE *f_src, long file_Len, long num_dm, double **uncomp_data, long time_ID)
{
  //long time_pass=0; //to mark whether the time_ID has been passed
  long index=0;

  long i=0;
  long j=0;
  long t=0;

  char src[LINE_LEN];
  char xc[LINE_LEN/10];

  src[0]='\0';
  fgets(src,LINE_LEN, f_src); //skip the first line about parameter names

  while (!feof(f_src) && (index<file_Len)) //index = 0, 1, ..., file_Len-1
    {
      src[0]='\0';	    
      fgets(src,LINE_LEN,f_src);
      i=0;
      //time_pass=0;
						
      //if (time_ID==-1)  //there is no time_ID
      //  {
          for (t=0;t<num_dm;t++) 
            {
              xc[0]='\0';
              j=0;
              while ((src[i]!='\0') && (src[i]!='\n') && (src[i]!=' ') && (src[i]!='\t'))
                {
                  xc[j]=src[i];
                  i++;
                  j++;
                }
		
              xc[j]='\0';	    
              i++;

              uncomp_data[index][t]=atof(xc);
            }	
       // }
      /*else
        {
          for (t=0;t<=num_dm;t++) //the time column needs to be skipped, so there are num_dm+1 columns
            {
              xc[0]='\0';
              j=0;
              while ((src[i]!='\0') && (src[i]!='\n') && (src[i]!=' ') && (src[i]!='\t'))
                {
                  xc[j]=src[i];
                  i++;
                  j++;
                }
		
              xc[j]='\0';	    
              i++;

              if (t==time_ID)
                {
                  time_pass=1;
                  continue;
                }
				
              if (time_pass)
                uncomp_data[index][t-1]=atof(xc);
              else
                uncomp_data[index][t]=atof(xc);
            }
        }*/ //commented by Yu Qian on Aug 31, 2010 as we do not want to check time_ID anymore        	
      index++;     	
      //fprintf(fout_ID,"%s",src);
    } //end of while
	
  if (DEBUG == 1)
    {
      printf("the last line of the source data is:\n");
      for (j=0;j<num_dm;j++)
        printf("%f ",uncomp_data[index-1][j]);
      printf("\n");
    }
}


/**************************************** Normalization ******************************************/
void tran(double **orig_data, long file_Len, long num_dm, long norm_used, double **matrix_to_cluster)
{
  long i=0;
  long j=0;

  double biggest=0;
  double smallest=MAX_VALUE;

  double *aver; //average of each column
  double *std; //standard deviation of each column

  aver=(double*)malloc(sizeof(double)*file_Len);
  memset(aver,0,sizeof(double)*file_Len);

  std=(double*)malloc(sizeof(double)*file_Len);
  memset(std,0,sizeof(double)*file_Len);	
		
  if (norm_used==2) //z-score normalization
    {
      for (j=0;j<num_dm;j++)
        {
          aver[j]=0;
          for (i=0;i<file_Len;i++)
            aver[j]=aver[j]+orig_data[i][j];
          aver[j]=aver[j]/(double)file_Len;

          std[j]=0;
          for (i=0;i<file_Len;i++)
            std[j]=std[j]+((orig_data[i][j]-aver[j])*(orig_data[i][j]-aver[j]));
          std[j]=sqrt(std[j]/(double)file_Len);
			
          for (i=0;i<file_Len;i++)
            matrix_to_cluster[i][j]=(orig_data[i][j]-aver[j])/std[j];  //z-score normalization
        }
    }

  if (norm_used==1) //0-1 min-max normalization
    {
      for (j=0;j<num_dm;j++)
        {
          biggest=0;
          smallest=MAX_VALUE;
          for (i=0;i<file_Len;i++)
            {
              if (orig_data[i][j]>biggest)
                biggest=orig_data[i][j];
              if (orig_data[i][j]<smallest)
                smallest=orig_data[i][j];
            }
			
          for (i=0;i<file_Len;i++)
            {
              if (biggest==smallest)
                matrix_to_cluster[i][j]=biggest;
              else
                matrix_to_cluster[i][j]=(orig_data[i][j]-smallest)/(biggest-smallest);
            }
        }
    }

  if (norm_used==0) //no normalization
    {
      for (i=0;i<file_Len;i++)
        for (j=0;j<num_dm;j++)
          matrix_to_cluster[i][j]=orig_data[i][j];
    }

}



/********************************************** RadixSort *******************************************/
/* Perform a radix sort using each dimension from the original data as a radix.
 * Outputs:
 * sorted_seq   -- a permutation vector mapping the ordered list onto the original data.
 *                  (sorted_seq[i] -> index in the original data of the ith element of the ordered list)
 * grid_ID      -- mapping between the original data and the "grids" (see below) found as a byproduct
 *                  of the sorting procedure.
 * num_nonempty -- the number of grids that occur in the data (= the number of distinct values assumed
 *                  by grid_ID)
 */

void radixsort_flock(long **position,long file_Len,long num_dm,long num_bin,long *sorted_seq,long *num_nonempty,long *grid_ID)
{
  long i=0;
  long length=0; 
  long start=0;
  long prev_ID=0;
  long curr_ID=0;
	
  long j=0;
  long t=0;
  long p=0;
  long loc=0;
  long temp=0;
  long equal=0;
	
  long *count; //count[i]=j means there are j numbers having value i at the processing digit
  long *index; //index[i]=j means the starting position of grid i is j
  long *cp; //current position
  long *mark; //mark[i]=0 means it is not an ending point of a part, 1 means it is (a "part" is a group of items with identical bins for all dimensions)
  long *seq; //temporary sequence

  count=(long*)malloc(sizeof(long)*num_bin);
  memset(count,0,sizeof(long)*num_bin);

  cp=(long*)malloc(sizeof(long)*num_bin);
  memset(cp,0,sizeof(long)*num_bin);

  index=(long*)malloc(sizeof(long)*num_bin); // initialized below

  seq=(long*)malloc(sizeof(long)*file_Len);
  memset(seq,0,sizeof(long)*file_Len);
	
  mark=(long*)malloc(sizeof(long)*file_Len);
  memset(mark,0,sizeof(long)*file_Len);
	
  for (i=0;i<file_Len;i++)
    {
      sorted_seq[i]=i;
      mark[i]=0;
      seq[i]=0;
    }
  for (i=0;i<num_bin;i++)
    {
      index[i]=0;
      cp[i]=0;
      count[i]=0;
    }

  for (j=0;j<num_dm;j++)
    {
      if (j==0) //compute the initial values of mark
        {
          for (i=0;i<file_Len;i++)
            count[position[i][j]]++; // initialize the count to the number of items in each bin of the 0th dimension

          index[0] = 0;
          for (i=0;i<num_bin-1;i++)
            {
              index[i+1]=index[i]+count[i];  //index[k]=x means k segment starts at x (in the ordered list)
              if ((index[i+1]>0) && (index[i+1]<=file_Len))
                {
                  mark[index[i+1]-1]=1; // Mark the end of the segment in the ordered list
                }
              else
                {
                  printf("out of myboundary for mark at index[i+1]-1.\n");
                }
            }
          mark[file_Len-1]=1;
			
          for (i=0;i<file_Len;i++)
            {
              /* Build a permutation vector for the partially ordered data.  Store the PV in sorted_seq */
              loc=position[i][j];
              temp=index[loc]+cp[loc]; //cp[i]=j means the offset from the starting position of grid i is j 
              sorted_seq[temp]=i;  //sorted_seq[i]=temp is also another way to sort
              cp[loc]++;
            }
        }
      else
        {
          //reset count, index, loc, temp, cp, start, and length
          length=0;
          loc=0;
          temp=0;
          start=0;
          for (p=0;p<num_bin;p++)
            {
              cp[p]=0;
              count[p]=0;
              index[p]=0;
            }

          for (i=0;i<file_Len;i++)
            {
              long iperm = sorted_seq[i]; // iperm allows us to traverse the data in sorted order.
              if (mark[i]!=1)
                {
                  /* Count the number of items in each bin of
                     dimension j, BUT we are going to reset at the end
                     of each "part".  Thus, the total effect is to do
                     a sort by bin on the jth dimension for each group
                     of data that has been identical for the
                     dimensions processed up to this point.  This is
                     the standard radix sort procedure, but doing it
                     this way saves us having to allocate buckets to
                     hold the data in each group of "identical-so-far"
                     elements. */
                  count[position[iperm][j]]++;  //count[position[i][j]]++;
                  length++;                     // This is the total length of the part, irrespective of the value of the jth component
                                                // (actually, less one, since we don't increment for the final element below)
                }
              if (mark[i]==1)
                {
                  //length++;
                  count[position[iperm][j]]++;//count[position[i][j]]++;  //the current point must be counted in
                  start=i-length; //this part starts from start to i: [start,i]
                  /* Now we sort on the jth radix, just like we did for the 0th above, but we restrict it to just the current part.
                     This would be a lot more clear if we broke this bit of code out into a separate function and processed recursively,
                     plus we could multi-thread over the parts.  (Hmmm...)
                  */
                  index[0] = start; // Let index give the offset within the whole ordered list.
                  for (t=0;t<num_bin-1;t++)
                    {
                      index[t+1]=index[t]+count[t];
						
                      if ((index[t+1]<=file_Len) && (index[t+1]>0))
                        {
                          mark[index[t+1]-1]=1; // update the part boundaries to include the differences in the current radix.
                        }
						
                    }
                  mark[i]=1;

                  /* Update the permutation vector for the current part (i.e., from start to i).  By the time we finish the loop over i
                     the PV will be completely updated for the partial ordering up to the current radix. */
                  for (t=start;t<=i;t++)
                    {
                      loc=position[sorted_seq[t]][j];//loc=position[t][j];
                      temp=index[loc]+cp[loc];
                      if ((temp<file_Len) && (temp>=0)) 
                        {
                          // seq is a temporary because we have to maintain the old PV until we have finished this step.
                          seq[temp]=sorted_seq[t];  //sorted_seq[i]=temp is also another way to sort
                          cp[loc]++;
                        }
                      else
                        {
                          printf("out of myboundary for seq at temp.\n");
                        }
                    }

                  for (t=start;t<=i;t++)
                    {
                      // copy the temporary back into sorted_seq.  sorted_seq is now updated for radix j up through
                      // entry i in the ordered list.
                      if ((t>=0) && (t<file_Len))
                        sorted_seq[t]=seq[t];
                      else
                        printf("out of myboundary for seq and sorted_seq at t.\n");
                    }
                  //reset count, index, seq, length, and cp
                  length=0;
                  loc=0;
                  temp=0;
                  for (p=0;p<num_bin;p++)
                    {
                      cp[p]=0;
                      count[p]=0;
                      index[p]=0;
                    }
                }
            }//end for i
        }//end else
    }//end for j

  /* sorted_seq[] now contains the ordered list for all radices.  mark[] gives the boundaries between groups of elements that are
     identical over all radices (= dimensions in the original data) (although it appears we aren't going to make use of this fact) */
  
  for (i=0;i<file_Len;i++)
    grid_ID[i]=0; //in case the initial value hasn't been assigned
  *num_nonempty=1; //starting from 1!	

  /* assign the "grid" identifiers for all of the data.  A grid will be what we were calling a "part" above.  We will number them
     serially and tag the *unordered* data with the grid IDs.  We will also count the number of populated grids (in general there will
     be many possible combinations of bin values that simply never occur) */
  
  for (i=1;i<file_Len;i++)
    {
      equal=1;
      prev_ID=sorted_seq[i-1];
      curr_ID=sorted_seq[i];
      for (j=0;j<num_dm;j++)
        {
          if (position[prev_ID][j]!=position[curr_ID][j])
            {	
              equal=0;  //not equal
              break;
            }
        }
		
      if (equal)
        {
          grid_ID[curr_ID]=grid_ID[prev_ID];
        }
      else
        {
          *num_nonempty=*num_nonempty+1;
          grid_ID[curr_ID]=grid_ID[prev_ID]+1;
        }
      //all_grid_vol[grid_ID[curr_ID]]++;
    }

  //free memory
  free(count);
  free(index);	
  free(cp);	
  free(seq);
  free(mark); 
  
}

/********************************************** Compute Position of Events ************************************************/
void compute_position(double **data_in, long file_Len, long num_dm, long num_bin, long **position, double *interval)
{
  /* What we are really doing here is binning the data, with the bins
     spanning the range of the data and number of bins = num_bin */
  long i=0;
  long j=0;

  double *small; //small[j] is the smallest value within dimension j
  double *big; //big[j] is the biggest value within dimension j
		
  small=(double*)malloc(sizeof(double)*num_dm);
  memset(small,0,sizeof(double)*num_dm);

  big=(double*)malloc(sizeof(double)*num_dm);
  memset(big,0,sizeof(double)*num_dm);
	
	
  for (j=0;j<num_dm;j++)
    {
      big[j]=MAX_VALUE*(-1);
      small[j]=MAX_VALUE;
      for (i=0;i<file_Len;i++)
        {
          if (data_in[i][j]>big[j])
            big[j]=data_in[i][j];

          if (data_in[i][j]<small[j])
            small[j]=data_in[i][j];
        }
		
      interval[j]=(big[j]-small[j])/(double)num_bin;	//interval is computed using the biggest value and smallest value instead of the channel limit
      /* XXX: I'm pretty sure the denominator of the fraction above should be num_bin-1. */
    }
    
  for (j=0;j<num_dm;j++)
  {
     for (i=0;i<file_Len;i++)
     {	
        if (data_in[i][j]>=big[j])
           position[i][j]=num_bin-1;
        else
        {
           position[i][j]=(long)((data_in[i][j]-small[j])/interval[j]); //position[i][j]=t means point i is at the t grid of dimensional j
           if ((position[i][j]>=num_bin) || (position[i][j]<0))
           {
               //printf("position mis-computed in density analysis!\n");
               //exit(0);
			   fprintf(stderr,"Incorrect input file format or input parameters (number of bins overflows)!\n"); //modified on July 23, 2010
				exit(0);
           }
        }
     }
  }


  free(small);
  free(big);
}

/********************************************** select_bin to select the number of bins **********************************/
//num_bin=select_bin(normalized_data, file_Len, num_dm, MIN_GRID, MAX_GRID, position, sorted_seq, all_grid_ID, &num_nonempty);
/* Determine the number of bins to use in each dimension.  Additionally sort the data elements according to the binned
 * values, and partition the data into "grids" with identical (binned) values.  We try progressively more bins until we
 * maximize a merit function, then return the results obtained using the optimal number of bins. 
 *
 * Outputs:
 * position     -- binned data values
 * sorted_seq   -- permutation vector mapping the ordered list to the original data
 * all_grid_ID  -- grid to which each data element was assigned.
 * num_nonempty -- number of distinct values assumed by all_grid_ID
 * interval     -- bin width for each data dimension
 * return value -- the number of bins selected.
 */

long select_bin(double **normalized_data, long file_Len, long num_dm, long min_grid, long max_grid, long **position, long *sorted_seq, 
                long *all_grid_ID, long *num_nonempty, double *interval, long user_num_bin)
{
 
  long num_bin=0;
  long select_num_bin=0;
  long m=0;
  long n=0;
	
  long i=0;
  long bin_scope=0;
  long temp_num_nonempty=0;

  long *temp_grid_ID;
  long *temp_sorted_seq;
  long **temp_position;

  //sorted_seq[i]=j means the event j ranks i

  double temp_index=0;
  double *bin_index;	
  double *temp_interval;
	
	
  temp_grid_ID=(long *)malloc(sizeof(long)*file_Len);
  memset(temp_grid_ID,0,sizeof(long)*file_Len);
	
  temp_sorted_seq=(long *)malloc(sizeof(long)*file_Len);
  memset(temp_sorted_seq,0,sizeof(long)*file_Len);

  temp_position=(long **)malloc(sizeof(long*)*file_Len);
  memset(temp_position,0,sizeof(long*)*file_Len);
  for (m=0;m<file_Len;m++)
    {
      temp_position[m]=(long*)malloc(sizeof(long)*num_dm);
      memset(temp_position[m],0,sizeof(long)*num_dm);
    }

  temp_interval=(double*)malloc(sizeof(double)*num_dm);
  memset(temp_interval,0,sizeof(double)*num_dm);

  bin_scope=max_grid-min_grid+1;
  bin_index=(double *)malloc(sizeof(double)*bin_scope);
  memset(bin_index,0,sizeof(double)*bin_scope);

  i=0;

  for (num_bin=min_grid;num_bin<=max_grid;num_bin++)
    {
      
	  temp_num_nonempty=0;
	   /* compute_position bins the data into num_bin bins.  Each
         dimension is binned independently.

         Outputs:
         temp_position[i][j] -- bin for the jth component of data element i.
         temp_interval[j]    -- bin-width for the jth component
      */
      compute_position(normalized_data, file_Len, num_dm, num_bin, temp_position, temp_interval);
      radixsort_flock(temp_position,file_Len,num_dm,num_bin,temp_sorted_seq,&temp_num_nonempty,temp_grid_ID);

      /* our figure of merit is the number of non-empty grids divided by number of bins per dimension.
         We declare victory when we have found a local maximum */
      bin_index[i]=((double)temp_num_nonempty)/((double)num_bin);
	  if ((double)(temp_num_nonempty)>=(double)(file_Len)*0.95)
		  break;
      if ((bin_index[i]<temp_index) && (user_num_bin==0))
         break;
	  if ((user_num_bin==num_bin-1) && (user_num_bin!=0))
		 break;
      
      /* Since we have accepted this trial bin, copy all the temporary results into
         the output buffers */
      memcpy(all_grid_ID,temp_grid_ID,sizeof(long)*file_Len);
      memcpy(sorted_seq,temp_sorted_seq,sizeof(long)*file_Len);
      memcpy(interval,temp_interval,sizeof(double)*num_dm);
		
      for (m=0;m<file_Len;m++)
        for (n=0;n<num_dm;n++)
          position[m][n]=temp_position[m][n];

      temp_index=bin_index[i];
      select_num_bin=num_bin;
      num_nonempty[0]=temp_num_nonempty;
      i++;
    }

   if ((select_num_bin<min_grid) || (select_num_bin>max_grid))
  {
    fprintf(stderr,"Number of events collected is too few in terms of number of markers used. The file should not be processed!\n"); //modified on Nov 04, 2010
	exit(0);
  }
	
  if (temp_index==0)
  {
	 fprintf(stderr,"Too many dimensions with too few events in the input file, or a too large number of bins used.\n"); //modified on July 23, 2010
	 exit(0);
  }

  free(temp_grid_ID);
  free(temp_sorted_seq);
  free(bin_index);
  free(temp_interval);
	
  for (m=0;m<file_Len;m++)
    free(temp_position[m]);
  free(temp_position);

  return select_num_bin; 
}

/************************************* Select dense grids **********************************/
// compute num_dense_grids, num_dense_events, dense_grid_reverse, and all_grid_vol
// den_cutoff=select_dense(file_Len, all_grid_ID, num_nonempty, &num_dense_grids, &num_dense_events, dense_grid_reverse);
/*
 * Prune away grids that are insufficiently "dense" (i.e., contain too few data items)
 *
 * Outputs:
 * num_dense_grids    -- number of dense grids
 * num_dense_events   -- total number of data items in all dense grids
 * dense_grid_reverse -- mapping from list of all grids to list of dense grids.
 * return value       -- density cutoff for separating dense from non-dense grids.
 */

int select_dense(long file_Len, long *all_grid_ID, long num_nonempty, long *num_dense_grids, long *num_dense_events, long *dense_grid_reverse, int den_t_event)
{
  

  long i=0;
  long vol_ID=0;
  long biggest_size=0; //biggest grid_size, to define grid_size_index
  long biggest_index=0;
  //long actual_threshold=0; //the actual threshold on grid_size, e.g., 1 implies 1 event in the grid
  //long num_remain=0; //number of remaining grids with different density thresholds
  long temp_num_dense_grids=0;
  long temp_num_dense_events=0;
	
  long *grid_size_index;
  long *all_grid_vol;
  long *grid_density_index;

  //double den_average=0;
 // double avr_index=0;
 

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // Compute all_grid_vol
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  all_grid_vol=(long *)malloc(sizeof(long)*num_nonempty);
  memset(all_grid_vol,0,sizeof(long)*num_nonempty);

  /* Grid "volume" is just the number of data contained in the grid. */
  for (i=0;i<file_Len;i++)
    {
      vol_ID=all_grid_ID[i]; //vol_ID=all_grid_ID[sorted_seq[i]];
      all_grid_vol[vol_ID]++;  //all_grid_vol[i]=j means grid i has j points
    }

 
	
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // Compute grid_size_index (histogram of grid sizes)
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  for (i=0;i<num_nonempty;i++)
    {
      if (biggest_size<all_grid_vol[i])
        {
          biggest_size=all_grid_vol[i];
        }
    }

   if (biggest_size<3)
  {
	 fprintf(stderr,"Too many dimensions with too few events in the input file, or a too large number of bins used.\n"); //modified on July 23, 2010
	 exit(0);
  }


  grid_size_index=(long*)malloc(sizeof(long)*biggest_size);
  memset(grid_size_index,0,sizeof(long)*biggest_size);
	
  for (i=0;i<num_nonempty;i++)
    {
      grid_size_index[all_grid_vol[i]-1]++; //grid_size_index[0]=5 means there are 5 grids having size 1
    }



  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // Compute den_cutoff
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	
  grid_density_index=(long *)malloc(sizeof(long)*(biggest_size-2));//from event 2 to biggest_size-1, i.e., from 0 to biggest_size-3
  memset(grid_density_index,0,sizeof(long)*(biggest_size-2));

  if (den_t_event==0)
  {
	  biggest_index=0;

	  for (i=2;i<biggest_size-1;i++) //the grid with 1 event will be skipped, i.e., grid_density_index[0] won't be defined
	  {
		  grid_density_index[i-1]=(grid_size_index[i-1]+grid_size_index[i+1]-2*grid_size_index[i]);
		  if (biggest_index<grid_density_index[i-1])
		  {
			biggest_index=grid_density_index[i-1];
			den_t_event=i+1;
		  }
	  }
  }

  if (den_t_event==0) //if biggest_size==3
  {
	 fprintf(stderr,"the densest hyperregion has only 3 events, smaller than the user-specified value. Therefore the density threshold is automatically changed to 3.\n");
	 den_t_event=3;
  }

  for (i=0;i<num_nonempty;i++)
	  if (all_grid_vol[i]>=den_t_event)
		temp_num_dense_grids++;

  if (temp_num_dense_grids<=1)
  {
	  //exit(0);
	  //printf("a too high density threshold is set! Please decrease your density threshold.\n");
	  fprintf(stderr,"a too high density threshold! Please decrease your density threshold.\n"); //modified on July 23, 2010
	  exit(0);
  }

   if (temp_num_dense_grids>=100000)
  {
	  //modified on July 23, 2010
	  //printf("a too low density threshold is set! Please increase your density threshold.\n");
	  //exit(0);
	  fprintf(stderr,"a too low density threshold! Please increase your density threshold.\n"); //modified on July 23, 2010
	  exit(0);
  }
  
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // Compute dense_grid_reverse
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  temp_num_dense_grids=0;
  temp_num_dense_events=0;
  for (i=0;i<num_nonempty;i++)
    {
      dense_grid_reverse[i]=-1;
		
      if (all_grid_vol[i]>=den_t_event) 
        {			
          dense_grid_reverse[i]=temp_num_dense_grids;  //dense_grid_reverse provides a mapping from all nonempty grids to dense grids.
          temp_num_dense_grids++;
          temp_num_dense_events+=all_grid_vol[i];						
        }
    }

  num_dense_grids[0]=temp_num_dense_grids;
  num_dense_events[0]=temp_num_dense_events;	

  free(grid_size_index);
  free(all_grid_vol);

  return den_t_event;
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Compute densegridID_To_gridevent and eventID_To_denseventID
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	grid_To_event(file_Len, dense_grid_reverse, all_grid_ID, eventID_To_denseventID, densegridID_To_gridevent);
/*
 * Filter the data so that only the data belonging to dense grids is left
 *
 * Output:
 * eventID_To_denseeventID   -- mapping from original event ID to ID in list containing only events contained in dense grids.
 * densegridID_To_gridevent  -- mapping from dense grids to prototype members of the grids.
 *
 */

void grid_To_event(long file_Len, long *dense_grid_reverse, long *all_grid_ID, long *eventID_To_denseventID, long *densegridID_To_gridevent)
{
  long i=0;
  long dense_grid_ID=0;
  long dense_event_ID=0;

  for (i=0;i<file_Len;i++)
    {
      dense_grid_ID=dense_grid_reverse[all_grid_ID[i]];
      eventID_To_denseventID[i]=-1;
      if (dense_grid_ID!=-1) //for point (i) belonging to dense grids 
        {			
          eventID_To_denseventID[i]=dense_event_ID;
          dense_event_ID++;
		
          if (densegridID_To_gridevent[dense_grid_ID]==-1) //for point i that hasn't been selected
            densegridID_To_gridevent[dense_grid_ID]=i; //densegridID_To_gridevent maps dense_grid_ID to its representative point			
        }		
    }
	
	
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Compute dense_grid_seq
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	generate_grid_seq(file_Len, num_dm, sorted_seq, num_dense_grids, densegridID_To_gridevent, position, dense_grid_rank, dense_grid_seq);
/* Construct a table of binned data values for each dense grid.
 *
 * Output:
 *
 * dense_grid_seq  -- table of binned data values for each dense grid (recall that all members of a grid have identical binned data values)
 */

void generate_grid_seq(long num_dm, long num_dense_grids, long *densegridID_To_gridevent, long **position, long **dense_grid_seq)
{  

  long i=0;
  long j=0;
  long ReventID=0; //representative event ID of the dense grid

  for (i=0;i<num_dense_grids;i++)
    {
      ReventID = densegridID_To_gridevent[i];

      for (j=0;j<num_dm;j++)
        dense_grid_seq[i][j]=position[ReventID][j];
    }
}

//compare two vectors
long compare_value(long num_dm, long *search_value, long *seq_value)
{
  long i=0;

  for (i=0;i<num_dm;i++)
    {
      if (search_value[i]<seq_value[i])
        return 1;
      if (search_value[i]>seq_value[i])
        return -1;
      if (search_value[i]==seq_value[i])
        continue;
    }
  return 0;
}

//binary search the dense_grid_seq to return the dense grid ID if it exists
long binary_search(long num_dense_grids, long num_dm, long *search_value, long **dense_grid_seq)
{

  long low=0;
  long high=0;
  long mid=0;
  long comp_result=0;
  long match=0;
  //long found=0;
	
  low = 0;
  high = num_dense_grids-1;
    
  while (low <= high) 
    {
      mid = (long)((low + high)/2);
	
      comp_result=compare_value(num_dm, search_value,dense_grid_seq[mid]);
	
		
      switch(comp_result)
        {
        case 1:
          high=mid-1;
          break;
        case -1:
          low=mid+1;
          break;
        case 0:
          match=1;
          break;
        }
      if (match==1)
        break;
    }
	


  if (match==1)
    {
      return mid;
    }
  else
    return -1;   
}


/********************************************** Computing Centers Using IDs **********************************************/

void ID2Center(double **data_in, long file_Len, long num_dm, long *eventID_To_denseventID, long num_clust, long *cluster_ID, double **population_center)
{
  long i=0;
  long j=0;
  long ID=0;
  long eventID=0;
  long *size_c;



  size_c=(long *)malloc(sizeof(long)*num_clust);
  memset(size_c,0,sizeof(long)*num_clust);

  for (i=0;i<num_clust;i++)
    for (j=0;j<num_dm;j++)
      population_center[i][j]=0;

  for (i=0;i<file_Len;i++)
    {
      eventID=eventID_To_denseventID[i];

      if (eventID!=-1) //only events in dense grids count
        {
          ID=cluster_ID[eventID];
			
          if (ID==-1)
            {
              //printf("ID==-1! in ID2Center\n");
              //exit(0);
			  fprintf(stderr,"Incorrect file format or input parameters (no dense regions found!)\n"); //modified on July 23, 2010
			  exit(0);
            }

          for (j=0;j<num_dm;j++)
            population_center[ID][j]=population_center[ID][j]+data_in[i][j];
			
          size_c[ID]++;				
        }
    }
	

  for (i=0;i<num_clust;i++)
    {
      for (j=0;j<num_dm;j++)
        if (size_c[i]!=0)
          population_center[i][j]=(population_center[i][j]/(double)(size_c[i]));
        else
		{
          population_center[i][j]=0;
		  printf("size_c[%ld]=0 in ID2center\n",i);
		}
    }

  free(size_c);

}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//  Compute Population Center with all events
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void ID2Center_all(double **data_in, long file_Len, long num_dm, long num_clust, long *cluster_ID, double **population_center)
{
  long i=0;
  long j=0;
  long ID=0;
  long *size_c;



  size_c=(long *)malloc(sizeof(long)*num_clust);
  memset(size_c,0,sizeof(long)*num_clust);

  for (i=0;i<num_clust;i++)
    for (j=0;j<num_dm;j++)
      population_center[i][j]=0;

	
  for (i=0;i<file_Len;i++)
    {
         ID=cluster_ID[i];
			
         if (ID==-1)
         {
            //printf("ID==-1! in ID2Center_all\n");
            //exit(0);
			fprintf(stderr,"Incorrect file format or input parameters (resulting in incorrect population IDs)\n"); //modified on July 23, 2010
			exit(0);
         }

         for (j=0;j<num_dm;j++)
           population_center[ID][j]=population_center[ID][j]+data_in[i][j];
			
         size_c[ID]++;        
    }
	
 
  for (i=0;i<num_clust;i++)
    {
      for (j=0;j<num_dm;j++)
        if (size_c[i]!=0)
          population_center[i][j]=(population_center[i][j]/(double)(size_c[i]));
        else
		{
          population_center[i][j]=0;
		  printf("size_c[%ld]=0 in ID2center_all\n",i);
		}
    }


  free(size_c);

}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Merge neighboring grids to clusters
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

long merge_grids(double **normalized_data, double *interval, long file_Len, long num_dm, long num_bin, long **position, long num_dense_grids, 
                 long *dense_grid_reverse, long **dense_grid_seq, long *eventID_To_denseventID, long *densegridID_To_gridevent, long *all_grid_ID,
                 long *cluster_ID, long *grid_ID, long *grid_clusterID)
{
  

  long i=0;
  long j=0;
  long t=0;
  long p=0;
  long num_clust=0;
  long ReventID=0;
  long denseID=0;
  long neighbor_ID=0;
  //long temp_grid=0;

  long *grid_value;
  long *search_value;
	
  long **graph_of_grid; //the graph constructed for dense grids: each dense grid is a graph node

  double real_dist=0;
  double **norm_grid_center;
	
  norm_grid_center=(double **)malloc(sizeof(double*)*num_dense_grids);
  memset(norm_grid_center,0,sizeof(double*)*num_dense_grids);
	
  for (i=0;i<num_dense_grids;i++)
    {
      norm_grid_center[i]=(double *)malloc(sizeof(double)*num_dm);
      memset(norm_grid_center[i],0,sizeof(double)*num_dm);
    }

  for (i=0;i<file_Len;i++)
    {
      denseID=eventID_To_denseventID[i];
      if (denseID!=-1) //only dense events can enter
        {
          grid_ID[denseID]=dense_grid_reverse[all_grid_ID[i]];
			
          if (grid_ID[denseID]==-1)
            {
              fprintf(stderr,"Incorrect input file format or input parameters (no dense region found)!\n");
              exit(0);
            }			
        }
    }

	
  /* Find centroid (in the normalized data) for each dense grid */
  ID2Center(normalized_data,file_Len,num_dm,eventID_To_denseventID,num_dense_grids,grid_ID,norm_grid_center);	
 
  //printf("pass the grid ID2 center\n");

	
  graph_of_grid=(long **)malloc(sizeof(long*)*num_dense_grids);
  memset(graph_of_grid,0,sizeof(long*)*num_dense_grids);
  for (i=0;i<num_dense_grids;i++)
    {
      graph_of_grid[i]=(long *)malloc(sizeof(long)*num_dm);
      memset(graph_of_grid[i],0,sizeof(long)*num_dm);		
		

      for (j=0;j<num_dm;j++)
        graph_of_grid[i][j]=-1;
    }	
	
  grid_value=(long *)malloc(sizeof(long)*num_dm);
  memset(grid_value,0,sizeof(long)*num_dm);

  search_value=(long *)malloc(sizeof(long)*num_dm);
  memset(search_value,0,sizeof(long)*num_dm);

  
  for (i=0;i<num_dense_grids;i++)
    {
      ReventID=densegridID_To_gridevent[i];

      for (j=0;j<num_dm;j++)
        {
          grid_value[j]=position[ReventID][j];
          
        }
     

      /* For each dimension, find the neighbor, if any, that is equal in all other dimensions and 1 greater in
         the chosen dimension.  If the neighbor's centroid is not too far away, add it to this grid's neighbor
         list. */
      for (t=0;t<num_dm;t++)
        {
          for (p=0;p<num_dm;p++)
            search_value[p]=grid_value[p];

          if (grid_value[t]==num_bin-1)
            continue;

          search_value[t]=grid_value[t]+1; //we only consider the neighbor at the bigger side

          neighbor_ID=binary_search(num_dense_grids, num_dm, search_value, dense_grid_seq);
			
          if (neighbor_ID!=-1) 
            {
              real_dist=norm_grid_center[i][t]-norm_grid_center[neighbor_ID][t];
	
              if (real_dist<0)
                real_dist=-real_dist;
				
              if (real_dist<2*interval[t]) //by using 2*interval, this switch is void
                graph_of_grid[i][t]=neighbor_ID;			
            }
        }
      grid_clusterID[i]=i; //initialize grid_clusterID
    } 
	
	
  //graph constructed
  //DFS-based search begins

  /* Use a depth-first search to construct a list of connected subgraphs (= "clusters").
     Note our graph as we have constructed it is a DAG, so we can use that to our advantage
     in our search. */
  //  num_clust=dfs(graph_of_grid,num_dense_grids,num_dm,grid_clusterID);
  num_clust=find_connected(graph_of_grid, num_dense_grids, num_dm, grid_clusterID);

	
  //computes grid_ID and cluster_ID
  for (i=0;i<file_Len;i++)
    {
      denseID=eventID_To_denseventID[i];
      if (denseID!=-1) //only dense events can enter
	  {
        cluster_ID[denseID]=grid_clusterID[grid_ID[denseID]];
		//if (cluster_ID[denseID]==-1)
		//	printf("catch you!\n");
	  }
    }
	
  free(search_value);
  free(grid_value);

  for (i=0;i<num_dense_grids;i++)
    {
      free(graph_of_grid[i]);
      free(norm_grid_center[i]);
    }
  free(graph_of_grid);
  free(norm_grid_center);

  return num_clust;
}

/********************************************* Merge Clusters to Populations *******************************************/

long merge_clusters(long num_clust, long num_dm, double *interval, double **cluster_center, long *cluster_populationID, long max_num_pop)
{
  long num_population=0;
  long temp_num_population=0;

  long i=0;
  long j=0;
  long t=0;
  long merge=0;
  long smid=0;
  long bgid=0;
  double merge_dist=1.1;

  long *map_ID;

  double diff=0;  

  map_ID=(long*)malloc(sizeof(long)*num_clust);
  memset(map_ID,0,sizeof(long)*num_clust);

  for (i=0;i<num_clust;i++)
  {
    	cluster_populationID[i]=i;
		map_ID[i]=i;
  }
    

  while (((num_population>max_num_pop) && (merge_dist<5.0)) || ((num_population<=1) && (merge_dist>0.1)))
  {
  
	  if (num_population<=1)
	  	  merge_dist=merge_dist-0.1;

	  if (num_population>max_num_pop)
          merge_dist=merge_dist+0.1;
	  
	    
	 for (i=0;i<num_clust;i++)
		cluster_populationID[i]=i;

    for (i=0;i<num_clust-1;i++)
    {
      for (j=i+1;j<num_clust;j++)
        {
          merge=1;
			
          for (t=0;t<num_dm;t++)
            {
              diff=cluster_center[i][t]-cluster_center[j][t];
				
              if (diff<0)
                diff=-diff;
              if (diff>(merge_dist*interval[t]))
                merge=0;
            }

          if ((merge) && (cluster_populationID[i]!=cluster_populationID[j]))
            {
              if (cluster_populationID[i]<cluster_populationID[j])  //they could not be equal
                {
                  smid = cluster_populationID[i];
                  bgid = cluster_populationID[j];
                }
              else
                {
                  smid = cluster_populationID[j];
                  bgid = cluster_populationID[i];
                }
              for (t=0;t<num_clust;t++)
                {
                  if (cluster_populationID[t]==bgid)
                    cluster_populationID[t]=smid;
                }
            }
        }
    }

  

  for (i=0;i<num_clust;i++)
    map_ID[i]=-1;

  for (i=0;i<num_clust;i++)
    map_ID[cluster_populationID[i]]=1;

  num_population=0;
  for (i=0;i<num_clust;i++)
    {
      if (map_ID[i]==1)
        {
          map_ID[i]=num_population;
          num_population++;
        }
    }

  if ((temp_num_population>max_num_pop) && (num_population==1))
	  break;
  else
	  temp_num_population=num_population;

  if (num_clust<=1)
	break;
  } //end while

  for (i=0;i<num_clust;i++)
    cluster_populationID[i]=map_ID[cluster_populationID[i]];
	
  free(map_ID);

  return num_population; 
}

///////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double kmeans(double **Matrix, long k, double kmean_term, long file_Len, long num_dm, long *shortest_id, double **center)
{
	
	long i=0;
	long j=0;
	long t=0;
	long random=0;
	long random1=0;
	long random2=0;
	long times=0;
	long dist_used=0; //0 is Euclidean and 1 is Pearson
	long random_init=0; //0: not use random seeds; 
	long real_Len=0;
	long skipped=0;
	
 	long *num;  //num[i]=t means the ith cluster has t points
	
	double vvv=1.0; // the biggest variation;
	double distance=0.0;
	double xv=0.0;
	double variation=0.0;

	double mean_dx=0;
	double mean_dy=0;
	double sum_var=0;
	double dx=0;
	double dy=0;
	double sd_x=0;
	double sd_y=0;	
	double diff=0;
	double distortion=0;
	double sd=0; //standard deviation
	double shortest_distance;

	double *temp_center;	
			
	double **sum;	
 
	temp_center = (double *)malloc(sizeof(double)*num_dm);
	memset(temp_center,0,sizeof(double)*num_dm);

	if (random_init)
	{
		for (i=0;i<k;i++)
		{	
			random1=rand()*rand();
			random2=abs((random1%5)+1);
			for (t=0;t<random2;t++)
				random2=random2*rand()+rand();
	
			random=abs(random2%file_Len);
			for (j=0;j<num_dm;j++)
				center[i][j]=Matrix[random][j];				
		}
	}


	num = (long *)malloc(sizeof(long)*k);
	memset(num,0,sizeof(long)*k);

	sum = (double **)malloc(sizeof(double*)*k);
	memset(sum,0,sizeof(double*)*k);
	for (i=0;i<k;i++)
	{
		sum[i] = (double *)malloc(sizeof(double)*num_dm);
		memset(sum[i],0,sizeof(double)*num_dm);
	}


	times=0;
	real_Len=0;

	while (((vvv>kmean_term) && (kmean_term<1)) || ((times<kmean_term) && (kmean_term>=1)))
	{
		for (i=0;i<k;i++)
		{
			num[i]=0;
			for (j=0;j<num_dm;j++)
				sum[i][j]=0.0;  
		}
		
		for (i=0;i<file_Len;i++)  //for each data point i, we compute the distance between Matrix[i] and center[j]
		{		
			skipped = 0;
			shortest_distance=MAX_VALUE;
			for (j=0;j<k;j++)  //for each center j
			{
				distance=0.0;
								
				if (dist_used==0)  //Euclidean distance
				{
					for (t=0;t<num_dm;t++) //for each dimension here num_dm is always 1 as we consider individual dimensions
					{
					
						diff=center[j][t]-Matrix[i][t];
					
						diff=diff*diff;  
						//this K-means is only used for dimension selection, and therefore for 1-dimensional data only, no need to use cube distance

						distance = distance+diff; //here we have a weight for each dimension
					}
				}
				else  //pearson correlation
				{
					mean_dx=0.0;
					mean_dy=0.0;
					sum_var=0.0;
					dx=0.0;
					dy=0.0;
					sd_x=0.0;
					sd_y=0.0;
					for (t=0;t<num_dm;t++)
					{
						mean_dx+=center[j][t];
						mean_dy+=Matrix[i][t];
					}
					mean_dx=mean_dx/(double)num_dm;
					mean_dy=mean_dy/(double)num_dm;
			
					for (t=0;t<num_dm;t++)
					{
						dx=center[j][t]-mean_dx;
						dy=Matrix[i][t]-mean_dy;
						sum_var+=dx*dy;
						sd_x+=dx*dx;
						sd_y+=dy*dy;
					}
					if (sqrt(sd_x*sd_y)==0)
						distance = 1.0;
					else
						distance = 1.0 - (sum_var/(sqrt(sd_x*sd_y))); // distance ranges from 0 to 2;
					//printf("distance=%f\n",distance);			
				}	//pearson correlation ends 				

				if ((distance<shortest_distance) && (skipped == 0))
				{
					shortest_distance=distance;						
					shortest_id[i]=j;  
				}			
			}//end for j
				real_Len++;
				num[shortest_id[i]]=num[shortest_id[i]]+1;
				for (t=0;t<num_dm;t++)
					sum[shortest_id[i]][t]=sum[shortest_id[i]][t]+Matrix[i][t];		
		}//end for i
	/* recompute the centers */
	//compute_mean(group);		
		vvv=0.0;
		for (j=0;j<k;j++)
		{
			memcpy(temp_center,center[j],sizeof(double)*num_dm);
			variation=0.0;
			if (num[j]!=0)
			{
				for (t=0;t<num_dm;t++)
				{
					center[j][t]=sum[j][t]/(double)num[j];
					xv=(temp_center[t]-center[j][t]);
					variation=variation+xv*xv;
				}
			}
			
			if (variation>vvv)
				vvv=variation;  //vvv is the biggest variation among the k clusters;			
		}
	//compute_variation;
		times++;
	} //end for while




	free(num);

	for (i=0;i<k;i++)
		free(sum[i]);
	free(sum);
	free(temp_center);


	return distortion;
		
}

//////////////////////////////////////////////////////
/*************************** Show *****************************/
void show(double **Matrix, long *cluster_id, long file_Len, long k, long num_dm, char *name_string)
{
	int situ1=0;
	int situ2=0;

	long i=0;
	long id=0;
	long j=0;
	long info_id=0;
	long nearest_id=0;
	long insert=0;
	long temp=0;
	long m=0;
	long n=0;
	long t=0;
	
	long *size_c;
	


	long **size_mybound_1;
	long **size_mybound_2;
	long **size_mybound_3;
	long **size_mybound_0;

	double interval=0.0;

	double *big;
	double *small;


	double **center;
	double **mybound;
	
	long **prof; //prof[i][j]=1 means population i is + at parameter j
	
	FILE *fpcnt_id; //proportion id
	//FILE *fcent_id; //center_id, i.e., centers of clusters within the original data
	FILE *fprof_id; //profile_id

	big=(double *)malloc(sizeof(double)*num_dm);
	memset(big,0,sizeof(double)*num_dm);

	small=(double *)malloc(sizeof(double)*num_dm);
	memset(small,0,sizeof(double)*num_dm);

	for (i=0;i<num_dm;i++)
	{
		big[i]=0.0;
		small[i]=(double)MAX_VALUE;
	}
	
	
	size_c=(long *)malloc(sizeof(long)*k);
	memset(size_c,0,sizeof(long)*k);

	center=(double**)malloc(sizeof(double*)*k);
	memset(center,0,sizeof(double*)*k);
	for (i=0;i<k;i++)
	{
		center[i]=(double*)malloc(sizeof(double)*num_dm);
		memset(center[i],0,sizeof(double)*num_dm);
	}

	mybound=(double**)malloc(sizeof(double*)*num_dm);
	memset(mybound,0,sizeof(double*)*num_dm);
	for (i=0;i<num_dm;i++) //there are 3 mybounds for 4 categories
	{
		mybound[i]=(double*)malloc(sizeof(double)*3);
		memset(mybound[i],0,sizeof(double)*3);
	}

	prof=(long **)malloc(sizeof(long*)*k);
	memset(prof,0,sizeof(long*)*k);
	for (i=0;i<k;i++)
	{
		prof[i]=(long *)malloc(sizeof(long)*num_dm);
		memset(prof[i],0,sizeof(long)*num_dm);
	}


	for (i=0;i<file_Len;i++)
	{
		id=cluster_id[i];
		for (j=0;j<num_dm;j++)
		{
			center[id][j]=center[id][j]+Matrix[i][j];
			if (big[j]<Matrix[i][j])
				big[j]=Matrix[i][j];
			if (small[j]>Matrix[i][j])
				small[j]=Matrix[i][j];
		}
		
		size_c[id]++;		
	}

	for (i=0;i<k;i++)
		for (j=0;j<num_dm;j++)
		{			
			if (size_c[i]!=0)
				center[i][j]=(center[i][j]/(double)(size_c[i]));
			else
				center[i][j]=0;	
		}

	for (j=0;j<num_dm;j++)
	{
		interval=((big[j]-small[j])/4.0);
		//printf("interval[%d] is %f\n",j,interval);
		for (i=0;i<3;i++)
			mybound[j][i]=small[j]+((double)(i+1)*interval);
	}
	

	size_mybound_0=(long **)malloc(sizeof(long*)*k);
	memset(size_mybound_0,0,sizeof(long*)*k);
	
	for (i=0;i<k;i++)
	{
		size_mybound_0[i]=(long*)malloc(sizeof(long)*num_dm);
		memset(size_mybound_0[i],0,sizeof(long)*num_dm);		
	}

	size_mybound_1=(long **)malloc(sizeof(long*)*k);
	memset(size_mybound_1,0,sizeof(long*)*k);
	
	for (i=0;i<k;i++)
	{
		size_mybound_1[i]=(long*)malloc(sizeof(long)*num_dm);
		memset(size_mybound_1[i],0,sizeof(long)*num_dm);		
	}

	size_mybound_2=(long **)malloc(sizeof(long*)*k);
	memset(size_mybound_2,0,sizeof(long*)*k);
	
	for (i=0;i<k;i++)
	{
		size_mybound_2[i]=(long*)malloc(sizeof(long)*num_dm);
		memset(size_mybound_2[i],0,sizeof(long)*num_dm);	
	}

	size_mybound_3=(long **)malloc(sizeof(long*)*k);
	memset(size_mybound_3,0,sizeof(long*)*k);

	for (i=0;i<k;i++)
	{
		size_mybound_3[i]=(long*)malloc(sizeof(long)*num_dm);
		memset(size_mybound_3[i],0,sizeof(long)*num_dm);
	}
	
	for (i=0;i<file_Len;i++)
		for (j=0;j<num_dm;j++)
			{
				if (Matrix[i][j]<mybound[j][0])// && ((Matrix[i][j]-small[j])>0)) //the smallest values excluded
					size_mybound_0[cluster_id[i]][j]++;
				else
				{
					if (Matrix[i][j]<mybound[j][1])
						size_mybound_1[cluster_id[i]][j]++;
					else
					{
						if (Matrix[i][j]<mybound[j][2])
							size_mybound_2[cluster_id[i]][j]++;
						else
							//if (Matrix[i][j]!=big[j]) //the biggest values excluded
								size_mybound_3[cluster_id[i]][j]++;
					}

				}
			}

	fprof_id=fopen("profile.txt","w");
	fprintf(fprof_id,"Population_ID\t");
	fprintf(fprof_id,"%s\n",name_string);
	
	for (i=0;i<k;i++)
	{
		fprintf(fprof_id,"%ld\t",i+1); //i changed to i+1 to start from 1 instead of 0: April 16, 2009
		for (j=0;j<num_dm;j++)
		{
			
			if (size_mybound_0[i][j]>size_mybound_1[i][j])
				situ1=0;
			else
				situ1=1;
			if (size_mybound_2[i][j]>size_mybound_3[i][j])
				situ2=2;
			else
				situ2=3;

			if ((situ1==0) && (situ2==2))
			{
				if (size_mybound_0[i][j]>size_mybound_2[i][j])
					prof[i][j]=0;
				else
					prof[i][j]=2;
			}
			if ((situ1==0) && (situ2==3))
			{
				if (size_mybound_0[i][j]>size_mybound_3[i][j])
					prof[i][j]=0;
				else
					prof[i][j]=3;
			}
			if ((situ1==1) && (situ2==2))
			{
				if (size_mybound_1[i][j]>size_mybound_2[i][j])
					prof[i][j]=1;
				else
					prof[i][j]=2;
			}
			if ((situ1==1) && (situ2==3))
			{
				if (size_mybound_1[i][j]>size_mybound_3[i][j])
					prof[i][j]=1;
				else
					prof[i][j]=3;
			}
			
			//begin to output profile
			if (j==num_dm-1)
			{
				if (prof[i][j]==0)
					fprintf(fprof_id,"1\n");
				if (prof[i][j]==1)
					fprintf(fprof_id,"2\n");
				if (prof[i][j]==2)
					fprintf(fprof_id,"3\n");
				if (prof[i][j]==3)
					fprintf(fprof_id,"4\n");
			}
			else
			{
				if (prof[i][j]==0)
					fprintf(fprof_id,"1\t");
				if (prof[i][j]==1)
					fprintf(fprof_id,"2\t");
				if (prof[i][j]==2)
					fprintf(fprof_id,"3\t");
				if (prof[i][j]==3)
					fprintf(fprof_id,"4\t");
			}
		}
	}
	fclose(fprof_id);

	///////////////////////////////////////////////////////////
	

	fpcnt_id=fopen("percentage.txt","w");
	fprintf(fpcnt_id,"Population_ID\tPercentage\n");

	for (t=0;t<k;t++)
	{
		fprintf(fpcnt_id,"%ld\t%.2f\n",t+1,(double)size_c[t]*100.0/(double)file_Len);	//t changed to t+1 to start from 1 instead of 0: April 16, 2009									
	}
	fclose(fpcnt_id);

	free(big);
	free(small);
	free(size_c);

	for (i=0;i<k;i++)
	{
		free(center[i]);
		free(prof[i]);
		free(size_mybound_0[i]);
		free(size_mybound_1[i]);
		free(size_mybound_2[i]);
		free(size_mybound_3[i]);
	}
	free(center);
	free(prof);
	free(size_mybound_0);
	free(size_mybound_1);
	free(size_mybound_2);
	free(size_mybound_3);

	for (i=0;i<num_dm;i++)
		free(mybound[i]);
	free(mybound);
	
}
///////////////////////////////////////////////////////////////////////////

double get_avg_dist(double *population_center, long smaller_pop_ID, long larger_pop_ID, long *population_ID, long num_population, long file_Len, 
				   long num_dm, double **normalized_data, long d1, long d2, long d3, long *size_c)
{
	long k=0;
	long temp=0;
	long i=0;
	long j=0;
	long t=0;
	long current_biggest=0;


	double total_dist=0.0;
	double distance=0.0;
	double dist=0.0;
	double biggest_distance=0.0;

	long *nearest_neighbors;
	double *nearest_distance;

	k=min(size_c[smaller_pop_ID],size_c[larger_pop_ID]);
	
	if (k>0)
	{
		k=(long)sqrt((double)k);
		if (num_dm<=3)
			k=(long)(2.7183*(double)k);  //e*k for low-dimensional space, added Nov. 4, 2010
		//printf("the k-nearest k is %d\n",k);
	}
	else
	{
		printf("error in k\n");
		exit(0);
	}

	nearest_neighbors=(long *)malloc(sizeof(long)*k);
	memset(nearest_neighbors,0,sizeof(long)*k);

	nearest_distance=(double *)malloc(sizeof(double)*k);
	memset(nearest_distance,0,sizeof(double)*k);

	temp=0;

	for (i=0;i<file_Len;i++)
	{
		if ((population_ID[i]==smaller_pop_ID) || (population_ID[i]==larger_pop_ID)) //the event belongs to the pair of populations
		{
			distance=0.0;
			for (t=0;t<num_dm;t++)
			{
				if ((t!=d1) && (t!=d2) && (t!=d3))
					continue;
				else
				{
					dist=population_center[t]-normalized_data[i][t];
					distance=distance+dist*dist;
				}
			}

			if (temp<k)
			{
				nearest_neighbors[temp]=i;
				nearest_distance[temp]=distance;
				temp++;
			}
			else
			{
			    biggest_distance=0.0;
				for (j=0;j<k;j++)
				{
					if (nearest_distance[j]>biggest_distance)
					{
						biggest_distance=nearest_distance[j];
						current_biggest=j;
					}
				}

				if (biggest_distance>distance)
				{
					nearest_neighbors[current_biggest]=i;
					nearest_distance[current_biggest]=distance;
				}
			}
		}
	}

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

	total_dist=total_dist/(double)k;

	free(nearest_distance);
	free(nearest_neighbors);
	
	return total_dist;
}


/******************************************************** Main Function **************************************************/
//for normalized data, there are five variables:
//cluster_ID
//population_center
//grid_clusterID
//grid_ID
//grid_Center

//the same five variables exist for the original data
//however, the three IDs (cluster_ID, grid_ID, grid_clusterID) don't change in either normalized or original data
//also, data + cluster_ID -> population_center
//data + grid_ID -> grid_Center

/* what is the final output */
//the final things we want are grid_Center in the original data and grid_clusterID //or population_center in the original data
//Sparse grids will not be considered when computing the two centroids (centroids of grids and centroids of clusters)

/*  what information should select_bin output */
//the size of all IDs are unknown to function main because we only consider the events in dense grids, and also the number of dense grids
//is unknown, therefore I must use a prescreening to output 
//how many bins I should use
//the number of dense grids
//total number of events in the dense grids

/* basic procedure of main function */ 
//1. read raw file and normalize the raw file
//2. select_bin
//3. allocate memory for eventID_To_denseventID, grid_clusterID, grid_ID, cluster_ID. 
//4. call select_dense and merge_grids with grid_clusterID, grid_ID, cluster_ID.
//5. release normalized data; allocate memory for grid_Center and population_center
//6. output grid_Center and population_center using ID2Center together with grid_clusterID //from IDs to centers

int main (int argc, char **argv)
{
  //inputs
  FILE *f_src; //source file pointer
 
  FILE *f_out; //coordinates
  FILE *f_cid; //population-ID of events
  FILE *f_ctr; //centroids of populations
  FILE *f_results; //coordinates file event and population column
  FILE *f_mfi; //added April 16, 2009 for mean fluorescence intensity
  FILE *f_parameters; //number of bins and density calculated by
                      //the algorithm. Used to update the database
  FILE *f_properties; //Properties file used by Image generation software

  //char tmpbuf[128];

  char para_name_string[LINE_LEN];

  long time_ID=-1;
  long num_bin=0; //the bins I will use on each dimension
  long num_pop=0;	
  long file_Len=0; //number of events		
  long num_dm=0;
  long num_clust=0;
  long num_dense_events=0;
  long num_dense_grids=0;
  long num_nonempty=0;
  long num_population=0;
  long num_real_pop=0;
  long keep_merge=0;
  long num_checked_range=0;

  //below are read from configuration file
  long i=0;
  long j=0;
  long p=0;
  long t=0;
  long q=0;
  long temp_i=0;
  long temp_j=0;
  long first_i=0;
  long first_j=0;
  long d1=0;
  long d2=0;
  long d3=0;
  long d_d=0;
  long max_num_pop=0; //upperbound of the number of populations that is equal to 2^d
  long index_id=0;
  long num_computed_population=0;
  long tot_size=0;

  int den_t_event=0;

  double distance=0.0;
  double dist=0.0;
  double current_smallest_dist=0;
  double max_d_dist=0.0;
  double Ei=0.0;
  double Ej=0.0;
  double E1=0.0;
  double E2=0.0;
  double E3=0.0;
  double Ep=0.0;
  double temp=0.0;
  double tmp=0.0;

  long *grid_clusterID; //(dense)gridID to clusterID
  long *grid_ID; //(dense)eventID to gridID
  long *cluster_ID; //(dense)eventID to clusterID
  long *eventID_To_denseventID; //eventID_To_denseventID[i]=k means event i is in a dense grid and its ID within dense events is k
  long *all_grid_ID; //gridID for all events
  long *densegridID_To_gridevent;
  long *sorted_seq;
  long *dense_grid_reverse;
  long *population_ID; //denseeventID to populationID
  long *cluster_populationID; //clusterID to populationID
  long *grid_populationID; //gridID to populationID
  long *all_population_ID; //populationID of event
  long *size_c;
  long *size_p_2;
  long *partit;
  long *size_p;
  long *temp_size_j;
  long *all_computed_pop_ID;


  long **position;
  long **dense_grid_seq;
  long **temp_population_ID;

  double *interval;
  double *center_1;
  double *center_2;
  double *center_3;
  double *aver;
  double *std;
  double *center_p_1;
  double *center_p_2;
  double *center_p_3;
	
  double **population_center; //population centroids in the raw/original data
  double **cluster_center; //cluster centroids in the raw/original data
  double **new_population_center;
  double **temp_population_center;
  double **min_pop;
  double **max_pop;

  double **input_data;
  double **normalized_data;

  double **real_pop_center;
  double **distance_pop;

  double ***ind_pop;
  double ***ind_pop_3;
		
  int min = 999999;
  int max = 0;

  printf( "Starting time:\t\t\t\t");
  fflush(stdout);
  system("/bin/date");

  /*
   * Windows Version
  _strtime( tmpbuf );
  printf( "Starting time:\t\t\t\t%s\n", tmpbuf );
  _strdate( tmpbuf );
  printf( "Starting date:\t\t\t\t%s\n", tmpbuf );
  */

  /////////////////////////////////////////////////////////////

  if ((argc!=2) && (argc!=3) && (argc!=4) && (argc!=5) && (argc!=6))
  {
      fprintf(stderr,"Incorrect number of input parameters!\n"); 
	  fprintf(stderr,"usage:\n");
	  fprintf(stderr,"basic mode: flock data_file\n");
	  fprintf(stderr,"advanced mode 0 (specify maximum # of pops): flock data_file max_num_pop\n");
	  fprintf(stderr,"advanced mode 1 (without # of pops): flock data_file num_bin density_index\n");
	  fprintf(stderr,"advanced mode 2 (specify # of pops): flock data_file num_bin density_index number_of_pop\n");
	  fprintf(stderr,"advanced mode 3 (specify both # of pops): flock data_file num_bin density_index number_of_pop max_num_pop\n");
      exit(0);
  }	

  f_src=fopen(argv[1],"r");

  if (argc==3)
  {
	max_num_pop=atoi(argv[2]);
	printf("num of maximum pop is %ld\n",max_num_pop);
  }
  
  if (argc==4)
  {
	num_bin=atoi(argv[2]);
	printf("num_bin is %ld\n",num_bin);

	den_t_event=atoi(argv[3]);
	printf("density_index is %d\n",den_t_event);
  }

  if (argc==5)
  {
	num_bin=atoi(argv[2]);
	printf("num_bin is %ld\n",num_bin);

	den_t_event=atoi(argv[3]);
	printf("density_index is %d\n",den_t_event);

	num_pop=atoi(argv[4]);
	printf("num of pop is %ld\n",num_pop);
  }

  if (argc==6)
  {
	num_bin=atoi(argv[2]);
	printf("num_bin is %ld\n",num_bin);

	den_t_event=atoi(argv[3]);
	printf("density_index is %d\n",den_t_event);

	num_pop=atoi(argv[4]);
	printf("num of pop is %ld\n",num_pop);

	max_num_pop=atoi(argv[5]);
	printf("num of pop is %ld\n",max_num_pop);
  }


  if (num_pop==1)
  {
	  printf("it is not allowed to specify only one population\n");
	  exit(0);
  }

    

  getfileinfo(f_src, &file_Len, &num_dm, para_name_string, &time_ID); //get the filelength, number of dimensions, and num/name of parameters

  
  if (max_num_pop==0)
  {
	  max_num_pop=(long)pow(2,num_dm);
		if (max_num_pop>MAX_POP_NUM)
			max_num_pop=MAX_POP_NUM;
  }

  /************************************************* Read the data *****************************************************/
	
  rewind(f_src); //reset the file pointer	

  input_data = (double **)malloc(sizeof(double*)*file_Len);
  memset(input_data,0,sizeof(double*)*file_Len);
  for (i=0;i<file_Len;i++)
  {
     input_data[i]=(double *)malloc(sizeof(double)*num_dm);
     memset(input_data[i],0,sizeof(double)*num_dm);
  }
	
  readsource(f_src, file_Len, num_dm, input_data, time_ID); //read the data;
	
  fclose(f_src);

  normalized_data=(double **)malloc(sizeof(double*)*file_Len);
  memset(normalized_data,0,sizeof(double*)*file_Len);
  for (i=0;i<file_Len;i++)
    {
      normalized_data[i]=(double *)malloc(sizeof(double)*num_dm);
      memset(normalized_data[i],0,sizeof(double)*num_dm);
    }
	
  tran(input_data, file_Len, num_dm, NORM_METHOD, normalized_data);
	

  position=(long **)malloc(sizeof(long*)*file_Len);
  memset(position,0,sizeof(long*)*file_Len);
  for (i=0;i<file_Len;i++)
    {
      position[i]=(long*)malloc(sizeof(long)*num_dm);
      memset(position[i],0,sizeof(long)*num_dm);
    }

  all_grid_ID=(long *)malloc(sizeof(long)*file_Len);
  memset(all_grid_ID,0,sizeof(long)*file_Len);

  sorted_seq=(long*)malloc(sizeof(long)*file_Len);
  memset(sorted_seq,0,sizeof(long)*file_Len);
	
  interval=(double*)malloc(sizeof(double)*num_dm);
  memset(interval,0,sizeof(double)*num_dm);

  /************************************************* select_bin *************************************************/
	
  if (num_bin>=1)  //num_bin has been selected by user
  {
  	select_bin(normalized_data, file_Len, num_dm, MIN_GRID, MAX_GRID, position, sorted_seq, all_grid_ID, &num_nonempty, interval,num_bin);
	printf("user set bin number is %ld\n",num_bin); 
  }
  else  //num_bin has not been selected by user
  {
	num_bin=select_bin(normalized_data, file_Len, num_dm, MIN_GRID, MAX_GRID, position, sorted_seq, all_grid_ID, &num_nonempty, interval,num_bin);
	printf("selected bin number is %ld\n",num_bin);    
  }
  printf("number of non-empty grids is %ld\n",num_nonempty);
		
 

  /* Although we return sorted_seq from select_bin(), we don't use it for anything, except possibly diagnostics */
  free(sorted_seq);
	
	
  dense_grid_reverse=(long*)malloc(sizeof(long)*num_nonempty);
  memset(dense_grid_reverse,0,sizeof(long)*num_nonempty);	

  /************************************************* select_dense **********************************************/

  if (den_t_event>=1) //den_t_event must be larger or equal to 2 if the user wants to set it
	den_t_event=select_dense(file_Len, all_grid_ID, num_nonempty, &num_dense_grids, &num_dense_events, dense_grid_reverse, den_t_event);
  else
  {
	den_t_event=select_dense(file_Len, all_grid_ID, num_nonempty, &num_dense_grids, &num_dense_events, dense_grid_reverse, den_t_event);
	printf("automated selected density threshold is %d\n",den_t_event);
  }		

  printf("Number of dense grids is %ld\n",num_dense_grids);

  densegridID_To_gridevent = (long *)malloc(sizeof(long)*num_dense_grids);
  memset(densegridID_To_gridevent,0,sizeof(long)*num_dense_grids);

  for (i=0;i<num_dense_grids;i++)
    densegridID_To_gridevent[i]=-1; //initialize all densegridID_To_gridevent values to -1
	

  eventID_To_denseventID=(long *)malloc(sizeof(long)*file_Len);
  memset(eventID_To_denseventID,0,sizeof(long)*file_Len);     //eventID_To_denseventID[i]=k means event i is in a dense grid and its ID within dense events is k


  grid_To_event(file_Len, dense_grid_reverse, all_grid_ID, eventID_To_denseventID, densegridID_To_gridevent);

	
  dense_grid_seq=(long **)malloc(sizeof(long*)*num_dense_grids);
  memset(dense_grid_seq,0,sizeof(long*)*num_dense_grids);
  for (i=0;i<num_dense_grids;i++)
    {
      dense_grid_seq[i]=(long *)malloc(sizeof(long)*num_dm);
      memset(dense_grid_seq[i],0,sizeof(long)*num_dm);
    }


  /* Look up the binned data values for each dense grid */
  generate_grid_seq(num_dm, num_dense_grids, densegridID_To_gridevent, position, dense_grid_seq); 	
	
	
  /************************************************* allocate memory *********************************************/
	
  grid_clusterID=(long *)malloc(sizeof(long)*num_dense_grids);
  memset(grid_clusterID,0,sizeof(long)*num_dense_grids);

  grid_ID=(long *)malloc(sizeof(long)*num_dense_events);
  memset(grid_ID,0,sizeof(long)*num_dense_events);

  cluster_ID=(long *)malloc(sizeof(long)*num_dense_events);
  memset(cluster_ID,0,sizeof(long)*num_dense_events);


  /*********************************************** merge_grids ***********************************************/
  //long merge_grids(long file_Len, long num_dm, long num_bin, long **position, long num_dense_grids, long *dense_grid_rank, long *dense_grid_reverse,
  //			 long **dense_grid_seq, long *eventID_To_denseventID, long *densegridID_To_gridevent, long *all_grid_ID,
  //			 long *cluster_ID, long *grid_ID, long *grid_clusterID)
	
  num_clust = merge_grids(normalized_data, interval, file_Len, num_dm, num_bin, position, num_dense_grids, dense_grid_reverse, dense_grid_seq, eventID_To_denseventID, densegridID_To_gridevent, all_grid_ID, cluster_ID, grid_ID, grid_clusterID);
	
  printf("computed number of groups is %ld\n",num_clust);	

	
  /************************************** release unnecessary memory and allocate memory and compute centers **********************************/
	
	
  for (i=0;i<file_Len;i++)
    free(position[i]);
  free(position);

  for (i=0;i<num_dense_grids;i++)
    free(dense_grid_seq[i]);
  free(dense_grid_seq);

  free(dense_grid_reverse);
	
  free(densegridID_To_gridevent);
  free(all_grid_ID);
	
  // cluster_center ////////////////////////////////////////////////////////////////////////////////////////////////////////
	
  cluster_center=(double **)malloc(sizeof(double*)*num_clust);
  memset(cluster_center,0,sizeof(double*)*num_clust);
  for (i=0;i<num_clust;i++)
  {
     cluster_center[i]=(double*)malloc(sizeof(double)*num_dm);
     memset(cluster_center[i],0,sizeof(double)*num_dm);
  }
	
  ID2Center(normalized_data,file_Len,num_dm,eventID_To_denseventID,num_clust,cluster_ID,cluster_center); //produce the centers with normalized data
	
  //printf("pass the first ID2center\n");

  /*** population_ID and grid_populationID **/
	
  cluster_populationID=(long*)malloc(sizeof(long)*num_clust);
  memset(cluster_populationID,0,sizeof(long)*num_clust);

  grid_populationID=(long*)malloc(sizeof(long)*num_dense_grids);
  memset(grid_populationID,0,sizeof(long)*num_dense_grids);

  population_ID=(long*)malloc(sizeof(long)*num_dense_events);
  memset(population_ID,0,sizeof(long)*num_dense_events);

  num_population = merge_clusters(num_clust, num_dm, interval, cluster_center, cluster_populationID, max_num_pop);

  

  for (i=0;i<num_clust;i++)
    free(cluster_center[i]);
  free(cluster_center);

  free(interval);
	
  for (i=0;i<num_dense_grids;i++)
    {
      grid_populationID[i]=cluster_populationID[grid_clusterID[i]];
    }

  for (i=0;i<num_dense_events;i++)
    {
      population_ID[i]=cluster_populationID[cluster_ID[i]];
    }

  printf("Number of merged groups before second-run partitioning is %ld\n",num_population);

  

  // population_center /////////////////////////////////////////////////////////////////////////////////////////////////////

 
  population_center=(double **)malloc(sizeof(double*)*num_population);
  memset(population_center,0,sizeof(double*)*num_population);
  for (i=0;i<num_population;i++)
  {
     population_center[i]=(double*)malloc(sizeof(double)*num_dm);
     memset(population_center[i],0,sizeof(double)*num_dm);
  }

  
	
  ID2Center(normalized_data,file_Len,num_dm,eventID_To_denseventID,num_population,population_ID,population_center); //produce population centers with normalized data
	

  // show ////////////////////////////////////////////////////////////////////////////////
  all_population_ID=(long*)malloc(sizeof(long)*file_Len);
  memset(all_population_ID,0,sizeof(long)*file_Len);

  kmeans(normalized_data, num_population, KMEANS_TERM, file_Len, num_dm, all_population_ID, population_center);

  for (i=0;i<num_population;i++)
	free(population_center[i]);

  free(population_center);
  
  ////// there needs to be another step to further partition the data
  
	
	  center_p_1=(double *)malloc(sizeof(double)*3);
	  memset(center_p_1,0,sizeof(double)*3);

	  center_p_2=(double *)malloc(sizeof(double)*3);
	  memset(center_p_2,0,sizeof(double)*3);

	  center_p_3=(double *)malloc(sizeof(double)*3);
	  memset(center_p_3,0,sizeof(double)*3);

	  size_p=(long *)malloc(sizeof(long)*num_population);
	  memset(size_p,0,sizeof(long)*num_population);
	
	  temp_size_j=(long *)malloc(sizeof(long)*num_population);
	  memset(temp_size_j,0,sizeof(long)*num_population);

	  all_computed_pop_ID=(long *)malloc(sizeof(long)*file_Len);
	  memset(all_computed_pop_ID,0,sizeof(long)*file_Len);

	  for (i=0;i<file_Len;i++)
	  {
		size_p[all_population_ID[i]]++;
		all_computed_pop_ID[i]=all_population_ID[i];
	  }

	  ind_pop=(double***)malloc(sizeof(double**)*num_population);
	  memset(ind_pop,0,sizeof(double**)*num_population);

	  ind_pop_3=(double***)malloc(sizeof(double**)*num_population);
	  memset(ind_pop_3,0,sizeof(double**)*num_population);

	  for (i=0;i<num_population;i++)
	  {
		  ind_pop[i]=(double**)malloc(sizeof(double*)*size_p[i]);
		  memset(ind_pop[i],0,sizeof(double*)*size_p[i]);

		  ind_pop_3[i]=(double**)malloc(sizeof(double*)*size_p[i]);
		  memset(ind_pop_3[i],0,sizeof(double*)*size_p[i]);

		  for (j=0;j<size_p[i];j++)
		  {
			  ind_pop[i][j]=(double*)malloc(sizeof(double)*num_dm);
			  memset(ind_pop[i][j],0,sizeof(double)*num_dm);

			  ind_pop_3[i][j]=(double*)malloc(sizeof(double)*3);
			  memset(ind_pop_3[i][j],0,sizeof(double)*3);
		  }
		  temp_size_j[i]=0;
	  }

	  for (i=0;i<file_Len;i++) //to generate ind_pop[i]
	  {
		  index_id=all_population_ID[i];
		  
		  j=temp_size_j[index_id];

		  for (t=0;t<num_dm;t++)
		  {
			ind_pop[index_id][j][t]=input_data[i][t];
		  }
		  temp_size_j[index_id]++;		
	  }

	  aver=(double*)malloc(sizeof(double)*num_dm);
	  memset(aver,0,sizeof(double)*num_dm);

	  std=(double*)malloc(sizeof(double)*num_dm);
	  memset(std,0,sizeof(double)*num_dm);

	  temp_population_center=(double**)malloc(sizeof(double*)*2);
	  memset(temp_population_center,0,sizeof(double*)*2);
	  for (i=0;i<2;i++)
	  {
		temp_population_center[i]=(double*)malloc(sizeof(double)*3);
		memset(temp_population_center[i],0,sizeof(double)*3);
	  }

	  size_p_2=(long*)malloc(sizeof(long)*2);
	  memset(size_p_2,0,sizeof(long)*2);

	  partit=(long*)malloc(sizeof(long)*num_population);
	  memset(partit,0,sizeof(long)*num_population);

	  num_computed_population=num_population;

	  temp_population_ID=(long**)malloc(sizeof(long*)*num_population);
	  memset(temp_population_ID,0,sizeof(long*)*num_population);

	  for (i=0;i<num_population;i++)
	  {
		temp_population_ID[i]=(long*)malloc(sizeof(long)*size_p[i]);
		memset(temp_population_ID[i],0,sizeof(long)*size_p[i]);
	  }
	
	  //printf("num_population=%d\n",num_population);
	  
	  for (i=0;i<num_population;i++) 
	  {
		  partit[i]=0;
		  tran(ind_pop[i], size_p[i], num_dm, 1, ind_pop[i]);//0-1 normalize ind_pop[i]
          
		  //find the 3 dimensions with the largest std for ind_pop[i]
		  d1=-1;
		  d2=-1;
		  d3=-1;

		  for (t=0;t<num_dm;t++)
		  {
			aver[t]=0;
			for (j=0;j<size_p[i];j++)
			{
				aver[t]=aver[t]+ind_pop[i][j][t];
			}
			aver[t]=aver[t]/(double)size_p[i];

			std[t]=0;
			for (j=0;j<size_p[i];j++)
			{
				std[t]=std[t]+((ind_pop[i][j][t]-aver[t])*(ind_pop[i][j][t]-aver[t]));
			}
			std[t]=sqrt(std[t]/(double)size_p[i]);
		  }

		  

		  for (j=0;j<3;j++)
		  {
			max_d_dist=0;
			for (t=0;t<num_dm;t++)
			{
				if ((t!=d1) && (t!=d2))
				{
					dist=std[t];
				}
				else
					dist=-1;

				if (dist>max_d_dist)
				{
					max_d_dist=dist;
					d_d=t;						
				}
			}

			if (j==0)
				d1=d_d;
			if (j==1)
				d2=d_d;
			if (j==2)
				d3=d_d;
		  }
		 
		 
		  for (t=0;t<size_p[i];t++)
		  {
				ind_pop_3[i][t][0]=ind_pop[i][t][d1];
				ind_pop_3[i][t][1]=ind_pop[i][t][d2];
				ind_pop_3[i][t][2]=ind_pop[i][t][d3];
		  }
		

		 temp_population_center[0][0]=(aver[d1])/2.0;
		 
		 temp_population_center[0][1]=aver[d2];
	     temp_population_center[0][2]=aver[d3];
		 
		 temp_population_center[1][0]=(aver[d1]+1.0)/2.0;

		 temp_population_center[1][1]=aver[d2];
		 temp_population_center[1][2]=aver[d3];
		 
		 

		 //run K-means with K=2
		 kmeans(ind_pop_3[i], 2, KMEANS_TERM, size_p[i], 3, temp_population_ID[i], temp_population_center);

		 for (j=0;j<2;j++)
			 size_p_2[j]=0;

		 for (j=0;j<size_p[i];j++)
			 size_p_2[temp_population_ID[i][j]]++;

		 for (j=0;j<2;j++)
			 if (size_p_2[j]==0)
				 printf("size_p_2=0 at i=%ld\n",i);


		 //check whether the 2 parts should be merged
		

		 Ei=get_avg_dist(temp_population_center[0], 0,1, temp_population_ID[i], 2,size_p[i], 3, ind_pop_3[i], 0,1,2,size_p_2);
		 Ej=get_avg_dist(temp_population_center[1], 0,1, temp_population_ID[i], 2,size_p[i], 3, ind_pop_3[i], 0,1,2,size_p_2);

		 for (t=0;t<3;t++)
		 {
			center_p_1[t]=temp_population_center[0][t]*0.25+temp_population_center[1][t]*0.75;
			center_p_2[t]=temp_population_center[0][t]*0.5+temp_population_center[1][t]*0.5;
			center_p_3[t]=temp_population_center[0][t]*0.75+temp_population_center[1][t]*0.25;
		 }
		
		 E1=get_avg_dist(center_p_1, 0,1, temp_population_ID[i], 2,size_p[i], 3, ind_pop_3[i], 0,1,2,size_p_2);
		
		 E2=get_avg_dist(center_p_2, 0,1, temp_population_ID[i], 2,size_p[i], 3, ind_pop_3[i], 0,1,2,size_p_2);

		 E3=get_avg_dist(center_p_3, 0,1, temp_population_ID[i], 2,size_p[i], 3, ind_pop_3[i], 0,1,2,size_p_2);

		 //printf("i=%d;E1=%f;E2=%f;E3=%f\n",i,E1,E2,E3);

		 if (E1<E2)
			Ep=E2;
		 else
			Ep=E1;
			
		 Ep=max(Ep,E3); //Ep is the most sparse area

		 

		 if ((Ep>Ei) && (Ep>Ej)) //the two parts should be partitioned
		 {
			 partit[i]=num_computed_population;
			 num_computed_population++;			 
		 }
		
	  } //end for (i=0;i<num_population;i++)

	  
	  for (i=0;i<num_population;i++)
		  temp_size_j[i]=0;

	  for (i=0;i<file_Len;i++)
	  {
		index_id=all_population_ID[i];
		if (partit[index_id]>0)
		{
			j=temp_size_j[index_id];
			if (temp_population_ID[index_id][j]==1)
				all_computed_pop_ID[i]=partit[index_id];
					
			temp_size_j[index_id]++;
		}
	  }

	  printf("Number of groups after partitioning is %ld\n", num_computed_population);

	
      free(size_p_2);
	  
	  free(aver);
	  free(std);
	  free(partit);

	  free(center_p_1);
	  free(center_p_2);
	  free(center_p_3);

	 

	  for (i=0;i<num_population;i++)
		  free(temp_population_ID[i]);
	  free(temp_population_ID);

	  for (i=0;i<num_population;i++)
	  {
		for (j=0;j<size_p[i];j++)
		{
			  free(ind_pop[i][j]);
			  free(ind_pop_3[i][j]);
		}
	  }

	 

	  for (i=0;i<num_population;i++)
	  {
		  free(ind_pop[i]);
		  free(ind_pop_3[i]);
	  }
	  free(ind_pop);
	  free(ind_pop_3);

	  
	 

	  for (i=0;i<2;i++)
		  free(temp_population_center[i]);

	  free(temp_population_center);

	  free(temp_size_j);
	  free(size_p);

	 

	  //update the IDs, Centers, and # of populations
	  for (i=0;i<file_Len;i++)
	  {
		  all_population_ID[i]=all_computed_pop_ID[i];
		  if (all_population_ID[i]>=num_computed_population)
			  printf("all_population_ID[%ld]=%ld\n",i,all_population_ID[i]);
	  }

	  num_population=num_computed_population;

	  free(all_computed_pop_ID);
    //end partitioning
  // since the num_population has changed, population_center needs to be redefined as below

  population_center=(double **)malloc(sizeof(double*)*num_population);
  memset(population_center,0,sizeof(double*)*num_population);
  for (i=0;i<num_population;i++)
  {
     population_center[i]=(double*)malloc(sizeof(double)*num_dm);
     memset(population_center[i],0,sizeof(double)*num_dm);
  }

  	  
  ID2Center_all(normalized_data,file_Len,num_dm,num_population,all_population_ID,population_center);
  


  ////// end of further partitioning

  ////////////////////////////////////////////////////////////
  //  to further merge populations to avoid overpartitioning
  //  Added June 20, 2010
  ////////////////////////////////////////////////////////////
  

  
  num_real_pop=num_population;
  keep_merge=1;

  if (num_pop>num_real_pop)
  {
		fprintf(stderr,"number of populations specified too large\n");
		exit(0);
  }

  center_1=(double *)malloc(sizeof(double)*num_dm);
  memset(center_1,0,sizeof(double)*num_dm);

  center_2=(double *)malloc(sizeof(double)*num_dm);
  memset(center_2,0,sizeof(double)*num_dm);

  center_3=(double *)malloc(sizeof(double)*num_dm);
  memset(center_3,0,sizeof(double)*num_dm);


    
  while (((num_real_pop>num_pop) && (num_pop!=0)) || ((keep_merge==1) && (num_pop==0) && (num_real_pop>2)))
  {
	  
	keep_merge=0;  //so, if no entering merge function, while will exit when num_pop==0

	real_pop_center=(double **)malloc(sizeof(double*)*num_real_pop);
    memset(real_pop_center,0,sizeof(double*)*num_real_pop);
    for (i=0;i<num_real_pop;i++)
    {
      real_pop_center[i]=(double*)malloc(sizeof(double)*num_dm);
      memset(real_pop_center[i],0,sizeof(double)*num_dm);
    }

	min_pop=(double **)malloc(sizeof(double*)*num_real_pop);
    memset(min_pop,0,sizeof(double*)*num_real_pop);
    for (i=0;i<num_real_pop;i++)
    {
      min_pop[i]=(double*)malloc(sizeof(double)*num_dm);
      memset(min_pop[i],0,sizeof(double)*num_dm);
    }

	max_pop=(double **)malloc(sizeof(double*)*num_real_pop);
    memset(max_pop,0,sizeof(double*)*num_real_pop);
    for (i=0;i<num_real_pop;i++)
    {
      max_pop[i]=(double*)malloc(sizeof(double)*num_dm);
      memset(max_pop[i],0,sizeof(double)*num_dm);
    }

	if (num_real_pop==num_population)
	{
		for (i=0;i<num_real_pop;i++)
			for (j=0;j<num_dm;j++)
				real_pop_center[i][j]=population_center[i][j];	
	}
	else
	{
		ID2Center_all(normalized_data,file_Len,num_dm,num_real_pop,all_population_ID,real_pop_center);
	}

	for (i=0;i<num_real_pop;i++)
	{
		for (j=0;j<num_dm;j++)
		{
			min_pop[i][j]=MAX_VALUE;
			
			max_pop[i][j]=0;
		}			
	}
	
	for (i=0;i<file_Len;i++)
	{
		index_id=all_population_ID[i];
		for (j=0;j<num_dm;j++)
		{
			if (normalized_data[i][j]<min_pop[index_id][j])
				min_pop[index_id][j]=normalized_data[i][j];
			
			if (normalized_data[i][j]>max_pop[index_id][j])
				max_pop[index_id][j]=normalized_data[i][j];
		}			
	}

/*	for (i=0;i<num_real_pop;i++)
	{
		for (j=0;j<num_dm;j++)
		{
			printf("min_pop is %f\t",min_pop[i][j]);
			//printf("max_pop is %f\n",max_pop[i][j]);
		}
		printf("\n");
	}
*/


	distance_pop=(double **)malloc(sizeof(double*)*num_real_pop);
    memset(distance_pop,0,sizeof(double*)*num_real_pop);
    for (i=0;i<num_real_pop;i++)
    {
      distance_pop[i]=(double*)malloc(sizeof(double)*num_real_pop);
      memset(distance_pop[i],0,sizeof(double)*num_real_pop);
    }

	for (i=0;i<num_real_pop-1;i++)
	{
		for (j=i+1;j<num_real_pop;j++)
		{
		  distance=0;
		
		  for (t=0;t<num_dm;t++)
		  {
			dist=real_pop_center[i][t]-real_pop_center[j][t];
			distance=distance+dist*dist;
		  }
		  distance_pop[i][j]=distance;
		  distance_pop[j][i]=distance;	      
		}
	}
 
	if (num_real_pop>2)
		num_checked_range=(long)((double)num_real_pop*(num_real_pop-1)/2.0); //to find the mergeable pair among the num_checked_range pairs of smallest distances
	else
		num_checked_range=1;

	size_c=(long *)malloc(sizeof(long)*num_real_pop);
    memset(size_c,0,sizeof(long)*num_real_pop);
	
	for (i=0;i<file_Len;i++)
	  size_c[all_population_ID[i]]++;

	for (p=0;p<num_checked_range;p++)  
	{
		current_smallest_dist=MAX_VALUE;

		for (i=0;i<num_real_pop-1;i++)
		{
			for (j=i+1;j<num_real_pop;j++)
			{
				if (distance_pop[i][j]<current_smallest_dist)
				{
					current_smallest_dist=distance_pop[i][j];
				
					temp_i=i;
					temp_j=j;
				}
			}
		}

		if (p==0)
		{
			first_i=temp_i;
			first_j=temp_j;
		}

		distance_pop[temp_i][temp_j]=MAX_VALUE+1; //make sure this pair won't be selected next time
		
		//normalize and calculate std for temp_i and temp_j
		
		//pick the largest 3 dimensions on standard deviation
		d1=-1;
		d2=-1;
		d3=-1;

		for (i=0;i<3;i++)
		{
			max_d_dist=0;
			for (j=0;j<num_dm;j++)
			{
				if ((j!=d1) && (j!=d2))
				{
					temp=min(min_pop[temp_i][j],min_pop[temp_j][j]);
					tmp=max(max_pop[temp_i][j],max_pop[temp_j][j]);
					if (tmp<=temp)
					{
						//printf("min_pop[%d][%d]=%f \t min_pop[%d][%d]=%f \t max_pop[%d][%d]=%f \t max_pop[%d][%d]=%f\n",temp_i,j,min_pop[temp_i][j],temp_j, j, min_pop[temp_j][j],temp_i,j,max_pop[temp_i][j],temp_j,j,max_pop[temp_j][j]);
						dist=0;
					}
					else
						dist=(real_pop_center[temp_i][j]-real_pop_center[temp_j][j])/(tmp-temp);
					if (dist<0)
						dist=-dist;
				}
				else
					dist=-1;

				if (dist>max_d_dist)
				{
					max_d_dist=dist;
					d_d=j;						
				}
			}

			if (i==0)
				d1=d_d;
			if (i==1)
				d2=d_d;
			if (i==2)
				d3=d_d;
		}

		//printf("d1=%d;d2=%d;d3=%d\n",d1,d2,d3);

		Ei=get_avg_dist(real_pop_center[temp_i], temp_i,temp_j, all_population_ID, num_real_pop,file_Len, num_dm, normalized_data, d1,d2,d3,size_c);
		Ej=get_avg_dist(real_pop_center[temp_j], temp_i,temp_j, all_population_ID, num_real_pop,file_Len, num_dm, normalized_data, d1,d2,d3,size_c);

		for (t=0;t<num_dm;t++)
		{
			center_1[t]=real_pop_center[temp_i][t]*0.25+real_pop_center[temp_j][t]*0.75;
			center_2[t]=real_pop_center[temp_i][t]*0.5+real_pop_center[temp_j][t]*0.5;
			center_3[t]=real_pop_center[temp_i][t]*0.75+real_pop_center[temp_j][t]*0.25;
		}
		
		E1=get_avg_dist(center_1, temp_i,temp_j, all_population_ID, num_real_pop,file_Len, num_dm, normalized_data, d1,d2,d3,size_c);
		
		E2=get_avg_dist(center_2, temp_i,temp_j, all_population_ID, num_real_pop,file_Len, num_dm, normalized_data, d1,d2,d3,size_c);

		E3=get_avg_dist(center_3, temp_i,temp_j, all_population_ID, num_real_pop,file_Len, num_dm, normalized_data, d1,d2,d3,size_c);

		if (E1<E2)
			Ep=E2;
		else
			Ep=E1;
			
		Ep=max(Ep,E3); //Ep is the most sparse area

		Ep=Ep*E_T;
		//printf("Ep=%f;Ei=%f;Ej=%f\n",Ep,Ei,Ej);

		if ((Ep<=Ei) || (Ep<=Ej))//if the most sparse area between i and j are still denser than one of them, temp_i and temp_j should be merged
		{
			keep_merge=1;
			break;		
		}//end if (Ep)
	} //end for p

	//printf("keep_merge=%d\n",keep_merge);

	for (i=0;i<num_real_pop;i++)
	{
		free(real_pop_center[i]);
		free(distance_pop[i]);
		free(min_pop[i]);
		free(max_pop[i]);
	}
	free(real_pop_center);
	free(min_pop);
	free(max_pop);
	free(distance_pop);
	free(size_c);

	//printf("temp_i=%d;temp_j=%d\n",temp_i,temp_j);

	if (keep_merge)  //found one within p loop
	{
		for (i=0;i<file_Len;i++)
		{
			if (all_population_ID[i]>temp_j)
			{
				all_population_ID[i]=all_population_ID[i]-1;
			}
			else 
			{
				if (all_population_ID[i]==temp_j)
				{
					all_population_ID[i]=temp_i;
				}
			}
		}

		num_real_pop--;
	}
	else
	{
		if ((num_pop!=0) && (num_pop<num_real_pop)) //not reach the specified num of pop
		{
			for (i=0;i<file_Len;i++)
			{
				if (all_population_ID[i]>first_j)
				{
					all_population_ID[i]=all_population_ID[i]-1;
				}
				else
				{
					if (all_population_ID[i]==first_j)
					{
						all_population_ID[i]=first_i;
					}
				}
			}

			num_real_pop--;

		

		}
	}
		//printf("num_real_pop is %d\n",num_real_pop);
	
  } //end of while

  

  free(center_1);
  free(center_2);
  free(center_3);
  printf("Final number of populations is %ld\n",num_real_pop);

  new_population_center=(double **)malloc(sizeof(double*)*num_real_pop);
  memset(new_population_center,0,sizeof(double*)*num_real_pop);
  for (i=0;i<num_real_pop;i++)
  {
      new_population_center[i]=(double*)malloc(sizeof(double)*num_dm);
      memset(new_population_center[i],0,sizeof(double)*num_dm);
  }
  
 
  ///////////////////////////////////////////
  //End of population mapping
  ///////////////////////////////////////////

  show(input_data, all_population_ID, file_Len, num_real_pop, num_dm, para_name_string);

  ID2Center_all(input_data,file_Len,num_dm,num_real_pop,all_population_ID,new_population_center);
  

  f_cid=fopen("population_id.txt","w");
  f_ctr=fopen("population_center.txt","w");
  f_out=fopen("coordinates.txt","w");
  f_results=fopen("flock_results.txt","w");

/*
  f_parameters=fopen("parameters.txt","w");
  fprintf(f_parameters,"Number_of_Bins\t%d\n",num_bin);
  fprintf(f_parameters,"Density\t%f\n",aver_index);
  fclose(f_parameters);
*/

  for (i=0;i<file_Len;i++)
	fprintf(f_cid,"%ld\n",all_population_ID[i]+1); //all_population_ID[i] changed to all_population_ID[i]+1 to start from 1 instead of 0: April 16, 2009

  /*
   * New to check for min/max to add to parameters.txt
   *
  */
  
  fprintf(f_out,"%s\n",para_name_string);
  //fprintf(f_results,"%s\tEvent\tPopulation\n",para_name_string);
  fprintf(f_results,"%s\tPopulation\n",para_name_string);
  for (i=0;i<file_Len;i++)
  {
	for (j=0;j<num_dm;j++)
	{
		if (input_data[i][j] < min) {
			min = (int)input_data[i][j];
		}
		if (input_data[i][j] > max) {
			max = (int)input_data[i][j];
		}
		if (j==num_dm-1)
		{
			fprintf(f_out,"%d\n",(int)input_data[i][j]);
			fprintf(f_results,"%d\t",(int)input_data[i][j]);
		}
		else
		{
			fprintf(f_out,"%d\t",(int)input_data[i][j]);
			fprintf(f_results,"%d\t",(int)input_data[i][j]);
		}
	}
	//fprintf(f_results,"%ld\t",i + 1);
	fprintf(f_results,"%ld\n",all_population_ID[i]+1); //all_population_ID[i] changed to all_population_ID[i]+1 to start from 1 instead of 0: April 16, 2009
  }

/*
  f_parameters=fopen("parameters.txt","w");
  fprintf(f_parameters,"Number_of_Bins\t%ld\n",num_bin);
  fprintf(f_parameters,"Density\t%d\n",den_t_event);
  fprintf(f_parameters,"Min\t%d\n",min);
  fprintf(f_parameters,"Max\t%d\n",max);
  fclose(f_parameters);
*/

  f_properties=fopen("fcs.properties","w");
  fprintf(f_properties,"Bins=%ld\n",num_bin);
  fprintf(f_properties,"Density=%d\n",den_t_event);
  fprintf(f_properties,"Min=%d\n",min);
  fprintf(f_properties,"Max=%d\n",max);
  fprintf(f_properties,"Populations=%ld\n",num_real_pop);
  fprintf(f_properties,"Events=%ld\n",file_Len);
  fprintf(f_properties,"Markers=%ld\n",num_dm);
  fclose(f_properties);

  for (i=0;i<num_real_pop;i++) {
	/* Add if we want to include population id in the output
	*/
	fprintf(f_ctr,"%ld\t",i+1);  //i changed to i+1 to start from 1 instead of 0: April 16, 2009

	for (j=0;j<num_dm;j++) {
		if (j==num_dm-1)
			fprintf(f_ctr,"%.0f\n",new_population_center[i][j]);
		else
			fprintf(f_ctr,"%.0f\t",new_population_center[i][j]);
	}
  }

  	//added April 16, 2009
	f_mfi=fopen("MFI.txt","w");

	for (i=0;i<num_real_pop;i++)
	{
		fprintf(f_mfi,"%ld\t",i+1);

		for (j=0;j<num_dm;j++)
		{
			if (j==num_dm-1)
				fprintf(f_mfi,"%.0f\n",new_population_center[i][j]);
			else
				fprintf(f_mfi,"%.0f\t",new_population_center[i][j]);
		}
	}
	fclose(f_mfi);

	//ended April 16, 2009
			
  fclose(f_cid);
  fclose(f_ctr);
  fclose(f_out);
  fclose(f_results);


  for (i=0;i<num_population;i++)
  	free(population_center[i]);
  
  free(population_center);
 
  for (i=0;i<num_real_pop;i++)
	  free(new_population_center[i]);
  
  free(new_population_center);

  for (i=0;i<file_Len;i++)
    free(normalized_data[i]);
  free(normalized_data);	
	
  free(grid_populationID);

  free(cluster_populationID);
  free(grid_clusterID);
  free(cluster_ID);

  for (i=0;i<file_Len;i++)
    free(input_data[i]);
  free(input_data);

  free(grid_ID);
  free(population_ID);
  free(all_population_ID);
  free(eventID_To_denseventID);
		
  ///////////////////////////////////////////////////////////
  printf("Ending time:\t\t\t\t");
  fflush(stdout);
  system("/bin/date");

  /*
   * Windows version
  _strtime( tmpbuf );
  printf( "Ending time:\t\t\t\t%s\n", tmpbuf );
  _strdate( tmpbuf );
  printf( "Ending date:\t\t\t\t%s\n", tmpbuf );
 */
  
  return 0;
}