//To compile: gcc -std=c99 -lm isrs.c isrs.exe 
//to use: ./isrs.exe <Ntol> <input data filename> <output cluster filename> 


#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <ctype.h>

//the following data structure is for the ordered sample
struct data_stru_def{
  long int Nevents;               //size of the sample
  double *X_poi;                      //array containing the ordered events data-sample
                                  //(e.g. the cumulative exposure at the time of the event for the gamma-ray data sample)
                                  //or the time of the event (if the exposure is flat)
  double *t_poi;                      //auxiliary info for plotting: array containing the time of the events (it differs from X when the exposure is not flat) 
  double X_start;                    // X start of observation
  double X_stop;                     // X stop of observation
  double t_start;                    // t start of observation
  double t_stop;                     // t stop of observation
};


typedef struct cluster_candidate {
  long int *level_id;
  long int *i_start;
  long int *i_stop;
  long int Nclusters;
} cluster_candidate_t;

typedef struct cluster {
  long int *level_id;
  double *chance_prob;
  double *t_start;
  double *t_stop;
  long int *Nele;
  double *len;
  double *elen;
  long int Nclusters;
} cluster_t;




//prototypes
void iSRS(struct data_stru_def *data_stru,long int N_tol,char *clusterfile,cluster_candidate_t *cluster_candidate_set,char *toutformat,char *xoutformat,bool use_cluster_buf);
long int decode_spacing_patch(int *spacing_patch,struct data_stru_def *data_stru,FILE *cfile,long int level_id,cluster_candidate_t *cluster_candidate_set,char *toutformat,char *xoutformat,bool use_cluster_buf);
void SRS(struct data_stru_def *data_stru,double Delta_thr,long int N_tol,int *spacing_patch);
void SRS_element(double Delta_thr,long int N_tol,long int ievt,struct data_stru_def *data_stru ,int *spacing_patch);
double get_distance(struct data_stru_def *data_stru,long int ifirst,long int ilast);
double get_max_distance_among_contiguous_events(struct data_stru_def *data_stru);
void fill_data_stru(struct data_stru_def *data_stru,char *datafilename,bool twocolumns_flag,bool differentialX_on_input_flag,bool binary_flag);
void fill_data_stru_onecolumn(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag);
void fill_data_stru_onecolumn_binary(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag);
void fill_data_stru_twocolumns(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag);
void fill_data_stru_twocolumns_binary(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag);
void save_cluster_candidate_info(struct data_stru_def *data_stru,long int i_cluster_candidate_start,long int i_cluster_candidate_stop,FILE *cfile,long int level_id,cluster_candidate_t *cluster_candidate_set,char *toutformat,char *xoutformat,bool use_cluster_buf);
FILE *init_cluster_candidate_file(char *clusterfile,struct data_stru_def *data_stru,char *toutformat,char *xoutformat,long int N_tol);
int fprintf_cluster(FILE *cfile,long int level_id,double t_start,double t_stop,long int Nele,double len,double elen,char *toutformat,char *xoutformat);
void close_cluster_candidate_file(FILE *cfile,cluster_candidate_t *cluster_candidate_set);
void init_cluster_candidate_set(struct data_stru_def *data_stru,cluster_candidate_t *cluster_candidate_set,bool use_cluster_buf);
bool cluster_candidate_is_new(long int i_cluster_candidate_start,long int i_cluster_candidate_stop,cluster_candidate_t *cluster_candidate_set,bool use_cluster_buf);
void free_data_stru(struct data_stru_def *data_stru,bool twocolumns_flag);
char *stolower(char *mystring);
void help(void);

