#include "scanstat.h"
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdio.h>
#include <sstream>
#include <stdlib.h>
#include<vector>
#include <algorithm>

using namespace std;

int main(int argc, char* argv[]) {  
  
  vector<double> Score_Data;
  vector<int> Position_Data;  
  string line;
  ifstream f;
  f.open(argv[5]);  
  getline(f, line);

  int i=0;

  while (line.length()){
    istringstream linestream(line);
    int val;
    linestream>>val;
    Position_Data.push_back(val);
    linestream>>val;
    Score_Data.push_back(val);
    line.clear();
    getline(f, line);
    i=i+1;
  }
  f.close();

  int _NbOcc_Data   = Score_Data.size();
  for (int i=0;i<(_NbOcc_Data-1);i++){
    if (Position_Data[i]>Position_Data[i+1]){
      cout << i << endl;
      cerr << "Data not in order" << endl;
      exit(2);
    }
  }

  vector<int> Position;
  vector<double> Score;

  int leftbound = atoi(argv[3]);
  int rightbound = atoi(argv[4]);
  
  for (int i=0;i<_NbOcc_Data;i++){
    if ( (Position_Data[i]>= leftbound) && (Position_Data[i]< rightbound)){
      Position.push_back(Position_Data[i]);
      Score.push_back(Score_Data[i]);         
    }
  }
  
  int _NbOcc   = Score.size();
  int minpos   = Position[0];

  for (int i=0;i<_NbOcc;i++){
    Position[i]= Position[i]-minpos+1-1;
  }

  int _lengthS = Position[_NbOcc-1]+1;
  int _u = atoi(argv[1]);
  scanstat Scan(_NbOcc, _lengthS,_u);

  cout.precision(6);  
  Scan.computeScans(Position, Score);

  Scan.computeQuantile(Scan._Mnu); 

  double alpha    = atof(argv[2]);
  double bound    = log(1-alpha*(1+0.0001));
  double quantile = Scan._Mnu;
  Scan.computeQuantile(quantile); 

  while (Scan._log1mPval > bound){	
    Scan.computeQuantile(quantile); 
    quantile = quantile/1.0005;
  }

  int j=0;
  while(j<(_lengthS-_u)){
	int h       = 0;
	int posmax  = 0;	
	int maxscan = Scan._WScan[j];
	while (Scan._WScan[j+h]>quantile){
		h++;
		if (Scan._WScan[j+h]>=maxscan){
			maxscan = Scan._WScan[j+h];
			posmax = h;
		}
      	}	
	if (h>0){
		cout << j+leftbound+_u/2 << "\t";
		cout << j+h+leftbound+_u/2 << "\t";   
		cout << j+posmax+leftbound+_u/2 << "\t";   
		cout << quantile << "\t";
		cout << maxscan << endl;
      		j = j+h+1;
    	} else{
		j++;
	}
  }
  return 0;
}

/* Author F. Picard franck.picard@univ-lyon1.fr */






