#ifndef TLIBHS_HPP
#define TLIBHS_HPP

#include "TLib.hpp"
#include "TMatrixLib.hpp"
#include <unordered_set>
#include <unordered_map>
#include<future>
#include<thread>

//#define NDEBUG
#include <assert.h>

using namespace std;


//const int NUMCLASSES=15;
/// for additive smoothing
const double alphaSmooth=0.5;
const int minNUMCLASSES=15;
const int dividend=2;

///May 2024: Add the population of reference for each block
struct Block {
    int refpop=0;
    int iniw=-1;
    int endw=0;

    Block(){ };
    Block(const int& refpop, const int& iniw, const int& endw) : refpop(refpop), iniw(iniw), endw(endw){ } // constructor for initialization an anonymous object
};


struct STAT{

    double stat;
    double pval;

};

struct ResultFst{

    string name; // locus name if any
    double pos; // real genome snp position  if any
    double fst;
    double pval;
    double qval;

    ResultFst(){ };
    ResultFst(const double& ps, const double& f,const double& pv,const double& qv, const string& id) : name(id), pos(ps), fst(f), pval(pv),qval(qv) { } // constructor for initialization of anonymous object

    double operator + (ResultFst s){
		return this->fst+s.fst;
	}

	double operator + (double s){
		return (this->fst+s);
	}

    double operator - (ResultFst s){
		return this->fst-s.fst;
	}

	double operator - (double s){
		return (this->fst-s);
	}

};

template <typename Tcol = vector<int> >
inline void HAC_freqs(const Tcol& hac,const int& L, vector<double>& freqs){

		size_t n= hac.size();
		double freq=0;

		for(int h=0;h<=L;++h){
            freq =count(hac.begin(), hac.end(), h); //
            if(freq>0.0)
                freqs.push_back(freq/n);
		}


}

// Compute HAC for each sequence:
//this version includes the candidate snp for the hac computation so that P1 belongs to (0, L-1) and P2 to (1,L)
template <typename RT= int, typename Tcol1 = deque<bool>,typename Tcol2 = vector<bool> >
inline RT computeHAC_incl(const Tcol1& sequence, const Tcol2& reference,const int& ini,const int& ends){

		RT conthac=RT(); // warrants zero initialization for fundamental types
		for(int h=ini;h<ends;++h){
                if(reference[h]^sequence[h]) ++conthac; // counts differences between reference and sequence
		}// close for h

		return conthac;
}


/// acr August 21
/// Compute partitions and hac values from [ini to end) for partitions 1 and 2 from a population sample and given a reference sequence:
/// The computation of HAC DO NOT count the candidate snp
template <typename TR=void, typename RT= double, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeHacs(const int& snp,const int& p, const vector<int>& samplesize, const Tcolcol& metapopsequences, const Tcol2& popreference,const int& ini,const int& ends,vector<RT>& Partition1, vector<RT>& Partition2, int& n1,int& n2){

        int shift = accumulate(samplesize.begin(),samplesize.begin()+p,0);
        /// acr Oct 21
        n1=0;n2=0;

        for(int smpl=0; smpl<samplesize[p]; ++smpl){ // for each haplotype

            if(popreference[snp]==metapopsequences[smpl+shift][snp]){ // first partition
                int conthac1=0;
                n1++; // count the elements in the partition
                for(int h=ini;h<ends;++h){
                    if(h==snp) // ignore the position of the candidate snp
                        continue;
                    if(popreference[h]^metapopsequences[smpl+shift][h]) ++conthac1; // compute the hac value of each haplotype in the sample
                }// close for h

                Partition1.push_back(conthac1);

            } // close if first partition
            else{ // second partition
                int conthac2=0;
                n2++;
                for(int h=ini;h<ends;++h){
                    if(h==snp)
                        continue;
                    if(popreference[h]^metapopsequences[smpl+shift][h]) ++conthac2; // compute the hac value of each haplotype in the sample
                }// close for h
                Partition2.push_back(conthac2);

            }

        } // close for samplesize

}

/// acr Jan 22
// For a given population. Compute hac values from [ini to end) for partitions 1 and 2 from a population sample and given a reference sequence:
/// The partitions are already performed. Computation of HAC DO NOT count the candidate snp
template <typename TR=void, typename RT= double, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeHacs(const int& snp,const Tcolcol& P1,const Tcolcol& P2, const Tcol2& popreference, const int& ini,const int& ends,vector<RT>& Hacs1, vector<RT>& Hacs2){

        ///
        int n1=P1.size(), n2=P2.size();

        if(!Hacs1.empty()){
            vector<RT>().swap(Hacs1);
        }
        if(!Hacs2.empty()){
            vector<RT>().swap(Hacs2);
        }
        /// First partition
        for(int smpl=0; smpl<n1; ++smpl){ // for each haplotype

            int conthac1=0;
            for(int h=ini;h<ends;++h){
                if(h==snp) // ignore the position of the candidate snp
                    continue;
                if(popreference[h]^P1[smpl][h]) ++conthac1; // compute the hac value of each haplotype in the sample
            }// close for h

            Hacs1.push_back(conthac1);
        }
        /// Second partition
        for(int smpl=0; smpl<n2; ++smpl){ // for each haplotype

            int conthac2=0;
            for(int h=ini;h<ends;++h){
                if(h==snp) // ignore the position of the candidate snp
                    continue;
                if(popreference[h]^P2[smpl][h]) ++conthac2; // compute the hac value of each haplotype in the sample
            }// close for h

            Hacs2.push_back(conthac2);
        }

}


/// acr Oct 2023
/// For a given sub population sample P from inisample to endsanmple, and a snp. Compute hac values from [ini to end) for a given reference sequence:
/// The hac are computed jointly for the two alleles of the snp. Computation of HAC DO NOT count the candidate snp
//template <typename TR=void, typename RT= double, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
//inline TR computePopHacs(const int& iniSample, const int& endSample,const int& snp,const Tcolcol& P,const Tcol2& popreference, const int& ini,const int& ends,vector<RT>& Hacs){
//
//        ///
//        //int n=P.size();
//
//        if(!Hacs.empty()){
//            vector<RT>().swap(Hacs);
//        }
//
//        /// Compute hac
//        for(int smpl=iniSample; smpl<endSample; ++smpl){ // for each haplotype
//
//            int conthac1=0;
//            for(int h=ini;h<ends;++h){
//                if(h==snp) // ignore the position of the candidate snp
//                    continue;
//                if(popreference[h]^P[smpl][h]) ++conthac1; // bitwise XOR, compute the hac value of each haplotype in the sample
//            }// close for h
//
//            Hacs.push_back(conthac1);
//        }
//
//}

/// acr April 2026 Controls added
/// For a given sub population sample P from inisample to endsanmple, and a snp. Compute hac values from [ini to end) for a given reference sequence:
/// The hac are computed jointly for the two alleles of the snp. Computation of HAC DO NOT count the candidate snp
template <typename TR=void, typename RT= double,
          typename Tcolcol = deque<deque<bool>>, typename Tcol2 = vector<bool> >
inline TR computePopHacs(const int& iniSample, const int& endSample, const int& snp,
                         const Tcolcol& P, const Tcol2& popreference,
                         const int& ini, const int& ends, vector<RT>& Hacs)
{
    if(!Hacs.empty()){
        vector<RT>().swap(Hacs);
    }

    if(P.empty()){
        cerr << "ERROR computePopHacs: empty data matrix P\n";
        exit(1);
    }

    if(iniSample < 0 || endSample < iniSample || endSample > (int)P.size()){
        cerr << "ERROR computePopHacs: bad sample range. "
             << "iniSample=" << iniSample
             << " endSample=" << endSample
             << " P.size()=" << P.size() << "\n";
        exit(1);
    }

    if(ini < 0 || ends < ini || ends > (int)popreference.size()){
        cerr << "ERROR computePopHacs: bad SNP window in reference. "
             << "ini=" << ini
             << " ends=" << ends
             << " popreference.size()=" << popreference.size()
             << " snp=" << snp << "\n";
        exit(1);
    }

    for(int smpl = iniSample; smpl < endSample; ++smpl){

        if(ends > (int)P[smpl].size()){
            cerr << "ERROR computePopHacs: bad SNP window in sample row. "
                 << "smpl=" << smpl
                 << " ini=" << ini
                 << " ends=" << ends
                 << " P[smpl].size()=" << P[smpl].size()
                 << " popreference.size()=" << popreference.size()
                 << " snp=" << snp << "\n";
            exit(1);
        }

        int conthac1 = 0;
        for(int h = ini; h < ends; ++h){
            if(h == snp)
                continue;
            if(popreference[h] ^ P[smpl][h]) ++conthac1;
        }

        Hacs.push_back(conthac1);
    }
}

/// acr nov 2023
/// compute the number of classes following the Diaconis rule:
///https://stats.stackexchange.com/questions/798/calculating-optimal-number-of-bins-in-a-histogram
/// https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule
/// class width w= 2×IQR×L^(−1/3)
/// K=(hmax-hmin)/w
template <typename TR=double, typename RT= double>
inline TR computeNumberofClasses(const int& pop, const vector<RT>& Hacs){

        ///

        if(Hacs.empty()){
            printval<int>(pop,"EmptyHac.err");
        }

        vector<RT> seq=Hacs;
        double denom=Hacs.size();
        denom=cbrt(denom); //cubic root
        double wide= 2*iqr(seq)/denom;

        int hmax=*max_element(Hacs.begin(),Hacs.end());
        int hmin=*min_element(Hacs.begin(),Hacs.end());

        int K= (hmax-hmin)/wide;

        return K;

}

/// acr Oct 2023
/// For each population hacs computed in computePopHacs make for each (P and Q) NUMCLASSES frequency classes S_i
/// such that each groups NUMCLASSES% of the hac values
/// IMPORTANT: The L values (window sizes-1) must be multiples of NUMCLASSES
/// SP_i = sum from j=u to U #hac[0]_(j-1)/n1 where i in [0,9] and u=1+L/NUMCLASSES*i and U=L/NUMCLASSES*(i+1)
/// SQ_i = sum from j=u to U #hac[1]_(j-1)/n2 where i in [0,9] and u=1+L/NUMCLASSES*i and U=L/NUMCLASSES*(i+1)
/// then for i from 0 to 8 Pi=SP_i and Qi=SQ_i
/// for i=9 Pi=S9+#h[0]_L/n1 and Qi=S9+#h[1]_L/n2
/// Store the classes for each pop y vectors P & Q then compute the Jeffreys divergence
/// For the J divergence Avoid 0's problem by additive smoothing with Jeffreys prior: https://www.wikiwand.com/en/Additive_smoothing
template <typename TR=double, typename RT= double, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR makeHacClassesandJ(const int& NUMCLASSES, const int& L, const vector<vector<RT>>& Hacs, vector<RT>& P, vector<RT>& Q){

        int n1=Hacs[0].size();
        int n2=Hacs[1].size();
        int u=1;
        int U=0;
        int Lfrac=L/NUMCLASSES;
        RT Sp=0.0,Sq=0.0;

        /// April 2026
        if(NUMCLASSES<=1){ return 0;}


        for(int c=0;c<NUMCLASSES;++c){
            u=1+Lfrac*c;
            U=Lfrac*(c+1);

            /// Compute P
            for(const auto&h: Hacs[0]){ // for each haplotype in pop1

                if(h>=(u-1) && h<=(U-1)){

                    Sp+=1;
                }

            }

            /// Compute Q
            for(const auto&h: Hacs[1]){ // for each haplotype in pop2

                if(h>=(u-1) && h<=(U-1)){

                    Sq+=1;
                }

            }
            /// Avoid 0's by additive smoothing (Laplace smoothing): https://www.wikiwand.com/en/Additive_smoothing
            if(c<NUMCLASSES-1){ /// not in the last class

                P[c]=(Sp+alphaSmooth)/(n1+alphaSmooth*NUMCLASSES);
                Sp=0;

                Q[c]=(Sq+alphaSmooth)/(n2+alphaSmooth*NUMCLASSES);
                Sq=0;
            }


        } // close for classes
        //include in the last class the hac values that correspond to h==L
        ///Update the last P class
        Sp+=count(Hacs[0].begin(),Hacs[0].end(),L);
        P[NUMCLASSES-1]=(Sp+alphaSmooth)/(n1+alphaSmooth*NUMCLASSES);
        ///Update the last Q class
        Sq+=count(Hacs[1].begin(),Hacs[1].end(),L);
        /// smooth
        Q[NUMCLASSES-1] = (Sq+alphaSmooth)/(n2+alphaSmooth*NUMCLASSES);

        /// *** Jeffreys divergence
        /*
        NOTE THAT "The Kullback–Leibler divergence is defined only if Q(i)=0 implies P(i)=0 for all i"
        https://math.stackexchange.com/questions/1228408/kullback-leibler-divergence-when-the-q-distribution-has-zero-values
        https://stackoverflow.com/questions/76566737/the-best-approach-to-preventing-the-value-of-the-kullback-leibler-divergence-bec
        https://www.wikiwand.com/en/Additive_smoothing
        */

        double n = n1+n2;
        n= (n1*n2)/n;
        double J=0;
        for(int c=0;c<NUMCLASSES;++c){
            if(P[c]!=Q[c])
                J+=(P[c]-Q[c])*log(P[c]/Q[c]);
        }
        //J=min(J,0.0);
        return n*J;

}