//************************************************
int main(int argc, char* argv[]){   
  //prepares clusters from a datasample, with  "N_{tol}" parameter 
  //two input: "N_{tol}" (see paper), and the name of the file with the data sample
  
  bool use_cluster_buf=true;
  bool twocolumns_flag=false;
  bool differentialX_on_input_flag=false;
  bool binary_flag=false;
  char toutformat[10];
  strcpy(toutformat,"%18.10lf");
  char xoutformat[10];
  strcpy(xoutformat,"%17.10lE");
  char myarg[1000];

  int iarg;
  for (iarg=0;iarg<argc;iarg++){
    strcpy(myarg,argv[iarg]);
    if(strcmp(stolower(myarg),"help") == 0){
      help();
      return 0;
    }
  }

  if(argc < 4){
    fprintf(stderr,"not enough parameters specified, bye\n");
    exit(-1);
  }
  long int N_tol=atol(argv[1]);
  char datafile[300];
  strcpy(datafile,argv[2]);
  char clusterfile[300];
  strcpy(clusterfile,argv[3]);

  iarg=4;
  bool xoutformat_changed=false;
  while (iarg<argc){
    strcpy(myarg,argv[iarg]);
    if(strcmp(stolower(myarg),"twocolumns") == 0){
      printf("reading two columns data\n");
      twocolumns_flag=true;
    }else if(strcmp(stolower(myarg),"no_rm_copies") == 0){
      printf("maintining cluster duplicates to save memory\n");
      use_cluster_buf=false;
    }else if(strcmp(stolower(myarg),"binary") == 0){
      printf("reading input file in binary mode\n");
      binary_flag=true;
    }else if(strcmp(stolower(myarg),"differentialx") == 0){
      printf("assuming differential X on input\n");
      differentialX_on_input_flag=true;
    }else if(strcmp(stolower(myarg),"t_format") == 0){
      if(iarg+1 < argc){
	strcpy(toutformat,argv[iarg+1]);
	printf("t format: %s\n",toutformat);
	if(!xoutformat_changed){
	  strcpy(xoutformat,argv[iarg+1]);
	}
	iarg++;
      }      
    }else if(strcmp(stolower(myarg),"x_format") == 0){
      if(iarg+1 < argc){
	strcpy(xoutformat,argv[iarg+1]);
	printf("X format: %s\n",xoutformat);
	iarg++;
	xoutformat_changed=true;
      }      
    }else{
      fprintf(stderr,"unknown parameter: %s, bye\n",argv[iarg]);
      exit(-1);
    }
    
    iarg++;
  }


  
  struct data_stru_def data_stru;
  fill_data_stru(&data_stru,datafile,twocolumns_flag, differentialX_on_input_flag,binary_flag);               //read sample, problem dependent
  
  cluster_candidate_t cluster_candidate_set;
  iSRS(&data_stru,N_tol,clusterfile,&cluster_candidate_set,toutformat,xoutformat, use_cluster_buf);                                  //performs iSRS clustering

  //printf("freeing datastru\n");
  free_data_stru(&data_stru,twocolumns_flag);
  
  return 0;
}

