#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include <sys/types.h>
#include <time.h>

#include <iostream>
#include <string>
#include <algorithm>
#include <iterator>
#include <fstream>
using namespace std;

#define MAX_ROOTPATH_LEN 100
#define MAX_FULLPATH_LEN 150
#define PATH_DATASET "PROC/DATASETS/"
#define PATH_NORMDATA "PROC/NORMDATA/"
#define MAX_NUM_SETS 100
#define DEBUG 1
#define SUBSET 10000

void quickSort(short unsigned int*, unsigned long long int);

int main(int argc, char *argv[]){

  char raw_fn[MAX_FULLPATH_LEN], 
    dataset_fn[MAX_FULLPATH_LEN],
    normdata_fn[MAX_FULLPATH_LEN],
    obj_count_fn[MAX_FULLPATH_LEN], 
    ROOT_PATH[MAX_ROOTPATH_LEN],
    curr_fn_index[5],
    curr_set_index[4],
    curr_obj_count[20],
    text[500],
    dataset_name[MAX_FULLPATH_LEN],
    error_string[200];

  int num_frames,
    num_sets,
    num_cols,
    set_index[MAX_NUM_SETS],
    curr_set,
    field_size,
    success,
    curr_frame,
    output_thresh,
    output_flag=0,
    i,
    j,
    temp;

  unsigned int random_index;

  //initialize random number generator
  time_t t1;
  (void) time(&t1);
  srand((unsigned int) t1);


  
  unsigned long long int *col_sums; // max val ullong=18446744073709551615ULL
  unsigned char *row_preserve;
  unsigned long long int num_rows = 0,
    curr_row = 0;
  short unsigned int bead_info[7];
  short unsigned int dataset_row[150];
  short unsigned int ** bead_vals;
  short unsigned int ** subset_bead_vals;
  float ** norm_bead_vals;
  float * base_sums_multiplier;


  FILE *infile,
    *outfile;

  //short unsigned int input_frame=0,
  short unsigned int input_frame=0,
    input_set=0,
    input_size=0,
    input_centroid_x=0,
    input_centroid_y=0,
    input_border=0,
    input_value=0,
    input_border_temp=0;


  //  short unsigned int input_frame[10];
  field_size = sizeof(unsigned short int);
  

  //****************** Parse command-line arguments **********************
  //
  strcpy(ROOT_PATH, argv[1]);
  strcpy(dataset_name, argv[2]);
  num_sets = atoi(argv[3]);
  num_cols = atoi(argv[3]) + 7;
  //num_sets = argc - 3;
  //
  //**********************************************************************

  if((col_sums = (unsigned long long int*) malloc(num_cols * sizeof(unsigned long long int)))==NULL){
    fprintf(stderr, "Could not allocate enough memory for col_sums list.\n");
    exit(42);
  }


  if((base_sums_multiplier = (float*) malloc(num_cols * sizeof(float)))==NULL){
    fprintf(stderr, "Could not allocate enough memory for base_sums_multiplier list.\n");
    exit(42);
  }
  
  if((bead_vals = (short unsigned int**) malloc(num_sets * sizeof(short unsigned int*)))==NULL){
    fprintf(stderr, "Could not allocate enough memory for bead_vals list.\n");
    exit(42);
  }
  
  if((norm_bead_vals = (float**) malloc(num_sets * sizeof(float*)))==NULL){
    fprintf(stderr, "Could not allocate enough memory for norm_bead_vals list.\n");
  }

  if((subset_bead_vals = (short unsigned int**) malloc(num_sets * sizeof(short unsigned int*)))==NULL){
    fprintf(stderr, "Could not allocate enough memory for subset_bead_vals list.\n");
    exit(42);
  }

  fprintf(stderr, "Allocating %iKB for subset_bead_vals array\n", (int) ( (SUBSET*num_sets*sizeof(short unsigned int)) / double(1024)));
  for(i=0; i< num_sets; i++){
    if((*(subset_bead_vals+i) = (short unsigned int*) malloc(SUBSET * sizeof(short unsigned int)))==NULL){
      fprintf(stderr, "Could not allocate enough memory for subset_bead_vals [%d] array.\n",i);
      exit(42);
    }
  }

  
  for(i=0; i<num_cols; i++){
    *(col_sums + i) = 0;
  }


  //********** Open input bead files and output dataset files ************
  //
#ifdef DEBUG
  fprintf(stderr, "Opening input bead files and output dataset file (%d sets).\n",num_sets);
#endif
  
  // Generate filename for current input (dataset) file
  strcpy(dataset_fn, "");
  strcat(dataset_fn, ROOT_PATH);
  strcat(dataset_fn, PATH_DATASET);
  strcat(dataset_fn, dataset_name);
  strcat(dataset_fn, ".dataset");
  
#ifdef DEBUG
  fprintf(stderr, "Opening file %s for input: ", dataset_fn);
#endif
  if((infile = fopen(dataset_fn, "r")) == NULL){
    fprintf(stderr, "ERROR OPENING INFILE.\n");
    exit(42);
  }
#ifdef DEBUG
  fprintf(stderr, "%s\n", strerror(errno));
#endif

  // Generate filename for current output (normdata) file
  strcpy(normdata_fn, "");
  strcat(normdata_fn, ROOT_PATH);
  strcat(normdata_fn, PATH_NORMDATA);
  strcat(normdata_fn, dataset_name);
  strcat(normdata_fn, ".normdata");
  
#ifdef DEBUG
  fprintf(stderr, "Opening file %s for output: ", normdata_fn);
#endif
  if((outfile = fopen(normdata_fn, "w")) == NULL){
    fprintf(stderr, "ERROR OPENING OUTFILE.\n");
    exit(42);
  }
#ifdef DEBUG
  fprintf(stderr, "%s\n", strerror(errno));
#endif


  curr_frame = 0;

  //run through the dataset file and determine the total number of rows (to initialize arrays)
  fprintf(stderr,"Counting number of rows in infile...");
  while(!feof(infile)){
    if((success=fread(&dataset_row, field_size, num_cols, infile))<1){
      //rintf(stderr, "ERROR counting number of rows in infile.\n"); 
      //exit(42);
    }
    num_rows++;
  }
  fprintf(stderr, "%ld\n", num_rows);


  /* Initialize data arrays
   */
  fprintf(stderr, "Allocating memory for %d cols by %ld rows of bead data...", num_sets, num_rows);
  for(i=0; i< num_sets; i++){
    if((*(bead_vals+i) = (short unsigned int*) malloc(num_rows * sizeof(short unsigned int)))==NULL){
      fprintf(stderr, "Could not allocate enough memory for bead_vals [%d] array.\n",i);
      exit(42);
    }
  }
  fprintf(stderr, " done.\n");

  fprintf(stderr, "Allocating memory for %d cols by %ld rows of normalized bead data...", num_sets, num_rows);
  for(i=0; i<num_sets; i++){
    if((*(norm_bead_vals+i) = (float*) malloc(num_rows * sizeof(float)))==NULL){
      fprintf(stderr, "Could not allocate enough memory for norm_bead_vals [%d] array.\n",i);
      exit(42);
    }

  }
  fprintf(stderr, " done.\n");
  
  fprintf(stderr, "Allocating memory for %ld rows of preserve data...", num_rows);
  if((row_preserve = (unsigned char*) malloc(num_rows * sizeof(unsigned char)))==NULL){
    fprintf(stderr, "Could not allocate enough memory for preserve [%ld] array.\n",num_rows);
    exit(42);
  }
  fprintf(stderr, " done.\n");
  

  //run through the dataset file and calculate sums of each column for normalization
  fprintf(stderr, "Loading data into memory...");
  rewind(infile);
  while(!feof(infile)){
    
    *(row_preserve + curr_row) = 1;
    
    if(((success=fread(&bead_info, field_size, 7, infile))<1) && (curr_row < num_rows-2)){
      fprintf(stderr, "ERROR reading bead_info from bead file at row %Ld", curr_row); 
      exit(42);
    }
    
    for(i=0; i< num_sets; i++){
      if(((success=fread(&input_value, field_size, 1, infile))<1)&&(curr_row < num_rows-2)){
	fprintf(stderr, "ERROR reading bead_value [%d] from bead file at row %Ld", i, curr_row); 
	exit(42);
      }
      //      *(col_sums + i) = *(col_sums + i) + input_value;
      *(*(bead_vals + i)+curr_row) = input_value;
      }

    // APPLY DATA FILTERS HERE TO SET row_preserve TO 0 IF WE WANT TO FILTER THE ROW OUT 
    curr_row++;
  }
  fprintf(stderr, " done.\n");
  
  fprintf(stderr, "Filling subset array for sorting... ");
  for(i=0;i<num_sets; i++){
    for(j=0; j<SUBSET; j++){
      random_index = (unsigned int) ( (rand()/(double)RAND_MAX) * (num_rows-1));
      if(random_index >= num_rows) fprintf(stderr, "index problem -> %u\n", random_index);
      *(*(subset_bead_vals + i) + j) = *(*(bead_vals + i) + random_index);
    }
  }
  fprintf(stderr, " done.\n");
  
  fprintf(stderr, "Sorting bead data in memory to determine percentile... ");
  for(i=0; i< num_sets; i++){
    quickSort(*(subset_bead_vals + i), SUBSET);
  }
  fprintf(stderr, " done.\n");

  /*
    for(i=0; i<num_sets; i++){
    fprintf(stderr, "%i\t%i\n", (int) (*(*(subset_bead_vals + i)+100)), (int) (*(*(subset_bead_vals + i)+SUBSET-1)));
  }
  */

  fprintf(stderr, "Subtracting first percentile from each value... ");
  for(i=0; i<num_sets; i++){
    for(j=0; j<num_rows; j++){
      temp = *(*(bead_vals + i)+j) - *(*(subset_bead_vals + i)+100);
      if(temp<0) temp = 0;
      *(*(bead_vals + i)+j) = (short unsigned int) temp;
    }
  }
  fprintf(stderr, " done.\n");

  fprintf(stderr, "Computing sums for all beads in each column...");
  for(i=0; i<num_sets; i++){
    for(j=0; j<num_rows; j++){
      *(col_sums + i) = *(col_sums + i) + *(*(bead_vals + i)+j);
    }
  }
  fprintf(stderr, " done.\n");

  fprintf(stderr, "Normalizing so that for a given base the sums of each cycle are equal...");
  unsigned long long int highest = 0;
  for(i=0; i<num_sets; i=i+4){

    highest = *(col_sums + i);
    if(*(col_sums + i + 1) > highest){
      highest = *(col_sums + i + 1);
    }
    if(*(col_sums + i + 2) > highest){
      highest = *(col_sums + i + 2);
    }
    if(*(col_sums + i + 3) > highest){
      highest = *(col_sums + i + 3);
    }

    *(base_sums_multiplier + i) = (float) (highest / (double)(*(col_sums + i)));
    *(base_sums_multiplier + i + 1) = (float) (highest / (double)(*(col_sums + i + 1)));
    *(base_sums_multiplier + i + 2) = (float) (highest / (double)(*(col_sums + i + 2)));
    *(base_sums_multiplier + i + 3) = (float) (highest / (double)(*(col_sums + i + 3)));
  }

  //all vals now in norm_bead_vals will sum to the same value for each base
  for(i=0; i<num_sets; i++){
    for(j=0; j<num_rows; j++){
      *(*(norm_bead_vals + i) + j) = *(base_sums_multiplier + i) * ( *(*(bead_vals + i) + j));
    }
  }

  for(i=0; i<100; i++){
    fprintf(stderr, "%d %d %f %f %f %f\n", *(*(bead_vals)+i), *(*(bead_vals+1)+i), *(*(norm_bead_vals)+i), *(base_sums_multiplier), *(*(norm_bead_vals+1)+i), *(base_sums_multiplier+1));
  }



  for(i=0; i< num_sets; i++){
    free(*(bead_vals + i));
    free(*(norm_bead_vals + i));
  }
  free(bead_vals);
  free(norm_bead_vals);
}    
  


  //run through the dataset file again and generate a list in memory of all values from the file normalized by the sums
  //include another boolean array with a "preserve" flag to know which rows need to be carried forward (apply thresholds here)
  //calculate M scores for each row by taking the max after normalization
  //calculate median of each cluster and re-assign M-scores
  //calculate CMS score for each bead (1 per 4-cycle)
  //apply CMS threshold
  //re-calculate median of each cluster and re-assign M-scores
  //re-calculate CMS score for each bead (1 per 4-cycle)
  //generate output matrix of chars



/* Quicksort implementation copied from:
http://www.personal.kent.edu/~rmuhamma/Algorithms/MyAlgorithms/Sorting/quickSort.htm
*/
void q_sort(short unsigned int numbers[], long long int left, long long int right){
  long long int pivot;
  long long int l_hold;
  long long int r_hold;

  l_hold = left;
  r_hold = right;
  pivot = numbers[left];
  while (left < right)
  {
    while ((numbers[right] >= pivot) && (left < right))
      right--;
    if (left != right)
    {
      numbers[left] = numbers[right];
      left++;
    }
    while ((numbers[left] <= pivot) && (left < right))
      left++;
    if (left != right)
    {
      numbers[right] = numbers[left];
      right--;
    }
  }
  numbers[left] = (short unsigned int) pivot;
  pivot = left;
  left = l_hold;
  right = r_hold;
  if (left < pivot)
    q_sort(numbers, left, pivot-1);
  if (right > pivot)
    q_sort(numbers, pivot+1, right);
}

void quickSort(short unsigned int numbers[], unsigned long long int array_size)
{
  q_sort(numbers, 0, (long long int) (array_size) - 1);
}