/// acr Jan 22
// Compute hac values from [ini to end) for a given sample P2 from different population samples and given a reference sequence
/// The partitions are already performed. Computation of HAC DO NOT count the candidate snp
/// The second partition P2 is accumulated for the different populations. The window size must be the same.
template <typename TR=void, typename RT= double, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeMetaHacc(const int& snp,const Tcolcol& P2, const Tcol2& popreference, const int& ini,const int& ends,vector<RT>& metaHacs2){

        ///
        int n2=P2.size();

        /// Second partition
        for(int smpl=0; smpl<n2; ++smpl){ // for each haplotype

            int conthac2=0;
            for(int h=ini;h<ends;++h){

                if(h==snp) // ignore the position of the candidate snp
                    continue;
                if(popreference[h]^P2[smpl][h]) ++conthac2; // compute the hac value of each haplotype in the sample
            }// close for h

            metaHacs2.push_back(conthac2);
        }

}

/// acr Jan 22
/// Separate the population sample in two partitions (with major allele or minor allele in the candidate snp)
/// This will allow to compute the haplotype blocks for a given partition
template <typename TR=void, typename RT= double, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR dataPartition(const int& snp,const int& p, const vector<int>& samplesize, const Tcolcol& metapopsequences, const Tcol2& popreference,Tcolcol& Partition1,Tcolcol& Partition2, int& n1, int& n2){

        int shift = accumulate(samplesize.begin(),samplesize.begin()+p,0);
        ///
        n1=0;n2=0;

        if(!Partition1.empty()){
            Tcolcol().swap(Partition1);
        }
        if(!Partition2.empty()){
            Tcolcol().swap(Partition2);
        }

        for(int smpl=0; smpl<samplesize[p]; ++smpl){ // for each haplotype

            if(popreference[snp]==metapopsequences[smpl+shift][snp]){ // first partition
                n1++; // count the elements in the partition

                Partition1.push_back(metapopsequences[smpl+shift]);

            } // close if first partition
            else{ // second partition
                n2++;
                Partition2.push_back(metapopsequences[smpl+shift]);

            }

        } // close for samplesize

}



/// acr July 21
// Compute Y (Cumulative test for Divergent Selection, Carvajal-Rodríguez 2021) for each of 4 partitions P1L, P1R, P2L, P2R:

template <typename TR=void, typename RT= int, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeY(const bool doPL,const bool doPR, const int& snp,const int& p, const vector<int>& samplesize, const Tcolcol& sequences, const Tcol2& reference,const int& iniL,const int& endL,const int& iniR,const int& endR,vector<int>& YPL1,vector<int>& YPR1,vector<int>& YPL2,vector<int>& YPR2, int& n1,int& n2){

        int shift = accumulate(samplesize.begin(),samplesize.begin()+p,0.0);

        /// Better separate left and right because ony could be empty
        /// Left side
        n1=0; n2=0; // sample size of the left side partitions

        if(doPL)
            for(int smpl=0; smpl<samplesize[p]; ++smpl){
//cout<<iniL<<"\t"<<endL<<endl;
                if(reference[snp]==sequences[smpl+shift][snp]){ // first partition
                    n1++;
                    int cl1=0; // counter
                    for(int h=iniL;h<endL;++h){

                            if(reference[h]^sequences[smpl+shift][h]){
                                YPL1[cl1]+=1; // counts differences between reference and sequence
                             }

                        cl1++;
                    }// close for h

                } // close if first partition
                else{ // second partition
                    n2++;
                    int cl2=0; // counter
                    for(int h=iniL;h<endL;++h){
                            if(reference[h]^sequences[smpl+shift][h]) YPL2[cl2]+=1; // counts differences between reference and sequence
                        cl2++;
                    }// close for h

                }

            } // close for samplesize

        /// Right side

        if(doPR)
            for(int smpl=0; smpl<samplesize[p]; ++smpl){

                if(reference[snp]==sequences[smpl+shift][snp]){ // first partition

                    int cr1=0; // counter
                    for(int h=iniR;h<endR;++h){
                            if(reference[h]^sequences[smpl+shift][h]) YPR1[cr1]+=1; // counts differences between reference and sequence
                        cr1++;
                    }// close for h


                } // close if first partition
                else{ // second partition

                    int cr2=0; // counter
                    for(int h=iniR;h<endR;++h){
                            if(reference[h]^sequences[smpl+shift][h]) YPR2[cr2]+=1; // counts differences between reference and sequence
                        cr2++;
                    }// close for h
                }

            } // close for samplesize

}

/// acr July 21
// Compute quantities Z and R (Cumulative test for Divergent Selection, Carvajal-Rodríguez 2021) for each of 4 partitions P1L, P1R, P2L, P2R:

template <typename TR=void, typename RT= int, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeZRd(const bool doPL,const bool doPR,const vector<int>& YPL1,const vector<int>& YPR1,const vector<int>& YPL2,const vector<int>& YPR2, double& RL1, double&RR1, double&RL2, double&RR2, vector<double>& ZL1, vector<double>& ZR1, vector<double>& ZL2,vector<double>& ZR2, vector<double>& dkL, vector<double>& dkR){

        // Marginal total of rows

        for(auto y: YPL1){
            RL1+=y;
        }
        for(auto y: YPR1){
            RR1+=y;
        }
        for(auto y: YPL2){
            RL2+=y;
        }
        for(auto y: YPR2){
            RR2+=y;
        }


        if(doPL){
            // Cummulative frequency of rows for the SNPs at the left
            size_t Ll = YPL1.size()-1;
            for(size_t k=0; k<(Ll); ++k){
                for(size_t i = 0; i<=k;++i){
                    ZL1[k]+= YPL1[i];
                    ZL2[k]+= YPL2[i];

                }
                dkL[k]= (k+1.0)/(Ll+1);
            } // close for k
        } //close if
        if(doPR){
            // Cummulative frequency of rows for the SNPs at the rifht
            size_t Lr = YPR1.size()-1;
            for(size_t k=0; k<(Lr); ++k){
                for(size_t i = 0; i<=k;++i){
                    ZR1[k]+= YPR1[i];
                    ZR2[k]+= YPR2[i];
                }
                dkR[k]= (k+1.0)/(Lr+1);
            }// close for k
        } // close if
}

// Cumulative chi square test (Hirotsu 1993)
template <typename TR=double, typename RT= int, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeTau(const int& N,const int& n,const double& R, const vector<double>& Z,  const vector<double>& dk, vector<double>&chiK){

    //
    double tau=0;
    size_t L = Z.size()+1;

    double R2 = N-R;
    vector<double>Z2(Z); //z2k=kn-z1k

    int c=1;
    for(auto& z2i:Z2){
        z2i=c*n-z2i;
        ++c;

    }

//cout<<N<<"\t"<<n<<"\t"<< L<<"\t"<<R<<endl;

    for(size_t k=0; k<(L-1); ++k){

        double t1 = Z[k]-R*dk[k];
        double d1 = R*dk[k];
        double t2 = Z2[k]-R2*dk[k];
        double d2 = R2*dk[k];
        double cdk = 1-dk[k];
        double t3 = (R-Z[k])-(R*cdk);
        double d3 =R-d1;
        double t4 =(R2-Z2[k])-(R2*cdk);
        double d4 = R2-d2;

        chiK[k] = (t1*t1)/d1 + (t2*t2)/d2 +(t3*t3)/d3 +(t4*t4)/d4;

        tau+= chiK[k];
    } // close for k

    return tau;

}

// Degrees of freedom for the cumulative chi square test (Hirotsu 1993)
template <typename TR=double, typename RT= int, typename Tcolcol = deque<deque<bool>>,typename Tcol2 = vector<bool> >
inline TR computeDF(const int& n,const int& L){

    //

    int N=n*L;
    vector<double>lambda(L-1,0);

    for(int k=1; k<L; ++k){
        double c=k*n;
        lambda[k-1]= c/(N-c);

    } // close for k

    double acl=0;
    for(int k=0; k<(L-2); ++k){
        double ac=0;
        for(int i=0; i<=k; ++i){
            ac+= lambda[i];
        }
        ac/=lambda[k+1];
        acl+=ac;
    }
    double d= 1 +(2*acl)/(L-1);
    double v=(L-1)/d;

    return v;

}


// compute a hacKL for a given window: haplotypes as defined by its hac value
template <typename T=double, typename Tcol=deque<bool>,  typename Tcolofcol = deque<Tcol> >
inline T hacKL(const Tcolofcol& data, const vector<vector<bool>>& references, const int numpops, const int samplesize, const int iniw, const int endw){
    // Allele (hac haplotype) distribution at both populations
    vector<vector<int>> HACdistrib(numpops, vector<int>(samplesize,0));

/// Compute hacLK distribution

    for (auto pop=0; pop<numpops; ++pop){

// Scan all over the sample to compute hac values

        for(int smpl=0; smpl<samplesize; ++smpl){

            HACdistrib[pop][smpl] = computeHAC_incl<double>(data[smpl], references[pop],iniw,endw);

        } // close for smpl

    } // close for numpops



    // Identify alleles
    vector<int> alleles;
    unordered_set<int> s; // perform operation like unique to get only non-repeated elements
    for (int i : HACdistrib[0])
        s.insert(i);
    for (int i : HACdistrib[1])
        s.insert(i);

    alleles.assign( s.begin(), s.end() );
    sort( alleles.begin(), alleles.end() );
    // number of alleles
    double nA = alleles.size();

    // compute the frequencies
    vector<double> f1(nA),f2(nA),meanf(nA);
    auto ca=0;
    for (int a: alleles){
        f1[ca] = count(HACdistrib[0].begin(), HACdistrib[0].end(), a) / nA;
        f2[ca] = count(HACdistrib[1].begin(), HACdistrib[1].end(), a) / nA;
        ++ca;
    }


    for(auto a=0; a < nA; ++a)
        meanf[a]= (f1[a]+f2[a])/2;

    T hacLK(0);

    for(auto a=0; a < nA; ++a){
        hacLK+= (f1[a]-meanf[a])*(f1[a]-meanf[a]);
        hacLK+= (f2[a]-meanf[a])*(f2[a]-meanf[a]);
    }

    return hacLK;
}

/// May 2024 the condition for block changed from Dprime or r² to Dprime AND r²
/// April 2024: This version receives precalculated D and r² vectors
/**
Get a list of haplotype blocks going each block from [iniw to endw). First element belongs to the block, but not the last.
An haplotype block is defined as a sequence of SNPs with significant pairwise r² (1,2)(2,3)...(x,x+1) or with D'>DBL (default DBL=0.5).
The minimum block size should be minwin if less the sequence of snps is ignored and the next block will begin in the last snp assayed.
I.e if the block x...x+n is less than minblock then next block begins in x+n.
Once we have a block x... x+n the center would be the integer ((x+n)-x)/2= n/2 (approximately because we dont matter if n is even or odd)
Return a vector of blocks where each block is a structure with iniw, endw.
From here it would be easy to set a candidate more or less at the center = (endw-iniw)/2
The computation is done for a given population

Enero 24: compute the desequilibrium for gamete with the same alleles as in the reference
**/
template <typename T = Block, typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol> >
inline void computeBlocks(const double& SL, const vector<double>& Dprime,const vector<double>& r2, const Tcolofcol& data, const int& minblock, const double & Dthreshold, const vector<int>& samplesizes, const int& pop, const vector<double>& freqspop, vector<T>& blocks){

    auto L = freqspop.size()-1;
    if(L<1){

        printval<int>(L,"computeBlocksError.log");
        exit(1);
    }

    int localblock= minblock;

    if((int)L<minblock){

        printval<int>(L,"Genome size " +NtoS(L) +" too short for block of size " + NtoS(minblock)+" so block is halved.log");
        //exit(1);
        localblock=minblock/2;
        if((int)L<localblock){
            printval<int>(L,"Genome size " +NtoS(L) +" too short for reduced block of size " + NtoS(localblock)+". Exiting the program.log");
            exit(localblock);
        }
    }

    if(!blocks.empty()){ // clean blocks if necessary

        vector<T>().swap(blocks); // why this is better than .clear(): https://www.techiedelight.com/delete-vector-free-memory-cpp/
    }

    int sampleshift = 0;
    for(int p=0; p<pop; ++p)
        sampleshift += samplesizes[p]; // for accesing the specific population in the metapopulation data


    /// Compute haplotype blocks

    bool block=false;
    double pval =1;
    size_t ini=0,bL =0; //blocklength


    for(size_t l=0; l < L; ++l){

        if(!block){ // set the beginning of the block
            ini=l;
            bL=1;
            block=true;
        }


        pval = getChiPval(r2[l]*samplesizes[pop],1);

        block = (pval<=SL) && (Dprime[l]>=Dthreshold);


        if(!block){
            if((int)bL<localblock){
                // next block in the next for iteration begins in the last assayed snp because we are comparing l with l+1
                continue;
            }
            else{ // add last block
                blocks.push_back(Block(pop,ini,l+1)); // remember block is defined as [ini,end) so we include current l becuase the prevous linkage Dl-1,l was true so we want to include
            }
        }
        else{
            ++bL; // increase block length
            if(l==L-1){ // the block reached the end recall L= size-1 so L-1=size-2
                if(l+2-ini >=(size_t)localblock)
                    blocks.push_back(Block(pop,ini,l+2)); //l+2=genomesize=L+1 the last position is l+1=L (because L=size-1)
            }
        }

    } //close for

}