//************************************************
void iSRS(struct data_stru_def *data_stru,long int N_tol,char *clusterfile,cluster_candidate_t *cluster_candidate_set,char *toutformat,char *xoutformat,bool use_cluster_buf){
  // the function "iSRS" performs the iterated SRS clustering procedure
  
  long int nlevels=400;
  double Delta_thr_max=get_max_distance_among_contiguous_events(data_stru);    // evaluates the max distance among contiguous events
  long int nclusters_withsize_larger_than_2;
  
  //data_stru contains the sample under investigation
  //
  //The array "spacing_patch" is used for the spacings patchwork.
  //The first element is for the spacing among the first and the second event of
  //the sample, the element "i" is for the spacing among the "ith" event and the "(i+1)th" event
  // To slightly ease the decoding,
  //it contains as much elements "Nevents"  as the size of the studied sample (instead of Nevents-1).
  //It will be filled with the value "1" at the index "i" if the element "i" of the data-set is an element of a cluster
  //But the last element of a cluster is the first "0" element after a sequence of "1".
  //NB: At the first element of the cluster corresponds the first "1" after a "0" 
  //example:
  //Two clusters are identified
  //      the 1st cluster is of size 6, and starts at element 4, ends at element 9
  //      the 2nd cluster is of size 4, and starts at element 13, ends at element 16
  //the array "spacing_patch" will be as follow:
  //spacing_patch= {0,0,0,0,1,1,1,1,1,0,0,0,0,1,1,1,0,0,0,0}
  //
  
  long int Nevents=data_stru->Nevents;               // events within the data set    
  int *spacing_patch=(int *)malloc(Nevents*sizeof(int));  //nb: the elements of spacing_patch are zeroed within patch()
  if(spacing_patch == NULL){
    fprintf(stderr,"unable to get memory for spacing vector\n");
    exit(-1);
  }

  //create the structure for the cluster set
  init_cluster_candidate_set(data_stru,cluster_candidate_set, use_cluster_buf);
  
  //create the file where to put the cluster set
  FILE *cfile=init_cluster_candidate_file(clusterfile,data_stru,toutformat,xoutformat,N_tol);

  //iteration on Delta_thr
  long int level;
  for(level=0;level<nlevels;level++){
    
    //compute Delta_thr value
    double index=(double)level/20;   //start from 0, 20 steps per decade
    double proexpo_level=pow(10.,index);
    double Delta_thr=Delta_thr_max/proexpo_level;

    // the function "SRS()" fills the "spacing_patch" array 
    SRS(data_stru,Delta_thr,N_tol,spacing_patch);
    
    //The return value of "decode_spacing_patch()" is the number of clusters found with size larger than 2    
    nclusters_withsize_larger_than_2=decode_spacing_patch(spacing_patch,data_stru,cfile,level,cluster_candidate_set,toutformat,xoutformat,use_cluster_buf);
    //if no clusters are found with size > 2, the iSRS procedure returns
    if(nclusters_withsize_larger_than_2 <= 0)goto after_for_patch_levels;
    
    
  }
 after_for_patch_levels:

  free(spacing_patch);  
  close_cluster_candidate_file(cfile,cluster_candidate_set);
  
  return;
}

//*****************************************************
long int decode_spacing_patch(int *spacing_patch,struct data_stru_def *data_stru,FILE *cfile,long int level_id,cluster_candidate_t *cluster_candidate_set,char *toutformat,char *xoutformat,bool use_cluster_buf){        
  //the function "decode_spacing_patch()" decodes the array "spacing_patch"
  //("spacing_patch" contains the clusters)
  // and returns the number of clusters with size > 2;

   long int nclusters_withsize_larger_than_2=0;  //this is the number of clusters found in "spacing_patch" with size > 2.
   int spacing_patch_value_old=0;
   int spacing_patch_value;
   long int i_cluster_candidate_start=-1;
   long int i_cluster_candidate_stop=-1;
   long int cluster_candidate_size=-1;

   long int ispacing_patch;    
   for(ispacing_patch=0;ispacing_patch<data_stru->Nevents;ispacing_patch++){
     spacing_patch_value=*(spacing_patch+ispacing_patch);
     if((spacing_patch_value_old == 0) && (spacing_patch_value == 1)){
        i_cluster_candidate_start=ispacing_patch;                                   //first element of the cluster
     }else if((spacing_patch_value_old == 1) && (spacing_patch_value == 0)){
       i_cluster_candidate_stop=ispacing_patch;                                    //last element of the cluster
       cluster_candidate_size=i_cluster_candidate_stop-i_cluster_candidate_start+1;               //cluster size
       if(cluster_candidate_size > 2){
         nclusters_withsize_larger_than_2++;                        //#clusters with size > 2 
       }
       if(cluster_candidate_is_new(i_cluster_candidate_start,i_cluster_candidate_stop,cluster_candidate_set,use_cluster_buf)){
	 save_cluster_candidate_info(data_stru,i_cluster_candidate_start,i_cluster_candidate_stop,cfile,level_id,cluster_candidate_set,toutformat,xoutformat,use_cluster_buf); // Here the user can put its own code to save cluster info.
       }
     }

     spacing_patch_value_old=spacing_patch_value;
   } //end loop on ispacing_patch

  return nclusters_withsize_larger_than_2;
}

//*****************************************************
void SRS(struct data_stru_def *data_stru,double Delta_thr,long int N_tol,int *spacing_patch){
  //the SRS function detects clusters within data_stru
  //according to eq. (2), with params Delta_thr and N_tol

  //clear "the spacing_patch" array.
  long int Nevents=data_stru->Nevents;
  long int ievt;
  for(ievt=0;ievt<Nevents;ievt++){  
    *(spacing_patch+ievt)=0;
  }

  for(ievt=0;ievt<Nevents;ievt++){
    SRS_element(Delta_thr,N_tol,ievt,data_stru,spacing_patch);    //verify cluster condistion (eq.2) for the event ievt

  }

  return;
}

