<commit_before>/**
 * Copyright (c) 2015, <NAME>, Quintelligence d.o.o. and contributors
 * All rights reserved.
 * 
 * This source code is licensed under the FreeBSD license found in the
 * LICENSE file in the root directory of this source tree.
 */

/////////////////////////////////////////////////
// Mathematical-Utilities
double TMath::E=2.71828182845904523536;
double TMath::Pi=3.14159265358979323846;
double TMath::LogOf2=log(double(2));

/////////////////////////////////////////////////
// Special-Functions
void TSpecFunc::GammaPSeries/*gser*/(
 double& gamser, const double& a, const double& x, double& gln){
  static const int ITMAX=100;
  static const double EPS=3.0e-7;
  int n;
  double sum, del, ap;

  gln=LnGamma(a);
  if (x <= 0.0){
    EAssert(x>=0); /*if (x < 0.0) nrerror("x less than 0 in routine gser");*/
    gamser=0.0;
    return;
  } else {
    ap=a;
    del=sum=1.0/a;
    for (n=1; n<=ITMAX; n++){
      ++ap;
      del *= x/ap;
      sum += del;
      if (fabs(del) < fabs(sum)*EPS){
        gamser=sum*exp(-x+a*log(x)-(gln));
        return;
      }
    }
    Fail; /*nrerror("a too large, ITMAX too small in routine gser");*/
    return;
  }
}

void TSpecFunc::GammaQContFrac/*gcf*/(
 double& gammcf, const double& a, const double& x, double& gln){
  static const int ITMAX=100;
  static const double EPS=3.0e-7;
  static const double  FPMIN=1.0e-30;
  int i;
  double an, b, c, d, del, h;

  gln=LnGamma(a);
  b=x+1.0-a;
  c=1.0/FPMIN;
  d=1.0/b;
  h=d;
  for (i=1;i<=ITMAX;i++){
    an = -i*(i-a);
    b += 2.0;
    d=an*d+b;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=b+an/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    del=d*c;
    h *= del;
    if (fabs(del-1.0) < EPS) break;
  }
  EAssert(i<=ITMAX);
  /*if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");*/
  gammcf=exp(-x+a*log(x)-(gln))*h;
}

double TSpecFunc::GammaQ/*gammq*/(const double& a, const double& x){
  EAssert((x>=0)&&(a>0));
  double gamser, gammcf, gln;
  if (x<(a+1.0)){
    GammaPSeries(gamser,a,x,gln);
    return 1.0-gamser;
  } else {
    GammaQContFrac(gammcf,a,x,gln);
    return gammcf;
  }
}

double TSpecFunc::LnGamma/*gammln*/(const double& xx){
  double x, y, tmp, ser;
  static double cof[6]={76.18009172947146,-86.50532032941677,
          24.01409824083091,-1.231739572450155,
          0.1208650973866179e-2,-0.5395239384953e-5};
  int j;

  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  return -tmp+log(2.5066282746310005*ser/x);
}

double TSpecFunc::LnComb(const int& n, const int& k){
  return LnGamma(n+1)-LnGamma(k+1)-LnGamma(n-k+1);
}

double TSpecFunc::BetaCf(const double& a, const double& b, const double& x){
  static const double MAXIT=200;
  static const double EPS=3.0e-7;
  static const double FPMIN=1.0e-30;
  int m,m2;
  double aa,c,d,del,h,qab,qam,qap;

  qab=a+b;
  qap=a+1.0;
  qam=a-1.0;
  c=1.0;
  d=1.0-qab*x/qap;
  if (fabs(d) < FPMIN) d=FPMIN;
  d=1.0/d;
  h=d;
  for (m=1;m<=MAXIT;m++) {
    m2=2*m;
    aa=m*(b-m)*x/((qam+m2)*(a+m2));
    d=1.0+aa*d;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    h *= d*c;
    aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
    d=1.0+aa*d;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    del=d*c;
    h *= del;
    if (fabs(del-1.0) < EPS) break;
  }
  if (m > MAXIT){Fail;}// a or b too big, or MAXIT too small in betacf
  return h;
}