/// May 2024 the condition for block changed from Dprime or r² to Dprime AND r²
/// ***** This April 24 version receives as argument two vectors with D' and r² values for every pair i,i+1
/// This is the version in use
/** Nov 21: Define window sizes from linkage disequilibrium blocks:
 For a given population. Compute D' and r²,
 the initial left position for the window [iniw] and the end at the right (endw) sides of a given candidate position depend of the length
 of the block with significant r² or D'>0.5.
 So that the window is [iniw,endw)
 Regarding D' and r² when the minor allele frequencies at two loci differ, there is nonlinearity in the dependence between the two loci
 (Smith 2020)
 D = x11 - p1q1 = freq gamete 11 - freqalleleLoc1*freqalleleLoc2
 Correlation between pairs of loci ij, r = D/sqrt(p_i(1-p_i)p_j(1-p_j))
 D' = D/Dmax and r

**/
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol> >
inline void computeWinSizesSeq(const double& SL, const int& L, const vector<double>& Dprime,const vector<double>& r2, const int& minblock, const double & Dthreshold, const vector<int>& samplesizes, const int& pop, const int& candidate, int& iniw, int& endw, int& conthits,double& aveDprime, double& aver2){

    if(L<2){

        printval<int>(L,"ShortGenomeComputeWinSizesSeq.log");
        exit(1);
    }

    assert(L-1 == (int)Dprime.size());
    assert(L-1== (int)r2.size());

    int sampleshift = 0;
    for(int p=0; p<pop; ++p)
        sampleshift += samplesizes[p];

    /// minimum windows sizes at left and right

    /// Define haplotype block
    /**
    When the minor allele frequencies at two loci differ, there is nonlinearity in the dependence between the two loci.
    Non linearity means that r² < D'² (Smith 2020)

    1.- Let x the index of the candidate
        left= x--1
        Compute r²(left,x) and check
        if nr²>= Chi²_sl or |D'|> Dthreshold then
            decrease --left and repeat until the condition doesnt hold.
        else end

        Similarly do with wsr.

    */
    conthits=0; /// May 2024
    /// block at the left
    if(candidate> 0){
        bool includelast=false;
        int left = candidate-1;
        double pval =1;

        do{

            pval = getChiPval(r2[left]*samplesizes[pop],1);

            includelast= (pval<=SL) && (Dprime[left]>=Dthreshold);


            if(includelast){ /// May 2024

                ++conthits; /// May 2024
                --left;
            }

        }while(includelast && left>0); /// May 2024

        if(includelast) iniw=0; // this can only be true under left=0 in this case we set the window beginning in 0
        else{


            iniw=min(left,candidate-minblock/2);
            iniw=max(0,iniw); // could happen if !includelast and genome length < minblock
        }


    } // close if candidate >0
    else{
        iniw=0; // here could be a conflict because iniw = candidate. The point is that there will be no computation at the lef side of the candidate
    }

    /// block at the right
    int Lminus1=L-1;
    if(candidate< Lminus1){ // recall L = genomesize
        bool includelast=false;
        int right = candidate+1; // L-1 is the last available position for the window
        double pval =1;

        do{

            pval = getChiPval(r2[right-1]*samplesizes[pop],1);
            includelast=(pval<=SL) && (Dprime[right-1]>=Dthreshold);
            if(includelast){
                ++right;
                ++conthits; /// May 2024

            }


        }while(includelast && right<L);

        if(includelast) endw=L;
        else{
            endw= max(right,candidate+minblock/2);
            endw=min(endw,L); // could happen if not includelast
        }



    } // close if candidate <Lminus
    else{
         endw = L;   //
    }


    ///May 2024
//    aveDprime=0;
//    aver2=0;
//    for(int i=iniw; i<endw-1; ++i){ // check mean values in the block
//        aveDprime+=Dprime[i];
//        aver2+=r2[i];
//
//    }
//    int denom=endw-1-iniw;
//    aveDprime/=denom;
//    aver2/=denom;

    /// April 2026
    if (iniw < 0 || endw > L || iniw >= endw) {
        cerr << "Bad block in computeWinSizesSeq: "
             << "candidate=" << candidate
             << " iniw=" << iniw
             << " endw=" << endw
             << " L=" << L << endl;
    }
    aveDprime = 0;
    aver2 = 0;

    int denom = endw - 1 - iniw;
    if (denom <= 0) {
        aveDprime = 0;
        aver2 = 0;
        return;
    }

    for (int i = iniw; i < endw - 1; ++i) {
        aveDprime += Dprime[i];
        aver2 += r2[i];
    }

    aveDprime /= denom;
    aver2 /= denom;
}


/// April 24:Compute D' and r² for every sequential SNP pair in a range of length L (D'(0,1), D'(1,2),...,D'(L-1,L)
/// As in the January 24 version the disequilibrium D and correlation r are computed for the alleles in the reference configuration
/// The value D'(x,x) = 1 and the number of LD values are L-1 from (0,1) to (L-2,L-1)
/// The computation is done for the given population pop
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol> >
inline void computeSeqDrsq(const vector<bool>& referencepop, const double& SL, const Tcolofcol& data, const vector<int>& samplesizes, const int& pop, const vector<double>& freqspop, vector<double>& Dprime, vector<double>& r2)
{
    auto L = freqspop.size();
    if(L<2){

        printval<int>(L,"computeSeqDrsq_L.err");
        exit(1);
    }


    if(!Dprime.empty()) vector<double>().swap(Dprime);
    Dprime.resize(L-1,0);
    if(!r2.empty()) vector<double>().swap(r2);
    r2.resize(L-1,0);

    int sampleshift = 0;

    for(int p=0; p<pop; ++p)
        sampleshift += samplesizes[p];


    double p1,q1,p2,q2,fA,fB;
    double AB=0, D=0,Dmax=0;

    for(size_t l = 0;l<L-1;++l){

        AB=0;
        for(auto s=0; s<samplesizes[pop]; ++s){
            if((data[s +sampleshift][l]==referencepop[l])&&(data[s +sampleshift][l+1]==referencepop[l+1])){
                ++AB; // count AB haplotypes where uppercase represenst the major allele in the reference population
            }
        }
        AB/=samplesizes[pop];
        p1 = freqspop[l]; //This is the frequency of the allele 1 (in the case of the ms format could not be the minor in the reference pop!)
        p2 = 1.0-p1;
        q1 = freqspop[l+1];
        q2= 1.0-q1;
        fA=referencepop[l]?p1:p2; /// if the reference in position l is 1 (1 is the major allele) then fA=p1 else p2. So this is correct whatever the allele in the reference (0 or 1) the comparison is made correctly
        fB=referencepop[l+1]?q1:q2;
        D= AB-fA*fB; // compute D = f(AB)-f(A)*f(B)
        // Dmax
        if(D>=0)
            Dmax=min(p1*q2,p2*q1);
        else
            Dmax=min(p1*q1,p2*q2);
        double pprod=p1*p2*q1*q2;
        if(pprod>0){ // only if alleles are segregating
            r2[l]= (D*D)/(pprod);
            Dprime[l]= D/ Dmax;
        }
        else{ //April 24 si no segregan NO consideramos inclusión en el bloque
            r2[l]=0;
            Dprime[l]=0;

        }

    } // close for l
}




/// Enero 24: Compute D' and r² for a candidate SNP x with every position in a range of length L (L values apart of x) so D'(0,x), D'(1,x),...,D'(x-1,x), 1, D'(x, x+1), D'(x,x+2)...D'(x,x+L-x)
/// ***NEW: The disequilibrium D and correlation r are computed for the alleles in the reference configuration
/// The value D'(x,x) = 1 and the number of LD values are L
/// The computation is done for the given population pop
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol> >
inline void computeDrsq(const vector<bool>& referencepop, const double& SL, const Tcolofcol& data, const vector<int>& samplesizes, const int& pop, const int& candidate, const vector<double>& freqspop, vector<double>& Dprime, vector<double>& r2)
{
    auto L = freqspop.size();
    if(L<2){

        printval<int>(L,"computeDrsq_L.err");
        exit(1);
    }

    if(candidate<0 || candidate>=(int)L){
        printval<int>(candidate,"computeDrsq_outofrange.err");
        exit(1);
    }

    if(!Dprime.empty()) vector<double>().swap(Dprime);
    Dprime.resize(L,0);
    if(!r2.empty()) vector<double>().swap(r2);
    r2.resize(L,0);

    Dprime[candidate]=1;
    r2[candidate]=1;

    int sampleshift = 0;

    for(int p=0; p<pop; ++p)
        sampleshift += samplesizes[p];


    double p1,q1,p2,q2,fA,fB;
    double AB=0, D=0,Dmax=0;

    for(size_t l = 0;l<L;++l){
        if((int)l!=candidate){
            AB=0;
            for(auto s=0; s<samplesizes[pop]; ++s){
                if((data[s +sampleshift][l]==referencepop[l])&&(data[s +sampleshift][candidate]==referencepop[candidate])){
                    ++AB; // count 00 haplotypes
                }
            }
            AB/=samplesizes[pop];
            p1 = freqspop[l];
            p2 = 1.0-p1;
            q1 = freqspop[candidate];
            q2= 1.0-q1;
            fA=referencepop[l]?p1:p2;
            fB=referencepop[candidate]?q1:q2;
            D= AB-fA*fB; // compute D = f(AB)-f(A)*f(B)
            // Dmax
            if(D>=0)
                Dmax=min(p1*q2,p2*q1);
            else
                Dmax=min(p1*q1,p2*q2);
            double pprod=p1*p2*q1*q2;
            if(pprod>0){ // only if alleles are segregating
                r2[l]= (D*D)/(pprod);
                Dprime[l]= D/ Dmax;
            }
            else{ //April 24 si no segregan NO consideramos inclusión en el bloque
                r2[l]=0;
                Dprime[l]=0;

            }

        } // close if l != candidate
    } // close for
}


// D = x11 - p1q1 = freq gamete 11 - freqalleleLoc1*freqalleleLoc2
//Compute correlation between adjacent pairs of loci, rho = D/sqrt(p(1-p)q(1-q))
// Return the vector of pairwise correlations
// Population level
/// Compute r' = D/Dmax. Because when the allele freqs differs between loci the range of r is smaller (Hedrick p. 532-533)

template <typename Tcol=deque<bool>,  typename Tcolofcol = deque<Tcol> >
inline void computeDpop(const Tcolofcol& data,const vector<int>& samplesizes, const int& pop, const vector<double>& freqspop, vector<double>& LDvalues){

    auto L = freqspop.size()-1;
    if(L<1){

        printval<int>(L,"computeDpop_Error.log");
        exit(1);
    }

    int sampleshift = 0;
    for(int p=0; p<pop; ++p)
        sampleshift += samplesizes[p];

    double p1,q1,p2,q2,Dmax;
    if(!LDvalues.empty()) vector<double>().swap(LDvalues);
    LDvalues.resize(L,0);
    for(auto s=0; s<samplesizes[pop]; ++s)
        for(auto loc=0; loc<L; ++loc)
            if(data[s + sampleshift][loc]&data[s + sampleshift][loc+1])
                ++LDvalues[loc]; // count 11 haplotypes
    transform(LDvalues.begin(),LDvalues.end(),LDvalues.begin(), [=](double element){return(element/samplesizes[pop]);} ); // freq of 11
    for(auto loc=0; loc<L; ++loc){
    /// Follows Hedrick p. 528 notation where he uses p for first locus and q for the second
        p1 = freqspop[loc];
        p2 = 1.0-p1;
        q1 = freqspop[loc+1];
        q2= 1.0-q1;
        LDvalues[loc]-= p1*q1; // compute D = f(AB)-f(A)*f(B)
        // Dmax
        if(LDvalues[loc]>=0)
            Dmax=min(p1*q2,p2*q1);
        else
            Dmax=min(p1*q1,p2*q2);
        //LDvalues[loc]/= sqrt(freqs[loc]*(1-freqs[loc])*freqs[loc+1]*(1-freqs[loc+1]));
        LDvalues[loc]/= Dmax;

    }
}
/// Compute D' = D/Dmax between adajacent pairs of loci. Because when the allele freqs differs between loci the range of r is smaller (Hedrick p. 532-533)
// D = x11 - p1q1 = freq gamete 11 - freqalleleLoc1*freqalleleLoc2
// Return the vector of pairwise correlations
// Metapopulation level
template <typename Tcol=deque<bool>,  typename Tcolofcol = deque<Tcol> >
inline void computeDmetapop(const Tcolofcol& data,const int& samplesize, const vector<double>& freqs, vector<double>& LDvalues){

    int L = freqs.size()-1;

    if(L<1){

        printval<int>(L,"computeDmetapop_Error.log");
        exit(1);
    }
    double p1,q1,p2,q2,Dmax;
    if(!LDvalues.empty()) vector<double>().swap(LDvalues);
    LDvalues.resize(L,0);

    for(auto s=0; s<samplesize; ++s)
        for(auto loc=0; loc<L; ++loc)
            if(data[s][loc]&data[s][loc+1]) // freq of gamete 11
                ++LDvalues[loc];
    transform(LDvalues.begin(),LDvalues.end(),LDvalues.begin(), [=](const double& element){return(element/samplesize);} );
    for(auto loc=0; loc<L; ++loc){
    /// Follows Hedrick p. 528 notation where he uses p for first locus and q for the second
        p1 = freqs[loc];
        p2 = 1.0-p1;
        q1 = freqs[loc+1];
        q2= 1.0-q1;
        LDvalues[loc]-= p1*q1;
        // Dmax
        if(LDvalues[loc]>=0)
            Dmax=min(p1*q2,p2*q1);
        else
            Dmax=min(p1*q1,p2*q2);
        //LDvalues[loc]/= sqrt(freqs[loc]*(1-freqs[loc])*freqs[loc+1]*(1-freqs[loc+1]));
        LDvalues[loc]/= Dmax;
    }

}

