#include "scanstat.h"
#include <fstream>
#include <list>
#include <cmath>
#include <cstdlib>
#include <memory>
#include <limits>
#include<vector>
#include <algorithm>


using namespace std;

void scanstat::computeQuantile(double y)
{  
  
  double lambda = double(_NbOcc) / double(_lengthS);
  double q      = double(_NbOcc) / double(_sumW);
  double p      = 1-q;
  double mu     = 1/q;
  double x      = y/_u;
  double a      = x/lambda;
  double z      = ((2*a*p+q)-sqrt(pow(2*a*p+q,2)-4*pow(a,2)*pow(p,2)))/(2*a*pow(p,2));
  double thetax = log(z);
  double alphax = lambda*z*q/(1-z*p);
  double M      = (1+p*z)/pow(1-p*z,2);
  double vux    = (1-p*z)/(1-p);
  double Ix     = -(alphax-lambda)+thetax*x;
  double pi     = 3.14159265358979323846;
  
  _log1mPval = -(_lengthS-_u)*vux*exp(-_u*Ix)*(x-lambda*mu)/sqrt(2*pi*_u*alphax*M);

}

void scanstat::computeScans(vector<int> Position,vector<double> Weights)
{  
  
  vector<double> WS(_lengthS, 0);  
  vector<double> Wcumsum(_lengthS, 0);

  vector<double> S(_lengthS, 0);  
  vector<double> cumsum(_lengthS, 0);

  for (int i = 0; i < _NbOcc ; i++){
    WS[Position[i]] = Weights[i];
    S[Position[i]]  = 1;
    _sumW += Weights[i];
  }  
 
  for (int j = 0; j < _lengthS ; j++){
    Wcumsum[j] = Wcumsum[j-1]+WS[j];
    cumsum[j]  = cumsum[j-1]+S[j];    
  }

  for (int j=0;j<(_lengthS-_u);j++){
    _WScan[j] = (Wcumsum[j+_u]-Wcumsum[j]);
    _Scan[j] = (cumsum[j+_u]-cumsum[j]);
  }
  
  _Mnu =*( max_element( _WScan.begin(), _WScan.end() )) ;


}


scanstat::scanstat(int nbocc, int lengthseq, int w)
{
  _NbOcc    = nbocc; 
  _lengthS  = lengthseq;
  _u        = w;
  _Scan.assign(_lengthS-_u,0);
  _WScan.assign(_lengthS-_u,0);

  int _Mnu          = 0;
  int _sumW         = 0;
  double _log1mPval = 0;


}

/** destructor **/
scanstat::~scanstat(){}

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