double TSpecFunc::BetaI(const double& a, const double& b, const double& x){
  double bt;

  if (x < 0.0 || x > 1.0){Fail;} // Bad x in routine betai
  if (x == 0.0 || x == 1.0) bt=0.0;
  else
    bt=exp(LnGamma(a+b)-LnGamma(a)-LnGamma(b)+a*log(x)+b*log(1.0-x));
  if (x < (a+1.0)/(a+b+2.0))
    return bt*BetaCf(a,b,x)/a;
  else
    return 1.0-bt*BetaCf(b,a,1.0-x)/b;
}

void TSpecFunc::LinearFit(
 const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // y = a + bx :: SigA (SigB): A's uncertainty; Chi2: std dev on all points
  int i;
  double t, sxoss, sx = 0.0, sy = 0.0, st2 = 0.0, ss, sigdat;

  A = B = SigA = SigB = Chi2 = 0.0;
  for (i = 0; i < XY.Len(); i++) {
    sx += XY[i].Val1;
    sy += XY[i].Val2;
  }
  ss = XY.Len();
  sxoss = sx / ss;
  for (i = 0; i <XY.Len(); i++) {
    t = XY[i].Val1 - sxoss;
    st2 += t*t;
    B += t * XY[i].Val2;
  }
  B /= st2;
  A = (sy - sx * B) / ss;
  SigA = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
  SigB = sqrt(1.0 / st2);
  for (i = 0; i < XY.Len(); i++)
    Chi2 += TMath::Sqr(XY[i].Val2 - A - B * XY[i].Val1);
  sigdat = sqrt(Chi2 / (XY.Len() - 2));
  SigA *= sigdat;
  SigB *= sigdat;

  // calculate R squared
  { double N = XY.Len(), sXY=0.0, sX=0.0, sY=0.0, sSqX=0.0, sSqY=0.0;
  for (int s =0; s < XY.Len(); s++) {
    sX += XY[s].Val1;  sY += XY[s].Val2;
    sXY += XY[s].Val1 * XY[s].Val2;
    sSqX += TMath::Sqr(XY[s].Val1);
    sSqY += TMath::Sqr(XY[s].Val2);
  }
  R2 = TMath::Sqr(N*sXY - sX*sY) / ((N*sSqX - sX*sX) * (N*sSqY - sY*sY)); }
  if (1.1 < R2 || -1.1 > R2) R2 = 0.0;
  if (_isnan(A) || ! _finite(A)) A = 0.0;
  if (_isnan(B) || ! _finite(B)) B = 0.0;
}

void TSpecFunc::PowerFit(const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // y = a x^b :: SigA (SigB): A's uncertainty; Chi2: std dev on all points
  // log fit
  double AA, BB;
  TFltPrV LogXY(XY.Len(), 0);
  for (int s = 0; s < XY.Len(); s++) {
    LogXY.Add(TFltPr(log((double)XY[s].Val1), log((double)XY[s].Val2)));
  }
  TSpecFunc::LinearFit(LogXY, AA, BB, SigA, SigB, Chi2, R2);
  A = exp(AA);  B = BB;
  if (_isnan(AA) || ! _finite(AA)) A = 0.0;
  if (_isnan(BB) || ! _finite(BB)) B = 0.0;
}

void TSpecFunc::LogFit(const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // Y = A + B*log(X)
  TFltPrV LogXY(XY.Len(), 0);
  for (int s = 0; s < XY.Len(); s++) {
    LogXY.Add(TFltPr(log((double)XY[s].Val1), XY[s].Val2));
  }
  TSpecFunc::LinearFit(LogXY, A, B, SigA, SigB, Chi2, R2);
}

void TSpecFunc::ExpFit(const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // Y = A * exp(B*X)
  TFltPrV XLogY(XY.Len(), 0);
  double AA, BB;
  for (int s = 0; s < XY.Len(); s++) {
    XLogY.Add(TFltPr(XY[s].Val1, log((double)XY[s].Val2)));
  }
  TSpecFunc::LinearFit(XLogY, AA, BB, SigA, SigB, Chi2, R2);
  A = exp(AA);
  B = BB;
}

double TSpecFunc::Entropy(const TIntV& ValV) {
  TFltV NewValV(ValV.Len());
  for (int i = 0; i < ValV.Len(); i++) { NewValV[i] = ValV[i]; }
  return Entropy(NewValV);
}