// Given a vector of pairwise rho computes between [iniw, endw) the mean of absolute rho values
// it also computes any k-order statistics or the quartile range
template <typename Tcol=deque<bool>,  typename Tcolofcol = deque<Tcol> >
inline void computeDstats(const int& iniw, const int& endw, const vector<double>& LDvalues, const size_t korder, double& kstat, double& absmean){

    int L = max(endw-iniw,1);
    if(L>LDvalues.size()|| L<4){
        printval<int>(L,"computerhostats_Error.log");
        exit(1);
    }

    vector<double>D(LDvalues);
    absmean=0;
    for(auto d=iniw; d <endw; ++d)
        absmean+=abs(D[d]);

    absmean/=D.size();

    // range interquartile
    if(korder>=L){
        int iq = iniw+0.75*L-1;
        nth_element(D.begin()+iniw,D.begin()+iq,D.begin()+endw);

        kstat = D[iq];
        iq = iniw+0.25*L-1;
        nth_element(D.begin()+iniw,D.begin()+iq,D.begin()+endw);

        kstat-=D[iq];

    }else{ //k-order statistic
        nth_element(D.begin()+iniw,D.begin()+iniw+korder,D.begin()+endw);
        kstat = D[korder];
    }

    //printcol(D);
    //mean = accumulate(LDvalues.begin(),LDvalues.end(),0.0)/LDvalues.size();

}

// IQR: interquartile range. The seq is partially sorted as done by nth_element
template <typename Tel=ResultFst,typename Tcol=vector<ResultFst>>
inline void iqrhs(Tcol & seq, Tel& Q1, Tel& Q3, double& IQR){

    auto const pQ1 = seq.size() / 4;
    auto const pQ2 = seq.size() / 2;
    auto const pQ3 = pQ1 + pQ2;
    std::nth_element(seq.begin(), seq.begin() + pQ1, seq.end(), [](const ResultFst& fst1,const ResultFst& fst2){return fst1.fst < fst2.fst;} );
    std::nth_element(seq.begin() + pQ1 + 1, seq.begin() + pQ3, seq.end(), [](const ResultFst& fst1,const ResultFst& fst2){return fst1.fst < fst2.fst;} );

    Q1 = seq[pQ1];
    Q3 = seq[pQ3];
    IQR = Q3-Q1;
//    cout<<Q3.fst<<endl;
//    cout<<Q1.fst<<endl;
//    cout<< (Q3-Q1)<< "\t"<<(Q1-Q3)<<endl; exit(0);

}

// Compute positive outliers > Q3 + 3IQR in a collection cdata of FST values
template <typename Tel=STAT,typename Tcoll=vector<Tel>>
inline void epos(const Tcoll & cdata,const vector<double> & cpos,const vector<string>&names, const double& sigma, vector<ResultFst>& eposet){

    vector<ResultFst>seq(cdata.size());


    for(size_t s=0; s < seq.size(); ++s){
        seq[s]=(ResultFst(cpos[s],cdata[s].stat,cdata[s].pval,1,names[s]));
//        seq[s].fst = cdata[s];
//        seq[s].pval=0;
    }

    if(!eposet.empty()) vector<ResultFst>().swap(eposet);

    eposet.reserve(.25*cdata.size());

    ResultFst Q1,Q3;
    double IQR;
    if(is_sorted(seq.begin(),seq.end(),[](const ResultFst& seq1,const ResultFst& seq2){return seq1.fst < seq2.fst;} )){

        Q1 = seq[seq.size() / 4];
        Q3 = seq[seq.size() / 4 + seq.size() / 2];
        IQR = Q3-Q1;
    }
    else
        iqrhs(seq,  Q1,  Q3, IQR);

    double lim = Q3.fst+ sigma*IQR;


// seq has all elements at the left of Q3 being <Q3 and so being < lim also min IQR is 1 so q3 is <lim
    int pq3=0;

    for(size_t p = pq3+1; p < seq.size(); ++p)
        if(seq[p].fst > lim){
            eposet.push_back(seq[p]);


        }

}



/// Compute k-percentile outliersin a collection cdata of FST values
/// argument k is percentile/100 i.e. for 90 percentile k=0.9
template <typename Tel=STAT,typename Tcoll=vector<Tel>>
inline void k_percentile(const size_t& genomesize,const double& k,const Tcoll & cdata,const vector<double> & cpos,const vector<string>&names, vector<ResultFst>& eposet){

    vector<ResultFst>seq(cdata.size());


    for(size_t s=0; s < seq.size(); ++s){
        seq[s]=(ResultFst(cpos[s],cdata[s].stat,cdata[s].pval,1,names[s]));
//        seq[s].fst = cdata[s];
//        seq[s].pval=0;
    }

    size_t percentil90=genomesize*(1.0-k);

    if(!eposet.empty()) vector<ResultFst>().swap(eposet);

    eposet.reserve(percentil90);
    //select the percentil90 highest FSTs
    std::nth_element(seq.begin(), seq.begin() + percentil90, seq.end(), [](const ResultFst& seq1, const ResultFst& seq2){return seq1.fst>seq2.fst;} );

    for(size_t p = 0; p < percentil90; ++p) //store this sites in eposet

        eposet.push_back(seq[p]);



}




/// K-means algorithm: modifies the object list leaving in it only the new detect outliers
template <typename T=ResultFst, typename Tcoll= vector<T>, typename Tcollofcoll = vector<Tcoll> >
inline void kmeans_outliers(Tcoll& objects, const vector<double> modes, const double& cutoff, const bool extreme)
{
    size_t k= modes.size();
    if(k<=1) {printval<int>(k,"K equal or less than 1 in kmeans","kmeans.log");exit(1);}
    if(objects.empty()) return;
    if(k>objects.size()) { // if more initial classes than objects
        //printval<int>(k,"K equal or greater than objects in kmeans","kmeans.log");
        k= objects.size();
        if(k==1){
            //double ru= (modes[0]+modes[1]) /2.0;
            double ru=cutoff;
            if(objects[0].fst<ru){
                objects.clear();
            }
            return;
        }

    }

    const double TINYDiff=1.0e-6;
    Tcollofcoll classes(k);

//    classes[0][0] = objects[0];
//    classes[k-1][0] = (objects.back());
//
//    for(auto clss=1; clss< k-1; ++clss){
//        classes[clss].push_back(objects[clss]);
//    }

    vector<double> centroids(k);
    centroids[0] = objects[0].fst; // if the objects collection is sorted this is the minimum fst
    centroids[k-1] = objects.back().fst; // if objects is sorted this is the maximum
    for(size_t clss=1; clss< k-1; ++clss) /// I use size_t because above we forze k > 1
        centroids[clss]=objects[clss].fst;

    bool end = false;
    int clase = -1;
    double maxd= modes.back()+1; // the last value of the modes vector should be a maximum
    double dist = maxd;
/// apply the k-means algorithm
    do{
        for(auto obj: objects){ // for each object-outlier
            for(size_t ctr=0; ctr< k; ++ctr){ // compute the distance between the object and the centroids
                double d = abs(obj.fst-centroids[ctr]);
                if( d< dist){
                    dist = d;
                    clase = ctr;
                    if(dist==0) break;
                }
            } // close for centroids
            classes[clase].push_back(obj);
            dist=maxd;
        } // close for objects

        // recompute centroids
        end =true;
        for(size_t cl=0; cl<classes.size(); ++cl){
            double z=0;
            for(auto x:classes[cl]){
                z+=x.fst;
            }
            z/=classes[cl].size();
            if(abs(centroids[cl]- z)>TINYDiff) {
                centroids[cl]=z;
                end=false;
            }

        }
        if(!end){ // clean classes for the next iteration
            for(size_t cl=0; cl<classes.size(); ++cl){
                classes[cl].clear();
            }
//            classes[0][0]= objects[0];
//            classes[k-1][0]= objects.back();
//            for(auto clss=1; clss< classes.size()-1; ++clss){
//                classes[clss][0]=(objects[clss]);
//            }
        }


    }while(!end);
    for(size_t cl=1; cl<classes.size(); ++cl)
        if(!classes[cl].empty() ) // sort inside each class in ascending order
            sort(classes[cl].begin(), classes[cl].end(), [](const ResultFst& fst1,const ResultFst& fst2){return fst1.fst < fst2.fst;} );

// Clean old outliers

    objects.clear();

/// This is the extra: Test if classes are new outliers
    double /*ave,sd,z=0,*/ ru{cutoff};
// Method 1: use cutoff as an optimal threshold

    if(extreme){ // only the elements in the highest class are compared
        if(classes[k-1].empty() ) ;
        else{
            //ru= (modes[k-2]+modes[k-1]) /2.0;
//cout<<k<<" "<<classes[k-1][0].fst<<endl;
            for(size_t el=0; el<classes[k-1].size(); ++el ) // copy every element with fst equal or higher than ru
                if(classes[k-1][el].fst >= ru){
                    copy(classes[k-1].cbegin()+el,classes[k-1].cend(), back_inserter( objects));
                    break;
                }
        }
    }
    else{ // compare each pair of classes
        for(size_t cl=0; cl<(k-1); ++cl)
            for(size_t cl2=cl+1; cl2<k; ++cl2){
                if(classes[cl2].empty() ) continue;
                ru= (modes[cl]+modes[cl2]) /2.0; // middle point between each pair of class means
                for(size_t el=0; el<classes[k-1].size(); ++el )
                    if(classes[cl2][el].fst >= ru){
                        copy(classes[cl2].cbegin(),classes[cl2].cend(), back_inserter( objects));
                        break;
                    }
            }
    }

// Method 2: Compute z scores in a class using the mean and sd of the other. Select the highest class with z scores >=3

//    if(!extreme){
//        if(objects.empty()){
//            if(classes[k-2].empty() ) sd=0;
//            else{
//                avesamplevar<double, vector<ResultFst> >(classes[k-2],ave,sd);
//                sd=sqrt(sd);
//            }
//            if(sd>0){
//                z = (classes[k-1][0].fst-ave) / sd;
//                if(z>sigma)
//                    copy(classes[k-1].cbegin(),classes[k-1].cend(), back_inserter( objects));
//            }
//
//        }
//
//    }
//    else{
//        if(objects.empty()){
//            if(classes[0].empty() ) sd=0;
//            else{
//                avesamplevar<double, vector<ResultFst> >(classes[0],ave,sd);
//                sd=sqrt(sd);
//            }
//            if(sd>0)
//                for(auto cl=1; cl<classes.size(); ++cl){
//
//                    z = (classes[cl][0].fst-ave) / sd;
//                    if(z>sigma)
//                        copy(classes[cl].cbegin(),classes[cl].cend(), back_inserter( objects));
//
//
//                }
//            else{
//                if(classes[1].empty() ) sd=0;
//                else{
//                    avesamplevar<double, vector<ResultFst> >(classes[1],ave,sd);
//                    sd=sqrt(sd);
//                }
//                if(sd>0){
//                    z = (classes[0][0].fst-ave) / sd;
//                    if(abs(z)>sigma)
//                        copy(classes[1].cbegin(),classes[1].cend(), back_inserter( objects));
//                    for(auto cl=2; cl<classes.size(); ++cl){
//
//                        z = (classes[cl][0].fst-ave) / sd;
//                        if(z>=3)
//                            copy(classes[cl].cbegin(),classes[cl].cend(), back_inserter( objects));
//                    }
//                }
//                else{
//                    if(k>2){
//                        int contk=2;
//                        do{
//                            if(classes[contk].empty() ) sd=0;
//                            else{
//                                avesamplevar<double, vector<ResultFst> >(classes[contk],ave,sd);
//                                sd=sqrt(sd);
//                            }
//
//                            if(sd>0){
//                                z = (classes[contk-1][0].fst-ave) / sd;
//                                if(abs(z)>sigma)
//                                    copy(classes[contk].cbegin(),classes[contk].cend(), back_inserter( objects));
//                                for(auto cl=contk+1; cl<k; ++cl){
//
//                                    z = (classes[cl][0].fst-ave) / sd;
//                                    if(z>sigma)
//                                        copy(classes[cl].cbegin(),classes[cl].cend(), back_inserter( objects));
//                                }
//                            }
//                            ++contk;
//                        }while( (sd<=0.0) && (contk<k) );
//
//                    }
//                } // close else
//
//            } // close else
//
//
//        } // close if objects is empty
//    } // close if no extreme


}

// compute the parameter n for the LBE location based estimator from (7) in Dalmasso et al 2005
template <typename T=float,typename Tint = int>
inline void calc_n_LBE(const T& var, const Tint& S, Tint& n){

	n=1; // initial guess, this is the minimum
	Tint nt{n};
	do{
        ++nt;
        T fnt{factorialNRC<T>(nt)};
        if( (factorialNRC<T>(2*nt) / (fnt*fnt) -1) <= (var*S) )
            n=nt;

	}while(n==nt);

}


// LBE location based estimator from Dalmasso et al 2005
template <typename T=double,typename Tcol=vector<T>, typename Tint = int>
inline T LBE(const Tcol& pvals,const Tint& n){
	T factorial{0};
	Tint pnumber=pvals.size();
	T p_i{0.0};
	for(auto i = 0; i< pnumber; i++){
		p_i+= pow(-log(1.0-pvals[i]), (T)n);

	}
//cout<<pvals.size()<<endl;
//cout<<pvals.back()<<endl; exit(9);
//cout<<pnumber<<endl;
//cout<<p_i<<endl;
	p_i/=pnumber;
//cout<<p_i<<endl;
	factorial = factorialNRC<T>(n);

	p_i/= factorial;
//cout<<p_i<<endl;

	return p_i;
}