//*******************************************************************************************
void SRS_element(double Delta_thr,long int N_tol,long int ievt,struct data_stru_def *data_stru ,int *spacing_patch){
  //the "SRS_element" function check eq. (2)
  //with params Delta_thr and N_tol but only for the event "ievt"
  
  double Delta;

  long int N_tol_step=1;
  long int N_tol_min=1;
  long int i_tol=N_tol;
  bool patch_found_flag=false;
  long int imin;
  long int imax;
  long int kk;
  long int iele;
  while ((i_tol >= N_tol_min) && (!patch_found_flag) ){
    if(ievt > i_tol-1){
      
      imin=ievt-i_tol;
      imax=ievt;
      kk=imax-imin;
      Delta=get_distance(data_stru,imin,imax);     //get distance among events imin and imax
      if(Delta <= Delta_thr*kk){                 // if eq. 2 is satisfied, put "1" in spacing_patch for all the events in between

	for(iele=imin;iele<imax;iele++){
	  *(spacing_patch+iele)=1;      
	}
	patch_found_flag=true;

      }
    }
    i_tol=i_tol-N_tol_step;  //for the next step 
  }

  return;
}

//**********************************
double get_distance(struct data_stru_def *data_stru,long int ifirst,long int ilast){
  
  //the function "get_distance" evaluates  the distance among the events "ifirst" and "ilast"  
  // get_distance depends on the data-set

  double X_last=*(data_stru->X_poi+ilast);
  double X_first=*(data_stru->X_poi+ifirst);
  double distance=X_last-X_first;

  return distance;
}

//********************************
double get_max_distance_among_contiguous_events(struct data_stru_def *data_stru){

  //the function "get_max_distance_among_contiguous_events" evaluate the max
  //distance among two contiguous events within the sample (it is called "\Delta^{max}" in the paper)

  double max_distance=0.;
  double distance;
  long int Nevents=data_stru->Nevents;
  long int ievt;
  for (ievt=0;ievt<Nevents-1;ievt++){
    distance=get_distance(data_stru,ievt,ievt+1);
    if(distance > max_distance){
      max_distance=distance;
    }
  }

  return max_distance;
}

//*****************************************************
void free_data_stru(struct data_stru_def *data_stru,bool twocolumns_flag){

  free(data_stru->X_poi);
  if(twocolumns_flag){
    free(data_stru->t_poi);
  }
    
  return;
}

//*****************************************************
void fill_data_stru(struct data_stru_def *data_stru,char *datafilename,bool twocolumns_flag,bool differentialX_on_input_flag,bool binary_flag){

  if(binary_flag){
    if(twocolumns_flag){
      printf("reading %s binary file in two columns mode\n",datafilename);
      fill_data_stru_twocolumns_binary(data_stru,datafilename, differentialX_on_input_flag);
      return;
    }
    printf("reading %s binary file in one column mode\n",datafilename);
    fill_data_stru_onecolumn_binary(data_stru,datafilename, differentialX_on_input_flag);
    return;
  }

  if(twocolumns_flag){
    printf("reading %s in two columns mode\n",datafilename);
    fill_data_stru_twocolumns(data_stru,datafilename, differentialX_on_input_flag);
    return;
  }
  printf("reading %s in one column mode\n",datafilename);
  fill_data_stru_onecolumn(data_stru,datafilename, differentialX_on_input_flag);
  return;

}