// Entropy of a distribution: ValV[i] contains the number of events i
double TSpecFunc::Entropy(const TFltV& ValV) {
  double Sum=0, Ent=0;
  for (int i = 0; i < ValV.Len(); i++) {
    const double& Val = ValV[i];
    if (Val > 0.0) { Ent -= Val * log(Val);  Sum += Val; }
  }
  if (Sum > 0.0) {
    Ent /= Sum;
    Ent += log(Sum);
    Ent /= TMath::LogOf2;
  } else return 1.0;
  return Ent;
}

double TSpecFunc::Entropy(const double& Prob) {
	if (Prob == 0 || Prob == 1) { return 0; }
	return -(Prob*TMath::Log2(Prob) + (1 - Prob)*TMath::Log2(1 - Prob));
}

void TSpecFunc::EntropyFracDim(const TIntV& ValV, TFltV& EntropyV) {
  TFltV NewValV(ValV.Len());
  for (int i = 0; i < ValV.Len(); i++) { 
    EAssert(ValV[i]==1 || ValV[i] == 0);
    NewValV[i] = ValV[i]; 
  }
  EntropyFracDim(NewValV, EntropyV);
}

// Entropy fractal dimension. Input is a vector {0,1}^n. 
// Where 0 means the event did not occur, and 1 means it occured.
// Works exactly as Mengzi Wang's code.
void TSpecFunc::EntropyFracDim(const TFltV& ValV, TFltV& EntropyV) {
  TFltV ValV1, ValV2;
  int Pow2 = 1;
  while (2*Pow2 <= ValV.Len()) { Pow2 *= 2; }
  ValV1.Gen(Pow2);
  for (int i = 0; i < Pow2; i++) { ValV1[i] = ValV[i]; 
    EAssert(ValV[i]==1.0 || ValV[i] == 0.0); }
  EntropyV.Clr();
  EntropyV.Add(Entropy(ValV1)); // 2^Pow2 windows
  while (ValV1.Len() > 2) {
    ValV2.Gen(ValV1.Len() / 2);
    for (int i = 0; i < ValV1.Len(); i++) {
      ValV2[i/2] += ValV1[i];
    }
    EntropyV.Add(Entropy(ValV2));
    ValV1.MoveFrom(ValV2);
  }
  EntropyV.Reverse();
}

// solves for p: B = p * log2(p) + (1-p) * log2(1-p)
double TSpecFunc::EntropyBias(const double& B){
  static TFltFltH BToP;
  if (BToP.Empty()) {
    for (double p = 0.5; p < 1.0; p +=0.0001) {
      double H = p * log(p) + (1.0-p) * log(1.0 - p);
      H = -H / log(2.0);
      BToP.AddDat(TMath::Round(H, 3), p);
    }
  }
  if (BToP.IsKey(TMath::Round(B, 3))) { return BToP.GetDat(TMath::Round(B, 3)); }
  else { return -1.0; }
}

// MLE of the power-law coefficient
double TSpecFunc::GetPowerCoef(const TFltV& XValV, double MinX) {
  for (int i = 0; MinX <= 0.0 && i < XValV.Len(); i++) { 
    MinX = XValV[i]; }
  EAssert(MinX > 0.0);
  double LnSum=0.0;
  for (int i = 0; i < XValV.Len(); i++) {
    if (XValV[i].Val < MinX) continue;
    LnSum += log(XValV[i] / MinX);
  }
  return 1.0 + double(XValV.Len()) / LnSum;
}

double TSpecFunc::GetPowerCoef(const TFltPrV& XValCntV, double MinX) {
  for (int i = 0; MinX <= 0.0 && i < XValCntV.Len(); i++) { 
    MinX = XValCntV[i].Val1; }
  EAssert(MinX > 0.0);
  double NSamples=0.0, LnSum=0.0;
  for (int i = 0; i < XValCntV.Len(); i++) {
    if (XValCntV[i].Val1() < MinX) continue;
    LnSum += XValCntV[i].Val2 * log(XValCntV[i].Val1 / MinX);
    NSamples += XValCntV[i].Val2;
  }
  return 1.0 + NSamples / LnSum;
}