// Method of Meinshausen and Rice 2006 to estimate the proportionof false null hypothesis: standard deviation-proportional bounding
// pvalues are copied and sorted if necessary
template <typename T=double,typename Tcol=vector<T>, typename Tint = int>
T sdpbound(const T& sdpbprecision, const Tcol& pvals0, const T& sl){
	T lambda, suplambda=0.0;
	T Fn;

	Tcol pvals{pvals0};

	T lowerb = 1.0/pvals.size();
	T upperb = 1.0 - lowerb;
	T delta, an, bn, betan,E;
	T indicator, n =pvals.size();

	if(!is_sorted(pvals.begin(),pvals.end()))
        sort(pvals.begin(),pvals.end());

// formulae (4) compute the empirical distribution of p-values

	for(T t=lowerb; t<upperb; t+=sdpbprecision){

		indicator=0;
		for(Tint p=0; p<n;++p)
			if(pvals[p]<=t) indicator++;
			else break; // because p-values are sorted

		Fn =indicator/n;

// compute bounding function delta 2.2
		delta= sqrt(t*(1-t));

// compute quantities on Lemma 1: an, bn and the Gumbel distribution E
		an= sqrt(2*n*log(log(n)));
		bn = 2* log(log(n)) + 0.5*log(log(log(n))) - .5*log(4*dPI);
		E = exp(-exp(sl-1));

// formulae (13) compute the sd-p b sequence
		betan = (1/E +bn) / an;

// Estimate the proportion of false null hypothesis as in (5) using the recommended interval 1/n, 1 - 1/n instead of (0,1) n being the number of tests
		lambda = (Fn - t - betan*delta) / (1-t);	// this is (13)
		if(lambda>suplambda) suplambda = lambda;

	} // close for t


	return suplambda;
}


/*
Folded Site Frequency Spectrum (FSFS):
For a deep of n=last-first sequences the frequency of each segregating site can go from 1 to n.
If the allelic state is unknown (folded SFS) we can distinguish the following categories:
1 or n-1, 2 or n-2.. n/2 or n - n/2. The observed distribution of this classes is the FSFS i.e. the
number of sites belonging to class 1, 2.. n/2. In general Si the number of sites in class i.
Result is an fsfs vector with dimension n/2 and the number of sites at each class stored. It also returns the number of segregating sites S
*/
template <typename T,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void FSFS(const Tcolofcol& sequences, size_t first, size_t last, size_t ini, size_t end, vector<double>& fsfs, size_t& S){

    auto exit_val=0;
    auto n = last - first; // number of sequences

    if(n<=1 || first<0){
        exit_val++;
        printval<int>(n,"FSFS_Error.log");
        exit(exit_val);
    }

    if( (sequences.size() < n) ){
        first = 0;
        last = sequences.size();
        n=sequences.size();
    }

    // Check the snp range [ini,end)
    if(sequences[first].size() < end){
        end =sequences[first].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exit_val+=2;
        printval<size_t>(ini,"FSFS_Error.log");
        exit(exit_val);
    }

    size_t WS = end-ini;

// no es necesario hacer clear porque si el nuevo tamaño acorta resize elimina lo restante, sino lo deja pero se sobrescribe despues
    fsfs.resize(n/2,0);
    vector<double> sites(WS,0.0);

    for(auto i=first; i< (last); ++i){

        for(auto s =0; s<WS; ++s)
            if(sequences[i][s+ini])
                ++sites[s]; // frequency of each segregating site s

    }
// Compute the number of segregating sites
    S=0;
    for_each(sites.begin(),sites.end(), [&](const double& site){if(site>0.0) ++S;});

    //fsfs[0] stores number of sites having value 1 or n-1 as frequency, fsfs[1] stores frequencies 2 and n-2, fsfs[i] freqs i+1 n-(i+1)

    for(auto i = 1; i<= fsfs.size(); ++i){
        fsfs[i-1] = count(sites.begin(),sites.end(),(i)) + count(sites.begin(),sites.end(),(n-i));
    }

    if(n%2==0) // n is even so the Kronecker delta is 1 just for the last element n/2 divide by two because  has been summed twice
        fsfs.back()/=2.0;

}

/// General unbiased estimator for the folded spectrum (formulae 14 in Achaz 2009)

template <typename T=double,typename Tcol=vector<T>>
inline T General_FS_Theta(const Tcol& fsfs, const Tcol& w, const size_t& n){

    auto exit_val=0;
    auto halfn = fsfs.size(); // number of sequences

    if(halfn<=1){
        exit_val++;
        printval<int>(halfn,"General_FS_Theta.log");
        exit(exit_val);
    }

    assert(fsfs.size() == n/2);
    assert(fsfs.size() == w.size());

//compute phi
    size_t cont=1;
    vector<T> phi(halfn,n);
    for_each(phi.begin(), phi.end(), [&](T& phi_i){phi_i/= (cont*(n-cont));++cont;});
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide by two because n/2 has been summed twice
        phi.back()/=2.0;

// Compute omega i.e. the difference between weights
    T sumw = accumulate(w.begin(),w.end(),0.0);

// This operation omega_i/ phi_i appears a lot of times so we may precompute it in a vector
    cont=0;
    T theta=0;

    for_each(w.begin(), w.end(), [&](const T& w_i){theta+=( (w_i/phi[cont]) * fsfs[cont] );++cont;});

    theta/= sumw;

    return theta;
}


// Compute ThetaPi using the general model (14) in Achaz09.
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol>>
inline double gfThetaPi(const Tcolofcol& sequences, size_t first, size_t last, const size_t ini, const size_t end, size_t& L){
    double gTpi=0;
    auto exit_val=0;
    auto n = last - first; // number of sequences

    if(n<=1 || first<0){
        exit_val++;
        printval<int>(n,"gfThetaPi_Error.log");
        exit(exit_val);
    }

    if( (sequences.size() < n) ){
        first = 0;
        last = sequences.size();
        n=sequences.size();

    }

    //if(ini>=end) return 0;

    L = end-ini;
    if(L > sequences[first].size() ){
        exit_val+=2;
        printval<size_t>(exit_val,"gfThetaPi_Error.log");
        exit(exit_val);
    }

// Compute folded spectrum
    vector<double> fsfs;
    FSFS<bool>(sequences, first, last, ini, end, fsfs,0);

/// The watterson theta corresponds to i(n-i)phi = n see Table 1 in Achaz09
    vector<double> w(fsfs.size(),n);

    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide last element by two because n/2 has been summed twice
        w.back()/=2.0;

    gTpi = General_FS_Theta(fsfs, w, n);

    return gTpi;
}

// Compute ThetaS using the general model (14) in Achaz09
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol>>
inline double gfThetaS(const Tcolofcol& sequences, size_t first, size_t last, const size_t ini, const size_t end, size_t& L){
    double gTS=0;
    auto exit_val=0;
    auto n = last - first; // number of sequences

    if(n<=1 || first<0){
        exit_val++;
        printval<int>(n,"gfThetaS_Error.log");
        exit(exit_val);
    }

    if( (sequences.size() < n) ){
        first = 0;
        last = sequences.size();
        n=sequences.size();

    }

    //if(ini>=end) return 0;

    L = end-ini;
    if(L > sequences[first].size() ){
        exit_val+=2;
        printval<size_t>(exit_val,"gfThetaS_Error.log");
        exit(exit_val);
    }

// Compute folded spectrum
    vector<double> fsfs;
    size_t S;
    FSFS<bool>(sequences, first, last, ini, end, fsfs,S);

/// The thetaS corresponds to phi see Table 1 in Achaz09
    vector<double> phi(fsfs.size(),n);
    size_t cont = 1;
    for_each(phi.begin(), phi.end(), [&](double& phi_i){phi_i/= (cont*(n-cont));++cont;});
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide last element by two because n/2 has been summed twice
        phi.back()/=2.0;

    gTS = General_FS_Theta(fsfs, phi, n);

    return gTS;
}

/// Folded spectrum neutrality tests (after Achaz 2009)

// Generic folded spectrum neutrality test following Achaz 2009 (eq 16)
// Given a pre-computed folded espectra fsfs and for a couple of weights vectors w1 and w2  performs a linear combination of the folded spectrum to perform a neutrality test T*
//The variable zeros indicate how many weights should be ignored. That means if zeros=0 all the spectra is considered
//if zeros=1 singletons are discarded and thus omega[0]=0 if zeros=2 omega[0]=omega[1]=0 and so on
template <typename T=double,typename Tcol=vector<T>>
inline T General_FS_Test(const Tcol& fsfs, const Tcol& weights1, const Tcol& weights2, const size_t& n, const size_t& S, const size_t zeros){

    auto exit_val=0;
    auto halfn = fsfs.size(); // number of sequences

    if(halfn<=1){
        exit_val++;
        printval<int>(halfn,"General_FS_Test.log");
        exit(exit_val);
    }

    assert(fsfs.size() == n/2);
    assert(fsfs.size() == weights1.size());
    assert(weights2.size() == weights1.size());

    vector<T> w1(weights1);
    vector<T> w2(weights2);

    vector<T> phi(halfn,n);



    for(auto i=0; i<halfn; ++i){
        phi[i]/= ((i+1)*(n-i-1));
    }

    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide by two because n/2 has been summed twice
        phi.back()/=2.0;

// Compute omega i.e. the difference between weights
    T sumw1 = accumulate(w1.begin(),w1.end(),0.0);
    T sumw2 = accumulate(w2.begin(),w2.end(),0.0);

    for_each(w1.begin(), w1.end(), [=](T& wi){wi/=sumw1;} );
    for_each(w2.begin(), w2.end(), [=](T& wi){wi/=sumw2;} );

    vector<T> omega(w1);
    size_t cont=0;
    for_each(omega.begin(), omega.end(), [&](T& omega_i){if(cont+1>zeros) omega_i-=w2[cont]; else omega_i=0;++cont;});

// This operation omega_i/ phi_i appears a lot of times so we may precompute it in a vector
    cont=0;
    vector<T> omegainvphi(omega);
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](T& omegainvphi_i){omegainvphi_i/=phi[cont];++cont;});

/// Computations for the denominator in (16)
    cont=0;
    T alphan=0;
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](const T& omegainvphi_i){alphan+=(omegainvphi_i*omega[cont]);++cont;});

// For computing betan* in (16) we need some auxiliary functions defined as lambdas (i will use notation like L_)
// From Fu95
    auto L_An = [](const size_t n) {T an =0.0;for(size_t i=1;i<n; ++i){an+= 1.0/i;} return an; }; // use as L_An(n)
    auto L_Bn = [=](const size_t n, const size_t i) {T bn = ( (2.0*n)/((n-i+1.0)*(n-i)))*(L_An(n+1) - L_An(i)) - 2.0/(n-i) ; return bn; };
    auto L_Sigmaii = [=](const size_t n, const size_t i){
        size_t nmed=n/2;
        if(i<nmed) return L_Bn(n,i+1);
        else if (i==nmed) return (2*(L_An(n) - L_An(i)) / (n-i) - 1.0/(i*i));
            else return (L_Bn(n,i) - 1.0/(i*i));
    };
    auto L_Sigmaij = [=](const size_t n, const size_t ci, const size_t cj){ // defined in Fu for i>j so we swap values i,j if necessary
        size_t imasj=ci+cj;
        size_t i=ci;
        size_t j = cj;
        // condition i > j
        if(i<j)
            swap(i,j);
        if(imasj<n) return (L_Bn(n,i+1) - L_Bn(n,i))/2.0;
        else
            if (imasj==n){
                T Ai = (L_An(n) - L_An(i)) / (n-i);
                T Aj = (L_An(n) - L_An(j)) / (n-j);
                T B = (L_Bn(n,i) + L_Bn(n,j+1)) / 2.0;

                return (Ai + Aj - B - 1.0/(i*j));
            }
            else
                return( (L_Bn(n,j) - L_Bn(n,j+1) )/2.0 - 1.0/(i*j) );

    };

    auto L_Rhoii = [=](const size_t n, const size_t i){
        T R = L_Sigmaii(n,i) + L_Sigmaii(n,n-i) + 2*L_Sigmaij(n,i,n-i) ;
        if(i==n-i) R/=4.0;
        return R;
    };
    auto L_Rhoij = [=](const size_t n, const size_t i, const size_t j){
        T R = L_Sigmaij(n,i,j) + L_Sigmaij(n,i,n-j) + L_Sigmaij(n,n-i,j) + L_Sigmaij(n,n-i,n-j) ;
        if((i==n-i) && (j==n-j) ) R/=4.0;
        else if(i==n-i) R/=2.0;
            else if (j==n-j) R/=2.0;

        return R;
    };

    cont=1;
    T betan=0.0;
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](const T& omegainvphi_i){betan+=(omegainvphi_i*omegainvphi_i*L_Rhoii(n,cont));++cont;});
    T bpart=0;
    for(auto i =0; i < (omegainvphi.size()-1); ++i)
        for(auto j =i+1; j < omegainvphi.size(); ++j){
            bpart+=omegainvphi[i]*omegainvphi[j]*L_Rhoij(n,i+1,j+1);
        }
//cout<<L_Sigmaii(n,1)<<" + " <<L_Sigmaii(n,n-1)<<" + "<<2*L_Sigmaij(n,1,n-1)<< " = "<<L_Rhoii(n,1) <<endl;
//cout<<betan<<" "<<L_Sigmaij(n,n-1,2)<<endl;
    betan += 2.0*bpart;

    T a1=1.0,a2=1.0; // (3) and (4) in Tajima 89
    for(auto i=2; i< n; ++i){
        a1+= 1.0/i;
        a2+= 1.0/(i*i);
    }