//*****************************************************
void fill_data_stru_twocolumns(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag){
  //fill the data_stru with the sample data contained in datafile
  // ... depends on the problem under investigation

  long int Nevents;
  long int iele0;

  FILE *datafile=fopen(datafilename,"r");
  if(datafile == NULL){
    fprintf(stderr,"error while opening datafile %s\n",datafilename);
    exit(-1);
  }
  char dummys[1000];
  fscanf(datafile,"%s %s %li\n",dummys,dummys,&Nevents);

  data_stru->Nevents=Nevents;
  //printf("Nevt: %li %li %li\n",Nevents,sizeof(double)*data_stru->Nevents,sizeof(double));
  data_stru->X_poi=(double *)malloc(sizeof(double)*data_stru->Nevents*2+10);
  if(data_stru->X_poi == NULL){
    fprintf(stderr,"unable to get memory for data\n");
    exit(-1);
  }
  data_stru->t_poi=(double *)malloc(sizeof(double)*data_stru->Nevents*2+10); 
  if(data_stru->t_poi == NULL){
    fprintf(stderr,"unable to get memory for data\n");
    exit(-1);
  }

  double time_appo;
  double expo_appo;

  fscanf(datafile,"%s %lf %lf\n",dummys,&expo_appo,&time_appo);
  data_stru->X_start=expo_appo;
  data_stru->t_start=time_appo;
  fscanf(datafile,"%s %lf %lf\n",dummys,&expo_appo,&time_appo);
  data_stru->X_stop=expo_appo;
  data_stru->t_stop=time_appo;

  long int iele=0;
  double cumulative_expo=0.;
  while(fscanf(datafile,"%lf %lf\n",&expo_appo,&time_appo)==2){
    if(differentialX_on_input_flag){
      //printf("exp: %lE %lE \n",expo_appo,cumulative_expo);
      *(data_stru->X_poi+iele)=cumulative_expo;
      cumulative_expo+=expo_appo;
    }else{
      *(data_stru->X_poi+iele)=expo_appo;
    }
    *(data_stru->t_poi+iele)=time_appo;
    iele++;
  }
  
  fclose(datafile);

  return;
}

//*****************************************************
void fill_data_stru_onecolumn(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag){
  //fill the data_stru with the sample data contained in datafile
  // ... depends on the problem under investigation

  long int Nevents;
  long int iele0;

  FILE *datafile=fopen(datafilename,"r");
  if(datafile == NULL){
    fprintf(stderr,"error while opening datafile %s\n",datafilename);
    exit(-1);
  }
  char dummys[1000];
  fscanf(datafile,"%s %s %li\n",dummys,dummys,&Nevents);
  data_stru->Nevents=Nevents;
  printf("Nevt: %li %li %li\n",Nevents,sizeof(double)*data_stru->Nevents,sizeof(double));
  data_stru->X_poi=(double *)malloc(sizeof(double)*data_stru->Nevents*2+10);
  if(data_stru->X_poi == NULL){
    fprintf(stderr,"unable to get memory for data\n");
    exit(-1);
  }
  data_stru->t_poi=data_stru->X_poi;


  double expo_appo;

  //printf("reading 1col\n");
  fscanf(datafile,"%s %lf\n",dummys,&expo_appo);
  data_stru->X_start=expo_appo;
  fscanf(datafile,"%s %lf\n",dummys,&expo_appo);
  data_stru->X_stop=expo_appo;
  data_stru->t_start=data_stru->X_start;
  data_stru->t_stop=data_stru->X_stop;
  printf("start:%lf, stop: %lf\n",data_stru->X_start,data_stru->X_stop);
  long int iele=0;
  double expo_appo_old=0.;
  while(fscanf(datafile,"%lf\n",&expo_appo)==1){

    if(differentialX_on_input_flag){
      *(data_stru->X_poi+iele)=expo_appo+expo_appo_old;
    }else{
      *(data_stru->X_poi+iele)=expo_appo;
    }
    expo_appo_old=expo_appo;
    iele++;
  }
  
  fclose(datafile);

  return;
}