double TSpecFunc::StudentCdf(const double& Val, const int& Df) {
	const double x = TFlt::Abs(Val) / TMath::Sqrt(Df);
	const double OnePlusSqX = 1 + x*x;

	const int Iters = Df / 2;
	double u = 1, s = 0;

	if ((Df & 1) == 0) { // if df is even
		for (int i = 1; i <= Iters; ++i) {
			s += u;
			u *= (1.0 - 1.0 / (2 * i)) / OnePlusSqX;
		}

		return 0.5 - 0.5 * (x / TMath::Sqrt(OnePlusSqX)) * s;
	}
	else { // if Df is odd
		for (int j = 2; j < Iters + 2; ++j) {
			s += u;
			u *= (1.0 - 1.0 / (2 * j - 1)) / OnePlusSqX;
		}

		return 0.5 - (x / OnePlusSqX * s + atan(x)) / TMath::Pi;
	}
}

double TSpecFunc::StudentCdf(const double& Val, const double& Mean,
		const double& Std, const int& Df) {
	return StudentCdf((Val - Mean) / sqrt(Std / (Df + 1)), Df);
}

/////////////////////////////////////////////////
// Statistical-Moments
TMom::TMom(const TFltV& _ValV):
  //WgtV(_ValV.Len(), 0), ValV(_ValV.Len(), 0),
  ValWgtV(_ValV.Len(), 0),
  SumW(), ValSumW(),
  UsableP(false), UnusableVal(-1),
  Mn(), Mx(),
  Mean(), Vari(), SDev(), SErr(),
  Median(), Quart1(), Quart3(),
  DecileV(), PercentileV(){
  for (int ValN=0; ValN<_ValV.Len(); ValN++){Add(_ValV[ValN], 1);}
  Def();
}