//From Tajima 89 p. 588
    T etheta = S / a1;
    T etheta2 = (S*(S-1)) / (a1*a1 + a2);

//cout<<etheta<<" "<<etheta2<<" "<<alphan<<" "<<betan<<endl;

/// The denominator in (16) is
    T test_denom = sqrt(alphan* etheta + betan*etheta2);


/// The Test is
    cont =0;
    T Tomtest = 0.0;
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](const T& omegainvphi_i){Tomtest+=(omegainvphi_i*fsfs[cont]);++cont;});

    Tomtest/=test_denom;

    return Tomtest;
}

// This second version of the Achaz09 test receives directly the omega vector instead of weights w1 and w2. This prove to be useful for computing optimal tests as in Ferreti et al 20101
template <typename T=double,typename Tcol=vector<T>>
inline T General_FS_Test(const Tcol& fsfs, const Tcol& gomega,const size_t& n, const size_t& S, const size_t zeros){

    auto exit_val=0;
    auto halfn = fsfs.size(); // number of sequences

    if(halfn<=1){
        exit_val++;
        printval<int>(halfn,"General_FS_Test.log");
        exit(exit_val);
    }

    assert(fsfs.size() == n/2);
    assert(fsfs.size() == gomega.size());

    vector<T> omega(gomega);
    vector<T> phi(halfn,n);

    for(auto i=0; i<halfn; ++i){
        phi[i]/= ((i+1)*(n-i-1));
    }
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide by two because n/2 has been summed twice
        phi.back()/=2.0;


    size_t cont=0; // correct of errors eliminating singletons, doubletons as indicated in the flag zeros
    for_each(omega.begin(), omega.end(), [&](T& omega_i){if(cont+1<=zeros) omega_i=0;++cont;});

// This operation omega_i/ phi_i appears a lot of times so we may precompute it in a vector
    cont=0;
    vector<T> omegainvphi(omega);
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](T& omegainvphi_i){omegainvphi_i/=phi[cont];++cont;});

/// Computations for the denominator in (16)
    cont=0;
    T alphan=0;
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](const T& omegainvphi_i){alphan+=(omegainvphi_i*omega[cont]);++cont;});

// For computing betan* in (16) we need some auxiliary functions defined as lambdas (i will use notation like L_)
// From Fu95
    auto L_An = [](const size_t n) {T an =0.0;for(size_t i=1;i<n; ++i){an+= 1.0/i;} return an; }; // use as L_An(n)
    auto L_Bn = [=](const size_t n, const size_t i) {T bn = ( (2.0*n)/((n-i+1.0)*(n-i)))*(L_An(n+1) - L_An(i)) - 2.0/(n-i) ; return bn; };
    auto L_Sigmaii = [=](const size_t n, const size_t i){
        size_t nmed=n/2;
        if(i<nmed) return L_Bn(n,i+1);
        else if (i==nmed) return (2*(L_An(n) - L_An(i)) / (n-i) - 1.0/(i*i));
            else return (L_Bn(n,i) - 1.0/(i*i));
    };
    auto L_Sigmaij = [=](const size_t n, const size_t ci, const size_t cj){ // defined in Fu for i>j so we swap values i,j if necessary
        size_t imasj=ci+cj;
        size_t i=ci;
        size_t j = cj;
        // condition i > j
        if(i<j)
            swap(i,j);
        if(imasj<n) return (L_Bn(n,i+1) - L_Bn(n,i))/2.0;
        else
            if (imasj==n){
                T Ai = (L_An(n) - L_An(i)) / (n-i);
                T Aj = (L_An(n) - L_An(j)) / (n-j);
                T B = (L_Bn(n,i) + L_Bn(n,j+1)) / 2.0;

                return (Ai + Aj - B - 1.0/(i*j));
            }
            else
                return( (L_Bn(n,j) - L_Bn(n,j+1) )/2.0 - 1.0/(i*j) );

    };

    auto L_Rhoii = [=](const size_t n, const size_t i){
        T R = L_Sigmaii(n,i) + L_Sigmaii(n,n-i) + 2*L_Sigmaij(n,i,n-i) ;
        if(i==n-i) R/=4.0;
        return R;
    };
    auto L_Rhoij = [=](const size_t n, const size_t i, const size_t j){
        T R = L_Sigmaij(n,i,j) + L_Sigmaij(n,i,n-j) + L_Sigmaij(n,n-i,j) + L_Sigmaij(n,n-i,n-j) ;
        if((i==n-i) && (j==n-j) ) R/=4.0;
        else if(i==n-i) R/=2.0;
            else if (j==n-j) R/=2.0;

        return R;
    };

    cont=1;
    T betan=0.0;
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](const T& omegainvphi_i){betan+=(omegainvphi_i*omegainvphi_i*L_Rhoii(n,cont));++cont;});
    T bpart=0;
    for(auto i =0; i < (omegainvphi.size()-1); ++i)
        for(auto j =i+1; j < omegainvphi.size(); ++j){
            bpart+=omegainvphi[i]*omegainvphi[j]*L_Rhoij(n,i+1,j+1);
        }
//cout<<L_Sigmaii(n,1)<<" + " <<L_Sigmaii(n,n-1)<<" + "<<2*L_Sigmaij(n,1,n-1)<< " = "<<L_Rhoii(n,1) <<endl;
//cout<<betan<<" "<<L_Sigmaij(n,n-1,2)<<endl;
    betan += 2.0*bpart;

    T a1=1.0,a2=1.0; // (3) and (4) in Tajima 89
    for(auto i=2; i< n; ++i){
        a1+= 1.0/i;
        a2+= 1.0/(i*i);
    }
//From Tajima 89 p. 588
    T etheta = S / a1;
    T etheta2 = (S*(S-1)) / (a1*a1 + a2);

//cout<<etheta<<" "<<etheta2<<" "<<alphan<<" "<<betan<<endl;

/// The denominator in (16) is
    T test_denom = sqrt(alphan* etheta + betan*etheta2);


/// The Test is
    cont =0;
    T Tomtest = 0.0;
    for_each(omegainvphi.begin(), omegainvphi.end(), [&](const T& omegainvphi_i){Tomtest+=(omegainvphi_i*fsfs[cont]);++cont;});

    Tomtest/=test_denom;

    return Tomtest;
}



// DTajima implemented via the general test in Achaz 2009. The variable zeros indicate how many weights should be ignored.
// That means if zeros=0 all the spectra is considered if zeros=1 singletons are discarded and thus omega[0]=0 if zeros=2 omega[0]=omega[1]=0 and so on
template <typename T=double,typename Tcol=vector<T>>
inline T gD_Test(const Tcol& fsfs, const size_t& n, const size_t& S, const size_t zeros){


/// The thetaS corresponds to phi see Table 1 in Achaz09
    vector<T> phi(fsfs.size(),n);
    size_t cont = 1;
    for_each(phi.begin(), phi.end(), [&](T& phi_i){phi_i/= (cont*(n-cont));++cont;});
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide last element by two because n/2 has been summed twice
        phi.back()/=2.0;

/// The watterson theta corresponds to i(n-i)phi = n see Table 1 in Achaz09
    vector<T> w(fsfs.size(),n);
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide last element by two because n/2 has been summed twice
        w.back()/=2.0;

    return General_FS_Test(fsfs, w,phi, n, S,zeros);
}

// Test for isolation with migration in Achaz 2009

template <typename T=double,typename Tcol=vector<T>>
inline T gIm_Test(const Tcol& fsfs, const size_t& n, const size_t& S, const double& p, const size_t zeros){

/// The thetaS corresponds to phi see Table 1 in Achaz09
    vector<T> w1(fsfs.size(),0.0);
    size_t cont = 1; // because the first position of sfsf stores the class for 1 or n-1

    for_each(w1.begin(), w1.end(), [&](T& w1_i){w1_i= (combinatory_tab(n,cont) * pow(p,(T)cont) * pow(1.0-p,(T)n-cont)  );++cont;});

/// The watterson theta corresponds to i(n-i)phi = n see Table 1 in Achaz09
    vector<T> w2(fsfs.size(),1.0);
    if(n%2==0){ // n is even so the Kronecker delta is 1 just in this case divide last element by two because n/2 has been summed twice
        w1.back()/=2.0;
        w2.back()/=2.0;
    }
    return General_FS_Test(fsfs, w1,w2, n, S,zeros);
}

// Test for an excess of intermediate frequencies. This one uses always p=0.5 which cancels out when computing omega in the general test so can be excluded from w1
template <typename T=double,typename Tcol=vector<T>>
inline T gIm_Test(const Tcol& fsfs, const size_t& n, const size_t& S, const size_t zeros){

/// The thetaS corresponds to phi see Table 1 in Achaz09
    vector<T> w1(fsfs.size(),0.0);
    size_t cont = 1; // because the first position of sfsf stores the class for 1 or n-1
    //double p=pow(.5,(double)n);
    for_each(w1.begin(), w1.end(), [&](T& w1_i){w1_i= (combinatory_tab(n,cont) );++cont;});

/// The watterson theta corresponds to i(n-i)phi = n see Table 1 in Achaz09
    vector<T> w2(fsfs.size(),1.0);

    if(n%2==0){ // n is even so the Kronecker delta is 1 just in this case divide last element by two because n/2 has been summed twice
        w1.back()/=2.0;
        w2.back()/=2.0;
    }
    return General_FS_Test(fsfs, w1,w2, n, S,zeros);
}


/// Optimal test following Ferreti et al 2010 (formulae 10 & 11 in Supplementary)
// efsfs is the expected folded spectrum under the scenario we desire to test
template <typename T=double,typename Tcol=vector<T>>
inline T Opt_Test(const Tcol& efsfs,const Tcol& ofsfs, const size_t& n, const size_t& S, const size_t zeros){

    auto exit_val=0;
    auto halfn = efsfs.size(); // number of sequences

    if(halfn<=1){
        exit_val++;
        printval<int>(halfn,"Opt_Test.log");
        exit(exit_val);
    }

    assert(efsfs.size() == n/2);

    T sumefsfs = accumulate(efsfs.begin(),efsfs.end(),0.0);

    Tcol Ef(efsfs);
    for_each(Ef.begin(), Ef.end(), [=](T& Ef_i){Ef_i/= (sumefsfs);});

    auto L_An = [](const size_t n) {T an =0.0;for(size_t i=1;i<n; ++i){an+= 1.0/i;} return an; }; // use as L_An(n)

    vector<T> phi(halfn,n);
    for(auto i=0; i<halfn; ++i){
        phi[i]/= ((i+1)*(n-i-1));
    }
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide by two because n/2 has been summed twice
        phi.back()/=2.0;

    for_each(phi.begin(), phi.end(), [=](T& phi_i){phi_i/= (L_An(n) );});

    Tcol omega(Ef);
    size_t cont = 0;
    for_each(omega.begin(), omega.end(), [&](T& omega_i){omega_i-=phi[cont]; ++cont;} );

    return General_FS_Test(ofsfs, omega, n, S,zeros);
}

/// Optimal test for detecting an excess of intermediate freqs

template <typename T=double,typename Tcol=vector<T>>
inline T Opt_IMTest(const Tcol& ofsfs, const size_t& n, const size_t& S, const size_t zeros){

    auto exit_val=0;
    auto halfn = ofsfs.size(); // number of sequences

    if(halfn<=1){
        exit_val++;
        printval<int>(halfn,"Opt_IMTest.log");
        exit(exit_val);
    }

    assert(ofsfs.size() == n/2);

    Tcol efsfs(ofsfs.size(),0.0);
    size_t cont=1;
    for_each(efsfs.begin(), efsfs.end(), [&](T& efsfs_i){efsfs_i= combinatory_tab(n,cont);++cont;});

    T sumefsfs = accumulate(efsfs.begin(),efsfs.end(),0.0);

// The number of segregating sites must be S also in the expected spectrum
// In fact not because the expected spectrum are passed only as relative weights
//    T factor = S/sumefsfs;
//    for_each(efsfs.begin(), efsfs.end(), [=](T& efsfs_i){efsfs_i*=factor;});
//    // Get the accumulation again that should have value S
//    sumefsfs = accumulate(efsfs.begin(),efsfs.end(),0.0);

    //if(S!=(size_t)sumefsfs )
        //cout<< sumefsfs<<" "<<S<<endl;
//    assert(S==(size_t)sumefsfs);

    Tcol Ef(efsfs);

    for_each(Ef.begin(), Ef.end(), [=](T& Ef_i){Ef_i/= (sumefsfs);});

    auto L_An = [](const size_t n) {T an =0.0;for(size_t i=1;i<n; ++i){an+= 1.0/i;} return an; }; // use as L_An(n)

    vector<T> phi(halfn,n);
    for(auto i=0; i<halfn; ++i){
        phi[i]/= ((i+1)*(n-i-1));
    }
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide by two because n/2 has been summed twice
        phi.back()/=2.0;

    for_each(phi.begin(), phi.end(), [=](T& phi_i){phi_i/= (L_An(n) );});

    Tcol omega(Ef);
    cont = 0;
    for_each(omega.begin(), omega.end(), [&](T& omega_i){omega_i-=phi[cont]; ++cont;} );

    return General_FS_Test(ofsfs, omega, n, S,zeros);
}