//**********************************************************************
void save_cluster_candidate_info(struct data_stru_def *data_stru,long int i_cluster_candidate_start,long int i_cluster_candidate_stop,FILE *cfile,long int level_id,cluster_candidate_t *cluster_candidate_set,char *toutformat,char *xoutformat,bool use_cluster_buf){

  //write cluster infos

  if(use_cluster_buf){
    long int iclu=cluster_candidate_set->Nclusters;
    *(cluster_candidate_set->level_id+iclu)=level_id;
    *(cluster_candidate_set->i_start+iclu)=i_cluster_candidate_start;
    *(cluster_candidate_set->i_stop+iclu)=i_cluster_candidate_stop;
  }
  cluster_candidate_set->Nclusters++;

  long int Nele=i_cluster_candidate_stop-i_cluster_candidate_start+1;
  double x_start=*((data_stru->X_poi)+i_cluster_candidate_start);
  double x_stop=*((data_stru->X_poi)+i_cluster_candidate_stop);
  double t_start=*((data_stru->t_poi)+i_cluster_candidate_start);
  double t_stop=*((data_stru->t_poi)+i_cluster_candidate_stop);
  double len=x_stop-x_start;
  double elen=(x_stop-x_start)*(1.+1./Nele);

  fprintf_cluster(cfile,level_id,t_start,t_stop,Nele,len,elen,toutformat,xoutformat);    
  return;
}

//****************************************************************************
FILE *init_cluster_candidate_file(char *clusterfile,struct data_stru_def *data_stru,char *toutformat,char *xoutformat,long int N_tol){

  //open cluster file and write down the root as first element  (with level id -1)
  
  FILE *cfile=fopen(clusterfile,"w");
  if(cfile == NULL){
    fprintf(stderr,"opening cluster output file\n");
    exit(-1);
  }


  long int nclusters_dummy=-1;
  fprintf(cfile,"Nclusters %8li\n",nclusters_dummy);
  fprintf(cfile,"N_tol %8li\n",N_tol);
  fprintf(cfile,"t_format %s\n",toutformat);
  fprintf(cfile,"x_format %s\n",xoutformat);
  //file columns:
  //level id of \Delta_{thr}   (-1 for root)
  //  t_start           
  //  t_stop
  //  cluster size
  //  cluster len
  //  cluster effective length    
  fprintf(cfile,"col1: \"level id\" (-1 is for the root)\n");
  fprintf(cfile,"col2: \"t start\"\n");
  fprintf(cfile,"col3: \"t stop\"\n");
  fprintf(cfile,"col4: \"cluster size\" (in X domain if \"twocolumns\" mode has been selected)\n");
  fprintf(cfile,"col5: \"cluster length\" (in X domain if \"twocolumns\" mode has been selected)\n");
  fprintf(cfile,"col6: \"effective cluster length\" (in X domain if \"twocolumns\" mode has been selected)\n");  
  fprintf(cfile,"col7: \"event density\" (in X domain if \"twocolumns\" mode has been selected)\n");  
  //write the root:
  long int level_id=-1;
  double t_start=data_stru->t_start;
  double t_stop=data_stru->t_stop;
  double len=data_stru->X_stop-data_stru->X_start;
  double elen=len;
  double Nele=data_stru->Nevents;
  fprintf_cluster(cfile,level_id,t_start,t_stop,Nele,len,elen,toutformat,xoutformat);

  return cfile;
}

//************************************
int fprintf_cluster(FILE *cfile,long int level_id,double t_start,double t_stop,long int Nele,double len,double elen,char *toutformat,char *xoutformat){
  char format[200];
  double rho=(double)Nele/elen;

  strcpy(format,"%8li ");
  strcat(format,toutformat);
  strcat(format," ");
  strcat(format,toutformat);
  strcat(format," %8li ");
  strcat(format,xoutformat);
  strcat(format," ");
  strcat(format,xoutformat);
  strcat(format," %10.3lE");
  strcat(format,"\n");

  int res=fprintf(cfile,format,level_id,t_start,t_stop,Nele,len,elen,rho);
  return res;
}

//*************************************
void close_cluster_candidate_file(FILE *cfile,cluster_candidate_t *cluster_candidate_set){
  
  //write ther true number of clusters at the head of the cluster file
  long int nclusters=cluster_candidate_set->Nclusters+1;   //+1 because the saved clusters will contains also the root (the cluster_candidate_set does not contain the root)
  fseek(cfile,0,SEEK_SET);
  fprintf(cfile,"Nclusters %8li\n",nclusters);
  fseek(cfile,0,SEEK_END);
  fclose(cfile);  
  return;
}