void TMom::Def(){
  EAssert(!DefP); DefP=true;
  UsableP=(SumW>0)&&(ValWgtV.Len()>0);
  if (UsableP){
    // Mn, Mx
    Mn=ValWgtV[0].Val1;
    Mx=ValWgtV[0].Val1;<commit_msg>fixed studentCdf function for 4 input parameters in xmath.cpp<commit_after>/**
 * Copyright (c) 2015, <NAME>, Quintelligence d.o.o. and contributors
 * All rights reserved.
 * 
 * This source code is licensed under the FreeBSD license found in the
 * LICENSE file in the root directory of this source tree.
 */

/////////////////////////////////////////////////
// Mathematical-Utilities
double TMath::E=2.71828182845904523536;
double TMath::Pi=3.14159265358979323846;
double TMath::LogOf2=log(double(2));

/////////////////////////////////////////////////
// Special-Functions
void TSpecFunc::GammaPSeries/*gser*/(
 double& gamser, const double& a, const double& x, double& gln){
  static const int ITMAX=100;
  static const double EPS=3.0e-7;
  int n;
  double sum, del, ap;

  gln=LnGamma(a);
  if (x <= 0.0){
    EAssert(x>=0); /*if (x < 0.0) nrerror("x less than 0 in routine gser");*/
    gamser=0.0;
    return;
  } else {
    ap=a;
    del=sum=1.0/a;
    for (n=1; n<=ITMAX; n++){
      ++ap;
      del *= x/ap;
      sum += del;
      if (fabs(del) < fabs(sum)*EPS){
        gamser=sum*exp(-x+a*log(x)-(gln));
        return;
      }
    }
    Fail; /*nrerror("a too large, ITMAX too small in routine gser");*/
    return;
  }
}

void TSpecFunc::GammaQContFrac/*gcf*/(
 double& gammcf, const double& a, const double& x, double& gln){
  static const int ITMAX=100;
  static const double EPS=3.0e-7;
  static const double  FPMIN=1.0e-30;
  int i;
  double an, b, c, d, del, h;

  gln=LnGamma(a);
  b=x+1.0-a;
  c=1.0/FPMIN;
  d=1.0/b;
  h=d;
  for (i=1;i<=ITMAX;i++){
    an = -i*(i-a);
    b += 2.0;
    d=an*d+b;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=b+an/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    del=d*c;
    h *= del;
    if (fabs(del-1.0) < EPS) break;
  }
  EAssert(i<=ITMAX);
  /*if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");*/
  gammcf=exp(-x+a*log(x)-(gln))*h;
}

double TSpecFunc::GammaQ/*gammq*/(const double& a, const double& x){
  EAssert((x>=0)&&(a>0));
  double gamser, gammcf, gln;
  if (x<(a+1.0)){
    GammaPSeries(gamser,a,x,gln);
    return 1.0-gamser;
  } else {
    GammaQContFrac(gammcf,a,x,gln);
    return gammcf;
  }
}

double TSpecFunc::LnGamma/*gammln*/(const double& xx){
  double x, y, tmp, ser;
  static double cof[6]={76.18009172947146,-86.50532032941677,
          24.01409824083091,-1.231739572450155,
          0.1208650973866179e-2,-0.5395239384953e-5};
  int j;

  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  return -tmp+log(2.5066282746310005*ser/x);
}

double TSpecFunc::LnComb(const int& n, const int& k){
  return LnGamma(n+1)-LnGamma(k+1)-LnGamma(n-k+1);
}

double TSpecFunc::BetaCf(const double& a, const double& b, const double& x){
  static const double MAXIT=200;
  static const double EPS=3.0e-7;
  static const double FPMIN=1.0e-30;
  int m,m2;
  double aa,c,d,del,h,qab,qam,qap;

  qab=a+b;
  qap=a+1.0;
  qam=a-1.0;
  c=1.0;
  d=1.0-qab*x/qap;
  if (fabs(d) < FPMIN) d=FPMIN;
  d=1.0/d;
  h=d;
  for (m=1;m<=MAXIT;m++) {
    m2=2*m;
    aa=m*(b-m)*x/((qam+m2)*(a+m2));
    d=1.0+aa*d;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    h *= d*c;
    aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
    d=1.0+aa*d;
    if (fabs(d) < FPMIN) d=FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < FPMIN) c=FPMIN;
    d=1.0/d;
    del=d*c;
    h *= del;
    if (fabs(del-1.0) < EPS) break;
  }
  if (m > MAXIT){Fail;}// a or b too big, or MAXIT too small in betacf
  return h;
}

double TSpecFunc::BetaI(const double& a, const double& b, const double& x){
  double bt;

  if (x < 0.0 || x > 1.0){Fail;} // Bad x in routine betai
  if (x == 0.0 || x == 1.0) bt=0.0;
  else
    bt=exp(LnGamma(a+b)-LnGamma(a)-LnGamma(b)+a*log(x)+b*log(1.0-x));
  if (x < (a+1.0)/(a+b+2.0))
    return bt*BetaCf(a,b,x)/a;
  else
    return 1.0-bt*BetaCf(b,a,1.0-x)/b;
}

void TSpecFunc::LinearFit(
 const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // y = a + bx :: SigA (SigB): A's uncertainty; Chi2: std dev on all points
  int i;
  double t, sxoss, sx = 0.0, sy = 0.0, st2 = 0.0, ss, sigdat;

  A = B = SigA = SigB = Chi2 = 0.0;
  for (i = 0; i < XY.Len(); i++) {
    sx += XY[i].Val1;
    sy += XY[i].Val2;
  }
  ss = XY.Len();
  sxoss = sx / ss;
  for (i = 0; i <XY.Len(); i++) {
    t = XY[i].Val1 - sxoss;
    st2 += t*t;
    B += t * XY[i].Val2;
  }
  B /= st2;
  A = (sy - sx * B) / ss;
  SigA = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
  SigB = sqrt(1.0 / st2);
  for (i = 0; i < XY.Len(); i++)
    Chi2 += TMath::Sqr(XY[i].Val2 - A - B * XY[i].Val1);
  sigdat = sqrt(Chi2 / (XY.Len() - 2));
  SigA *= sigdat;
  SigB *= sigdat;

  // calculate R squared
  { double N = XY.Len(), sXY=0.0, sX=0.0, sY=0.0, sSqX=0.0, sSqY=0.0;
  for (int s =0; s < XY.Len(); s++) {
    sX += XY[s].Val1;  sY += XY[s].Val2;
    sXY += XY[s].Val1 * XY[s].Val2;
    sSqX += TMath::Sqr(XY[s].Val1);
    sSqY += TMath::Sqr(XY[s].Val2);
  }
  R2 = TMath::Sqr(N*sXY - sX*sY) / ((N*sSqX - sX*sX) * (N*sSqY - sY*sY)); }
  if (1.1 < R2 || -1.1 > R2) R2 = 0.0;
  if (_isnan(A) || ! _finite(A)) A = 0.0;
  if (_isnan(B) || ! _finite(B)) B = 0.0;
}

void TSpecFunc::PowerFit(const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // y = a x^b :: SigA (SigB): A's uncertainty; Chi2: std dev on all points
  // log fit
  double AA, BB;
  TFltPrV LogXY(XY.Len(), 0);
  for (int s = 0; s < XY.Len(); s++) {
    LogXY.Add(TFltPr(log((double)XY[s].Val1), log((double)XY[s].Val2)));
  }
  TSpecFunc::LinearFit(LogXY, AA, BB, SigA, SigB, Chi2, R2);
  A = exp(AA);  B = BB;
  if (_isnan(AA) || ! _finite(AA)) A = 0.0;
  if (_isnan(BB) || ! _finite(BB)) B = 0.0;
}

void TSpecFunc::LogFit(const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // Y = A + B*log(X)
  TFltPrV LogXY(XY.Len(), 0);
  for (int s = 0; s < XY.Len(); s++) {
    LogXY.Add(TFltPr(log((double)XY[s].Val1), XY[s].Val2));
  }
  TSpecFunc::LinearFit(LogXY, A, B, SigA, SigB, Chi2, R2);
}

void TSpecFunc::ExpFit(const TVec<TFltPr>& XY, double& A, double& B,
 double& SigA, double& SigB, double& Chi2, double& R2) {
  // Y = A * exp(B*X)
  TFltPrV XLogY(XY.Len(), 0);
  double AA, BB;
  for (int s = 0; s < XY.Len(); s++) {
    XLogY.Add(TFltPr(XY[s].Val1, log((double)XY[s].Val2)));
  }
  TSpecFunc::LinearFit(XLogY, AA, BB, SigA, SigB, Chi2, R2);
  A = exp(AA);
  B = BB;
}

double TSpecFunc::Entropy(const TIntV& ValV) {
  TFltV NewValV(ValV.Len());
  for (int i = 0; i < ValV.Len(); i++) { NewValV[i] = ValV[i]; }
  return Entropy(NewValV);
}

// Entropy of a distribution: ValV[i] contains the number of events i
double TSpecFunc::Entropy(const TFltV& ValV) {
  double Sum=0, Ent=0;
  for (int i = 0; i < ValV.Len(); i++) {
    const double& Val = ValV[i];
    if (Val > 0.0) { Ent -= Val * log(Val);  Sum += Val; }
  }
  if (Sum > 0.0) {
    Ent /= Sum;
    Ent += log(Sum);
    Ent /= TMath::LogOf2;
  } else return 1.0;
  return Ent;
}

double TSpecFunc::Entropy(const double& Prob) {
	if (Prob == 0 || Prob == 1) { return 0; }
	return -(Prob*TMath::Log2(Prob) + (1 - Prob)*TMath::Log2(1 - Prob));
}

void TSpecFunc::EntropyFracDim(const TIntV& ValV, TFltV& EntropyV) {
  TFltV NewValV(ValV.Len());
  for (int i = 0; i < ValV.Len(); i++) { 
    EAssert(ValV[i]==1 || ValV[i] == 0);
    NewValV[i] = ValV[i]; 
  }
  EntropyFracDim(NewValV, EntropyV);
}

// Entropy fractal dimension. Input is a vector {0,1}^n. 
// Where 0 means the event did not occur, and 1 means it occured.
// Works exactly as Mengzi Wang's code.
void TSpecFunc::EntropyFracDim(const TFltV& ValV, TFltV& EntropyV) {
  TFltV ValV1, ValV2;
  int Pow2 = 1;
  while (2*Pow2 <= ValV.Len()) { Pow2 *= 2; }
  ValV1.Gen(Pow2);
  for (int i = 0; i < Pow2; i++) { ValV1[i] = ValV[i]; 
    EAssert(ValV[i]==1.0 || ValV[i] == 0.0); }
  EntropyV.Clr();
  EntropyV.Add(Entropy(ValV1)); // 2^Pow2 windows
  while (ValV1.Len() > 2) {
    ValV2.Gen(ValV1.Len() / 2);
    for (int i = 0; i < ValV1.Len(); i++) {
      ValV2[i/2] += ValV1[i];
    }
    EntropyV.Add(Entropy(ValV2));
    ValV1.MoveFrom(ValV2);
  }
  EntropyV.Reverse();
}

// solves for p: B = p * log2(p) + (1-p) * log2(1-p)
double TSpecFunc::EntropyBias(const double& B){
  static TFltFltH BToP;
  if (BToP.Empty()) {
    for (double p = 0.5; p < 1.0; p +=0.0001) {
      double H = p * log(p) + (1.0-p) * log(1.0 - p);
      H = -H / log(2.0);
      BToP.AddDat(TMath::Round(H, 3), p);
    }
  }
  if (BToP.IsKey(TMath::Round(B, 3))) { return BToP.GetDat(TMath::Round(B, 3)); }
  else { return -1.0; }
}

// MLE of the power-law coefficient
double TSpecFunc::GetPowerCoef(const TFltV& XValV, double MinX) {
  for (int i = 0; MinX <= 0.0 && i < XValV.Len(); i++) { 
    MinX = XValV[i]; }
  EAssert(MinX > 0.0);
  double LnSum=0.0;
  for (int i = 0; i < XValV.Len(); i++) {
    if (XValV[i].Val < MinX) continue;
    LnSum += log(XValV[i] / MinX);
  }
  return 1.0 + double(XValV.Len()) / LnSum;
}

double TSpecFunc::GetPowerCoef(const TFltPrV& XValCntV, double MinX) {
  for (int i = 0; MinX <= 0.0 && i < XValCntV.Len(); i++) { 
    MinX = XValCntV[i].Val1; }
  EAssert(MinX > 0.0);
  double NSamples=0.0, LnSum=0.0;
  for (int i = 0; i < XValCntV.Len(); i++) {
    if (XValCntV[i].Val1() < MinX) continue;
    LnSum += XValCntV[i].Val2 * log(XValCntV[i].Val1 / MinX);
    NSamples += XValCntV[i].Val2;
  }
  return 1.0 + NSamples / LnSum;
}

double TSpecFunc::StudentCdf(const double& Val, const int& Df) {
	const double x = TFlt::Abs(Val) / TMath::Sqrt(Df);
	const double OnePlusSqX = 1 + x*x;

	const int Iters = Df / 2;
	double u = 1, s = 0;

	if ((Df & 1) == 0) { // if df is even
		for (int i = 1; i <= Iters; ++i) {
			s += u;
			u *= (1.0 - 1.0 / (2 * i)) / OnePlusSqX;
		}

		return 0.5 - 0.5 * (x / TMath::Sqrt(OnePlusSqX)) * s;
	}
	else { // if Df is odd
		for (int j = 2; j < Iters + 2; ++j) {
			s += u;
			u *= (1.0 - 1.0 / (2 * j - 1)) / OnePlusSqX;
		}

		return 0.5 - (x / OnePlusSqX * s + atan(x)) / TMath::Pi;
	}
}

double TSpecFunc::StudentCdf(const double& Val, const double& Mean,
		const double& Std, const int& Df) {
	return StudentCdf((Val - Mean) / (Std / sqrt(Df + 1)), Df);
}

/////////////////////////////////////////////////
// Statistical-Moments
TMom::TMom(const TFltV& _ValV):
  //WgtV(_ValV.Len(), 0), ValV(_ValV.Len(), 0),
  ValWgtV(_ValV.Len(), 0),
  SumW(), ValSumW(),
  UsableP(false), UnusableVal(-1),
  Mn(), Mx(),
  Mean(), Vari(), SDev(), SErr(),
  Median(), Quart1(), Quart3(),
  DecileV(), PercentileV(){
  for (int ValN=0; ValN<_ValV.Len(); ValN++){Add(_ValV[ValN], 1);}
  Def();
}

void TMom::Def(){
  EAssert(!DefP); DefP=true;
  UsableP=(SumW>0)&&(ValWgtV.Len()>0);
  if (UsableP){
    // Mn, Mx
    Mn=ValWgtV[0].Val1;
    Mx=ValWgtV[0].Val1;
{'fraction_non_alphanumeric': 0.2232817394508937, 'fraction_numerical': 0.049705177814630554, 'mean_word_length': 2.78205574912892, 'pattern_counts': {'":': 0, '<': 105, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6504302', 'n_tokens_mistral': 10765, 'n_tokens_neox': 10118, 'n_words': 2319}