text: /**
   * Copyright (c) 2015, Jozef Stefan Institute, 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.
   */
  
  #include "bd.h"
  
  /////////////////////////////////////////////////
  // Mathematics-Utilities
  class TMath{
  public:
    static double E;
    static double Pi;
    static double LogOf2;
  
    static double Abs(const double& x) { return TFlt::Abs(x); }
    static int Abs(const int& x) { return TInt::Abs(x); }
    static double Inv(const double& x){IAssert(x!=0.0); return (1.0/x);}
    static double Sqr(const double& x){return x*x;}
    static double Sqrt(const double& x){IAssert(!(x<0.0)); return sqrt(x);}
    static double Log(const double& Val){return log(Val);}
    static double Log2(const double& Val){return log(Val)/LogOf2;}
    static double Round(const double& Val){
      return (Val>0)?floor(Val+0.5):ceil(Val-0.5);}
    static double Round(const double & Val, int Decs){
      const double pwr=pow(10.0, Decs); return Round(Val * pwr) / pwr;}
    static int Fac(const int& Val){
      if (Val<=1){return 1;} else {return Val*Fac(Val-1);}}
    static int Choose(const int& N, const int& K){ // binomial coefficient
      return Fac(N)/(Fac(K)*Fac(N-K)); }
    static uint Pow2(const int& pow){return uint(1u<<pow);}
    static double Power(const double& Base, const double& Exponent){
      return exp(log(Base)*Exponent);}
  
    template <typename T>
    static int Sign(const T& Val){return Val<0?-1:(Val>0?1:0);}
  
    template <class T>
    static const T& Mx(const T& LVal, const T& RVal) {
      return LVal > RVal ? LVal : RVal;}
  
    template <class T>
    static const T& Mn(const T& LVal, const T& RVal){
      return LVal < RVal ? LVal : RVal;}
  
    template <class T>
    static const T& Mx(const T& Val1, const T& Val2, const T& Val3) {
      if (Val1 > Val2) {
        if (Val1 > Val3) return Val1;
        else return Val3;
      } else {
        if (Val2 > Val3) return Val2;
        else return Val3;
      }
    }
  
    template <class T>
    static const T& Mn(const T& Val1, const T& Val2, const T& Val3) {
      if(Val1 < Val2) {
        if (Val1 < Val3) return Val1;
        else return Val3;
      } else {
        if (Val2 < Val3) return Val2;
        else return Val3;
      }
    }
  
    template <class T>
    static const T& Median(const T& Val1, const T& Val2, const T& Val3) {
      if (Val1 < Val2) {
        if (Val2 < Val3) return Val2;
        else if (Val3 < Val1) return Val1;
        else return Val3;
      } else {
        if (Val1 < Val3) return Val1;
        else if (Val3 < Val2)  return Val2;
        else return Val3;
      }
    }
  
    template <class T>
    static const T& InRange(const T& Val, const T& Mn, const T& Mx) {
      IAssert(Mn <= Mx); return Val < Mn ? Mn : (Val > Mx ? Mx : Val);}
  
    template <class T>
    static bool IsInRange(const T& Val, const T& Mn, const T& Mx) {
      IAssert(Mn <= Mx); return Val >= Mn && Val <= Mx;}
  
    template <class T>
    static bool IsInEps(const T& Val, const T& Eps) {
      return Val >= -Eps && Val <= Eps;}
  };
  
  /////////////////////////////////////////////////
  // Special-Functions
  class TSpecFunc{
  public:
    static void GammaPSeries/*gser*/(
     double& gamser, const double& a, const double& x, double& gln);
    static void GammaQContFrac/*gcf*/(
     double& gammcf, const double& a, const double& x, double& gln);
    static double GammaQ/*gammq*/(const double& a, const double& x);
    static double LnGamma/*gammln*/(const double& xx);
    static double BetaCf/*betacf*/(
     const double& a, const double& b, const double& x);
    static double BetaI(const double& a, const double& b, const double& x);
  
    static void LinearFit( // Y = A + B*X
     const TVec<TFltPr>& XY, double& A, double& B,
     double& SigA, double& SigB, double& Chi2, double& R2);
    /// Y = A + B*X, uses basic linear algebra operations
    static void LinearFit( // Y = A + B*X
  	  const TFltV& X, const TFltV& Y, double& A, double& B);
    static void PowerFit( // Y = A * X^B
     const TVec<TFltPr>& XY, double& A, double& B,
     double& SigA, double& SigB, double& Chi2, double& R2);
    static void LogFit( // Y = A + B*log(X)
     const TVec<TFltPr>& XY, double& A, double& B,
     double& SigA, double& SigB, double& Chi2, double& R2);
    static void ExpFit( // Y = A * exp(B*X)
     const TVec<TFltPr>& XY, double& A, double& B,
     double& SigA, double& SigB, double& Chi2, double& R2);
  public:
    static double LnComb(const int& n, const int& k);
  public:
    static double Entropy(const TIntV& ValV);
    static double Entropy(const TFltV& ValV);
    static double Entropy(const double& Prob);
    static void EntropyFracDim(const TIntV& ValV, TFltV& EntropyV);
    static void EntropyFracDim(const TFltV& ValV, TFltV& EntropyV);
  public:
    static double EntropyBias(const double& B); // solves for p: B = p*log2(p)+(1-p)log2(1-p)
    //MLE of the power-law coefficient
    static double GetPowerCoef(const TFltV& XValV, double MinX=-1.0); // values (sampled from the distribution)
    static double GetPowerCoef(const TFltPrV& XValCntV, double MinX=-1.0); // (value, count) pairs
  
    // returns the value of the cumulative distribution function of the value specified be Val
    // with Df degrees of freedom
    static double StudentCdf(const double& Val, const int& Df);
    // returns the value of the cumulative distribution function for a value 'Val' sampled from
    // the normal distribution N(Mean, Std) with Df degrees of freedom
    static double StudentCdf(const double& Val, const double& Mean,
  	  const double& Std, const int& Df);
  };
  
  /////////////////////////////////////////////////
  // Statistical-Moments
  ClassTPV(TMom, PMom, TMomV)//{
  private:
    TBool DefP;
    TFltPrV ValWgtV;
    TFlt SumW, ValSumW;
    TInt Vals;
    TBool UsableP;
    TFlt UnusableVal;
    TFlt Mn, Mx;
    TFlt Mean, Vari, SDev, SErr;
    TFlt Median, Quart1, Quart3;
    TFlt Mode;
    TFltV DecileV; // 0=min 1=1.decile, ..., 9=9.decile, 10=max
    TFltV PercentileV; // 0=min 1=1.percentile, ..., 9=9.percentile, 10=max
  public:
    TMom():
      DefP(false), ValWgtV(),
      SumW(), ValSumW(), Vals(),
      UsableP(false), UnusableVal(-1),
      Mn(), Mx(),
      Mean(), Vari(), SDev(), SErr(),
      Median(), Quart1(), Quart3(), Mode(),
      DecileV(), PercentileV(){}
    TMom(const TMom& Mom):
      DefP(Mom.DefP), ValWgtV(Mom.ValWgtV),
      SumW(Mom.SumW), ValSumW(Mom.ValSumW), Vals(Mom.Vals),
      UsableP(Mom.UsableP), UnusableVal(Mom.UnusableVal),
      Mn(Mom.Mn), Mx(Mom.Mx),
      Mean(Mom.Mean), Vari(Mom.Vari), SDev(Mom.SDev), SErr(Mom.SErr),
      Median(Mom.Median), Quart1(Mom.Quart1), Quart3(Mom.Quart3), Mode(Mom.Mode),
      DecileV(Mom.DecileV), PercentileV(Mom.PercentileV){}
    static PMom New(){return PMom(new TMom());}
    static void NewV(TMomV& MomV, const int& Moms){
      MomV.Gen(Moms); for (int MomN=0; MomN<Moms; MomN++){MomV[MomN]=New();}}
    static void NewVV(TVVec<PMom>& MomVV, const int& XMoms, const int& YMoms){
      MomVV.Gen(XMoms, YMoms);
      for (int XMomN=0; XMomN<XMoms; XMomN++){
        for (int YMomN=0; YMomN<YMoms; YMomN++){
          MomVV.At(XMomN, YMomN)=New();}}}
    TMom(const TFltV& _ValV);
    static PMom New(const TFltV& ValV){
      return PMom(new TMom(ValV));}
    TMom(TSIn& SIn):
      DefP(SIn),
      ValWgtV(SIn),
      SumW(SIn), ValSumW(SIn), Vals(SIn),
      UsableP(SIn), UnusableVal(SIn),
      Mn(SIn), Mx(SIn),
      Mean(SIn), Vari(SIn), SDev(SIn), SErr(SIn),
      Median(SIn), Quart1(SIn), Quart3(SIn), Mode(SIn),
      DecileV(SIn), PercentileV(SIn){}
    static PMom Load(TSIn& SIn){return new TMom(SIn);}
    void Save(TSOut& SOut) const {
      DefP.Save(SOut);
      ValWgtV.Save(SOut);
      SumW.Save(SOut); ValSumW.Save(SOut); Vals.Save(SOut);
      UsableP.Save(SOut); UnusableVal.Save(SOut);
      Mn.Save(SOut); Mx.Save(SOut);
      Mean.Save(SOut); Vari.Save(SOut); SDev.Save(SOut); SErr.Save(SOut);
      Median.Save(SOut); Quart1.Save(SOut); Quart3.Save(SOut); Mode.Save(SOut);
      DecileV.Save(SOut); PercentileV.Save(SOut);}
  
    TMom& operator=(const TMom& Mom){
      Assert(!DefP); DefP=Mom.DefP;
      ValWgtV=Mom.ValWgtV;
      SumW=Mom.SumW; ValSumW=Mom.ValSumW; Vals=Mom.Vals;
      UsableP=Mom.UsableP; UnusableVal=Mom.UnusableVal;
      Mn=Mom.Mn; Mx=Mom.Mx;
      Mean=Mom.Mean; Vari=Mom.Vari; SDev=Mom.SDev; SErr=Mom.SErr;
      Median=Mom.Median; Quart1=Mom.Quart1; Quart3=Mom.Quart3; Mode=Mom.Mode;
      DecileV=Mom.DecileV; PercentileV=Mom.PercentileV;
      return *this;}
    bool operator==(const TMom& Mom) const {
      return Vals==Mom.Vals;}
    bool operator<(const TMom& Mom) const {
      return Vals<Mom.Vals;}
  
    // define
    void Def();
    static void DefV(TMomV& MomV){
      for (int MomN=0; MomN<MomV.Len(); MomN++){MomV[MomN]->Def();}}
    static void DefVV(TVVec<PMom>& MomVV){
      for (int XMomN=0; XMomN<MomVV.GetXDim(); XMomN++){
        for (int YMomN=0; YMomN<MomVV.GetYDim(); YMomN++){
          MomVV.At(XMomN, YMomN)->Def();}}}
    bool IsDef() const {return DefP;}
  
    // values
    void Add(const TFlt& Val, const TFlt& Wgt=1){Assert(!DefP);
      ValWgtV.Add(TFltPr(Val, Wgt)); SumW+=Wgt; ValSumW+=Wgt*Val; Vals++;}
    double GetWgt() const {return SumW;}
    int GetVals() const {return Vals;}
    TFlt GetVal(const int& ValN) const {IAssert(!IsDef()); return ValWgtV[ValN].Val1;}
    //const TFltV& GetValV() const {IAssert(!IsDef()); return ValV;} //J:
  
    // usability
    bool IsUsable() const {Assert(DefP); return UsableP;}
    static bool IsUsableV(const TMomV& MomV){
      for (int MomN=0; MomN<MomV.Len(); MomN++){
        if (!MomV[MomN]->IsUsable()){return false;}}
      return true;}
    static bool IsUsableVV(const TVVec<PMom>& MomVV){
      for (int XMomN=0; XMomN<MomVV.GetXDim(); XMomN++){
        for (int YMomN=0; YMomN<MomVV.GetYDim(); YMomN++){
          if (!MomVV.At(XMomN, YMomN)->IsUsable()){return false;}}}
      return true;}
  
    // moments
    double GetMn() const {Assert(DefP&&UsableP); return Mn;}
    double GetMx() const {Assert(DefP&&UsableP); return Mx;}
    double GetExtent() const {Assert(DefP&&UsableP); return Mx-Mn;}
    double GetMean() const {Assert(DefP&&UsableP); return Mean;}
    double GetVari() const {Assert(DefP&&UsableP); return Vari;}
    double GetSDev() const {Assert(DefP&&UsableP); return SDev;}
    double GetSErr() const {Assert(DefP&&UsableP); return SErr;}
    double GetMedian() const {Assert(DefP&&UsableP); return Median;}
    double GetQuart1() const {Assert(DefP&&UsableP); return Quart1;}
    double GetQuart3() const {Assert(DefP&&UsableP); return Quart3;}
    double GetMode() const {Assert(DefP&&UsableP); return Mode;}
    double GetDecile(const int& DecileN) const {
      Assert(DefP&&UsableP); return DecileV[DecileN];}
    double GetPercentile(const int& PercentileN) const {
      Assert(DefP&&UsableP); return PercentileV[PercentileN];}
    double GetByNm(const TStr& MomNm) const;
    TStr GetStrByNm(const TStr& MomNm, char* FmtStr=NULL) const;
  
    // strings
    TStr GetStr(const char& SepCh=' ', const char& DelimCh=':',
     const bool& DecileP=true, const bool& PercentileP=true, const TStr& FmtStr="%g") const;
    static TStr GetNmVStr(const TStr& VarPfx,
     const char& SepCh='\t', const bool& DecileP=true, const bool& PercentileP=true);
    TStr GetValVStr(const char& SepCh='\t', const bool& DecileP=true, const bool& PercentileP=true) const;
  };
  typedef TVVec<PMom> TMomVV;
  typedef THash<TInt, PMom> TIntMomH;
  typedef THash<TInt, TMomV> TIntMomVH;
  typedef THash<TInt, TMomVV> TIntMomVVH;
  
  /////////////////////////////////////////////////
  // Correlation
  ClassTP(TCorr, PCorr)//{
  private:
    int ValVLen;
    double CorrCf;
    double CorrCfPrb;
    double FisherZ;
  public:
    TCorr(){}
    TCorr(const TFltV& ValV1, const TFltV& ValV2);
    static PCorr New(const TFltV& ValV1, const TFltV& ValV2){
      return PCorr(new TCorr(ValV1, ValV2));}
    TCorr(TSIn&){Fail;}
    static PCorr Load(TSIn& SIn){return new TCorr(SIn);}
    void Save(TSOut&){Fail;}
  
    TCorr& operator=(const TCorr&){Fail; return *this;}
  
    double GetCorrCf() const {return CorrCf;}
    double GetCorrCfPrb() const {return CorrCfPrb;}
  
    TStr GetStr() const;
  };
  
  /////////////////////////////////////////////////
  // Statistical Tests
  class TStatTest {
  private:
    static void AveVar(const TFltV& ValV, double& Ave, double& Var);
    static double KsProb(const double& Alam);
  public:
    static void ChiSquareOne(
     const TFltV& ObservedBinV, const TFltV& ExpectedBinV,
     double& ChiSquareVal, double& SignificancePrb);
    static void ChiSquareTwo(
     const TFltV& ObservedBin1V, const TFltV& ObservedBin2V,
     double& ChiSquareVal, double& SignificancePrb);
  
    static void TTest(
     const TFltV& ValV1, const TFltV& ValV2, double& TTestVal, double& TTestPrb);
  
    // Kolmogorov-Smirnov (are two distributions the same)
    static void KsTest(const TFltV& ValV1, const TFltV& ValV2, double& DStat, double& PVal);
    static void KsTest(const TFltPrV& ValCntV1, const TFltPrV& ValCntV2, double& DStat, double& PVal);
  };
  
  /////////////////////////////////////////////////
  // Combinations
  ClassTP(TComb, PComb)//{
  public:
    int Items;
    int Order;
    int CombN;
    TIntV ItemV;
  public:
    TComb(): Items(-1), Order(-1), CombN(-1), ItemV(){}
    TComb(const int& _Items, const int& _Order):
      Items(_Items), Order(_Order), CombN(0), ItemV(){
      IAssert((Order>0)&&(Order<=Items));}
    static PComb New(const int& Items, const int& Order){
      return PComb(new TComb(Items, Order));}
    ~TComb(){}
    TComb(TSIn&){Fail;}
    static PComb Load(TSIn& SIn){return new TComb(SIn);}
    void Save(TSOut&){Fail;}
  
    TComb& operator=(const TComb&){Fail; return *this;}
  
    bool GetNext();
    TIntV& GetItemV(){return ItemV;}
    int GetCombN() const {return CombN;}
    int GetCombs() const;
    void Wr();
  };
  
  /////////////////////////////////////////////////
  // Linear-Regression
  ClassTP(TLinReg, PLinReg)//{
  public:
    TFltVV XVV;
    TFltV YV;
    TFltV SigV;
    int Recs, Vars;
    TFltVV CovarVV; // 1 based
    TFltV CfV; // 1 based
    double ChiSq;
    void GetXV(const int RecN, TFltV& VarV) const {
      VarV.Gen(Vars+1);
      for (int VarN=0; VarN<Vars; VarN++){VarV[VarN+1]=XVV.At(RecN-1, VarN);}
    }
    double GetY(const int RecN) const {return YV[RecN-1];}
    double GetSig(const int RecN) const {return SigV[RecN-1];}
    void NR_covsrt(TFltVV& CovarVV, const int& Vars, const TIntV& ia, const int& mfit);
    void NR_gaussj(TFltVV& a, const int& n, TFltVV& b, const int& m);
    void NR_lfit();
  public:
    TLinReg(){}
    static PLinReg New(
     const TFltVV& XVV, const TFltV& YV, const TFltV& SigV=TFltV());
    ~TLinReg(){}
    TLinReg(TSIn&){Fail;}
    static PLinReg Load(TSIn& SIn){return new TLinReg(SIn);}
    void Save(TSOut&){Fail;}
  
    TLinReg& operator=(const TLinReg&){Fail; return *this;}
  
    int GetRecs() const {return Recs;}
    int GetVars() const {return Vars;}
  
    double GetCf(const int& VarN) const {return CfV[VarN+1];}
    double GetCfUncer(const int& VarN) const {
      return sqrt(double(CovarVV.At(VarN+1, VarN+1)));}
    double GetCovar(const int& VarN1, const int& VarN2) const {
      return CovarVV.At(VarN1+1, VarN2+1);}
  
    double GetChiSq() const {return ChiSq;}
  
    static double LinInterp(const double& x1, const double& y1,
     const double& x2, const double& y2, const double& AtX) _CMPWARN{
      if (x1 == x2) return (y1+y2)/2.0;
      const double k = (y2 - y1) / (x2 - x1);
      return k*(AtX - x1) + y1;
    }
  
    void Wr() const;
  };
  
  /////////////////////////////////////////////////
  // Singular-Value-Decomposition
  ClassTP(TSvd, PSvd)//{
  public:
    TFltVV XVV;
    TFltV YV;
    TFltV SigV;
    int Recs, Vars;
    TFltVV CovarVV; // 1 based
    TFltV CfV; // 1 based
    double ChiSq;
    void GetXV(const int RecN, TFltV& VarV) const {
      VarV.Gen(Vars+1);
      for (int VarN=0; VarN<Vars; VarN++){VarV[VarN+1]=XVV.At(RecN-1, VarN);}
    }
    double GetY(const int RecN) const {return YV[RecN-1];}
    double GetSig(const int RecN) const {return SigV[RecN-1];}
    static double NR_SIGN(double a, double b){return b >= 0.0 ? fabs(a) : -fabs(a);}
    static double NR_FMAX(double maxarg1, double maxarg2){
      return maxarg1 > maxarg2 ? maxarg1 : maxarg2;}
    static int NR_IMIN(int iminarg1, int iminarg2){
      return iminarg1 < iminarg2 ? iminarg1 : iminarg2;}
    static double NR_pythag(double a, double b);
    static void NR_svdcmp(TFltVV& a, int m, int n, TFltV& w, TFltVV& v);
    void NR_svbksb(
     TFltVV& u, TFltV& w, TFltVV& v, int m, int n, TFltV& b, TFltV& x);
    void NR_svdvar(TFltVV& v, int ma, TFltV& w, TFltVV& cvm);
    void NR_svdfit();
  public:
    TSvd(){}
    static PSvd New(
     const TFltVV& XVV, const TFltV& YV, const TFltV& SigV=TFltV());
    ~TSvd(){}
    TSvd(TSIn&){Fail;}
    static PSvd Load(TSIn& SIn){return new TSvd(SIn);}
    void Save(TSOut&){Fail;}
  
    TSvd& operator=(const TSvd&){Fail; return *this;}
  
    int GetRecs() const {return Recs;}
    int GetVars() const {return Vars;}
  
    double GetCf(const int& VarN) const {return CfV[VarN+1];}
    double GetCfUncer(const int& VarN) const {
      return sqrt(double(CovarVV.At(VarN+1, VarN+1)));}
    double GetCovar(const int& VarN1, const int& VarN2) const {
      return CovarVV.At(VarN1+1, VarN2+1);}
  
    double GetChiSq() const {return ChiSq;}
  
    void GetCfV(TFltV& _CfV);
    void GetCfUncerV(TFltV& CfUncerV);
  
    static void Svd(const TFltVV& InMtx, TFltVV& LSingV, TFltV& SingValV, TFltVV& RSingV);
    static void Svd1Based(const TFltVV& InMtx1, TFltVV& LSingV, TFltV& SingValV, TFltVV& RSingV);
  
    void Wr() const;
  };
  
  /////////////////////////////////////////////////
  // Histogram
  class THist {
  private:
      TFlt MnVal;
      TFlt MxVal;
      TIntV BucketV;
      TFlt BucketSize;
      TInt Vals;
  public:
      THist() {}
      THist(const double& _MnVal, const double& _MxVal, const int& Buckets) :
          MnVal(_MnVal), MxVal(_MxVal), BucketV(Buckets) {
          BucketSize = (MxVal == MnVal) ? 1.0 : (1.01 * double(MxVal - MnVal) / double(Buckets));
      }
  
      void Add(const double& Val, const bool& OnlyInP);
  
      int GetVals() const { return Vals; }
      int GetBuckets() const { return BucketV.Len(); }
      double GetBucketMn(const int& BucketN) const { return MnVal + BucketN * BucketSize; }
      double GetBucketMx(const int& BucketN) const { return MnVal + (BucketN + 1) * BucketSize; }
      int GetBucketVal(const int& BucketN) const { return BucketV[BucketN]; }
      double GetBucketValPerc(const int& BucketN) const {
          return (Vals > 0) ? (double(BucketV[BucketN]) / double(Vals)) : 0.0;
      }
  
      void SaveStat(const TStr& ValNm, TSOut& FOut) const;
      void SaveTxt(const TStr& ValNm, const TStr& FNm) const {
          TFOut FOut(FNm); SaveStat(ValNm, FOut);
      }
  };
  
  /////////////////////////////////////////////////
  // Statistics
  class TStatFun {
  public:
      static void ChiSquare(const TFltV& OutValVX, const TFltV& OutValVY, const TInt& Df,
          TFlt& Chi2, TFlt& P);
  };
  