//****************************************
void init_cluster_candidate_set(struct data_stru_def *data_stru,cluster_candidate_t *cluster_candidate_set,bool use_cluster_buf){
  //init the cluster set (which does not include root)
  
  cluster_candidate_set->Nclusters=0;  
  if(use_cluster_buf){
    printf("clu candidate set size %li\n",data_stru->Nevents-1);
    cluster_candidate_set->level_id=(long int  *)malloc((data_stru->Nevents-1)*sizeof(long int));
    if(cluster_candidate_set->i_start == NULL){
      fprintf(stderr,"unable to get memory in init_cluster_candidate_set\n");
      exit(-1);
    }
    cluster_candidate_set->i_start=(long int  *)malloc((data_stru->Nevents-1)*sizeof(long int));
    if(cluster_candidate_set->i_start == NULL){
      fprintf(stderr,"unable to get memory in init_cluster_candidate_set\n");
      exit(-1);
    }
    cluster_candidate_set->i_stop=(long int  *)malloc((data_stru->Nevents-1)*sizeof(long int));
    if(cluster_candidate_set->i_stop == NULL){
      fprintf(stderr,"unable to get memory in init_cluster_candidate_set\n");
      exit(-1);
    }
  }

   return;
}

//*****************************************
bool cluster_candidate_is_new(long int i_cluster_candidate_start,long int i_cluster_candidate_stop,cluster_candidate_t *cluster_candidate_set,bool use_cluster_buf){

  if(!use_cluster_buf){
    return true;
  }


  long int i_cluster_candidate_start_this;
  long int i_cluster_candidate_stop_this;
  long int iclu;

  for(iclu=0;iclu<cluster_candidate_set->Nclusters;iclu++){

    i_cluster_candidate_start_this=*(cluster_candidate_set->i_start+iclu);
    if(i_cluster_candidate_start == i_cluster_candidate_start_this){

      i_cluster_candidate_stop_this=*(cluster_candidate_set->i_stop+iclu);
      if(i_cluster_candidate_stop == i_cluster_candidate_stop_this){
	return false;
      }

    }

  }
  
  return true;
}

//*************************************
char *stolower(char *str){
  long int i;
  for(i = 0; str[i]; i++){
    str[i] = tolower(str[i]);
  }
  return str;
}

//************************************
void help(){
  printf("isrs usage:\n");
  printf("./isrs <N tol> <input event list> <output cluster list>\n");
  printf("\n");
  printf("parameters:\n");
  printf("  <N tol>: tolerance parameter\n");
  printf("  <input event list>: contains the original set of events, see example for its format\n");
  printf("  <output cluster list>: contains the list of built clusters\n");
  printf("\n");
  printf("the input event list could be the list of times of occurrence of the events\n");
  printf("or it could contain the list of time and the cumulative exposure at the time\n");
  printf("of the occurrence of the events.\n");
  printf("In this second case, the input event file has two columns\n");
  printf("(i.e., \"time\" and \"cumulative exposure\")\n");
  printf("\n");
  printf("if the <input event list> is in two columns format (there is a time and an expo column)\n");
  printf("the optional flag \"twocolumns\" must be set:\n");
  printf("e.g., ./isrs <N tol> <input event list> <output cluster list> twocolumns\n");
  printf("\n");
  printf("optional parameters:\n");
  printf("    t_format <t format value>\n");
  printf("    x_format <x format value>\n");
  printf("    no_rm_copies\n");
  printf("\n");
  printf("\"t_format\" can be set to specify the format of the\n");
  printf("boundaries of the clusters (in time domain).\n");
  printf("Format must be specified in c style (e.g., %%13.8lf or %%15.2lE)\n");
  printf("\n");
  printf("\"x_format\" can be set to specify the format of the \"cluster length\" and\n");
  printf("of \"effective cluster length\"  (in the domain of exposure)\n");
  printf("Format must be specified in c style (e.g., %%13.8lf or %%15.2lE)");
  printf("\n");
  printf("\"no_rm_copies\" can be set to avoid excessive memory consumption, at the price of\n");
  printf("maintaining cluster duplicates in the output file\n");
  printf("\n");
  return;
}