/// Optimal test for detecting an excess of freqs of value k via a poisson with parameter k
/// E.g. for an excess of doubletons k=2
template <typename T=double,typename Tcol=vector<T>>
inline T Opt_kTest(const Tcol& ofsfs, const size_t& n, const size_t& k, const size_t& S, const size_t zeros){

    auto exit_val=0;
    auto halfn = ofsfs.size(); // number of sequences

    if(halfn<=1){
        exit_val++;
        printval<int>(halfn,"Opt_kTest.log");
        exit(exit_val);
    }

    assert(ofsfs.size() == n/2);


    std::default_random_engine generator;
    std::poisson_distribution<int> poisdistrib((double)k);

    Tcol efsfs(ofsfs.size(),0.0);

    //for_each(efsfs.begin(), efsfs.end(), [=](T& efsfs_i){efsfs_i= (T)number;});
    for(auto i=0;i < efsfs.size(); ++i){
        int number = poisdistrib(generator);
        efsfs[i]=number;
    }

    T sumefsfs = accumulate(efsfs.begin(),efsfs.end(),0.0);

// The number of segregating sites must be S also in the expected spectrum
// In fact not because the expected spectrum are passed only as relative weights
//    T factor = S/sumefsfs;
//    for_each(efsfs.begin(), efsfs.end(), [=](T& efsfs_i){efsfs_i*=factor;});
//    // Get the accumulation again that should have value S
//    sumefsfs = accumulate(efsfs.begin(),efsfs.end(),0.0);
    Tcol Ef(efsfs);

    for_each(Ef.begin(), Ef.end(), [=](T& Ef_i){Ef_i/= (sumefsfs);});

    auto L_An = [](const size_t n) {T an =0.0;for(size_t i=1;i<n; ++i){an+= 1.0/i;} return an; }; // use as L_An(n)

    vector<T> phi(halfn,n);
    for(auto i=0; i<halfn; ++i){
        phi[i]/= ((i+1)*(n-i-1));
    }
    if(n%2==0) // n is even so the Kronecker delta is 1 just in this case divide by two because n/2 has been summed twice
        phi.back()/=2.0;

    for_each(phi.begin(), phi.end(), [=](T& phi_i){phi_i/= (L_An(n) );});

    Tcol omega(Ef);
    size_t cont = 0;
    for_each(omega.begin(), omega.end(), [&](T& omega_i){omega_i-=phi[cont]; ++cont;} );

    return General_FS_Test(ofsfs, omega, n, S,zeros);
}


// ThetaS from Hedrick 2005 pg 433: It returns thetaS and also as argument the variance Vd from 8.7c in Hedrick 3rd ed  p. 434
// This variance is for the unfolded spectrum. The variance for the folded is computed in the general test.
template <typename T,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline double ThetaS(const Tcolofcol& sequences, size_t first, size_t last, const size_t ini, const size_t end, double& Vd){

    double S=0;
    auto exit_val=0;
    auto n = last - first; // number of sequences

    if(n<=1 || first<0){
        exit_val++;
        printval<int>(n,"Theta_S_Error.log");
        exit(exit_val);
    }

    if( (sequences.size() < n) ){
        first = 0;
        last = sequences.size();
        n=sequences.size();

    }
    //if(ini>=end) return 0;

    size_t WS = end-ini; // if the whole set of sequences is considered (n=sequences.size) then the number of segregating sites S = WS
    if(WS > sequences[first].size() ){
        exit_val+=2;
        printval<size_t>(exit_val,"Theta_S_Error.log");
        exit(exit_val);
    }

    if(n==sequences.size())
        S=WS;
    else{
        for(auto site =ini; site<end; ++site)
            for(auto i=first; i< (last-1); ++i)
    // Count up number of snps
                if(sequences[i][site]!=sequences[i+1][site]){
                    ++S;

                    break;
                }
    }

    double a1=1.0,a2=1.0;
    for(auto i=2; i< n; ++i){
        a1+= 1.0/i;
        a2+= 1.0/(i*i);
    }

    double c = (2.0*(n*n + n + 3))/ (9.0*n*(n-1));
    c-= (n+2)/(a1*n);
    c+= a2/(a1*a1);

    double inva1 = 1.0 / a1;

    double e1 = ((n+1.0)/(3.0*(n-1)) - inva1)*inva1;
    double e2 = c / (a1*a1 + a2);

    Vd = e1*S + e2*S*(S-1);

    return S/a1;

}


// As proposed in Ferretti et al 2013
// ThetaPi = Pi/L: Pairwise nucleotide difference (PI): average number of differences between pairs of sequences
// The function returns Pi (the average number of differences not by site) and as an argument L
template <typename T,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline double ThetaPi(const Tcolofcol& sequences, size_t first, size_t last, const size_t ini, const size_t end, size_t& L){

    double pnd=0;
    auto exit_val=0;
    auto n = last - first; // number of sequences

    if(n<=1 || first<0){
        exit_val++;
        printval<int>(n,"Theta_Pi_Error.log");
        exit(exit_val);
    }

    if( (sequences.size() < n) ){
        first = 0;
        last = sequences.size();
        n=sequences.size();

    }

    //if(ini>=end) return 0;

    L = end-ini;
    if(L > sequences[first].size() ){
        exit_val+=2;
        printval<size_t>(exit_val,"Theta_Pi_Error.log");
        exit(exit_val);
    }

    for(auto i=first; i< (last-1); ++i)
        for(auto j=i+1; j < last; ++j)
    // sum up number of differences between pairs of sequences
            for(auto site =ini; site<end; ++site)
                if(sequences[i][site]!=sequences[j][site])
                    ++pnd;

    // the mean pairwise sumed up WS sites
    pnd /= (n*(n-1)/2);

    return pnd;
}

// Partially specialization for deque bool and use binary operation ^
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol>>
inline double ThetaPi(const Tcolofcol& sequences, size_t first, size_t last, const size_t ini, const size_t end, size_t& L){
    double pnd=0;
    auto exit_val=0;
    auto n = last - first; // number of sequences

    if(n<=1 || first<0){
        exit_val++;
        printval<int>(n,"Theta_Pi_Error.log");
        exit(exit_val);
    }

    if( (sequences.size() < n) ){
        first = 0;
        last = sequences.size();
        n=sequences.size();

    }

    //if(ini>=end) return 0;

    L = end-ini;
    if(L > sequences[first].size() ){
        exit_val+=2;
        printval<size_t>(exit_val,"Theta_Pi_Error.log");
        exit(exit_val);
    }

    for(auto i=first; i< (last-1); ++i)
        for(auto j=i+1; j < last; ++j)
    // sum up number of differences between pairs of sequences
            for(auto site =ini; site<end; ++site)
                if(sequences[i][site]^sequences[j][site])
                    ++pnd;

    // the mean pairwise sumed up L sites
    pnd /= (n*(n-1)/2);


    return pnd;
}


// As proposed in Ferretti et al 2013
// ThetaPia = Pia/L: Pairwise nucleotide difference (PIa): average number of differences between pairs of sequences from different populations
template <typename T, typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline double ThetaPia(const Tcolofcol& metapop, const size_t p1, const size_t p2, const vector<int>& samplesizes,const size_t ini, const size_t end){

    size_t WS = end-ini;
    if(WS > metapop[p1].size() ){
        size_t val=1;
        printval<size_t>(WS,"ThetaPia_Error.log");
        exit(val);
    }

    double pnd=0;
    int shiftp1{accumulate(samplesizes.begin(),samplesizes.begin()+p1,0)}; // the sum of elements in populations previous to p1
    int shiftp2{accumulate(samplesizes.begin(),samplesizes.begin()+p2,0)}; // the sum of elements in populations previous to p2

    for(auto seqP1=0; seqP1< samplesizes[p1]; ++seqP1)
        for(auto seqP2=0; seqP2< samplesizes[p2]; ++seqP2)
// sum up number of differences between pairs different population sequences
            for(auto site =ini; site<end; ++site)
                if(metapop[seqP1 + shiftp1][site]!=metapop[seqP2 + shiftp2][site])
                    ++pnd;
    // the mean
    pnd /= (samplesizes[p1]*samplesizes[p2]);

    pnd /= WS;

    return pnd;
}

// Partial specialization for deque<bool>
template <typename Tcol=deque<bool>, typename Tcolofcol = deque<Tcol>>
inline double ThetaPia(const Tcolofcol& metapop, const size_t p1, const size_t p2, const vector<int>& samplesizes,const size_t ini, const size_t end){

    size_t WS = end-ini;
    if(WS > metapop[p1].size() ){
        size_t val=1;
        printval<size_t>(WS,"ThetaPia_Error.log");
        exit(val);
    }

    double pnd=0;
    int shiftp1{accumulate(samplesizes.begin(),samplesizes.begin()+p1,0)}; // the sum of elements in populations previous to p1
    int shiftp2{accumulate(samplesizes.begin(),samplesizes.begin()+p2,0)}; // the sum of elements in populations previous to p2

    for(auto seqP1=0; seqP1< samplesizes[p1]; ++seqP1)
        for(auto seqP2=0; seqP2< samplesizes[p2]; ++seqP2)
// sum up number of differences between pairs different population sequences
            for(auto site =ini; site<end; ++site)
                if(metapop[seqP1 + shiftp1][site]^metapop[seqP2 + shiftp2][site])
                    ++pnd;
    // the mean
    pnd /= (samplesizes[p1]*samplesizes[p2]);

    pnd /= WS;

    return pnd;
}


// This function follows Hs estimate as proposed in Ferretti et al 2013
template <typename T, typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline double Hs(const Tcolofcol& metapop,const size_t numpops, const vector<int>& samplesizes, const size_t ini, const size_t end){

    //size_t numpops = metapop.size()/ samplesize;
    size_t L=0;
    double hs=0;
    int first{0};
    int last{0};
    for(size_t p = 0; p < numpops; ++p){
        last+=samplesizes[p];
        hs+=(ThetaPi(metapop, first, last, ini, end,L)/L );
        first+=samplesizes[p];
    }

    hs/= numpops;

    return hs;

}

// This function follows Ht estimate as proposed in Ferretti et al 2013
template <typename T,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline double Ht(const Tcolofcol& metapop, const size_t numpops, const double& Hs, const vector<int>& samplesizes, const size_t ini, const size_t end){

    double ht=0;


    for(size_t p1 = 0; p1 < (numpops-1); ++p1)
        for(size_t p2 = (p1+1); p2 < numpops; ++p2)
            ht+= ThetaPia(metapop,p1,p2, samplesizes, ini, end);

    ht*=2.0;
    ht/=(numpops*numpops);
    ht+= (Hs/numpops);

    return ht;
}

// This function follows Fst estimate as proposed in Ferretti et al 2013 for NGS data
//It computes the mean fst from [ini,end) if end-ini=1 just compute the fst for the ini position
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline double FST(const Tcolofcol& metapop, const vector<int>& samplesizes, size_t ini, size_t end){

    auto numpops = samplesizes.size();

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FST_Error.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FST.log");
        return 0;
    }

// Check the sequence length (number of snps) at each population is the same
    int acum{samplesizes[0]};
    for(size_t p2=1; p2< numpops; ++p2){
        if(metapop[0].size()!=metapop[acum].size()){
            exitval+=2;
            printval<size_t>(exitval,"FST_Error.log");
            exit(exitval);
        }
        acum+=samplesizes[p2];
    }
// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FST_Error.log");
        exit(exitval);
    }

    double hs = Hs<T>(metapop,numpops, samplesizes, ini, end);
    double ht = Ht<T>(metapop, numpops, hs, samplesizes, ini, end);

    if(ht>0)
        return (max((1.0 - hs/ht),0.0));
        //return (1.0 - hs/ht);
    else{
//return 0;
        exitval+=4;
        printval<size_t>(metapop.size(),"FST_Error.log");
        printval<size_t>(numpops,"FST_Error.log");
        printval<size_t>(ini,"FST_Error.log");
        printval<size_t>(end,"FST_Error.log");
        printval<double>(hs,"FST_Error.log");
        printval<double>(ht,"FST_Error.log");

        //Tcolofcol tdata;
        //transpose<bool, deque<bool>>(metapop,tdata);
        //printcol(tdata[ini],"FST_Error.log");
        exit(exitval);
    }

}


// Compute FST distribution at each site
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void FSTdistrib(const Tcolofcol& metapop, const vector<int>& samplesizes, size_t ini, size_t end, vector<double>& fstdistrib){

    auto numpops = samplesizes.size();

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FSTd_Error.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FSTd.log");
        return;
    }

// Check the sequence length (number of snps) at each population is the same
    int acum{samplesizes[0]};
    for(auto p2=1; p2< numpops; ++p2){
        if(metapop[0].size()!=metapop[acum].size()){
            exitval+=2;
            printval<size_t>(exitval,"FSTd_Error.log");
            exit(exitval);
        }
        acum+=samplesizes[p2];
    }
// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FSTd_Error.log");
        exit(exitval);
    }

    int L = end -ini;
    if(L> fstdistrib.size()){
        exitval+=4;
        printval<size_t>(ini,"FSTd_Error.log");
        exit(exitval);
    }
    double hs=0,ht=0;
    //vector<double> fst(L,0);

    for(auto i=0; i<L; ++i){
        hs = Hs<T>(metapop,numpops, samplesizes, i, i+1);
        ht = Ht<T>(metapop, numpops, hs, samplesizes, i, i+1);
        if(ht>0){

            fstdistrib[i]=max((1.0 - hs/ht),0.0);
        }
         else{
            fstdistrib[i]=0;


//            exitval+=4;
//            printval<size_t>(metapop.size(),"FST_Error.log");
//            printval<size_t>(numpops,"FST_Error.log");
//            printval<size_t>(ini,"FST_Error.log");
//            printval<size_t>(end,"FST_Error.log");
//            printval<double>(hs,"FST_Error.log");
//            printval<double>(ht,"FST_Error.log");
//
//            //Tcolofcol tdata;
//            //transpose<bool, deque<bool>>(metapop,tdata);
//            //printcol(tdata[ini],"FST_Error.log");
//            exit(exitval);
        }

    }

}