//*****************************************************
void fill_data_stru_twocolumns_binary(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag){
  //fill the data_stru with the sample data contained in datafile
  // ... depends on the problem under investigation

  long int Nevents;
  long int iele0;

  FILE *datafile=fopen(datafilename,"rb");
  if(datafile == NULL){
    fprintf(stderr,"error while opening datafile %s\n",datafilename);
    exit(-1);
  }
  char dummys[1000];
  fread(&Nevents,sizeof(long int),1,datafile);

  data_stru->Nevents=Nevents;
  printf("Nevt: %li %li %li\n",Nevents,sizeof(double)*data_stru->Nevents,sizeof(double));
  data_stru->X_poi=(double *)malloc(sizeof(double)*data_stru->Nevents*2+10);
  if(data_stru->X_poi == NULL){
    fprintf(stderr,"unable to get memory for data\n");
    exit(-1);
  }
  data_stru->t_poi=(double *)malloc(sizeof(double)*data_stru->Nevents*2+10); 
  if(data_stru->t_poi == NULL){
    fprintf(stderr,"unable to get memory for data\n");
    exit(-1);
  }

  double time_appo;
  double expo_appo;

  fread(&expo_appo,sizeof(double),1,datafile);
  fread(&time_appo,sizeof(double),1,datafile);
  data_stru->X_start=expo_appo;
  data_stru->t_start=time_appo;

  fread(&expo_appo,sizeof(double),1,datafile);
  fread(&time_appo,sizeof(double),1,datafile);
  data_stru->X_stop=expo_appo;
  data_stru->t_stop=time_appo;

  long int iele=0;
  double cumulative_expo=0.;
  while(1 == fread(&expo_appo,sizeof(double),1,datafile)){
    fread(&time_appo,sizeof(double),1,datafile);
    if(differentialX_on_input_flag){
      //printf("exp: %lE %lE \n",expo_appo,cumulative_expo);
      *(data_stru->X_poi+iele)=cumulative_expo;
      cumulative_expo+=expo_appo;
    }else{
      *(data_stru->X_poi+iele)=expo_appo;
    }
    *(data_stru->t_poi+iele)=time_appo;
    iele++;
  }
  
  fclose(datafile);

  return;
}


//*****************************************************
void fill_data_stru_onecolumn_binary(struct data_stru_def *data_stru,char *datafilename,bool differentialX_on_input_flag){
  //fill the data_stru with the sample data contained in datafile
  // ... depends on the problem under investigation

  long int Nevents;
  long int iele0;

  FILE *datafile=fopen(datafilename,"rb");
  if(datafile == NULL){
    fprintf(stderr,"error while opening datafile %s\n",datafilename);
    exit(-1);
  }
  char dummys[1000];
  fread(&Nevents,sizeof(long int),1,datafile);

  data_stru->Nevents=Nevents;

  data_stru->X_poi=(double *)malloc(sizeof(double)*data_stru->Nevents*2+10);
  if(data_stru->X_poi == NULL){
    fprintf(stderr,"unable to get memory for data\n");
    exit(-1);
  }
  data_stru->t_poi=data_stru->X_poi;

  double expo_appo;

  fread(&expo_appo,sizeof(double),1,datafile);
  data_stru->X_start=expo_appo;

  fread(&expo_appo,sizeof(double),1,datafile);
  data_stru->X_stop=expo_appo;

  data_stru->t_start=data_stru->X_start;
  data_stru->t_stop=data_stru->X_stop;
  printf("start:%lf, stop: %lf\n",data_stru->X_start,data_stru->X_stop);


  long int iele=0;
  double cumulative_expo=0.;
  while(1 == fread(&expo_appo,sizeof(double),1,datafile)){
    if(differentialX_on_input_flag){
      //printf("exp: %lE %lE \n",expo_appo,cumulative_expo);
      *(data_stru->X_poi+iele)=cumulative_expo;
      cumulative_expo+=expo_appo;
    }else{
      *(data_stru->X_poi+iele)=expo_appo;
    }
    iele++;
  }
  
  fclose(datafile);

  return;
}