// Compute FST distribution at each site and the average (sfttot) between ini and end
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void FSTdistribandTot(const Tcolofcol& metapop, const vector<int>& samplesizes, size_t ini, size_t end, vector<STAT>& fstdistrib, double& fsttot){

    auto numpops = samplesizes.size();

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FSTdT_Error1.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FSTdT.log");
        return;
    }

// Check the sequence length (number of snps) at each population is the same
    int acum{samplesizes[0]};
    for(size_t p2=1; p2< numpops; ++p2){
        if(metapop[0].size()!=metapop[acum].size()){
            exitval+=2;
            string msg="L1 = "+NtoS(metapop[0].size())+"; L2 = "+NtoS(metapop[acum].size())+"\n";
            printval<size_t>(exitval,msg,"FSTdT_Error2.log");
            exit(exitval);
        }
        acum+=samplesizes[p2];
    }
// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FSTdT_Error3.log");
        exit(exitval);
    }

    int L = end -ini;
    if(L> (int)fstdistrib.size()){
        exitval+=4;
        printval<size_t>(L-fstdistrib.size(),"FSTdT_Error4.log");
        exit(exitval);
    }
    double hs=0,ht=0;
    double avehs{0},aveht{0};
    //vector<double> fst(L,0);

    for(auto i=0; i<L; ++i){
        hs = Hs<T>(metapop,numpops, samplesizes, i, i+1);
        ht = Ht<T>(metapop, numpops, hs, samplesizes, i, i+1);
        if(ht>0){

            fstdistrib[i].stat=max((1.0 - hs/ht),0.0);
            //fstdistrib[i].stat=((1.0 - hs/ht));
        }
         else{
            fstdistrib[i].stat=0;


//            exitval+=4;
//            printval<size_t>(metapop.size(),"FST_Error.log");
//            printval<size_t>(numpops,"FST_Error.log");
//            printval<size_t>(ini,"FST_Error.log");
//            printval<size_t>(end,"FST_Error.log");
//            printval<double>(hs,"FST_Error.log");
//            printval<double>(ht,"FST_Error.log");
//
//            //Tcolofcol tdata;
//            //transpose<bool, deque<bool>>(metapop,tdata);
//            //printcol(tdata[ini],"FST_Error.log");
//            exit(exitval);
        }
        avehs+=hs;
        aveht+=ht;
    }
    fsttot = max((1.0 - avehs/aveht),0.0);
}

// Compute FST distribution at each site including mean and sample variance
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void FSTdistrib(const Tcolofcol& metapop, const vector<int>& samplesizes, size_t ini, size_t end, vector<double>& fstdistrib, double& mean, double& variance){

    auto numpops = samplesizes.size();

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FSTd_Error.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FSTd.log");
        return;
    }

// Check the sequence length (number of snps) at each population is the same
    int acum{samplesizes[0]};
    for(auto p2=1; p2< numpops; ++p2){
        if(metapop[0].size()!=metapop[acum].size()){
            exitval+=2;
            printval<size_t>(exitval,"FSTd_Error.log");
            exit(exitval);
        }
        acum+=samplesizes[p2];
    }// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FSTd_Error.log");
        exit(exitval);
    }

    int L = end -ini;
    if(L> fstdistrib.size()){
        exitval+=4;
        printval<size_t>(ini,"FSTd_Error.log");
        exit(exitval);
    }
    double hs=0,ht=0;
    //vector<double> fst(L,0);
    mean=0;
    variance=0;
    for(auto i=ini; i<end; ++i){
        hs = Hs<T>(metapop,numpops, samplesizes, i, i+1);
        ht = Ht<T>(metapop, numpops, hs, samplesizes, i, i+1);
        if(ht>0){

            fstdistrib[i]=max((1.0 - hs/ht),0.0);
        }
         else{
            fstdistrib[i]=0;

        }
        mean+=fstdistrib[i];
        variance+=(fstdistrib[i]*fstdistrib[i]);

    } // close for
    double s2 =mean*mean;
    mean/=L;
    s2/=L;
    variance = max((variance -s2) / (L-1),0.0); // sample variance

}

// Compute robust location parameters of an FST distribution
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void robavevarFST(const Tcolofcol& metapop, const size_t samplesize, size_t ini, size_t end,vector<double>& fstdistrib, double& locfst, double& scalfst){

    auto numpops = metapop.size()/ samplesize;

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FST_Error.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FST.log");
        return;
    }

// Check the sequence length (number of snps) at each population is the same
    for(auto p2=1; p2< numpops; ++p2)
        if(metapop[0].size()!=metapop[p2*samplesize].size()){
            exitval+=2;
            printval<size_t>(exitval,"FST_Error.log");
            exit(exitval);
        }
// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FST_Error.log");
        exit(exitval);
    }

    int L = end -ini;
    if(L> fstdistrib.size()){
        exitval+=4;
        printval<size_t>(ini,"FST_Error.log");
        exit(exitval);
    }
    double hs=0,ht=0;
    //vector<double> fst(L,0);

    for(auto i=0; i<L; ++i){
        hs = Hs<T>(metapop,numpops, samplesize, i, i+1);
        ht = Ht<T>(metapop, numpops, hs, samplesize, i, i+1);
        if(ht>0){

            fstdistrib[i]=max((1.0 - hs/ht),0.0);
        }
         else{
            fstdistrib[i]=0;


//            exitval+=4;
//            printval<size_t>(metapop.size(),"FST_Error.log");
//            printval<size_t>(numpops,"FST_Error.log");
//            printval<size_t>(ini,"FST_Error.log");
//            printval<size_t>(end,"FST_Error.log");
//            printval<double>(hs,"FST_Error.log");
//            printval<double>(ht,"FST_Error.log");
//
//            //Tcolofcol tdata;
//            //transpose<bool, deque<bool>>(metapop,tdata);
//            //printcol(tdata[ini],"FST_Error.log");
//            exit(exitval);
        }

    }

    locfst = HLe(fstdistrib);
    scalfst = Qn<double>(fstdistrib);

}




// Compute robust location and scale parameters of an FST distribution
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void robavevarFST(const Tcolofcol& metapop, const size_t samplesize, size_t ini, size_t end, double& locfst, double& scalfst){

    auto numpops = metapop.size()/ samplesize;

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FST_Error.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FST.log");
        return;
    }

// Check the sequence length (number of snps) at each population is the same
    for(auto p2=1; p2< numpops; ++p2)
        if(metapop[0].size()!=metapop[p2*samplesize].size()){
            exitval+=2;
            printval<size_t>(exitval,"FST_Error.log");
            exit(exitval);
        }
// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FST_Error.log");
        exit(exitval);
    }

    int L = end -ini;

    double hs=0,ht=0;
    vector<double> fst(L,0);

    for(auto i=0; i<L; ++i){
        hs = Hs<T>(metapop,numpops, samplesize, i, i+1);
        ht = Ht<T>(metapop, numpops, hs, samplesize, i, i+1);
        if(ht>0){

            fst[i]=max((1.0 - hs/ht),0.0);
        }
         else{
            fst[i]=0;

//            exitval+=4;
//            printval<size_t>(metapop.size(),"FST_Error.log");
//            printval<size_t>(numpops,"FST_Error.log");
//            printval<size_t>(ini,"FST_Error.log");
//            printval<size_t>(end,"FST_Error.log");
//            printval<double>(hs,"FST_Error.log");
//            printval<double>(ht,"FST_Error.log");
//
//            //Tcolofcol tdata;
//            //transpose<bool, deque<bool>>(metapop,tdata);
//            //printcol(tdata[ini],"FST_Error.log");
//            exit(exitval);
        }

    }

    locfst = HLe(fst);
    //scalfst = MAD<double>(fst, 1.4826);
    scalfst = Qn<double>(fst);
}


// Compute robust scale parameters of an FST distribution
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void robvarFST(const Tcolofcol& metapop, const size_t samplesize, size_t ini, size_t end, double& scalfst){

    auto numpops = metapop.size()/ samplesize;

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"FST_Error.log");
        exit(exitval);
    }
    if(numpops==1){
        printval<size_t>(metapop.size(),"FST.log");
        return;
    }

// Check the sequence length (number of snps) at each population is the same
    for(auto p2=1; p2< numpops; ++p2)
        if(metapop[0].size()!=metapop[p2*samplesize].size()){
            exitval+=2;
            printval<size_t>(exitval,"FST_Error.log");
            exit(exitval);
        }
// Check the snp range [ini,end) for which Fst will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"FST_Error.log");
        exit(exitval);
    }

    int L = end -ini;

    double hs=0,ht=0;
    vector<double> fst(L,0);

    for(auto i=0; i<L; ++i){
        hs = Hs<T>(metapop,numpops, samplesize, i, i+1);
        ht = Ht<T>(metapop, numpops, hs, samplesize, i, i+1);
        if(ht>0){

            fst[i]=max((1.0 - hs/ht),0.0);
        }
         else{
            fst[i]=0;

//            exitval+=4;
//            printval<size_t>(metapop.size(),"FST_Error.log");
//            printval<size_t>(numpops,"FST_Error.log");
//            printval<size_t>(ini,"FST_Error.log");
//            printval<size_t>(end,"FST_Error.log");
//            printval<double>(hs,"FST_Error.log");
//            printval<double>(ht,"FST_Error.log");
//
//            //Tcolofcol tdata;
//            //transpose<bool, deque<bool>>(metapop,tdata);
//            //printcol(tdata[ini],"FST_Error.log");
//            exit(exitval);
        }

    }


    //scalfst = MAD<double>(fst, 1.4826);
    scalfst = Qn<double>(fst);
}



// Compute vector of unfolded TajimaD's for the different populations in a metapopulation. Algorithm from Hedrick 3rd Ed. p. 434
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void DTajima(const Tcolofcol& metapop, const size_t samplesize, size_t ini, size_t end, vector<double>& DT, size_t& L){

    auto numpops = metapop.size()/ samplesize;

    auto exitval=0;
    if(numpops<1){
        exitval++;
        printval<size_t>(metapop.size(),"D_Error.log");
        exit(exitval);
    }

// Check the sequence length (number of snps) at each population is the same
    for(auto p2=1; p2< numpops; ++p2)
        if(metapop[0].size()!=metapop[p2*samplesize].size()){
            exitval+=2;
            printval<size_t>(exitval,"D_Error.log");
            exit(exitval);
        }

    DT.resize(numpops);

// Check the snp range [ini,end) for which D will be computed
    if(metapop[0].size() < end){
        end =metapop[0].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"D_Error.log");
        exit(exitval);
    }

    L=0;
    double k=0,ThS=0,Vd=0,D=0;

    for(auto p = 0; p < numpops; ++p){

        size_t first = p*samplesize;
        size_t last = (p+1)*samplesize;

        k= ThetaPi(metapop, first, last, ini, end,L) ;
        ThS = ThetaS<bool>(metapop, first, last, ini, end, Vd);

        D =  (k - ThS);

        if(Vd>0) // variance for the unfolded spectrum
            D /= sqrt(Vd);
            //return (1.0 - hs/ht);
        else{
            exitval+=4;
            printval<size_t>(metapop.size(),"D_Error.log");
            printval<size_t>(p,"D_Error.log");
            printval<size_t>(ini,"D_Error.log");
            printval<size_t>(end,"D_Error.log");
            printval<double>(D,"D_Error.log");
            printval<double>(Vd,"D_Error.log");

            Tcolofcol tdata;
            transpose<bool, deque<bool>>(metapop,tdata);
            //printcol(tdata[ini],"FST_Error.log");
            exit(exitval);
        }

        DT[p]=D;

    } // close for pops
}

// Compute vector of unfolded TajimaD's for the different populations in a metapopulation. Algorithm from Hedrick 3rd Ed. p. 434
template <typename T=bool,typename Tcol=deque<T>, typename Tcolofcol = deque<Tcol>>
inline void DTajima1Pop(const Tcolofcol& metapop, const size_t pop, const size_t samplesize, size_t ini, size_t end, double& DT, size_t& L){

    auto exitval=0;

    size_t inipop = pop*samplesize;
// Check the snp range [ini,end) for which D will be computed
    if(metapop[inipop].size() < end){
        end =metapop[inipop].size();
    }
    if(ini< 0)
        ini=0;

    if(ini>=end){
        exitval+=3;
        printval<size_t>(ini,"D1p_Error.log");
        exit(exitval);
    }

    L=0;
    double k=0,ThS=0,Vd=0,D=0;

    size_t first = inipop;
    size_t last = (pop+1)*samplesize;

    k= ThetaPi(metapop, first, last, ini, end,L) ;
    ThS = ThetaS<bool>(metapop, first, last, ini, end, Vd);

    D =  (k - ThS);
    if(Vd>0) // variance for the unfolded spectrum
        D /= sqrt(Vd);
        //return (1.0 - hs/ht);
    else{
        exitval+=4;
        printval<size_t>(metapop.size(),"D1p_Error.log");
        printval<size_t>(pop,"D1p_Error.log");
        printval<size_t>(ini,"D1p_Error.log");
        printval<size_t>(end,"D1p_Error.log");
        printval<double>(D,"D1p_Error.log");
        printval<double>(Vd,"D1p_Error.log");

        Tcolofcol tdata;
        transpose<bool, deque<bool>>(metapop,tdata);
        //printcol(tdata[ini],"FST_Error.log");
        exit(exitval);
    }

    DT=D;


}

# endif
