/** *************************************************************************************
 *
 * PROGRAM:
 *
 * 	iHDSel: Software for estimating selection from informed haplotype allelic patterns
 *
 * 	https://acraaj.webs.uvigo.es/iHDSel.html
 *
 * AUTHOR:
 *
 * 	Antonio Carvajal Rodriguez (acraaj@uvigo.es) October 2023 for v 0
 *
 * 	Area de Genetica. Universidad de Vigo.
 *
 * 	Group of Genetics of Populations and Cytogenetics (XB2)
 *
 * *************************************************************************************
 *
 * DISCLAIMER:
 * This program is free software; you can redistribute it and/or modify it under the terms of
 * the GNU General Public License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version. This program is
 * distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
 * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details. You should have received a copy of the GNU
 * General Public License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. See the file
 * "gpl.html" under the "license" directory.
 *
 * *************************************************************************************
 *
 * VERSION HISTORY:
 *
 * iHDivSel (April 2026):
 *  0.5.3:
 *  0.5.2: Minor changes to improve speed (May 2024)
 *  0.4: The computation of the outlier-centered blocks was slightly corrected to compute the blocks in the same way as in the useblcoks case (April 2024)
 *  0.3: Correct a bug so that the haplotypes blocks could be unusually larger if some fixed SNPs happened in the refrence population
 *  0.2:
 *  - The SNPs at the beginning and end of each block are given in the output file
 *  - Under the -useblocks 0 option if mean Fst is >0.5 there is not outlier test (EOS) performed but the 90-percentile are used as centers
 * of haplotpe blocks
 *
 *
 * iHDivSel (November 2023):
 *  0: this version evolves from HACDivSel v1.6 program. It eliminates nvd computations and substitute them by Jeffreys (Jhac) divergence.
 *  - VCF input format added
 *  - nvd (normalized variance difference) method changed to an information-based statistic with chi-squared distribution.
 *  - window sizes are automatically computed by estimating haplotype blocks
 *
 * ****************************
 * Old HacDivSel:
 *  0.7.1: this version improves the nvd test fast option that becomes the default when the number of true SNPs is above 100,000
 *  1.4: Added the input format from FastPhase
 * 1.5: Added the possibility of study candidate(s) SNP(s)
 * 1.6:
 *  - the input value of Minor Allele Frequency maf changed to be a per population relative value (now by default 0.05) previously was an absolute metapopulation value (by default 4)
 *  - nvd (normalized variance difference) method changed to a robust variance homogeneity test.
 *  - the test can now be applied to each population
 *  - window sizes are computed by estimating haplotype blocks

 * ************************************************************************************* */

#include "SNPSelTools.h"
#include "RandomFuncs.h"

/// ************ GLOBAL CONSTANTS **********

const double defaultMAF{0.01}; //
const size_t defaultseed = SEED;

/// acr July 21: test mode for testing code purposes only
const bool TESTMODE=false;

// These are upper bounds for reserving memory
const int MEM_SAMPLESIZE = 100;
const int MEM_NUMPOPS = 5;
const int MEM_NUMSNPS = 1E6;
const double sdpbprecision=0.001;

// Iterations for resampling in the computation of nvd p-values when there are many SNPs
const size_t rsmp_lowiter =1E2;
const size_t rsmp_mediumiter=2E2;
const size_t rsmp_iter=5E2;

/// verbosity
bool verbose=false;

/// acr Nov 21 Threshold for LD block depending if sweep is recent or old
bool old=false;
double DBT =0.4;

//const int resampnumiter = 1;

// ************ TYPES DEFINITION **********
struct CRange{ // Store the initial and final position of a candidates set
    int inisite{0};
    int endsite{0};
};

struct VarTest{
    vector<double> v1;
    vector<double> v2;
    vector<double> F;
    vector<double> Fpval;
};


/// ************************************ GLOBAL VARIABLES **********************************

    bool flag_therearecommons=true;
    size_t resampnumiter = rsmp_iter; //5E2 by default
	typedef chrono::high_resolution_clock myclock;
    myclock::time_point beginning = myclock::now();

// Initialize the seed for the random generator
    auto seed = std::chrono::duration_cast<std::chrono::seconds>(beginning.time_since_epoch()).count();

//auto seed2 = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();


/// ************************************ FUNCTION PROTOTYPES **********

template<typename TBcol=vector<bool>>
TBcol computereference(const deque < deque<bool> >&datapop,const int& Pref, const vector<int>& samplesizes);

template<typename TBcol=vector<bool>>
vector<TBcol> computereferencepops(const deque < deque<bool> >&datapop, const vector<int>& samplesizes);

template<typename TR=void>
TR computeSNPFreqinPop(const deque< deque<bool> >&data,const int& pop, const vector<int>&samplesizes, vector<double>& snpfreqs);

template<typename TR=void>
TR computemeanSNPFreqinMetaPop(const deque< deque<bool> >&data, vector<double>& snpfreqs);

template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void copycommonwithMAF(const vector<int>& samplesizes,const Tcolofcol& positions,const Tcolofcol& freqs, Tcol& commons, const double& maf);

template<typename T=double>
void updateMAFData(deque< deque<bool> >&data,const int& pop, const vector<int>&samplesizes, const T& MAF,
               vector<T>& snpfreqs, vector<T>& snppos, const vector<T>& commons, const bool onlycommons);
template<typename T=double>
void deleteNonCommonData(deque< deque<bool> >&data,const int& pop, const vector<int>&samplesizes, vector<T>& snpfreqs, vector<T>& snppos, const vector<T>& commons);

template<typename T=double>
bool dataFusion(const bool shared,const deque< deque<bool> >&data, deque< deque<bool> >& newdata, const vector<vector<T>>& snpositions, vector<T>& metapopsnpositions, const vector<int>&samplesizes, int& metapopsamplesize);

template<typename T=double>
bool dataFusion(const deque< deque<bool> >&data, deque< deque<bool> >& newdata, const vector<vector<T>>& snpositions,
            vector<T>& metapopsnpositions, const vector<int>&samplesizes, int& metapopsamplesize);

template<typename TR=void>
TR bufferAnalyzedSNPs(const vector<string>& Popnames, const vector<string>& rsids, const vector<vector<double>> snpositions,const vector<vector<double>> snpfreqs, string& bufferdata);

template<typename T=string, typename TR=void>
TR writeAnalyzedSNPs(const vector<T>& outnames, const vector<T>&buffers);

template<typename T=string, typename TR=void>
TR writeJtest(const vector<bool>& highmeanFST,const vector<int>&totsnps,const vector<double>& fstmean,const size_t& totruns,const string& foutname,const vector<string>& names,const vector <vector<ResultIndexes> > indexesfiles,
              vector <vector<ResultFst> > fstfiles, const vector<int>& contSLsig,const double& SL,
              double& SegSites,vector<HapMapTable>& snpstable, const size_t& totalmaxcandidates,
              const int& numpops, const vector<vector<int>>& samplesizePops);

template<typename T=string, typename TR=void>
TR writeVartest(const string& foutname,const vector<string>& names,const vector <vector<ResultIndexes> > indexesfiles,
              vector <vector<ResultFst> > fstfiles, const vector<int>& contSLsig,const double& SL,
              double& SegSites,
              const vector<int>& numoutliers, const vector<int>& totsnps,const vector<int>& numindepsnps,const double& iniMAF,
              const vector<double>& fstmax,const vector<double>& fstmean,const vector<double>& pi0,const vector<double>& LBqval,
              vector<HapMapTable>& snpstable, const size_t& totalmaxcandidates,
              const int& numpops, const vector<vector<vector<int>>>& samplesizeP1, const vector<vector<vector<int>>>& samplesizeP2,
              const vector<vector<VarTest>>&VarTests, const vector<vector<vector<bool>>>& skip, const vector<size_t>& testtype);

template<typename TR=void>
TR writeFSTs(const string& foutname,const vector<string>& names,const vector <vector<ResultFst> > fstfiles0,
              const vector <vector<ResultFst> > fstfiles,const double& SL,
              const vector<int>& totsnps, const double& iniMAF);



/// ************ FUNCTIONS **********

template<typename T=double>
void writestats(const string& outname, const vector<T>& ave, const vector<T>& sde,const int& lastwinsize, const int& winredfactor){

    ofstream canalsalida(outname.c_str(),ios::app);
    canalsalida<<"WS\tAve\tSDE"<<endl;
    double ws= lastwinsize;
    for(int wi=ave.size()-1; wi>=0; --wi){
        canalsalida<<ws<<"\t"<<ave[wi]<<"\t"<<sde[wi]<<endl;
        if(wi==(int)ave.size()-1) ws+=.5;
        ws*= winredfactor;
    }
    canalsalida.close();
}

template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void doIndexNullMetapop(const Tcol& meansnpfreqs,const size_t numpops,const size_t nummetapops,
            const vector<int>& samplesizes, vector<ResultIndexes>& problemvalue,
            const size_t metapopindex, const size_t numiter, const Tcol&fstdistrib, const T& SL){

//
// This is just for one metapopulation

    if(nummetapops<(metapopindex+1)){
        string msg = "\nERROR in doIndexNullMetapop. Metapop index is too high "+ NtoS(metapopindex)+"\n";
        cout<<msg<<endl;
        SNPSelTools::texterr(1,msg,true);
    }

    size_t metapopsize= accumulate(samplesizes.begin(),samplesizes.end(),0);

    // Check that the original data was grouped i.e. number of snps has been settled to be the same among populations
    //assert(selectivesites.size()==numpops);

    // The haplotype length
    size_t L = meansnpfreqs.size();

    assert(fstdistrib.size()==L);
    // Data
    deque<deque<bool>> datanull(metapopsize,deque<bool>(L,0));

    double U1;
    struct simPval{
        T pval;
        T qval;
        int pos;
        simPval():pval{0},qval{1}, pos{-1}{}
        simPval(int p):pval{0},qval{1}, pos{p}{}
    //bool operator==(const simPval& rhs) const {return this->pval==rhs.pval;}
    bool operator<(const simPval& rhs) const {return this->pval<rhs.pval;}
    bool operator<=(const simPval& rhs) const {return this->pval<=rhs.pval;}

    };


    // If correlation dependence exits among snps this is already reflected in the frequencies that will be used to generate the samples
    vector<simPval>pvals;
    pvals.reserve(L);

    for(size_t l=0; l<L; ++l)
        pvals.push_back(simPval(l)); // each pval is associated to a specific position

    Tcol simfreqs(L,0); // To store the simulated frequency of the L sites in the metapop
    // Because stat[i] = fstsel[i] - meanfst
    T meanfst = fstdistrib[problemvalue[0].site] - problemvalue[0].fstindex.stat;
    // Repeat numiter times
    cout<<"Beginning empirical neutral resampling..."<<endl;
    for(size_t iter=0; iter<numiter; ++iter){


    /// Generate the metapop data. Note that at some sites population frequencies could be zero if we are working with non-common snps

        for(size_t smpl = 0; smpl < metapopsize; ++smpl){
     // For each site i generate 1 with probability snpfreqs[i] or 0 otherwise
            for(size_t site=0; site < L; ++site){
                U1=uniform(gen);

                if( U1 < meansnpfreqs[site]){
                    datanull[smpl][site]=1;
                    ++simfreqs[site];
                }
                else
                    datanull[smpl][site]=0;
            }

        }//close for sample

    /// Compute Fst distribution

        Tcol simfstdistrib(L,0);
        T simmeanfst{0};
        FSTdistribandTot(datanull, samplesizes, 0, datanull[0].size(), simfstdistrib,simmeanfst);
    /// Accumulate occurrences for pvalues
        for(size_t snp=0;snp< L; ++snp){
            if( (simfstdistrib[snp]-simmeanfst)>= (fstdistrib[snp]-meanfst)){
                pvals[snp].pval+=1;
            }
        }

        //double fstSel=0;
        //double simfstindex = 0;


    // Reset frequencies to zero at each iteration
        for(auto& f:simfreqs){

            f=0;
        }
        if(verbose) cout<<"iteration "+NtoS(iter)+" done\n";
    } // Close empirical neutral resampling iterations
    cout<<"End empirical neutral resampling..."<<endl;

    /// Compute p-values

    for(auto& item:pvals)
        item.pval/=numiter;

    /// sort pvalues in ascending order. This can be done without adding a comparator because operator < is overloaded in the struct simPval definition (see above)
    //sort(pvals.begin(),pvals.end(), [](const simPval& p1, const simPval& p2){return p1.pval<p2.pval;} );
    sort(pvals.begin(),pvals.end());
    /// Compute Pi0
    T pi0{1.0 - sdpbound<double, vector<simPval>>(sdpbprecision, pvals, SL)};

    if(pi0<=0){
        string msgerr = "ERROR: pi_0 = ";
        std::ostringstream strs;
        strs << pi0;
        msgerr += strs.str();
        msgerr += " using SDB method";
        cout<<msgerr<<endl;
        SNPSelTools::texterr(2,msgerr,true);
    }

    /// Compute q values

    double minq{1.0};

    for(int t=(pvals.size()-1);t>=0;--t){

        pvals[t].qval = pi0*L*pvals[t].pval /((t+1));
        if(minq<pvals[t].qval) pvals[t].qval=minq; // because the qvalue is the minimum pFDR that can occur when rejecting a statistic with value t for the set of nested significance regions (Storey 2003)
        else minq=pvals[t].qval;
    //cout<<pvals[t].pval<<"\t"<<(t+1)<<"\t"<<pvals[t].pos<<endl;
    }

    //cout<<pvals[problemvalue[problemvalue.size()-1].site].pos<<endl;
    //cout<<pvals[problemvalue[problemvalue.size()-1].site].pval<<endl;
    //cout<<pvals[problemvalue[problemvalue.size()-1].site].qval<<endl;

    /// Store the q-values for the sites identified as selective
    // sort the pvalues now by position in ascending order
    sort(pvals.begin(),pvals.end(), [](const simPval& p1, const simPval& p2){return p1.pos<p2.pos;} );

    for(size_t candidate=0; candidate<problemvalue.size(); ++candidate){

        problemvalue[candidate].fstindex.pval=pvals[problemvalue[candidate].site].pval;
        problemvalue[candidate].fstindex.qval=pvals[problemvalue[candidate].site].qval;

    }
    //cout<<problemvalue[0].site<<"\t"<<pvals[problemvalue[0].site].pval<<"\t"<<problemvalue[0].fstindex.pval<<endl;
    //cout<<pi0<<"\t"<<L<<"\t"<<problemvalue[0].fstindex.pval<<"\t"<<problemvalue[0].fstindex.qval<<endl;exit(0);

    //    for(auto& candidates: problemvalue)
    //        candidates.fstindex.pval/= numiter;

    // clean data
    for(auto& dn: datanull){
        //dn.clear();
        deque<bool>().swap(dn);
    }
    //datanull.clear();
    deque<deque<bool>>().swap(datanull);

} // close doindexnullmetapop

/// ACR June 21: doindex null without computing q-values
template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void doIndexNullMetapop2(const Tcol& meansnpfreqs,const size_t numpops,const size_t nummetapops,
            const vector<int>& samplesizes, vector<ResultIndexes>& problemvalue,
            const size_t metapopindex, const size_t numiter, const Tcol&fstdistrib, const T& SL){

//
// This is just for one metapopulation

    if(nummetapops<(metapopindex+1)){
        string msg = "\nERROR in doIndexNullMetapop2. Metapop index is too high "+ NtoS(metapopindex)+"\n";
        cout<<msg<<endl;
        SNPSelTools::texterr(1,msg,true);
    }

    size_t metapopsize= accumulate(samplesizes.begin(),samplesizes.end(),0);

// Check that the original data was grouped i.e. number of snps has been settled to be the same among populations
//assert(selectivesites.size()==numpops);

    // The haplotype length
    size_t L = meansnpfreqs.size();

    assert(fstdistrib.size()==L);
// Data
    deque<deque<bool>> datanull(metapopsize,deque<bool>(L,0));

    double U1;
    struct simPval{
        T pval;
        int pos;
        simPval():pval{0}, pos{-1}{}
        simPval(int p):pval{0}, pos{p}{}

    };

  // If correlation dependence exits among snps this is already reflected in the frequencies that will be used to generate the samples
    vector<simPval>pvals;
    pvals.reserve(L);

    for(size_t l=0; l<L; ++l)
        pvals.push_back(simPval(l)); // each pval is associated to a specific position

    T meanfst = fstdistrib[problemvalue[0].site] - problemvalue[0].fstindex.stat; // because stat[candidate] = fst[candidate]-fstmean so fstmean=fst[candiate]-stat[candidate] y fst mean is the same whatever the candidate
// Repeat numiter times
cout<<"Beginning empirical neutral resampling..."<<endl;
    for(size_t iter=0; iter<numiter; ++iter){


/// Generate the metapop data. Note that at some sites population frequencies could be zero if we are working with non-common snps

        for(size_t smpl = 0; smpl < metapopsize; ++smpl){
     // For each site i generate 1 with probability snpfreqs[i] or 0 otherwise
            for(size_t site=0; site < L; ++site){
                U1=uniform(gen);

                if( U1 < meansnpfreqs[site]){
                    datanull[smpl][site]=1;
                }
                else
                    datanull[smpl][site]=0;
            }

        }//close for sample

/// Compute Fst distribution

        Tcol simfstdistrib(L,0);
        T simmeanfst{0};
        FSTdistribandTot(datanull, samplesizes, 0, datanull[0].size(), simfstdistrib,simmeanfst);
/// Accumulate occurrences for pvalues
        for(size_t snp=0;snp< L; ++snp){
            if( (simfstdistrib[snp]-simmeanfst)>= (fstdistrib[snp]-meanfst)){
                pvals[snp].pval+=1;
            }
        }

        if(verbose) cout<<"iteration "+NtoS(iter)+" done\n";
    } // Close empirical neutral resampling iterations

        if(verbose) cout<<"End empirical neutral resampling..."<<endl;

/// Compute p-values

    for(auto& item:pvals)
        item.pval/=numiter;

    for(size_t candidate=0; candidate<problemvalue.size(); ++candidate){
        if(verbose){
            cout<<"Candidate "<<problemvalue[candidate].site<<":"<<endl;
            cout<<"Fst p-val "<<problemvalue[candidate].fstindex.pval<<"."<<endl;
        }
        problemvalue[candidate].fstindex.pval=pvals[problemvalue[candidate].site].pval;

    }

    // clean data
    for(auto& dn: datanull){
        //dn.clear();
        deque<bool>().swap(dn);
    }
    //datanull.clear();
    deque<deque<bool>>().swap(datanull);

} // close doIndexNullMetapop2



/// close do_nvdnull

/// acr June 21 Generate population sequence file in a given format and write to hard disk
///meansnpfreqs shoud be a vector of SNP frequencies

template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void GenPop(const Tcol& meansnpfreqs,const size_t numpops, const vector<int>& samplesizes, const string& format,const string& foutname){


    size_t metapopsize= accumulate(samplesizes.begin(),samplesizes.end(),0);


    // The haplotype length
    size_t L = meansnpfreqs.size();

// Data
    deque<deque<bool>> datanull(metapopsize,deque<bool>(L,0));

    double U1;

    Tcol simfreqs(L,0); // To store the simulated frequency of the L sites in the metapop

/// Generate the metapop data. Note that at some sites population frequencies could be zero if we are working with non-common snps

    for(size_t smpl = 0; smpl < metapopsize; ++smpl){
 // For each site i generate 1 with probability snpfreqs[i] or 0 otherwise
        for(size_t site=0; site < L; ++site){
            U1=uniform(gen);

            if( U1 < meansnpfreqs[site]){
                datanull[smpl][site]=1;
                ++simfreqs[site];
            }
            else
                datanull[smpl][site]=0;
        }

    }//close for sample


   // Reset frequencies to zero at each iteration
        for(auto& f:simfreqs){

            f=0;
        }


cout<<"Writing the file..."<<endl;


    for(auto& dn: datanull){
        //dn.clear();
        deque<bool>().swap(dn);
    }
    //datanull.clear();
    deque<deque<bool>>().swap(datanull);

}


/// end GenPop



/// delete some snps from the whole data matrix delete from data and the positions vector.
//The frequencies should be computed after this function
// The number of deleted snps depends on the step.
//Let the total number of snps be L then int(L/steps) SNPs are deleted beginning in the first SNP (the 0 position)
// if step=1 all the snps are deleted. If step =L only the first SNP is deleted. If step>L no snp is deleted.

template<typename T=double>
void deleteSNPs(deque< deque<bool> >&data,vector<T>& metapopsnpositions, const int& step){

    // 1.- Update Data, frequencies and position vectors
    auto L=metapopsnpositions.size();
    auto totindivs=data.size();
    int delsnp{0};
    if(step<=L)
        do{
            for(int indiv=0; indiv<totindivs;++indiv){

                data[indiv].erase(data[indiv].begin()+delsnp);

            }
            metapopsnpositions.erase(metapopsnpositions.begin()+delsnp);
            --L;
            // update the next snp to delete
            delsnp+=(step-1); // minus 1 because each time one position has been already deleted

        }while(delsnp<L);
}

// Change the phase of some SNPs separated by step in randomly selected gametes
template<typename T=double>
void changeSNPs(deque< deque<bool> >&data,const int& step){

    auto totindivs=data.size();
    auto L=data[0].size();
    int snp{0};
    if(step<=L)
        do{

            // Randomly choose an individual and flip the value
            int indiv=uniform(gen)*totindivs;

            data[indiv][snp]=!data[indiv][snp];

            // update the next snp to flip
            snp+=(step); // minus 1 because each time one position has been already deleted

        }while(snp<L);
}

// Change the phase of some fraction of randomly selected SNPs in every gamete
template<typename T=double>
void changeSNPsbygam(deque< deque<bool> >&data,const T& fr){

    auto totindivs=data.size();
    auto L=data[0].size();
    int errors = L*fr;
    int snp{0};
    if(fr<=1)
    for(int indiv=0; indiv<totindivs;++indiv){
        for(int e=0; e<errors; ++e){
            // Randomly choose an individual and flip the value
            int snp=uniform(gen)*L;

            data[indiv][snp]=!data[indiv][snp];

        }
    }
}


/// ************************************ MAIN PROGRAM **********************************

int main(int argc, char** argv){
	cout<<endl;
	cout<< "---------------------------------------------------------------------\n"<<endl;
	cout<< " iHDSel: Program to identify patterns of natural selection using information-based statistics\n"<<endl;
	cout<< " by A Carvajal-Rodriguez			\n"<<endl;
	cout<< " https://acraaj.webs.uvigo.es/iHDSel.html			\n"<<endl;
	cout<< "---------------------------------------------------------------------\n"<<endl;
	//cout<< " ***********************************************************\n"<<endl;
	cout<<endl;

const float ver = 0.3;
std::remove("iHDSel_ERROR.log"); /// 28 abril 2026. Borramos logs de error antiguos

/// ********************************************* PARAMETERS
//int testintmax = numeric_limits<int>::max();
//cout << testintmax<<endl; exit(0);

	string tag=""; //FEB 2024
	bool store=false; //FEB 2024
	//bool doEOS=true; //FEB 2024
	bool doEOS=false; //April 2026
	bool nogaps=false; //FEB 2024
	bool userfile=false,fixseed=false;
	bool skiprun=false;
	bool vcfformat=false, msformat=true,bsformat=false, fastaformat=false, gpformat=false, hapmapformat=false, plformat{false},plnumeric{false},fastPhaseformat{false};
	bool soloEOS{false},EPOhighest{false};
	uint_fast8_t diploid=1;
	bool onepopbyfile{false}; // this identifies if the two populations are in separate files (true)
    bool writersidInfo{false};
    bool onlyrsidinfo{false};
    size_t totalmaxcandidates{0}; // this is useful when analyzing several files with different selective sites
	size_t totruns{0};
	int numpops=2,nummetapops=1, /*inisamplesize,*/samplesize=50,samplesize2=0,metapopsamplesize=0, NumOfFiles=1;
	int ws,winsize,fixwinsize=0; // default window size

	int iniw,endw,center;
	int MaxWinSize=1000,minWinSize=11, minblocksize=minWinSize;

	/// Dec 21
	int popforHB=-1;
	int Pref=0;
	/// May 2024
	bool autorefpop=true;
	/// Jan 22 OB threshold
	//int OBtrh=40; // over a value of 30-40 violation of normality should not cause major problems (Ghasemi and Zahediasl 2012)

	float start,secs;
	double SegSites{0};

	double numdesvs{1.25};

	double SL=0.05/*,correctedSL=0*/;
	string dir,mensaje="",path="", inputname="", inputname2="",defaultFname="iHDSF",outheader="iHDS_out";
    string format="MS";

    // acr July 21
    vector < vector<bool> >popreferences;
    /// May 2024
    vector<bool> referencehap;
    /// snps freqs store the frequency of 1's, which for most formats are the minor allele but for the ms format could be the major,
    vector < vector<double> > snpfreqs; // Dimensions are: pop, frequencies. So, snpfreqs size will give the number of pops
    vector < vector<double> > snpositions; // Dimensions are: pop, positions. So, snpositions size will give the number of pops
    vector < vector<double> > metapopsnpfreqs;
    vector<int> numsnps;

    vector<int> samplesizes;

// Vectors for the different runs or files
	vector<string> namefiles;
	vector <vector<ResultIndexes> > indexesfiles;
    vector <vector<ResultFst> > fstfiles;
    vector <vector<ResultFst> > fstfiles0;
// vectors for selection analysis

    vector<double> Dprime,r2, DprimeB,r2B; //May 2024 DprimeB and r2B to compute for the second population
    vector<vector<double>> Dp_runs,r2_runs;

    /// MAF variables:
    /** Minimum Allele Frequency (default is defaultMAF=0.01).


    **/
    double iniMAF{defaultMAF}, popMAF{defaultMAF};

    ///acr June 21
    vector<double> selectives; // to store candidate snps positions. Use double because the positions could be relative as in ms e.g. position 0.5. so if ms format selective 0.5 indicate not the position named as 0.5 but to use the position in the middle of the set of snps
    vector<int> selindexes;

    vector<double> candidates; //store as candidates the position of outliers selected in epos
    /// acr August 21: variances for partition 1 and 2
    //vector<VarTest>VarFTest(numpops);
    //vector<vector<int>>n1(numpops),n2(numpops); // second dimension is for each selindex
    //vector<int>metan2;

    /// Oct 2023
    vector<double>JHac;

    int NUMCLASSES=0;
    bool useblocks=true;
    //bool nooutliers=false;

    /// Oct 21 Vectors for the different runs or files
    //vector<vector<VarTest>>runsVarFTest;
    //vector<vector<vector<int>>>runsn1,runsn2;
    ///Dec 21
     vector<vector<bool>>skip(numpops); // second dimension is for each selindex
     vector<vector<vector<bool>>>runsskip; //dim1: num runs; dim2: pop; dim3: selective site
    /// January 22
    vector<size_t>testtype; // 1:Bartlett 2:Brown-Forsythe 3:F-test 4:Levene 5:LeveneHSK2 6:O'Brian


// **********
try{
    start = clock();

    // Set the current local date
    time_t rawtime;
    struct tm * timeinfo;

    time ( &rawtime );
    timeinfo = localtime ( &rawtime );

    cout<<asctime(timeinfo)<<endl;
    cout<<"------------------------------------------------\n"<<endl;
    // Store for the log file
    mensaje=asctime (timeinfo);
	mensaje+=
	"---------------------------------------------------------------------\n";
	mensaje+=" iHDSel: Program to identify patterns of natural selection using information-based statistics\n";
	mensaje+=" by A. Carvajal-Rodríguez			\n";
	mensaje+=" https://acraaj.webs.uvigo.es/iHDSel.html			\n";
	mensaje+="---------------------------------------------------------------------\n";


/// **** Arguments in the command line

    SNPSelTools::ParseArgs(argc, argv,
                tag,
                nogaps,
                store,
                old,
                path,
                format,
                onepopbyfile,
                inputname,
                inputname2,
				outheader,
                userfile,
                NumOfFiles,
                popforHB,
                selectives,
                fixwinsize,
                samplesize,
                samplesize2,
                minWinSize,MaxWinSize,iniMAF,SL,
                doEOS,soloEOS,EPOhighest,numdesvs,fixseed,useblocks,verbose, writersidInfo,onlyrsidinfo,ver);

    /// May 2024
    if(fixwinsize){ /// ONly if windows size was defined to be not 0
        if(fixwinsize>MaxWinSize || fixwinsize<minWinSize){

            string msg = "\nERROR: \nDefined window size of "+NtoS(fixwinsize)+ " is out of limits: ["+NtoS(minWinSize)+", "+NtoS(MaxWinSize)+"]\n";
            cout<<msg<<endl;
            SNPSelTools::texterr(useblocks+doEOS, msg,true);

        }
    }

    /// Enero 24
    if(old){DBT=2;} // if block is old then only check r and exlucde D (never can be greater than 1 so cannot be > DBT)

    ///Feb 24
    if(soloEOS){
            doEOS=true;
            string msg = "\nWARNING:\nsoloEOS is enabled; setting doEOS to true.\n";
            cout<<msg<<endl;
            SNPSelTools::texterr(useblocks+doEOS, msg,false);

    } /// if soloEOS is activated then activate doEOS


    ///May24 if not useblocks then we required doEOS
    if(!useblocks && !doEOS){
        string msg = "\nWARNING: \ndOption with outliers (without automatic blocks) is enabled; setting doEOS to true.\n";
		cout<<msg<<endl;
		SNPSelTools::texterr(useblocks+doEOS, msg,false);
    }

    //cout<<DBT<<endl; exit(1);

    /// acr June 21
    //bool usercandidate = false; // this will mark if the candidates are user-defined or not
    popMAF = iniMAF;
    /// Dec 21: User-defined selective sites
    if(!selectives.empty()){

        if(selectives[0]<0){ /// There is the -1 flag for no user defined selective sites. So this is no selectives
            vector<double>().swap(selectives);
        }
        else{ ///April 2024: User-defined selective sites
            if(useblocks && fixwinsize){ /// if fixwinsize there is not necessary to calculate blocks
                useblocks=0;
                mensaje+="\nUseblocks option desactivated when there is user defined window size.\n";
                cout<<"\nUseblocks option desactivated when there is user defined window size.\n"<<endl;
            }

            if(!useblocks && !fixwinsize){ /// There are selective sites and if no predefined window size it is required to useblcoks
                useblocks=1; /// If there are user candidates the compute the automatic blocks for them
                cout<<"\nUseblocks is mandatory when there are user defined sites. It was changed to 1.\n"<<endl;
                 mensaje+="\nUseblocks is mandatory when there are user defined sites. It was changed to 1.\n";
            }
        }
    }

    if(selectives.empty()){ ///April 2024: This is not an else to enter here if the selectives argument was -1 in the previous if
    ///April 2024:
        if(useblocks && fixwinsize){ /// if fixwinsize there is not necessary to calculate blocks
            useblocks=0;
            cout<<"\nUseblocks option desactivated when there is user defined window size.\n"<<endl;
            mensaje+="\nUseblocks option desactivated when there is user defined window size.\n";
        }
    }
    minblocksize=minWinSize;

    /// acr Oct 21 if desire all available snps as window size should be
    /// MaxWinSize=0 & fixwinsize=0
    /// Nov 21 if desire window size depending on linkage fixwinsize=0
    if(fixwinsize){if(MaxWinSize<=0){MaxWinSize=fixwinsize;}}


    /// Dec 21
    if(popforHB>=numpops){
		string msg = "\nERROR: \npopforHB index: "+ NtoS(popforHB)+" is out of range. It should be lower than the number of populations.\n";
		cout<<msg<<endl;
		SNPSelTools::texterr(popforHB, msg,true);

    }

    transform(format.begin(),format.end(), format.begin(), (int(*)(int)) toupper);

/// ************************ DATA FORMATS FLAGS *************************************

  // default format is ms
    if(format.find("FASTA")!=string::npos){
        msformat=false;
        bsformat=false;
        fastaformat=true;
        vcfformat=false;
        hapmapformat=false;
        fastPhaseformat=false;
        gpformat=false;
        plformat=false;

        if(!inputname2.empty()){onepopbyfile=true;}
    }
    else
    if(format.find("HM")!=string::npos){ // hapmap phased format
        msformat=false;
        bsformat=false;
        fastaformat=false;
        vcfformat=false;
        hapmapformat=true;
        fastPhaseformat=false;
        gpformat=false;
        plformat=false;
    }
    else
    if(format.find("FP")!=string::npos){ // fastphase phased format
        msformat=false;
        bsformat=false;
        fastaformat=false;
        vcfformat=false;
        hapmapformat=false;
        fastPhaseformat=true;
        gpformat=false;
        plformat=false;

        onepopbyfile=true;
    }
    else
    if(format.find("VCF")!=string::npos){ // VCF phased format
        msformat=false;
        bsformat=false;
        fastaformat=false;
        vcfformat=true;
        hapmapformat=false;
        fastPhaseformat=false;
        gpformat=false;
        plformat=false;

        onepopbyfile=true;
    }
    else
    if(format.find("BS")!=string::npos){
        msformat=false;
        bsformat=true;
        fastaformat=false;
        vcfformat=false;
        hapmapformat=false;
        fastPhaseformat=false;
        gpformat=false;
        plformat=false;

        soloEOS=true; // there are no haplotypes under this format
        if(!doEOS){std::cout << "Detected BayeScan format: setting doEOS = true." << std::endl;}
        doEOS=true;

    }
    else
    if(format.find("GP")!=string::npos){
        msformat=false;
        bsformat=false;
        fastaformat=false;
        vcfformat=false;
        hapmapformat=false;
        fastPhaseformat=false;
        gpformat=true;
        plformat=false;

        soloEOS=true; // there are no haplotypes under this format
        if(!doEOS){std::cout << "Detected Genepop format: setting doEOS = true." << std::endl;}
        doEOS=true;
    }

    else
    if(format.find("PL2")!=string::npos){
        msformat=false;
        bsformat=false;
        fastaformat=false;
        vcfformat=false;
        hapmapformat=false;
        fastPhaseformat=false;
        gpformat=false;
        plformat=true;
        plnumeric=true;
        NumOfFiles=1;
        soloEOS=true; // there are no haplotypes under this format
        if(!doEOS){std::cout << "Detected Plink format: setting doEOS = true." << std::endl;}
        doEOS=true;
    }
    else
    if(format.find("PL")!=string::npos){
        msformat=false;
        bsformat=false;
        fastaformat=false;
        vcfformat=false;
        hapmapformat=false;
        fastPhaseformat=false;
        gpformat=false;
        plformat=true;

        NumOfFiles=1;
        soloEOS=true; // there are no haplotypes under this format
        if(!doEOS){std::cout << "Detected Plink format: setting doEOS = true." << std::endl;}
        doEOS=true;
    }

    // if fixwinsize==0 we will set automatically the winsize, but if not 0 then check:
    if(fixwinsize && (fixwinsize<minWinSize)){ // default minWinSize is 11
		string msg = "\nERROR: \nfixwinsize is lower than minWinSize: "+ NtoS(minWinSize)+". You should define a higher value or change minWinSize\n";
		cout<<msg<<endl;
		SNPSelTools::texterr(fixwinsize, msg,true);
	}

	if(MaxWinSize>0 && MaxWinSize<minWinSize){
		string msg = "\nERROR: \nMaxWinSize is lower than minWinSize: "+ NtoS(minWinSize)+". You should define a higher value\n";
		cout<<msg<<endl;
		SNPSelTools::texterr(MaxWinSize, msg,true);
	}

    if(fixseed) {

        seed = defaultseed;
        mensaje+="\nWARNING: Fixed seed mode is on. Fixed seed is "+NtoS(seed);
    }
    else{
        mensaje+="\nInitial seed "+NtoS(seed)+"\n";
    }

   // Initialize seed in the engine
    gen.seed(seed);

    nummetapops=1; /// BY NOW ONLY 1 METAPOP ALLOWED
    numsnps.resize(nummetapops); // this depends on the metapop flag
    metapopsnpfreqs.resize(nummetapops);

/// **** Output directory
    string dirname="iHDSel_Results";
    portablemkdir(dirname);
	dir=dirname+"/"+outheader+"_";
/*
    if(samplesize2>0)
        inisamplesize=samplesize+samplesize2;
    else
        inisamplesize=samplesize;
*/
    if(onepopbyfile){
        userfile=true;
        NumOfFiles=1;
        if(inputname==inputname2){
            string msg = "\nERROR: \nEqual file names under onepopbyfile option\n";
            cout<<msg<<endl;
            SNPSelTools::texterr(0, msg,true);
        }
    }
    else inputname2=inputname;

    if(userfile && NumOfFiles>1)
        defaultFname=inputname;


/// **** Output file name

    string foutnamemax="",outname="";

    if(soloEOS){
        outname=dir+"EOS_";
    }

    indexesfiles.reserve(NumOfFiles); // the size could be lower is some run is skipped
    Dp_runs.reserve(NumOfFiles);
    fstfiles.reserve(NumOfFiles);
    fstfiles0.reserve(NumOfFiles);
    /// acr Oct 21
    //runsVarFTest.reserve(NumOfFiles);
    //runsn1.reserve(NumOfFiles);
    //runsn2.reserve(NumOfFiles);

    ///Dec 21
    runsskip.reserve(NumOfFiles);
    /// Jan 22
    testtype.reserve(numpops*NumOfFiles);
    /// Oct 2023
    vector<vector<int>>runsamplesizes(NumOfFiles,vector<int>(2,0));

//    vector<int>minsamplesizes;
//    minsamplesizes.reserve(NumOfFiles); // min sample size found at each analyzed file

    vector<bool> highmeanFST; //April 2024
    highmeanFST.reserve(NumOfFiles);

    ///May 2024
    bool liberal=false;

//    vector<int>numoutliers;
    vector<int> totsnps;
    //vector<int> numindepsnps;
//    numoutliers.reserve(NumOfFiles);
    totsnps.reserve(NumOfFiles);
    //numindepsnps.reserve(NumOfFiles);
    vector<int>contSLsig;
    contSLsig.reserve(NumOfFiles);
    vector<double>fstmax;
    fstmax.reserve(NumOfFiles);
    vector<double>fstmean;
    fstmean.reserve(NumOfFiles);

//    vector<double>lbqval;
//    lbqval.reserve(NumOfFiles);

    vector<HapMapTable> snpstable;
    if(hapmapformat) snpstable.resize(NumOfFiles);
    vector<string> bufferdata,outfilesbufdat;

    if(onlyrsidinfo || writersidInfo) {bufferdata.resize(NumOfFiles);outfilesbufdat.resize(NumOfFiles);}

    snpfreqs.resize(numpops); // This will store the snp frequencies at each population
    snpositions.resize(numpops);

    popreferences.resize(numpops); // acr July 21

    if(verbose){
        if(popforHB<0)
            cout<<"\nIf sample sizes are equal the reference for the analysis is computed from population "<<(Pref+1)<<endl;
        else
            cout<<"\nIf sample sizes are equal the reference for the analysis is computed from population "<<(popforHB+1)<<endl;
    }

/// ***************************************************** FOR EACH FILE

int refpopforHB=popforHB; /// May 2024 because if there are more than one file to analyze popforHB first time could be -1 (not user-defined) but second time is no longer -1

cout<<"\nOpening and parsing file(s) "<<endl;
for(int run=0;run<NumOfFiles;++run){

    popforHB=refpopforHB; ///May 2024 if it had the -1 value here that value will be restored
    //usercandidate = false;
    flag_therearecommons=true;

    if(fixwinsize) {winsize = fixwinsize; } // updated here because can change in some runs

    skiprun=false;

	if(!userfile) inputname =defaultFname;
	if((NumOfFiles>1))
        inputname =defaultFname +NtoS(run)+".txt";

/// **********************************************************************************************************************************
/// *************************************** Begin PARSING the input files *************************************************
/// **********************************************************************************************************************************

/// ********** READ DATA UNDER DIFFERENT FORMATS
 	deque < deque<bool> > data; // This stores the raw data
	/*
	// data will be a matrix of bools with variable row length (number of segregating sites) in different pops
	// and as much rows as numpops x sample size


	*/
    deque<deque<bool>> metapopdata; // Same structure as data. To manage each population data as a whole. Cleared and charged in datafusion

	vector<double> metapopsnpositions, commons;
    vector<STAT>fstdistrib;
    vector<string> filelocnames;

    /// Nov 21 haplotype blocks
    /**
    Block stores the beginning and end of the block and is defined in TLibHACSel as

    struct Block {
        int iniw=-1;
        int endw=0;
        Block(){ };
        Block(const int& iniw, const int& endw) : iniw(iniw), endw(endw){ } // constructor for initialization an anonymous object
    };

    */
    int numblocks=0;
    vector<Block> blocks;

    double fsttot{0},maxfst{1}, minHs{0}, maxHt{1};

///
    if(verbose) cout<<"... Identifying file format ..."<<endl;
    if(msformat){ // if !samplesize2 this ms format has the same sample size in every population
        if(samplesize2==0){
            if(!SNPSelTools::readmsdata(path,inputname,samplesize,data,snpositions) ) {

                if(NumOfFiles<2){
                    string msg = "\nERROR: \nFail when reading "+ inputname+" under MS format\n";
                    cout<<msg<<endl;
                    SNPSelTools::texterr(1,msg,true);
                }
                else {
                    string msg = "\nERROR: \nFail when reading "+ inputname+" under MS format\n";
                    cout<<msg<<endl;
                    mensaje+=msg;

                    continue;
                }
            }
            else {

                mensaje+= "\niHDivSel of file "+inputname+":\n";

                samplesizes.resize(numpops,samplesize);

                mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                if(verbose){
                    cout << "\nRead file "+inputname<<endl;
                    cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                }
                // This is only necessary when reading the ms format
                for(int p=0; p<numpops;++p){ ///May 2024 Compute freq of 1's the problem is that under the ms format the 1 is not necessarily the minor
                     computeSNPFreqinPop(data,p, samplesizes, snpfreqs[p]);
                }
            /// FST upper-bound from Carvajal-Rodriguez 2017 (Appendix)
    //            double samplesizequotient{0.0};
    //            for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);});
                // because under fasta format sample size is constant

                minHs = 2*(1-popMAF)*popMAF*samplesize/(samplesize-1.0);
                maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);



            }
        }
        else{ // samplesize2 >0: under this ms format sample sizes can be different
            if(!SNPSelTools::readmspopsdata(path,inputname,samplesize+samplesize2,samplesize,samplesize2,data,snpositions) ) {

                if(NumOfFiles<2){
                    string msg = "\nERROR: \nFail when reading "+ inputname+" under MS format\n";
                    cout<<msg<<endl;
                    SNPSelTools::texterr(1,msg,true);
                }
                else {
                    string msg = "\nERROR: \nFail when reading "+ inputname+" under MS format\n";
                    cout<<msg<<endl;
                    mensaje+=msg;

                    continue;
                }
            }
            else {

                mensaje+= "\niHDivSel of file "+inputname+":\n";
                samplesizes.resize(numpops,samplesize);
                samplesizes[1]=samplesize2;

                mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";

                if(verbose){
                    cout << "\nRead file "+inputname<<endl;
                    cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                }

                if(samplesize!=samplesize2) // choose as ref that having the highest sample size
                    //Pref = distance(samplesizes.begin(),max_element(samplesizes.begin(),samplesizes.end()));
                    Pref= samplesizes[0]>=samplesizes[1]? 0:1; ///May 2024

                // This is only necessary when reading the ms format
                for(int p=0; p<numpops;++p){
                    computeSNPFreqinPop(data,p, samplesizes, snpfreqs[p]); ///Compute freq of 1's
                }

            /// FST upper-bound from Carvajal-Rodriguez 2016 (Appendix)

            double samplesizequotient{0.0};
            for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
            minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
            maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);


            }
        }
    } // close ms format
    else
    if(bsformat){ // the sample size is user defined but some locus may have lower due to missing alleles, that will only affect the computation of frequency of the alleles in that loci
        vector<vector<vector<float>>> alfreqs; // loci, pops, allelic freqs
        vector<vector<float>>majoralfreqs;

        if(!SNPSelTools::readBSSNPdata(path,inputname,numpops, alfreqs,majoralfreqs,diploid, filelocnames)){

            if(NumOfFiles<2){
                string msg = "\nERROR: \nFail when reading "+ inputname+" under BS format\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(1,msg,true);
            }
            else {
                string msg = "\nERROR: \nFail when reading "+ inputname+" under BS format\n";
                cout<<msg<<endl;
                mensaje+=msg;

                continue;
            }
        }
        else {
            mensaje+= "\niHDivSel of file "+inputname+":\n";

            samplesizes.resize(numpops,samplesize);

            mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
            if(verbose){
                cout << "\nRead file "+inputname<<endl;
                cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
            }

            // there are no haplotypes under this format

// The loci are the same in every population but can be monomorphic at some population so delete sites with major allele freq above 1- popmaf
          bool deleteloc{false};
          auto nloc = alfreqs.size();
          auto npops = alfreqs[0].size();
          vector<string> templocnames;
          templocnames.reserve(nloc);
          vector<vector<double>> Hs(npops);
          vector<double> Ht;
          Ht.reserve(nloc);

          for(auto& h:Hs)
            h.reserve(nloc);
// cout<<nloc<<endl;

          vector<uint_fast32_t> numal(numpops,0);

          for(size_t l=0; l<nloc; ++l){ /// for each locus

            for(size_t p=0; p<npops; ++p){
                numal[p] = accumulate(begin(alfreqs[l][p]),end(alfreqs[l][p]),0); // num of copies of non-null alleles for this locus in pop p

//cout<<alfreqs[l][p][0]<<"\t"<<alfreqs[l][p][1]<<"\t"<<majoralfreqs[l][p]<<"\t"<<(int)numal[p]<<"\t"<<(1-popmaf)<<endl;

                if( (majoralfreqs[l][p]/numal[p] > (1-popMAF)) ){
                   deleteloc=true;


                    break;
                }
                else{
                    snpfreqs[p].push_back(1.0 - majoralfreqs[l][p]/numal[p]); ///For bayescan input format: store frequency of minor allele
                    snpositions[p].push_back(l);
                }
    //if(l==(nloc-2))
    //    cout<<locnames[run][l]<<"\t"<<p<<"\t"<<majoralfreqs[l][p]/numal[p]<<"\t"<<(1-popmaf)<<endl;

            }

//exit(1);
            if(!deleteloc){ // compute Hs and Ht for this locus
                templocnames.push_back(filelocnames[l]);
                float h2{0.0},ht2{0.0},total{0};
                // Hs
                for(size_t p=0; p<npops; ++p){

                    total+=numal[p]; // total num of non-null alleles for this locus in all populations
                    for(auto al:alfreqs[l][p] ){
                        al/=numal[p];
                        h2+= al*al;
                    }
                    Hs[p].push_back(1-h2);
                    h2=0;

                } // close for p
                // Ht
                // Count the frequency of each allele in the metapopulation
                size_t maxnumal{0}; // if for the locus l the number of alleles at each population is different then get the maximum
                for_each(begin(alfreqs[l]),end(alfreqs[l]), [&](vector<float> element){maxnumal= max(maxnumal,element.size());});
                vector<double> metapopfreqal(maxnumal,0);
                for(size_t p=0; p<npops; ++p){
                    for(size_t al=0; al < alfreqs[l][p].size(); ++al)

                        metapopfreqal[al]+= alfreqs[l][p][al]; // this sum the copies of a given allele in the metapop

                } // close for p
                for(auto& freqal: metapopfreqal){
                    freqal/=total; // this is the metapop allele freq for each allele
                    ht2+= freqal*freqal;
                }


                Ht.push_back(1-ht2);

            } // close if not deleteloc
            deleteloc=false;
          } // close for l
        /// uptdate the list of locus names

          vector<string>().swap(filelocnames);
          filelocnames= templocnames;
          vector<string>().swap(templocnames);


            for(size_t p=0; p<npops; ++p)
                assert(Hs[p].size() == Ht.size());
            assert(filelocnames.size() == Ht.size());

            numsnps[nummetapops-1]=Ht.size(); // by now just one metatpop

            metapopsnpositions.reserve(numsnps[nummetapops-1]);
    //cout<<numsnps[0]<<endl;
        /// Construct gstdistrib and compute gsttot (Hedrick chapter 9 p.491)

            fstdistrib.resize(Ht.size());
            for(size_t l=0; l< Ht.size(); ++l){

                metapopsnpositions.push_back((float)l/numsnps[0]) ; // by the way fullfil this vector with uniform positions

                float fst{0},hs{0};
                for(size_t p=0; p<npops; ++p)
                    hs+=Hs[p][l]; //
                hs/=npops;
                fst=max(1.0- hs/Ht[l],0.0); // avoid negative gst
                fstdistrib[l].stat=fst;
                fsttot+=fst;

    //if(l==2) cout<< fst<<"\t"<<hs<<"\t"<<Ht[l]<<endl;
            } // close for loci
            fsttot/=Ht.size();
        /// GST upper-bound from Carvajal-Rodriguez 2016 (Appendix)
            minHs = 1.0 - (1-popMAF)*(1-popMAF) - popMAF*popMAF;
            maxHt = minHs + ((numpops-1.0)/numpops)*((1.0-popMAF)*(1.0-popMAF) + popMAF*popMAF);
            maxfst = 1.0 -minHs/maxHt;


        }
    } // close bs format
    else
    if(fastaformat){
        //snpfreqs store the frequency of 1's which are the minor alleles in the reference pop (and could be the minor or major in the other)
        if(verbose){
            cout<<"Fasta file"<<endl;
        }
        if(onepopbyfile){
            //samplesizes.resize(numpops);
            if(!SNPSelTools::readfastadata(nogaps,path,inputname,inputname2, samplesizes,data, snpfreqs,snpositions,Pref)){ //snpfreqs store the frequency of 1's

                if(NumOfFiles<2){
                    string msg = "\nERROR: \nFail when reading "+ inputname+" under fasta format\n";
                    cout<<msg<<endl;
                    SNPSelTools::texterr(1,msg,true);
                }
                else {
                    string msg = "\nERROR: \nFail when reading "+ inputname+" under fasta format\n";
                    cout<<msg<<endl;
                    mensaje+=msg;

                    continue;
                }
            }
            else {
                mensaje+= "\niHDivSel of file "+inputname+":\n";
                mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                if(verbose){
                    cout << "\nRead file "+inputname<<endl;
                    cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                }

       /// FST upper-bound from Carvajal-Rodriguez 2016 (Appendix pag 17)
            double samplesizequotient{0.0};
            for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
            minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
            maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);
            }

        } // close if onepopbyfile
        else{
            if(!tag.empty()){
                tag=ToLowerCase(tag);

                if(tag=="rnd0"){
                    if(!SNPSelTools::readfastadatarand(store,nogaps,path,inputname, samplesizes,data, snpfreqs,snpositions,Pref)){ //snpfreqs store the frequency of 1's

                        if(NumOfFiles<2){
                            string msg = "\nERROR: \nFail when reading "+ inputname+" under fasta format\n";
                            cout<<msg<<endl;
                            SNPSelTools::texterr(1,msg,true);
                        }
                        else {
                            string msg = "\nERROR: \nFail when reading "+ inputname+" under fasta format\n";
                            cout<<msg<<endl;
                            mensaje+=msg;

                            continue;
                        }
                    }
                    else {
                        /// May 2024 comment this that at best do nothing and on the contrary is wrong:
                        ///samplesizes.resize(numpops,samplesize);

                        mensaje+= "\niHDivSel of randomized file "+inputname+":\n";
                        mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                        if(verbose){
                            cout << "\nRead file "+inputname<<endl;
                            cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                        }

       /// FST upper-bound from Carvajal-Rodriguez 2016 (Appendix pag 17)
                        double samplesizequotient{0.0};
                        for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
                        minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
                        maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);
                    }
                }
                else{ //tag is normal not for randomization
                    if(!SNPSelTools::readfastadatatag(tag,store,nogaps,path,inputname, samplesizes,data, snpfreqs,snpositions,Pref)){ //snpfreqs store the frequency of 1's

                        if(NumOfFiles<2){
                            string msg = "\nERROR: \nFail when reading "+ inputname+" under fasta format\n";
                            cout<<msg<<endl;
                            SNPSelTools::texterr(1,msg,true);
                        }
                        else {
                            string msg = "\nERROR: \nFail when reading "+ inputname+" under fasta format\n";
                            cout<<msg<<endl;
                            mensaje+=msg;

                            continue;
                        }
                    }
                    else {

                        mensaje+= "\niHDivSel of file "+inputname+":\n";
                        mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                        if(verbose){
                            cout << "\nRead file "+inputname<<endl;
                            cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                        }

                        double samplesizequotient{0.0};
                        for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
                        minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
                        maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);

                    }

                } //close else so tag is not rnd0

            } // close tagged pop
            else{ // pop is not tagged so the sequences are assumed to be the first samplesize ones from one population and the following samplesize ones to the other
                if(!SNPSelTools::readfastadata(nogaps,path,inputname, samplesizes,data, snpfreqs,snpositions)){ //snpfreqs store the frequency of 1's

                    if(NumOfFiles<2){
                        string msg = "\nERROR:\nFail when reading "+ inputname+" under fasta format\n";
                        cout<<msg<<endl;
                        SNPSelTools::texterr(1,msg,true);
                    }
                    else {
                        string msg = "\nERROR:\nFail when reading "+ inputname+" under fasta format\n";
                        cout<<msg<<endl;
                        mensaje+=msg;

                        continue;
                    }
                }
                else {


                    mensaje+= "\niHDivSel of file "+inputname+":\n";
                    mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                    if(verbose){
                        cout << "\nRead file "+inputname<<endl;
                        cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                    }

                    double samplesizequotient{0.0};
                    for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
                    minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
                    maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);

                }

            } // pop is not tagged


        } // close else to if onepopbyfile
    }//close if fasta format
    else
    if(vcfformat){// for this format snpfreqs store the frequency of the non-reference allele which are noted as 1's
        //if(!SNPSelTools::readphasedVCF(path,inputname,inputname2, samplesizes,data, snpfreqs,snpositions, Pref))
        auto t0 = chrono::high_resolution_clock::now();
        if(verbose){
            cout<<"Reading vcf files"<<endl;


        }
        if(!SNPSelTools::readphasedVCF_intersec(path,inputname,inputname2, samplesizes,data, snpfreqs,snpositions, Pref)){
            if(NumOfFiles<2){
                string msg = "\nERROR: \nFail when reading "+ inputname+ " or "+inputname2 + " under VCF phased format\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(1,msg,true);
            }
            else {
                string msg = "\nERROR: \nFail when reading "+ inputname+ " or "+inputname2+" under VCF phased format\n";
                cout<<msg<<endl;
                mensaje+=msg;

                continue;
            }
        }
        else {

            mensaje+= "\niHDivSel of files "+inputname+" and "+inputname2+ "\n";
            mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
            if(verbose){
                cout<< "\nRead files "+inputname+" and "+inputname2<<endl;
                cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                auto t1 = chrono::high_resolution_clock::now();
                cout << "Time required: "
                     << chrono::duration<double>(t1 - t0).count()
                     << " s" << endl;
            }

            if(NumOfFiles<2)
                defaultFname=inputname;
            inputname += " vs "+inputname2;


   /// FST upper-bound from Carvajal-Rodriguez 2016 (Appendix)
             double samplesizequotient{0.0};
            for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
            minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
            maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);
        }
    }
    else
    if(hapmapformat){ // hapmap phased format: it allows different sample sizes at different populations
        // if sample sizes are different pref is fixed to the pop with the highest sample size
        //snpfreqs store the frequency of 1's which are the minor alleles in the reference pop (and could be the minor or major in the other)
        if(!SNPSelTools::readhmphased(path,inputname,inputname2,onepopbyfile, samplesizes,data, snpfreqs,snpositions,snpstable[run],Pref)){

            if(NumOfFiles<2){
                string msg = "\nERROR: \nFail when reading "+ inputname+" under hapmap phased format\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(1,msg,true);
            }
            else {
                string msg = "\nERROR: \nFail when reading "+ inputname+" under hapmap phased format\n";
                cout<<msg<<endl;
                mensaje+=msg;

                continue;
            }
        }
        else {
            if(onepopbyfile){
                mensaje+= "\niHDivSel of files "+inputname+" and "+inputname2+ "\n";
                mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                if(verbose){
                    cout<< "\nRead files "+inputname+" and "+inputname2<<endl;
                    cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                }

                if(NumOfFiles<2)
                    defaultFname=inputname;
                inputname += " vs "+inputname2;
            }
            else{
                mensaje+= "\niHDivSel of file "+inputname+":\n";
                mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                if(verbose){
                    cout << "\nRead file "+inputname<<endl;
                    cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
                }
            }

   /// FST upper-bound from Carvajal-Rodriguez 2016 (Appendix)
             double samplesizequotient{0.0};
            for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note 1.0 is key for obtaining a non-integer division that in other case would be always 1
            minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
            maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);
//cout<<mensaje<<endl;
//for(auto& dat:data)
//printcol(dat);
//
//
//  for ( auto& rsid:snpstable[run])
//    std::cout << " " << rsid.second << "\t" << snpstable[run][rsid.first]<<endl;


        }
    }
    else
    if(fastPhaseformat){
        bool firsthaplo=true;
        //snpfreqs store the frequency of 1's which are the minor alleles in the reference pop (and could be the minor or major in the other)
        if(!SNPSelTools::readfastPHASE(path,inputname,inputname2,samplesizes,data, snpfreqs,snpositions, Pref, firsthaplo)){

            if(NumOfFiles<2){
                string msg = "\nERROR: \nFail when reading "+ inputname+" under fastPhase format\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(1,msg,true);
            }
            else {
                string msg = "\nERROR: \nFail when reading "+ inputname+" under fastPhase format\n";
                cout<<msg<<endl;
                mensaje+=msg;

                continue;
            }
        }
        else {

            mensaje+= "\niHDivSel of files "+inputname+" and "+inputname2+ "\n";
            mensaje+= "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
            if(verbose){
                cout<< "\nRead files "+inputname+" and "+inputname2<<endl;
                cout<< "\nSample sizes are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
            }

            if(NumOfFiles<2)
                defaultFname=inputname;
            inputname += " vs "+inputname2;

   /// FST upper-bound from Carvajal-Rodriguez 2016 (Appendix)

            double samplesizequotient{0.0};
            for_each(samplesizes.begin(),samplesizes.end(), [&](const int& element){samplesizequotient+= element/(element-1.0);}); // Note here float 1.0 is key for obtaining a non-integer division that in other case would be always 1
            minHs = 2*(1-popMAF)*popMAF*samplesizequotient/(numpops);
            maxfst = 1.0 - numpops*minHs/(minHs + numpops -1.0);

        }

    }
    else
    if(gpformat){ // it allows different sample sizes at different populations

        vector<vector<vector<float>>> alfreqs; // loci, pops, allelic freqs
        vector<vector<float>>majoralfreqs;
        vector<vector<vector<uint_fast32_t>>> alelids;
        vector<uint_fast32_t> uniquelist;

        if(!SNPSelTools::readgenpopdata(path,inputname, samplesizes,alfreqs,majoralfreqs,alelids, diploid,filelocnames)){
            string msg = "\nERROR: \nFail when reading "+ inputname+" under genepop format\n";
            cout<<msg<<endl;
            SNPSelTools::texterr(1,msg,true);
        }
        else{
            mensaje+= "\nSample sizes in " +inputname+ " are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
            if(verbose)
                cout<<"\nSample sizes in " +inputname+ " are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
        }

      // The loci are the same in every population but can be monomorphic at some population so delete sites with major allele freq above 1- popmaf
      bool deleteloc{false};
      auto nloc = alfreqs.size();
      auto npops = alfreqs[0].size();
      vector<string> templocnames;
      templocnames.reserve(nloc);
      vector<vector<double>> Hs(npops);
      vector<double> Ht;
      Ht.reserve(nloc);

      for(auto& h:Hs)
        h.reserve(nloc);

 //cout<<nloc<<endl;
      vector<uint_fast32_t> numal(numpops,0);
      for(size_t l=0; l<nloc; ++l){
if(verbose)
    cout<<"\nLOCUS "<<l+1<<endl;
        for(size_t p=0; p<npops; ++p){
            numal[p] = accumulate(begin(alfreqs[l][p]),end(alfreqs[l][p]),0); // num of non-null alleles for this locus in pop p

            if( (majoralfreqs[l][p]/numal[p] > (1-popMAF)) || majoralfreqs[l][p] == MISSINGVAL ){
               deleteloc=true;

                break;
            }
            else{
                snpositions[p].push_back(l);
            }

        } //close for p

        if(!deleteloc){ // compute Hs and Ht for this locus
            templocnames.push_back(filelocnames[l]);
            float h2{0.0},tal{0},ht2{0.0},total{0};
            // Hs
            for(size_t p=0; p<npops; ++p){

                total+=numal[p]; // num of non-null alleles for this locus in all populations
                for(auto al:alfreqs[l][p] ){
                    al/=numal[p];
                    h2+= al*al;
                }
                Hs[p].push_back(1.0-h2);
                h2=0;
            } // close for p
            // Ht
            do_unique_list(alelids[l], size_t(0), uniquelist);
            // Count the frequency of each allele in the metapopulation
            for(auto al: uniquelist){
                for(size_t p=0; p<npops; ++p){
                    // get the position of al in the given population
                    auto iterpos=find(alelids[l][p].begin(),alelids[l][p].end(), al);
                    if(iterpos!=alelids[l][p].end())

                        tal+= alfreqs[l][p][distance(alelids[l][p].begin(),iterpos)]; // this sum the copies of a given allele in the metapop

                } // close for p
                tal/=total; // this is the metapop allele freq
                ht2+= tal*tal;
                tal=0;
            } // close for each allele in uniquelist

            Ht.push_back(1-ht2);
            if(verbose){
                cout<<"read."<<endl;
            }

        } // close if not deleteloc
        else{

            if(verbose){

                cout<<"deleted."<<endl;
            }

        }
        deleteloc=false;
      } // close for l
    /// uptdate the list of locus names
      vector<string>().swap(filelocnames);
      filelocnames= templocnames;
      vector<string>().swap(templocnames);

//printcol(filelocnames);

        for(size_t p=0; p<npops; ++p)
            assert(Hs[p].size() == Ht.size());

        assert(filelocnames.size() == Ht.size());

        numsnps[nummetapops-1]=Ht.size(); // by now just one metatpop

        metapopsnpositions.reserve(numsnps[nummetapops-1]);
//cout<<numsnps[0]<<endl;
    /// Construct gstdistrib and compute gsttot (Hedrick chapter 9 p.491)

        fstdistrib.resize(Ht.size());
        for(size_t l=0; l< Ht.size(); ++l){

            metapopsnpositions.push_back((float)l/numsnps[0]) ; // by the way fullfil this vector with uniform positions

            float fst{0},hs{0};
            for(size_t p=0; p<npops; ++p)
                hs+=Hs[p][l]; //
            hs/=npops;
            fst=max(1.0- hs/Ht[l],0.0); // avoid negative gst
            fstdistrib[l].stat=fst;
            fsttot+=fst;

//if(l==2) cout<< fst<<"\t"<<hs<<"\t"<<Ht[l]<<endl;
        } // close for loci
//printcol(fstdistrib);
        fsttot/=Ht.size();

    /// GST upper-bound from Carvajal-Rodriguez 2016 (Appendix)
        minHs = 1.0 - (1-popMAF)*(1-popMAF) - popMAF*popMAF;
        maxHt = minHs + ((numpops-1.0)/numpops)*((1.0-popMAF)*(1.0-popMAF) + popMAF*popMAF);
        maxfst = 1.0 -minHs/maxHt;

    } // close if gpformat
    else
    if(plformat){

        inputname2+=".ped";
        inputname+=".map";

        vector<vector<vector<float>>> alfreqs; // loci, pops, allelic freqs
        vector<vector<float>>majoralfreqs;
        vector<vector<vector<uint_fast32_t>>> alelids;
        vector<uint_fast32_t> uniquelist;

        uint_fast32_t nullval='N';

//string msg = "\n By 28 June 2016- SORRY: PLINK format should be ready in a couple of days.\n";
//cout<<msg<<endl;
//SNPSelTools::texterr(99, msg,true);
        if(plnumeric){
            if(!SNPSelTools::readplinkdata2(path,inputname,inputname2, numpops,  snpositions,  samplesizes,alfreqs,majoralfreqs, alelids, diploid, filelocnames)){
                string msg = "\nERROR: \nFail when reading "+ inputname2+" under plink numeric format\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(1,msg,true);
            }
            nullval=0;
        }
        else
        if(!SNPSelTools::readplinkdata(path,inputname,inputname2, numpops,  snpositions,  samplesizes,alfreqs,majoralfreqs, alelids, diploid, filelocnames)){
                string msg = "\nERROR: \nFail when reading "+ inputname2+" under plink nucleotide format\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(1,msg,true);
        }

        mensaje+= "\niHDivSel of files "+inputname+" and "+inputname2+ ":\n";
        mensaje+= "\nSample sizes in " +inputname2+ " are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
        if(verbose){
                cout<< "\nRead files "+inputname+" and "+inputname2<<endl;
                cout<< "\nSample sizes in " +inputname2+ " are "+NtoS(samplesizes[0])+" and "+NtoS(samplesizes[1])+ "\n";
        }
        inputname += " and "+inputname2;

    // The loci are the same in every population but can be monomorphci at some or both populations these sites are deleted the same as the sites with major allele freq above 1- popmaf
      bool deleteloc{false};
      auto nloc = alfreqs.size();
      auto npops = alfreqs[0].size();
assert((int)npops==numpops);
      vector<string> templocnames;
      templocnames.reserve(nloc);
      vector<vector<double>> Hs(npops);
      vector<double> Ht;
      Ht.reserve(nloc);

      for(auto& h:Hs)
        h.reserve(nloc);
//cout<<nloc<<endl;
      vector<uint_fast32_t> numal(numpops,0);

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

        for(size_t p=0; p<npops; ++p){
            numal[p] = accumulate(begin(alfreqs[l][p]),end(alfreqs[l][p]),0); // num of non-null alleles for this locus in pop p

            if( (majoralfreqs[l][p]/numal[p] > (1-popMAF)) || majoralfreqs[l][p] == MISSINGVAL ){
               deleteloc=true;
//cout<<l<<"\t"<<p<<"\t"<<majoralfreqs[l][p]/numal[p]<<"\t"<<(1-popmaf)<<endl;
                break;
            }
            else{
                snpfreqs[p].push_back(1.0 - majoralfreqs[l][p]/numal[p]); //the frequency of the minor allele
            }
//if(l==(nloc-2))
//    cout<<locnames[run][l]<<"\t"<<p<<"\t"<<majoralfreqs[l][p]/numal[p]<<"\t"<<(1-popmaf)<<endl;
        }

        if(!deleteloc){ // compute Hs and Ht for this locus
            templocnames.push_back(filelocnames[l]);
            metapopsnpositions.push_back(snpositions[0][l]); //snpositions[0] and snpositions[1] are the same
            float h2{0.0},tal{0},ht2{0.0},total{0};
            // Hs
            for(size_t p=0; p<npops; ++p){

                total+=numal[p]; // num of non-null alleles for this locus in all populations
                for(auto al:alfreqs[l][p] ){ // pass by value because we continue to use the absolute frequency later
                    al/=numal[p];
                    h2+= al*al;
                }
                Hs[p].push_back(1-h2);
                h2=0;

            } // close for p
            // Ht
            do_unique_list<uint_fast32_t>(alelids[l], nullval, uniquelist);
    // Count the frequency of each allele in the metapopulation
            for(auto al: uniquelist){
                for(size_t p=0; p<npops; ++p){
                // get the position of al in the given population
                    auto iterpos=find(alelids[l][p].begin(),alelids[l][p].end(), al);
                    if(iterpos!=alelids[l][p].end())
                        tal+= alfreqs[l][p][distance(alelids[l][p].begin(),iterpos)]; // this sum the number of copies of a given allele in the metapop

                } // close for p
                tal/=total; // this is the metapop allele freq
                ht2+= tal*tal;
                tal=0;
            } // close for each allele in uniquelist

            Ht.push_back(1-ht2);

        } // close if not deleteloc
        deleteloc=false;
      } // close for l
    /// uptdate the list of locus names

      vector<string>().swap(filelocnames);
      filelocnames= templocnames;
      vector<string>().swap(templocnames);

for(size_t p=0; p<npops; ++p)
    assert(Hs[p].size() == Ht.size());

assert(filelocnames.size() == Ht.size());

        numsnps[nummetapops-1]=Ht.size(); // by now just one metatpop
assert(filelocnames.size() == metapopsnpositions.size());
//cout<<numsnps[0]<<endl;
    /// Construct gstdistrib and compute gsttot (Hedrick chapter 9 p.491)

        fstdistrib.resize(Ht.size());
        for(size_t l=0; l< Ht.size(); ++l){

            float fst{0},hs{0};
            for(size_t p=0; p<npops; ++p)
                hs+=Hs[p][l]; //
            hs/=npops;
            fst=max(1.0- hs/Ht[l],0.0); // avoid negative gst
            fstdistrib[l].stat=fst;
            fsttot+=fst;

//if(l==2) cout<< fst<<"\t"<<hs<<"\t"<<Ht[l]<<endl;
        } // close for loci
//cout<<fsttot<<" "<<Ht.size()<<endl;
//printcol(fstdistrib);
        fsttot/=Ht.size();


    /// GST upper-bound from Carvajal-Rodriguez 2016 (Appendix)
        minHs = 1.0 - (1-popMAF)*(1-popMAF) - popMAF*popMAF;
        maxHt = minHs + ((numpops-1.0)/numpops)*((1.0-popMAF)*(1.0-popMAF) + popMAF*popMAF);
        maxfst = 1.0 -minHs/maxHt;
//printcol(locnames[run]);

    } // close pl format

/// *******************************************************************************************************************
/// *************************************** END OF THE PART FOR READING FORMATS **************************************** ///
/// *******************************************************************************************************************
    size_t shorternumsnps=infinito<size_t>();
    size_t largernumsnps=snpositions[0].size();


    cout<<"\nFile "<<inputname<<" parsed"<<endl;

    if(!soloEOS){
        if(((size_t)Pref >= snpositions.size()) ){
            string msg = "\nERROR: In file" +inputname+ "\nReference population value is out of range "+ NtoS(Pref+1)+" > "+ NtoS(snpositions.size()) +"\n";
            cout<<msg<<endl;
            SNPSelTools::texterr(2,msg,true);
        }

        if((size_t)numpops != snpositions.size()){
            string msg = "\nERROR: In file" +inputname+ "\nThe number of populations "+ NtoS(snpositions.size())+" does not coincide with the expected "+ NtoS(numpops) +"\n";
            cout<<msg<<endl;
            SNPSelTools::texterr(3,msg,true);
        }

    // Compute the minimum number of SNPs which is the largest possible number of common snps in any population sample
        for_each(snpositions.begin(),snpositions.end(), [&](vector<double> element){shorternumsnps= min(shorternumsnps,element.size());});
        commons.reserve(shorternumsnps);
        if(largernumsnps==shorternumsnps)
            largernumsnps=snpositions[1].size();

    }

/// **** Manage population structure

    if(verbose) cout<<"Managing population structure"<<endl;


    int* iniwin = new int[numpops];
    int* endwin = new int[numpops];

/// *******************************************************************************************************************
/// *************************************** ELIMINATE LOW FREQUENCY DATA **************************************** ///
/// *******************************************************************************************************************
/**
SNP frequency is required to be in the interval [1-maf, maf] in the metapopulation.
CHANGES from October 2023

NEW for the iHDivSel version:
    The onlycommonsSNP variable was eliminated

    1.- AS before only common SNPs are considered (sites present in data from both populations)
    2.- But for each population now fixed sites are allowed
    3.- Sites almost fixed for the metapopulation will be eliminated:
        freq allele <maf in pop 1 and 2 or freq allele > 1-maf in pop 1 and 2

**/
    /// Dec 21
    bool blocksdetected = (fixwinsize==0)&&(MaxWinSize>0);
    bool seldetection=true;
    size_t genomesize=0;

    /// Oct 2023. Here to avoid cross initiallizaton with the goto
    vector<ResultIndexes> indexes;
    bool haploformat = (msformat || fastaformat || hapmapformat || fastPhaseformat) && (metapopsnpfreqs[0].size()>1) && (numsnps[0]>=minWinSize) ;
    //double maxMeff{static_cast<double>( numsnps[0])}; // for using as total number of independent test in the multitest correction

    if(!gpformat && !bsformat && !plformat){ /// ***** if phased haplotype formats (not genpop, bayescan or plink)
            /// this should be the same that !soloEOS ? not because we can have fastaformat and the user set soloEOS=True so we do not get in here in this case but we want to get in even if soloEOS=True

        if(verbose) cout<< "Management of low frequency polymorphisms."<<endl;

    //
    //October 2023: Mark for elimination SNPs that do not reach a minimum frequency in the metapopulation
    // so in commons are stored any population SNP that has a frequency > maf in the sum of both populations

            copycommonwithMAF(samplesizes,snpositions,snpfreqs, commons, popMAF);

            // Check that there are at least 1 common snp

            flag_therearecommons = commons.size();
            if(!flag_therearecommons){
                mensaje+="\n From a total of "+ NtoS(snpositions[1].size());
                mensaje+="\nThere are not common SNPs at frequency >= "+NtoS(popMAF)+ "in the file "+ inputname +"\n";
                if(verbose){
                    cout<<"From "+NtoS(snpositions[0].size()) +" SNPs"<<endl;
                    cout<<"\nThere are not common SNPs at frequency >= "+NtoS(popMAF)+ "in the file "+ inputname +"\n";
                    cout<<"The total number of SNPs is "+NtoS(shorternumsnps)<<endl;
                    cout<<"The number of commons SNPs with a frequency higher than maf is "<<NtoS(commons.size())<<endl;
                    string msg="A number of "+NtoS(shorternumsnps-commons.size()) +" SNPs at frequency lower than "+NtoS(popMAF) +" has been discarded:"+"\n";
                    cout<<msg;

                    cout<<"This run will be skipped"<<endl;
                }

            /// New v1.6 and October 2023
                mensaje+="A number of "+NtoS(shorternumsnps-commons.size()) +" SNPs at frequency lower than "+NtoS(popMAF) +" has been discarded:"+"\n";
                skiprun=true;
                /// Enero 2024 descuenta esta réplica
                --totruns;


                goto NOSNPS;
            }
            else{

                auto t0 = chrono::high_resolution_clock::now(); /// 28 abril 2026
                if(verbose){
                    cout<<NtoS(commons.size())<< " common SNPs passed. Deleting non-common ones..."<<endl;


                }
                /// Enero 2024 cuenta esta réplica
                ++totruns;
                for(int p=0; p<numpops;++p){

                    deleteNonCommonData(data,p, samplesizes, snpfreqs[p], snpositions[p], commons); // paralelizable ?

                    assert(snpfreqs[p].size()==snpositions[p].size());

                    if(verbose){cout<<"Population "<<p<<" filtered"<<endl;}

                }

                if(verbose){

                    cout<<"The total number of SNPs is "+NtoS(shorternumsnps)<<endl;
                    cout<<"The number of SNPs with a frequency higher than maf is "<<NtoS(commons.size())<<endl;
                    string msg="A number of "+NtoS(shorternumsnps-commons.size()) +" SNPs at frequency lower than "+NtoS(popMAF) +" has been discarded:"+"\n";
                    cout<<msg;
                    mensaje+=msg;
                    auto t1 = chrono::high_resolution_clock::now();
                    cout << "Time required: "
                         << chrono::duration<double>(t1 - t0).count()
                         << " s" << endl;

                }
                /// Oct 2023

                assert(snpositions[0].size()==snpositions[1].size());

            }// close else



        // this crashes if positions are repeated
        //cout<<snpositions[0].size()<<" "<<commons.size()<<endl; exit(7);
        if(snpositions[0].size()!= commons.size()){
            int repeatedpos=-1;
            int p1 = count_distinct_els(snpositions[0], MISSINGVAL,repeatedpos);
            int p2 = count_distinct_els(snpositions[1], MISSINGVAL,repeatedpos);
            int repeatedels = max(snpfreqs[0].size() - p1, snpfreqs[1].size() - p2);
            if(repeatedels){
                string msg = "\nERROR:" +NtoS(snpositions[0].size())+ " is different from "+ NtoS(commons.size())+". This may be caused by " + NtoS(repeatedels)+ " repeated positions in your files. Please check position "+NtoS(repeatedpos) +"\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(4,msg,true);
            }
            else{
                string msg = "\nERROR:" +NtoS(snpositions[0].size())+ "positions is different number that "+ NtoS(commons.size())+" common SNPs. Please check your files." +"\n";
                cout<<msg<<endl;
                SNPSelTools::texterr(5,msg,true);
            }

        }
        //assert(snpositions[0].size()==commons.size());

        if(verbose) cout<< "\n... Updating metapopulation information ...\n"<<endl;


    //Original data is copied in metapopdata. The new data will store the data of all populations as if they were one.

        dataFusion(data, metapopdata, snpositions, metapopsnpositions, samplesizes,metapopsamplesize);



    assert((int)metapopdata.size()==metapopsamplesize);
    assert(metapopsnpositions.size() == metapopdata[0].size());
 ///Oct 2023
    assert(metapopsnpositions.size() == snpositions[0].size());

    assert(commons.size()==metapopsnpositions.size());

        genomesize = metapopsnpositions.size();

        if(hapmapformat){

            for(size_t l=0; l< commons.size(); ++l)
                filelocnames.push_back(snpstable[run][commons[l]]);
        }
        else{


            for(size_t l=0; l< commons.size(); ++l)
                filelocnames.push_back(NtoS(metapopsnpositions[l],outputprecision));


        }

        for(int mp=0; mp<nummetapops;++mp){
            metapopsamplesize/=nummetapops; // nummetapos is 1 by now
            //computeSNPFreqinPop(metapopdata,mp, metapopsamplesize, metapopsnpfreqs[mp]);
            computemeanSNPFreqinMetaPop(metapopdata, metapopsnpfreqs[mp]);
        }

        for(int p=0;p<nummetapops;++p){
            numsnps[p] = metapopsnpositions.size();
            if(numsnps[p]>=1E4) // if there are more than 10000 snps reduce the number of iterations (by default to 100)
                resampnumiter=rsmp_lowiter;
            else
            if(numsnps[p]>=5E3)
                resampnumiter=rsmp_mediumiter; // if there are more than 5000 snps reduce the number of iterations (by def to 500)
            }
            fstdistrib.resize(numsnps[0]);
    } // close if not gpformat && not bsformat && not plformat
    else /// ***** if non-phased formats i.e. genpop or bayescan or plink
        assert(numsnps[0] == (int)filelocnames.size());



if(writersidInfo || onlyrsidinfo){

    vector<string> Popnames(numpops);
    if(inputname==inputname2){
        Popnames[0]="1";
        Popnames[1]="2";
    }
    else{
        Popnames[0]=defaultFname;
        Popnames[1]=inputname2;
    }
    bufferAnalyzedSNPs(Popnames, filelocnames, snpositions,snpfreqs, bufferdata[run]);

    if(onlyrsidinfo){
        if(verbose){
            SNPSelTools::texterr(0,mensaje,false);
            cout<< "\nThe SNPs information table for files"<<Popnames[0]<<" and "<<Popnames[1] <<" was done"<<endl;
        }
        if(NumOfFiles==1){
            outfilesbufdat[0]=dir+".txt";
            writeAnalyzedSNPs(outfilesbufdat, bufferdata);
        }
        return 0;
    }

}

/// **** PERFORM ANALYSIS

    runsamplesizes[run]=samplesizes;


	mensaje+="Number of pops is "+NtoS(numpops);
	//double maxMeff{static_cast<double>( numsnps[0])}; // for using as total number of independent test in the multitest correction

/// Haplotype formats
    if(verbose) cout<< "... Preparing iHDivSel analysis ..."<<endl;

	//bool haploformat = (msformat || fastaformat || hapmapformat || fastPhaseformat) && (metapopsnpfreqs[0].size()>1) && (numsnps[0]>=minWinSize) ;

    if(!haploformat && metapopsnpfreqs[0].size()<=1) {
            skiprun=true; /// Enero 2024 descuenta esta réplica
            --totruns;
    }


  /// Oct 2023. EOS is performed before to get putative candidates for the haplotype analysis

/// ****************************************************************************************************************************
/// ************************************* EOS **********************************************************************************
/// ****************************************************************************************************************************

/// **** PERFORM Fst COMPUTATIONS if there are SNPs to do so

    if(numsnps[0]>=minWinSize){ // if there are enough common SNPs

    vector<ResultFst> eposet;
    vector<ResultFst> outliers;
    if(doEOS){ //Feb 2024

            if(verbose) cout<< "\n... Beginning EOS analysis ...\n"<<endl;

            //numindepsnps.push_back(maxMeff);
            double maxdev{((1.0/3)*maxfst)}; // maximum standard deviation given an asymmetric unimodal (DharmadhikariandKhumar89)
            //double varupb{(maxdev*maxdev)}; // The upper variance bound DharmadhikariandKhumar89 note that the range [a,b] is [0,maxfsst]

            ResultFst tempfst;

    /// FST tot
            if(!gpformat && !bsformat && !plformat){ // because in those formats fstdistrib is already computed
    /// fst distribution
                FSTdistribandTot(metapopdata, samplesizes, 0, metapopdata[0].size(),fstdistrib,fsttot);
                if(verbose) cout<< "\nComputed Fst distribution (n= "<<fstdistrib.size()<<"). Fst mean = "<<fsttot<<endl<<endl;
    assert(metapopdata[0].size()==fstdistrib.size());
    //April 2024
    assert(genomesize==fstdistrib.size());
            }


    /// Get the whole distribution of pvals for computing the q-values of the EOS set
            vector<double> pvals(fstdistrib.size(),1);
            double stat{0};
            int contSLsigs{0};

        /// ********************* LK test
   ///April 2024 check fst distribution is not biased towards 1 because the populations are too much differentiated
            if(fsttot<=0.5){
                ///April 2024
                highmeanFST.push_back(false);

                for(size_t ifst=0;ifst<fstdistrib.size(); ++ifst){

                    stat = (numpops-1)*fstdistrib[ifst].stat;

                    if(fsttot>0)
                        stat/=fsttot;

                    if(stat>0)
                        pvals[ifst] = getChiPval(stat,(numpops-1));

                    fstdistrib[ifst].pval=pvals[ifst];

                    if(pvals[ifst]==1){
                        pvals[ifst] = quasione;

                    }
                    else
                    if(pvals[ifst] <= SL)
                        ++contSLsigs;


                }
                if(verbose) cout<< "\nPerformed Lk test. #* = "<<contSLsigs<<endl;

        /// get Positive Outliers:
        ///**NOTE** that fstdistrib, pvals and metapopsnpositions must be correlative (e.g. we cannot sort pvals without the others)
        /// store the positions of the outliers also in the vector candidates for non consecutive outliers

                epos(fstdistrib, metapopsnpositions,filelocnames, 3, eposet); /// Compute positive outliers > Q3 + 3IQR

                size_t numoutliers=eposet.size();  ///Store the positive outliers > Q3 + 3IQR

                fstmax.push_back(maxfst);

                if(verbose) cout<< "\nLooking for positive outliers. Found "<<eposet.size()<<endl;

                //double piD{1.0},LBq{0};

                if(!eposet.empty()){

                // are there outliers inside the outliers? use k-mean variant with k=2 the cutting point between modes would be  the mean plus the upp-bound of the standard deviation fsttot+maxdev
                    vector<double>modes{0,maxfst};

        //cout<<eposet.size()<<endl;
                    bool extreme=true;

           /// sort eposet by the fst values in ascending order (which should sort the pvalues in the reverse order)
                    sort(eposet.begin(),eposet.end(),[](const ResultFst& epo1,const ResultFst& epo2){return epo1.fst < epo2.fst;} );
            ///The whole set of outliers
                    outliers = eposet;
            /// also sort the whole set of pvalues as they are needed for sdpbound and we are interested also in those that are <= SL
                    sort(pvals.begin(),pvals.end());

                    if(contSLsigs){
                        stat = (numpops-1)*maxfst;
                        if(fsttot>0)
                            stat/=fsttot;

                    }


          /// ******************** Compute the extreme outlier set

                    double cutoff= (fsttot + numdesvs*maxdev);

                    ///it modifies eposet leaving only the eos outliers: outliers contains all the outliers and eposet only the eos outliers
                    kmeans_outliers(eposet, modes,cutoff,extreme);

                    outliers.resize(numoutliers- eposet.size()); // conserve the sorted fsts of the first class. The non EPO outliers

                /// Nov 2023 ******* Store outliers positions for using as candidates
                    if(!candidates.empty()) //clean the vector for each file
                        //vector<double>().swap(candidates);
                        candidates.clear();

                    double lastpval=1;// to compare if two candidates are close for choosing the more significant

                    for(const auto& eo:eposet){
                        if(candidates.empty()){ // first candidate
                            candidates.push_back(eo.pos); // store the positions of non consecutive extreme outliers
                            lastpval=eo.pval;
                        }
                        else if(fabs(candidates.back() - eo.pos)>=1 ){ /// May 2024
                                candidates.push_back(eo.pos); // store the positions of non consecutive extreme outliers
                                lastpval=eo.pval;
                            }
                            else{ // candidates are too close choose the more significant
                                if(eo.pval<lastpval){ // if more significant substitute by the new one, if not, do nothing
                                    candidates.back()=eo.pos;
                                    lastpval=eo.pval;
                                }
                            }
                    } // close for eposet
                    /// Now store the normal outliers
                    for(const auto& o:outliers){
                        if(candidates.empty()){ // if there were no eos this is the first candidate
                            candidates.push_back(o.pos); // store the positions of non consecutive extreme outliers
                            lastpval=o.pval;
                        }
                        else if(fabs(candidates.back() - o.pos)>=1 ){ /// May 2024
                                candidates.push_back(o.pos); // store the positions of non consecutive extreme outliers
                                lastpval=o.pval;
                            }
                            else{ // candidates are too close choose the conservatively the less significant because this are not extreme outliers
                                if(o.pval>lastpval){ // if less significant, be conservative and substitute by the new one, if not, do nothing
                                    candidates.back()=o.pos;
                                    lastpval=o.pval;
                                }
                            }
                    } // close for eposet

                    //nooutliers=false;
                } // close if not eposet empty
                else{
                    //nooutliers=true;
                    if(!candidates.empty())
                        vector<double>().swap(candidates);

                }
         /// Store for output

                //lbqval.push_back(LBq);
                contSLsig.push_back(contSLsigs);
                if(verbose) cout<< "\nEOS analysis completed"<<endl;
            }// close if fsttot<0.5
            else{ ///April 2024 just get the highest 10% Fst values
                if(verbose) cout<< "\nFST mean too high (>0.5). Sites with Fst in the 90th percentile are considered as candidates if option -useblocks 0 was chosen."<<endl;
                ///April 2024
                highmeanFST.push_back(true);
                k_percentile(genomesize,0.9,fstdistrib, metapopsnpositions,filelocnames, eposet);

                /// April 2024 ******* Store 95 percentile positions for using as candidates
                    if(!candidates.empty()) //clean the vector for each file
                        //vector<double>().swap(candidates);
                        candidates.clear();

                    for(const auto& eo:eposet){
                        if(candidates.empty()){ // first candidate
                            candidates.push_back(eo.pos); // store the positions of non consecutive extreme outliers
                        }
                        else if(fabs(candidates.back() - eo.pos)> 1 ){ ///May 2024
                            candidates.push_back(eo.pos); // store the positions of non consecutive extreme outliers

                        }
                    } // close for eposet


                    /// Clean the eposet because there were no true outliers they are only used as candidates because they high FST
                    vector<ResultFst>().swap(eposet);

            } //close else


        } // close if doEOS
/// ****************************** EOS ENDED

/// ************************************************** CANDIDATES FOR HAPLOTYPE ANALYSES *******************************************************************
   if(!soloEOS){


/// ACR June 21: Now update the candidate selective positions if any
/// First update the reference population if not user-defined
    if(refpopforHB<0){ // not user defined
        popforHB=Pref;
    } // close if popforHB
    else{
        autorefpop=false;
    }

/// Search the index in metapopsnpositions for the values in selective. The current index metapopsnpositions of the selective site
/// corresponds to its position in the metapopdata vector used for the computations
// the function get_indx(const Col& v,const T& value) Search value within sorted v (USE Tfindindx if v is non-sorted)
// and returns the index if v[index]==value or -1 if value is not in v

        /// Enero 24 calcularlo aquí en vez de 2583
        /// acr: July 21 Compute the reference (the haplotype carrying major alleles) for each population.
        popreferences = computereferencepops(metapopdata,samplesizes);

        //cout<<std::accumulate(popreferences[popforHB].begin(), popreferences[popforHB].end(), 0)<<endl;
        //assert(std::accumulate(popreferences[popforHB].begin(), popreferences[popforHB].end(), 0)==0);

     /// Oct 2023: the reference for using in the K-L method is popreferences[popforHB] that correspond the one with the highest sample size (if not user-defined)

    /// April 2024 Compute Dprime and r2 for every SNP pair just one time, here.
    computeSeqDrsq(popreferences[popforHB], SL,metapopdata, samplesizes, popforHB, snpfreqs[popforHB], Dprime, r2);

    /// May 2024 compute Dprime and r2 for the other population only if the user has not fixed the reference population
    if(autorefpop)
        computeSeqDrsq(popreferences[1-popforHB], SL,metapopdata, samplesizes, 1-popforHB, snpfreqs[1-popforHB], DprimeB, r2B);

    /// May 2024
    int conthits=0,accconthits=0;
    double aveDprime=0,aver2=0,aveDprimeB=0,aver2B=0;

    if(!selectives.empty()){ /// If there are user def selective positions then update the positions after cleaning any user-defined if it is non-shared and non-maf SNP

        /// v0.5.3 April 2026
        bool select_all = false;
        for (const auto& sel : selectives) {
            if (std::isinf(sel)) { // el valor inf o -inf indican que se usarán todos los sitios como candidatos
                select_all = true;
                break;
            }
        }

        if(select_all){

            selindexes.resize(metapopsnpositions.size());
            std::iota(selindexes.begin(), selindexes.end(), 0);

        }
        else{
            if(msformat){ // positions are relative
                int site=-1;
                for(auto sel: selectives){
                    int si = get_indx<double>(metapopsnpositions,sel);

                    if(si>=0){
                        selindexes.push_back(si);
                    }
                    else{ // if it is not a shared position change by a close one if the distribution is uniform
                        site=metapopsnpositions.size()*sel; // because positions are relative
                        selindexes.push_back(site);
                        if(verbose){
                            cout<< "\nCandidate changed to position "<<site <<endl;
                            mensaje+="\nBecause proposed candidate is not shared the site is changed to position "+NtoS(metapopsnpositions[site])+"\n";
                        }
                    }
                }
            } // close if ms format
            else{

                for(auto sel: selectives){ /// *** it is assumed metapopsnpspositions is sorted (and it is)
                    int si = get_indx<double>(metapopsnpositions,sel);
                    if(si>=0){
                        selindexes.push_back(si); // charge sel indexes with the updated metapopulation candidate positions
                    }
                    else{ // if it is not a shared position change by a close one
                        auto minmaxpos = minmax_element(begin(metapopsnpositions),end(metapopsnpositions));
                        int half = (*minmaxpos.second - *minmaxpos.first)/2;
                        int pos = sel;
                        bool upstream=true;
                        if(pos> half)upstream=false;

                        do{

                            if(upstream){
                                if(pos<*minmaxpos.second) // if not the maximum position increase
                                    ++pos;
                                else break;
                            }
                            else{
                                if(pos>*minmaxpos.first) // if not the minimum position decrease
                                    --pos;
                                else break;
                            }

                            si = get_indx<double>(metapopsnpositions,pos);

                        }while(si<0);

                        if(si>=0){
                            selindexes.push_back(si);
                            if(verbose){
                                cout<< "\nCandidate changed to position "<<si <<endl;
                                mensaje+="\nBecause proposed candidate is not shared the site is changed to position "+NtoS(metapopsnpositions[si])+"\n";
                            }
                        }

                    }
                }
            } // close else ms format
        } // close else select_all

        if(selindexes.empty()){ //the candidate SNPs were not shared between populations or were not at MAF
                string m = "iHDivSel ERROR:\nCannot found between-population-shared SNP with frequency above MAF and close to any of the proposed candidate SNPs for studying selection.";
                SNPSelTools::texterr(1,m,true);
                cout<< "\nCannot found between-population-shared SNP with frequency above MAF and close to any of the proposed candidate SNPs for studying selection." <<endl;
        }
        else{ ///Enero 2024 The selective candidates are present in the data

                /// May 2024. Do only if not fixwinsize
                if(!fixwinsize){
                    auto t0 = chrono::high_resolution_clock::now();
                    if(verbose){
                            cout<< "Going into estimating window sizes for each snp"<<endl;

                    }

                    for(auto snp : selindexes){

                    /// Estimate the window size for the whole sample of the pop popforHB
                    /// for each snp the window size will be the same for both populations
                    ///April 2024
                        computeWinSizesSeq(SL, genomesize, Dprime,r2, minblocksize, DBT, samplesizes, popforHB, snp, iniw,endw,conthits,aveDprime,aver2);
                        ///May 2024
                        if(endw-iniw>=minblocksize){
                            blocks.push_back(Block(popforHB,iniw,endw));
                            accconthits+=conthits;
                        }
                        else{///May 2024
                            if(autorefpop)
                                computeWinSizesSeq(SL, genomesize, DprimeB,r2B, minblocksize, DBT, samplesizes, 1-popforHB, snp, iniw,endw,conthits,aveDprimeB,aver2B);
                            if(endw-iniw>=minblocksize){
                                blocks.push_back(Block(1-popforHB,iniw,endw));
                                accconthits+=conthits;
                            }
                            else{ // April 2026 si seguimos sin bloque lo forzamos de tamaño 1 y lo comprobamos al calcular el estadístico

                                blocks.push_back(Block(popforHB,snp,snp+1));

                            }

                        }


                    }// close for selindexes

                    if(verbose){
                            cout<< "Window sizes estimated for each snp"<<endl;
                            auto t1 = chrono::high_resolution_clock::now();
                            cout << "Time required: "
                                 << chrono::duration<double>(t1 - t0).count()
                                 << " s" << endl;
                    }


                    blocksdetected=!blocks.empty();
                    if(blocksdetected){

                        numblocks=blocks.size();

                        if (int(selindexes.size()) > numblocks) {
                            cerr << "ERROR: number of blocks (" << blocks.size()
                                 << ") does not match number of candidates (" << selindexes.size()
                                 << ")." << endl;
                            exit(1);
                        }

                    }
                    else{
                        if(verbose){
                            cout<<"No blocks detected in run "<<run;
                            cout<<". The HAC method cannot detect selection without haplotype blocks."<< endl;
                        }

                        seldetection=false;
                        skiprun=true;
                        cout<<"\n***NO SELECTION DETECTED IN THIS FILE.***\n"<< endl;
                        mensaje+="\n\n***NO SELECTION DETECTED IN THIS FILE.***\n";
                        goto NOSNPS;
                    }

                } // closing if not fixwinsize

        }//close else is selindexes empty
    } // close if not selectives.empty()
    else{ /// No user-defined selective sites:
        /// Oct 2023 charge selindexes with the positions of candidate outliers
        if(!useblocks){
            for(auto sel: candidates){ /// *** it is assumed metapopsnpspositions is sorted (and it is)
                    //int si = get_indx<double>(metapopsnpositions,sel);
                    selindexes.push_back(get_indx<double>(metapopsnpositions,sel)); // charge sel indexes with the updated metapopulation candidate positions

            }
            if(selindexes.empty()){ // there were no outliers neither mean Fst>0.5
                if(verbose){
                    cout<<"No outliers detected in run "<<run<<" so blocks centered on outliers cannot be made."<<endl;
                    cout<<"The JHac method cannot detect selection without haplotype blocks."<< endl;
                }

                seldetection=false;
                skiprun=true;
                cout<<"\n***NO SELECTION DETECTED IN THIS FILE.***\n"<< endl;
                goto NOSNPS;

            }// close if selindexes.empty()

        } // close if !useblocks
        //if(useblocks || selindexes.empty())
        else { // use blocks
            ///April 2024
            computeBlocks(SL, Dprime,r2, metapopdata, minblocksize, DBT, samplesizes, popforHB, snpfreqs[popforHB], blocks); //r² significant or DBT > 0.5 is the LD threshold used for considering block
            blocksdetected=!blocks.empty();
            if(blocksdetected){
                numblocks=blocks.size();
                for(int s=0; s< numblocks;++s){
                    selindexes.push_back((blocks[s].iniw+blocks[s].endw)/2);// middle point between ini and end
                }
            }
            else{
                /// First else: May 2024 check blocks in the other population
                if(autorefpop){
                    computeBlocks(SL, DprimeB,r2B, metapopdata, minblocksize, DBT, samplesizes, 1-popforHB, snpfreqs[1-popforHB], blocks);
                    blocksdetected=!blocks.empty();
                }

                if(blocksdetected){
                    numblocks=blocks.size();
                    for(int s=0; s< numblocks;++s){
                        selindexes.push_back((blocks[s].iniw+blocks[s].endw)/2);// middle point between ini and end
                    }


                }
                else{ /// 2nd else no blocks detected with the minimum window size minWinSize so define 1 block in the middle of the genome with minWinSize and skip the run

                    if(verbose){
                            cout<<"No blocks detected in run "<<run;
                            cout<<". The HAC method cannot detect selection without haplotype blocks."<< endl;
                            mensaje+="\n\n******* No blocks detected in file "+inputname+". The HAC method cannot detect selection without haplotype blocks.*******\n\n";
                    }

    //                selindexes.push_back(genomesize/2); // middle point
    //                blocks.push_back(Block(0,minWinSize));
    //                numblocks=1;
                    seldetection=false;
                    skiprun=true;
                    cout<<"\n***NO SELECTION DETECTED IN THIS FILE.***\n"<< endl;
                    goto NOSNPS;
                } // close second else
            } // close first else

        } // close if selindexes empty
    } // close no user-defined selective sites

    //metan2.resize(selindexes.size(),0);

    JHac.resize(selindexes.size(),0);


    int ini{0};
    /// Oct 2023
    vector<vector<double>> popHACS(numpops);

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

        /// Dec 21
        skip[p].resize(selindexes.size(),false); // if the sample of partition 2 is too low the analysis will be skipped and this should be account in the output
        if(!seldetection){skip[p][0]=true;}

        assert(data[ini].size() == snpfreqs[p].size());
        assert(data[ini].size() == snpositions[p].size());
		mensaje+="\nNumber of true SNPs shared in pop "+NtoS(p+1) +" is "+ NtoS(snpfreqs[p].size());
        ini+=samplesizes[p];

        /// Oct 2023: prepare data for the new analyses
        popHACS[p].reserve(samplesizes[p]);
    }

// *** Now we should set a window size to study the SNPs and compute the HACS
// *** We may use the whole set of SNPS winsize=numsnps[p] or subsets
/// acr June 21: however if usercandidate is true we will use a specific snp position...
/// acr July 21: ...to test the cluster effect around it


/// ****************************************** HAPLOTYPE ANALYSES *******************************************************************

/// ********** PERFORM D and haplotype-based COMPUTATIONS

    int metapop{0};


/// Compute initial window sizes for the metapopulation and each population case
/// Oct 2023: FORCE WINDOW SIZES TO BE MULTIPLES OF 10 TO EASILY CONTRUCT FREQUENCY CLASSES FOR THE NEW METHODS
    if(fixwinsize){// the window size is defined by the user
        if(numsnps[metapop]<minWinSize){
            string msg="\nNot enough snps (" + NtoS(numsnps[metapop]) + ") in metapopulation "+NtoS(metapop)+" in file "+inputname;
            cout<<msg<<endl;
            mensaje+= msg;
            SNPSelTools::texterr(numsnps[metapop],msg,false);
            skiprun=true; // skip the variance analysis
            /// Enero 2024 descuenta esta réplica
            --totruns;
        }
        else{
            if(numsnps[metapop]<winsize)
            {


                string msg="\nFixed Window size is: " + NtoS(winsize) +".\nIn metapopulation "+NtoS(metapop+1)+" only "+NtoS(numsnps[metapop])+" true SNPs in the sample." " The window size is changed for file "+inputname+". New window size is: ";


                /// Oct 2023 force NUMCLASSES to be half of the window size and...
                winsize=numsnps[metapop];
                NUMCLASSES=winsize/dividend; // dividend is =2
                if(NUMCLASSES<minNUMCLASSES) NUMCLASSES=winsize; // minNUMCLASSES is 15
                /// ... force window size to be multiple of NUMCLASSES +1
                int r= numsnps[metapop]%NUMCLASSES;
                winsize-=r;
                if(winsize<numsnps[metapop]) /// acr nov 2024 add the if to void core dump if r =0 and winsize al ready equal to numsnps[metapop] before adding 1
                    winsize+=1;

                mensaje+=msg;
                SNPSelTools::texterr(winsize,msg,false);

            } // close if the window is higher than the number of available snps
            /// /// Oct 2023 force NUMCLASSES to be half of the window size and...force window size to be multiple of NUMCLASSES +1
            else{

                NUMCLASSES=winsize/dividend;
                if(NUMCLASSES<minNUMCLASSES) NUMCLASSES=winsize;
                int r=winsize%NUMCLASSES;
                (r<=(NUMCLASSES/2))? winsize-=r:winsize+=(NUMCLASSES-r);
                if(winsize<numsnps[metapop]) /// acr nov 2024 add the if to void core dump if r =0 and winsize al ready equal to numsnps[metapop] before adding 1
                    winsize+=1;

            }// close else
            ///May 2024
            referencehap=popreferences[popforHB];
        } //close else addded in enero 2024
    } // Close if fixwinsize
    else{ // no window size predefined we will compute it depending on linkage

        if(numsnps[metapop]<minWinSize){
            string msg="\nNot enough snps (" + NtoS(numsnps[metapop]) + ") in metapopulation "+NtoS(metapop)+" in file "+inputname;
            cout<<msg<<endl;
            mensaje+= msg;
            SNPSelTools::texterr(numsnps[metapop],msg,false);
            skiprun=true; // skip the analysis
            /// Enero 2024 descuenta esta réplica
            --totruns;
        }
        else{
            if(MaxWinSize==0){ // If MaxWinSize is undefined (value 0) we use all the snps for defining the size
               winsize=numsnps[metapop]; // The winsize is set to the maximum possible
               ///May 2024
               if(winsize>MaxWinSize) winsize=MaxWinSize;
                //else // MaxWinSize <0: will use haplotype blocks
                 /// /// Oct 2023 force NUMCLASSES to be half of the window size and force window size to be multiple of NUMCLASSES +1
                NUMCLASSES=winsize/dividend;
                if(NUMCLASSES<minNUMCLASSES) NUMCLASSES=winsize;
                int r= numsnps[metapop]%NUMCLASSES;
                winsize-=r;
                if(winsize<MaxWinSize) /// acr nov 2024 add the if to void core dump if r =0 and winsize already equal to MaxWinSize before adding 1
                    winsize+=1;
            }
            else{ ///April 2024
                winsize=minblocksize;
            }
        } // close else addded in enero 2024
    }

    if(!skiprun){

         // For computing the average number of segregating sites through all analyzed files (runs)
assert((size_t)numsnps[0]==metapopsnpositions.size());

        totsnps.push_back(numsnps[0]);

        SegSites+=(metapopsnpositions.size());

        /// ***** Window overlapping

        ws=winsize;

        // Compute the mamixum number of canditates among all files
        size_t maxnumcandidates = selindexes.size();
        totalmaxcandidates= max(totalmaxcandidates,maxnumcandidates);


        if(verbose){
            cout<<"Population references computed: "<<popreferences[0].size()<<","<<popreferences[1].size()<<endl;
        }

        if(verbose){cout<<"Number of candidates to analize "<<maxnumcandidates<<"."<<endl;}

        center=(ws-1)/2;

    /// ***********************************************************************************************************
    ///     ******************************* J Divergence Selection Test **************************************
    /// ***********************************************************************************************************

        if(verbose) cout<< "\n... Beginning J analysis ...\n"<<endl;
        // for the desired candidate snp
        int ss=0; // candidate counter
        ResultIndexes tempIndex;
        //tempIndex.varTest.resize(numpops,0);
        tempIndex.iniW.resize(numpops,-1);
        tempIndex.endW.resize(numpops,-1);

        /// 1.- for each candidate SNP
        for(auto snp : selindexes){ // each element snp is the position in the metapopdata.

    ///


            iniw=snp-center; // center =(ws-1)/2
            /*
            e.g. snp=22, ws= 23, center=(23-1)/2= 11, iniw=22-11=11
            end2 = ini+ws = 11+23=34
            endw>23, endw=23, iniw = 23-23 =0
            */
            endw=iniw+ws;
            if(iniw<0){iniw=0;endw=ws;}
            else if(endw>(numsnps[metapop])){endw=numsnps[metapop]; iniw=numsnps[metapop]-ws;}

            if(verbose){
                cout<<"\nAnalizing CANDIDATE "<<ss+1<<" which corresponds to SNP "<<snp<<" in the metapopulation."<<endl;
            }

//cout<< numsnps[metapop] << " - "<<ws<<endl;
            /// Recompute window size if not fixwinsize and not the whole genome as windowsize
            if(!fixwinsize && MaxWinSize>0){

                if(!useblocks && !candidates.empty()){ /// blocks were not computed because there are outliers

                /// Oct 2023 compute Dprime and r2 using the whole reference population popforHB instead of a partition for each p
                /// Oct 23 estimate the window size for the whole sample of the pop popforHB
                /// for each snp the window size will be the same for both populations
                    ///April 2024

                    computeWinSizesSeq(SL, genomesize, Dprime,r2, minblocksize, DBT, samplesizes, popforHB, snp, iniw,endw,conthits,aveDprime,aver2);
                    referencehap=popreferences[popforHB]; ///May 2024 also update the reference haplotype
                    accconthits+=conthits;

                    if(endw-iniw<minblocksize){

                    /// May 2024 recompute for the other population
                        if(autorefpop){
                            computeWinSizesSeq(SL, genomesize, DprimeB,r2B, minblocksize, DBT, samplesizes, 1-popforHB, snp, iniw,endw,conthits,aveDprimeB,aver2B);
                            referencehap=popreferences[1-popforHB]; ///May 2024 also update the reference haplotype
                        }
                        if(endw-iniw<minblocksize){
                            endw = (iniw+minblocksize);
                            if(endw>(int)genomesize){
                                string msg="\niHDS ERROR: window end " + NtoS(endw) +" out of the genome range and changed to "+NtoS(genomesize)+".\n";
                                endw=genomesize;
                                iniw=genomesize-minblocksize;
                                SNPSelTools::texterr(iniw,msg,false);
                            }
                        }
                    }

                } // close if  not useblocks && candidates not empty
                else{ // useblocks or candidates empty
                    iniw = blocks[ss].iniw;
                    endw = blocks[ss].endw;
                    referencehap=popreferences[blocks[ss].refpop];

                }

                /// Update window size
                /// Oct 2023 force NUMCLASSES=ws/2 and ws to be multiple of NUMCLASSES +1 so that the hac values has a maximum L that is multiple of NUMCLASSES
                ws =endw-iniw;
                ///May 2024
                if(ws>MaxWinSize) ws=MaxWinSize;
                NUMCLASSES=ws/dividend;
                if(NUMCLASSES<minNUMCLASSES) NUMCLASSES=ws; // April 2026 ws es 1 si definí el bloque con iniw=snp endw=snp+1 entonces NUMCLASSES queda a 1 pq minNUMCLASSES está definido a 15
                int r= ws%NUMCLASSES;
                (r<=(NUMCLASSES/2))? ws-=r:ws+=(NUMCLASSES-r); //round lower or upper
                ws+=1;

            } // close if( !fixwinsize && MaxWinSize>0)


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

                /// store window size for each population
                iniwin[p]=iniw;
                endwin[p]=endw;

                /// Make the partitions for each population

                //dataPartition(snp,p, samplesizes, metapopdata, popreferences[p],Partition1[p],Partition2[p], n1[p][ss],n2[p][ss]);
                //metan2[ss]+=n2[p][ss];

                /// ***** Compute Hacs for each population

                //int iniSample=p*samplesizes[0];
                /// April 2026
                int iniSample = accumulate(samplesizes.begin(), samplesizes.begin() + p, 0);

                computePopHacs(iniSample, iniSample+samplesizes[p],snp,metapopdata,referencehap, iniw,endw,popHACS[p]);

            assert(popHACS[p].size()==(size_t)samplesizes[p]);


                if(verbose){
                    cout<<"Population "<<p <<". Window range "<<iniw<<"-"<<endw<<"."<<endl;
                }

    /// acr nov 23: Compute the number of classes following the Diaconis_rule https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule
                //K[p] = computeNumberofClasses(p,popHACS[p]);

            } // close for pops


            tempIndex.nclasses=NUMCLASSES;

            /// Oct 2023 Hac frequency classes
            vector<double>P(NUMCLASSES,0),Q(NUMCLASSES,0);

            /// Oct 2023 Compute the JHac divergence
            /** May 2024
            for the case with outliers as candidates the block is computed with ComputeWinszieSeq function
            in such case may there are not blocks but the values of iniw and endw are put to minblocksize
            when the jac test is passed to such "blocks" sometimes they are significant but they are not true blocks
            and may be are false positives. So in such cases we consider to pass the test only if such "blocks" hold at least some
            minimal conditions. For example we know that the parts of the block may not hold the condition of Dprime>DBT or r2 significant
            but if we compute the average DBT and the averager2 of the block we can check if the sum of the averages > 0
            which will mean that there were some pair of values that were connected one way or another

            if the case if a highmean FST should be calculate always or let the same condition aveDprime+aver2??

            **/
            //if(useblocks || fixwinsize ||(aveDprime>0)||(aveDprimeB>0)|| (liberal&&highmeanFST[run])){ ///May 2024
            if(useblocks || fixwinsize ||(aveDprime>0)||(aveDprimeB>0) || (liberal&&highmeanFST[run])){ ///May 2024 liberal is set to false by default in line 1034



                //if(NUMCLASSES>=minNUMCLASSES) //5 April 2026
                if(NUMCLASSES>1){ // 12 April 2026
                    JHac[ss]=makeHacClassesandJ(NUMCLASSES,ws-1, popHACS, P, Q);
                    tempIndex.JHac.stat=JHac[ss];
                    tempIndex.JHac.pval=getChiPval(JHac[ss],NUMCLASSES-1);
                }
                else{
                    JHac[ss]=0;
                    tempIndex.JHac.stat=JHac[ss];
                    tempIndex.JHac.pval=1.0;
                }

            }
            else{ ///this is the case of !blocks i.e. blocks centered in outliers and for this snp there were no block

                tempIndex.JHac.stat=0;
                tempIndex.JHac.pval=1;

            }

            /// acr April 2024
            tempIndex.wide=endw-iniw;
            for(int p=0; p< numpops; ++p){
            /// We have already the window size for the population with the highest sample size: popforHB.
            /// Now compute the hac values for partition P2. Accumulate the HAC values between populations.
            /// Set the same window size for every population
            ///Enero 2024 comenta el computeMetaHacc que ya no se usa
            //computeMetaHacc(snp,Partition2[p], popreferences[p], iniwin[popforHB],endwin[popforHB],metaHACdistri2);
                tempIndex.iniW[p] = metapopsnpositions[iniwin[popforHB]]; ///May 2024 this is ok because iniwin[] holds the same value for both populations (see line 2860)
                tempIndex.endW[p] = metapopsnpositions[endwin[popforHB]-1];
            } // close for p


        ///Enero 2024 comenta todo el bucle de cálculo de HACdistri1 y 2 que ya no se usa

            //vector<deque < deque<bool> >>().swap(Partition1);
            //vector<deque < deque<bool> >>().swap(Partition2);

/// Enero 2024

            tempIndex.candidate=ss+1; // candidate is de index of the selective site +1
            tempIndex.site=selindexes[ss]; // site is the number of the SNP in the data (i.e. SNP number 3 in a range from 0 to metapopdata.size()-1)
            tempIndex.pos = metapopsnpositions[selindexes[ss]]; // real genomic snp position as indexed in the metapopsnppositions vector

            indexes.push_back(tempIndex);

            ss++;
        } // close for the snp candidates in selindexes

        /// Dec 21
        runsskip.push_back(skip);

    } // close !skiprun

    } // close if not soloEOS


/// ************************** Fst for HAPLOTYPES:
        if((!skiprun) && (!soloEOS)){
 // Fst in selective positions
            double fstsel = 0.0;

///ACR May16 warning. WE use as expectation the weighted (by the sample size) mean of the two population frequencies which is equivalent to the pooled frequency
            //computemeanSNPFreqinMetaPop(metapopdata, newsnpfreqs);
// Oct 2023
if(verbose)

    cout<< "\nStoring Fst index p-values"<<endl; /// Nov 21

            for(size_t candidate=0; candidate<indexes.size(); ++candidate){

                //fstsel=max(FST(metapopdata, samplesizes,indexes[candidate].site ,indexes[candidate].site+1),0.0);
                fstsel = fstdistrib[indexes[candidate].site].stat;

                indexes[candidate].fst= fstsel;

                indexes[candidate].fstindex.stat= (fstsel - fsttot);


                /// Oct 2023: store the LK test value instead of doing the resampling
                //indexes[candidate].fstindex.pval=0;
                indexes[candidate].fstindex.pval=fstdistrib[indexes[candidate].site].pval;
            }

            ///April 2024: update here the mean Fst because to store only when the run has selective blocks and will be written to the output file
            fstmean.push_back(fsttot);

/// **** Oct 2023 Sort indexes
// sort the indexes in descending order from higher to lower JHac value
            //sort(indexes.begin(), indexes.end(), [](const ResultIndexes& candidate1,const ResultIndexes& candidate2){return candidate1.JHac.stat > candidate2.JHac.stat;} );
/// FROM April 2024 se recupera el orden natural por ello se comenta el sort anterior y se añade uno que ordena por posición
/// ello es necesario porque si los candidatos se eligieron a partir de outliers el orden de posición no está garantizado
/// el Fst mean no se ve afectado porque es el mismo para todos los candidatos
            sort(indexes.begin(), indexes.end(), [](const ResultIndexes& candidate1,const ResultIndexes& candidate2){return candidate1.pos < candidate2.pos;} );


        /// Resample using mean metapopulation frequencies

            //size_t resamplings = resampnumiter;
            // here the p-values are computed and the ones corresponding to candidates stored in indexes[candidate].pval

            //doIndexNullMetapop2(metapopsnpfreqs[0],numpops,nummetapops,samplesizes,indexes,0, resamplings,fstdistrib,SL); // without q-value estimation

            indexesfiles.push_back(indexes);

        } // close if no skiprun && no soloEOS

        fstfiles.push_back(eposet); /// extreme outliers set
        fstfiles0.push_back(outliers); /// non extreme outlier set
        //locnames.push_back(filelocnames);

        namefiles.push_back((inputname));

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

        if(!outliers.empty()){
            vector<ResultFst>().swap(outliers);
        }


        vector<ResultIndexes>().swap(indexes);



    } // Close if numsnp>minWinSize for fst
    else{ // too few SNPs
        string msg = "\nIn file "+inputname+ ": too few SNPS "+ NtoS(numsnps[0])+" < "+NtoS(minWinSize)+"\n";
		cout<<msg<<endl;
		SNPSelTools::texterr(numsnps[0],msg,false);




    } //SNPs less than min win size


    //minsamplesizes.push_back(*min_element(samplesizes.begin(),samplesizes.end())); // write the minimum samplesize analyzed in the latest file


    /// New Oct21
    if(verbose){cout<<"selection indexes dimension "<<selindexes.size()<<endl;}
    if(!selindexes.empty())
        //vector<int>().swap(selindexes);
        selindexes.clear();


    NOSNPS:


    if(verbose){cout<<"Population references "<<popreferences.size()<<endl;}
    for(auto& refp:popreferences){
        //vector<bool>().swap(refp);
        refp.clear();
    }

    if(verbose){cout<<"Data size to clean "<<data.size()<<endl;}

    for(auto& dat:data){
        //deque<bool>().swap(dat);
        dat.clear();
    }
	//deque < deque<bool> >().swap(data);
	data.clear();

    if(verbose){cout<<"SNP positions dimension "<<snpositions.size()<<endl;}
    for(auto& snp:snpositions){
        //vector<double>().swap(snp);
        snp.clear();
    }

	//vector < vector<double> >().swap(snpositions);

	if(verbose){cout<<"SNP freqs dimension "<<snpfreqs.size()<<endl;}
    for(auto& snp:snpfreqs){

        //vector<double>().swap(snp);
        snp.clear();
    }


    if(verbose){cout<<"Metapop SNP positions dimension "<<metapopsnpositions.size()<<endl;}

	if(!metapopsnpositions.empty())
        vector<double>().swap(metapopsnpositions);
        //metapopsnpositions.clear();

    if(verbose){cout<<"commons dimension "<<commons.size()<<endl;}
	if(!commons.empty())
        vector<double>().swap(commons);

    if(verbose){cout<<"Fst distribution dimension "<<fstdistrib.size()<<endl;}
    if(!fstdistrib.empty())
        //vector<STAT>().swap(fstdistrib);
        fstdistrib.clear();

    /// Nov 21
    vector<Block>().swap(blocks);
    delete[] iniwin;
    delete[] endwin;

    if(verbose){cout<<"Memory freed"<<endl;}

    // Update random generator
    if(!fixseed){
        //seed = (myclock::now() - beginning).count();
        seed = myclock::now().time_since_epoch().count();
        gen.seed(seed);

        mensaje+="\nNew seed "+NtoS(seed)+ "\n";

        if(verbose) cout<<"\nNew seed "+NtoS(seed)<<endl;
    }



    cout<<"\nAnalysis of file "<<inputname<<" completed"<<endl;
	cout<<"\n------------------------------------------------\n"<<endl;

} // close for files

/// ********* OUTPUT

    if(!soloEOS){

        if(autorefpop){
            outname=dir+"AutoRefPop"+NtoS(popforHB+1)+"_MAF_"+NtoS(iniMAF);
        }
        else
            outname=dir+"RefPop"+NtoS(popforHB+1)+"_MAF_"+NtoS(iniMAF);

        if(NumOfFiles>1)
            outname+=("_NF"+NtoS(NumOfFiles));
        if(old)
            outname+="_OldSweep";
        if(useblocks)
            outname+="_Blocks";

        if(nogaps)
            outname+="_NoGaps";

        if(tag=="rnd0")
            outname+="_RND";
        if(fixwinsize>0){
            outname+="_W"+NtoS(fixwinsize);

        }

    }
    else{
        outname+="_MAF"+NtoS(iniMAF);
    }

    foutnamemax= outname+".tsv";

    cout<<"WRITING IHDivSel OUTPUT FILE " << foutnamemax<<endl;

    //writeVartest(foutnamemax,namefiles,indexesfiles,fstfiles, contSLsig,SL,SegSites,numoutliers, totsnps,numindepsnps,
    //          iniMAF,fstmax,fstmean,pi0,lbqval,snpstable, totalmaxcandidates,numpops, runsn1,runsn2,runsVarFTest,runsskip,testtype);
    //fst files correspond to the eposet data
    /// Enero 2024: incluyo la variable totruns que en el caso de useblocks cuenta también las réplicas en las que no
    /// haya bloques descartándose solo al principio aquellas que no tienen suficientes snps
    writeJtest(highmeanFST,totsnps,fstmean,totruns,foutnamemax,namefiles,indexesfiles,fstfiles, contSLsig,SL,SegSites,snpstable, totalmaxcandidates,numpops,runsamplesizes);
    if(doEOS){
        foutnamemax= outname+"_FSTs.txt";
        writeFSTs(foutnamemax,namefiles,fstfiles0,fstfiles,SL,totsnps,iniMAF);

    }


// **********

secs = (double)(clock() - start) / CLOCKS_PER_SEC;
cout<<"\n\niHDivSel FINISHED. ";
mensaje+="\n\niHDivSel FINISHED. ";
if(secs>60.0){
	secs/=60.0;
	string msg="Time processing: "+NtoS(secs)+" minutes";
	cout<<msg<<endl;
	//SNPSelTools::texterr(msg,false);
	mensaje+=msg;
}
else{
	string msg="Time processing: "+NtoS(secs)+" seconds";
	cout<<msg<<endl;
	//SNPSelTools::texterr(msg,false);
	mensaje+=msg;

}
mensaje+="\n---------------------------------------------------------------------\n";

if(verbose) SNPSelTools::texterr(0,mensaje,false);


} catch(const std::overflow_error& e) {
    cout<< e.what()<<endl;
} catch (const std::exception& ex) {
    // ...
    cout<<"\nAn exception occurred during iHDivSel execution"<<endl;
} catch (...) {
    // ...
    cout<<"\nUnexpected problem occurred during iHDivSel execution"<<endl;
}

return 0;

} // ********** CLOSE MAIN
// using
//Note that data could include information from more than one pop and different sample sizes so we need to pass a vector of samplesizes and the pop index to access the correct sequences
//Compute the frquency of 1's (be minor or major)
template<typename TR>
TR computeSNPFreqinPop(const deque< deque<bool> >&data,const int& pop, const vector<int>&samplesizes, vector<double>& snpfreqs){

    int ini = accumulate(samplesizes.begin(),samplesizes.begin()+pop,0); // the sum of elements in populations previous to pop
    auto numsnps = data[ini].size();
    if(!snpfreqs.empty()){

        vector<double>().swap(snpfreqs);
    }
    snpfreqs.resize(numsnps,0);
    // Compute the SNPs freq
    for(auto smpl=ini; smpl<(ini+samplesizes[pop]); ++smpl){
        for(size_t snp=0; snp<numsnps; ++snp)
            if(data[smpl][snp]){ // This is the frequency of the derived allele because has an 1
                snpfreqs[snp]+=1.0;
            }
    }

    for(auto& snp: snpfreqs)
        snp/=samplesizes[pop];
}

// This computes a vector with mean frequencies in the metapopulation which results as considering all data as coming from one population.
template<typename TR>
TR computemeanSNPFreqinMetaPop(const deque< deque<bool> >&data, vector<double>& snpfreqs){


    auto numsnps = data[0].size();
    if(!snpfreqs.empty())
        vector<double>().swap(snpfreqs);
    snpfreqs.resize(numsnps,0);

    // Compute the SNPs freq
    for(size_t smpl=0; smpl<(data.size()); ++smpl){
        for(size_t snp=0; snp<numsnps; ++snp)
            if(data[smpl][snp]){ // This is the frequency of the allele 1
                snpfreqs[snp]+=1.0;
            }

        }

    for(auto& snp: snpfreqs)
        snp/=(data.size());

}


/// This procedure only update the non-common SNPs just checking the frequencies of the non-commons are between 1-maf and maf

template<typename T>
void updateMAFData(deque< deque<bool> >&data,const int& pop, const vector<int>&samplesizes, const T& MAF,
                   vector<T>& snpfreqs, vector<T>& snppos, const vector<T>& commons, const bool onlycommons){

    auto cMAF = 1.0 -MAF;
    int ini = accumulate(samplesizes.begin(),samplesizes.begin()+pop,0); // the sum of elements in populations previous to pop
    // 1.- Update Data, frequencies and position vectors
    for(int snp=0; snp<(int)snpfreqs.size(); ++snp){ //

        if( !binary_search(commons.begin(), commons.end(), snppos[snp]) ){ // If this snp is not a common one

            if( (snpfreqs[snp]>(cMAF)) || (snpfreqs[snp]<(MAF))){ // delete if it is non-polymorphic frequency

                for(int d=(ini);d<(ini+samplesizes[pop]);++d)
                    data[d].erase(data[d].begin()+snp); // delete position deletepositions1[np] of data[a]

                snpfreqs.erase(snpfreqs.begin()+snp);
                snppos.erase(snppos.begin()+snp);
                --snp;
            }

        }

    }
}
/// Oct 2023
///Mark for elimination SNPs that do not reach a minimum frequency in the metapopulation
///April 2024 clarification:
/**
Common here means that the SNP position appears in both samples. It could happen that a position x is fixed with allele 0 in Pop1
and fixed with allele 1 in Pop2. Even in this case is considered common if the positions appears in both samples.
Regarding the maf, the frequency of any allele of x in the metasample will be (1+0)/2=0.5 which is > maf so the position is maintained


**/
// From different collections i.e. a collection of collections copy the common elements in a new collection
// Copy only those having a minimum of maf freq considering jointly all the collections
// if used with maf=0 then copy all common elements
// first dimension of positions refers to the populations
// freqs first dimension is the population and the second is the frequency of 1's (which are the minor alleles in the reference population)
/// IMPORTANT: freqs for the haplotype formats refers to the same allele in both populations

/// The vectors of commons is sorted before return
template <typename T, typename Tcol, typename Tcolofcol>
inline void copycommonwithMAF(const vector<int>& samplesizes, const Tcolofcol& positions,const Tcolofcol& freqs, Tcol& commons, const double& maf){

    int n = accumulate(samplesizes.begin(),samplesizes.end(),0);
    bool consider=true;

    auto minel=infinito<size_t>(); // minel is the minimum number of elements among the collections so is the maximum possible number of common elements

    // Get the size of the lowest collection and store in minel
    for(const auto& sizepos: positions)
        if(minel>sizepos.size()) // number of positions for the corresponding population
            minel =sizepos.size();

    if(!commons.empty())

         Tcol().swap(commons);

    int numel=0,element2=-1;

// Search elements present in all collections

    int numpops=positions.size();

    for(const auto& element: positions[0]){ // each element in the first population

        double n1= freqs[0][numel]*samplesizes[0];

        for(int p2=1;p2<numpops;++p2){ // compare with elements in other populations

         element2 = Tfindindx(positions[p2], element);
         consider= (element2>=0);
         if(!consider){ // the element is not shared by all populations

            break; // break the inner for
         }
         /// if the element is common but not reaches maf frequency in the metapopulation

         double n2= freqs[p2][element2]*samplesizes[p2];
         if( ((n1+n2)/n < maf) )
            consider=false;
        }

         if(consider){ // if it is a common snp
//            if(fabs(freqs[0][numel] - freqs[p2][element2])>0.5) // the snp seems to be differentiated between both populations
//                candidatepos.push_back(element);

            if(!binary_search(commons.begin(),commons.end(),element)){
                commons.push_back(element); // add the position of the common snp


            }


        } // close if consider

        consider=true;


        ++numel;

        if( (commons.size()) ==minel) break; // cannot exists more common snps



    } // close for elements

    if(!is_sorted(commons.begin(),commons.end()))
        sort(commons.begin(), commons.end()); // sort positions

}

/// This procedure just delete the SNPs corresponding to positions not included in the vector commons.

//template<typename T>
//void deleteNonCommonData(deque< deque<bool> >&data,const int& pop, const vector<int>&samplesizes,
//                   vector<T>& snpfreqs, vector<T>& snppos, const vector<T>& commons){
//
//    int ini = accumulate(samplesizes.begin(),samplesizes.begin()+pop,0); // the sum of elements in populations previous to pop
//    // 1.- Update Data, frequencies and position vectors
//    for(int snp=0; snp<(int)snpfreqs.size(); ++snp){ //
//        /// binary_search requires commons to be sorted
//        if(!binary_search(commons.begin(), commons.end(), snppos[snp]) ){ // Delete if this snp is not a common one
//
//                for(int d=(ini);d<(ini+samplesizes[pop]);++d)
//                    data[d].erase(data[d].begin()+snp); // delete the column of this snp in the population pop
//
//                snpfreqs.erase(snpfreqs.begin()+snp);
//                snppos.erase(snppos.begin()+snp);
//                --snp;
//        }
//
//    }
//}


///More efficient version that reconstruct the lists instead of in place deleting
template<typename T>
void deleteNonCommonData(
    deque< deque<bool> >& data,
    const int& pop,
    const vector<int>& samplesizes,
    vector<T>& snpfreqs,
    vector<T>& snppos,
    const vector<T>& commons
){
    int ini = accumulate(samplesizes.begin(), samplesizes.begin() + pop, 0);
    int end = ini + samplesizes[pop];

    vector<T> new_freqs;
    vector<T> new_pos;

    new_freqs.reserve(snpfreqs.size());
    new_pos.reserve(snppos.size());

    vector<int> keep_indices;
    keep_indices.reserve(snppos.size());

    for (int snp = 0; snp < (int)snppos.size(); ++snp) {
        if (binary_search(commons.begin(), commons.end(), snppos[snp])) {
            keep_indices.push_back(snp);
            new_freqs.push_back(snpfreqs[snp]);
            new_pos.push_back(snppos[snp]);
        }
    }

    for (int d = ini; d < end; ++d) {
        deque<bool> new_row;
        new_row.resize(keep_indices.size());

        for (int j = 0; j < (int)keep_indices.size(); ++j) {
            new_row[j] = data[d][keep_indices[j]];
        }

        data[d].swap(new_row);
    }

    snpfreqs.swap(new_freqs);
    snppos.swap(new_pos);
}



// Original data is copied in newdata. The number of populations correspond to the dimension of snpositions.
// The new data will store the data of all populations as if they were one.
// It will consider only common snps if indicated or the whole set otherwise.
// The snps positions for the whole metapop will also be updated in the metapopsnppositions vector
// Finally the samplesizes will be updated in metapopsamplesize
/// This version assumes that if we only consider the shared snps the non-shared has been already deleted from data
template<typename T>
bool dataFusion(const bool shared, const deque< deque<bool> >&data, deque< deque<bool> >& newdata, const vector<vector<T>>& snpositions,
                vector<T>& metapopsnpositions, const vector<int>&samplesizes, int& metapopsamplesize){

        auto numpops=snpositions.size();

        if(data.empty() || numpops<1 ){
            string msg = "Empty data in dataFusion. It seems there are not shared SNPs in your data.";
            cout<<msg<<endl;
            SNPSelTools::texterr(data.size(),msg,true);
        }
       // Clean newdata if necessary
        if(!newdata.empty()){

            for(auto& seq: newdata){

                deque<bool>().swap(seq);
            }

            deque< deque<bool> >().swap(newdata);
        }

        if(!metapopsnpositions.empty()){

            vector<double>().swap(metapopsnpositions);
        }

        if(numpops==1){
            newdata=data;
            metapopsnpositions = snpositions[0];

        }
        else{

            newdata = data;

            size_t largernumsnps=0;
            // count the maximum number of snps in any of the populations if they have different numbers
            for_each(snpositions.begin(),snpositions.end(), [&](vector<double> element){largernumsnps= max(largernumsnps,element.size());});


            metapopsnpositions.reserve(numpops*largernumsnps);
        // For each population copy the positions to the metatpopsnpositions vector
            for(const auto& element:snpositions){ /// note that each element here is a vector of positions
                copy(element.cbegin(),element.cend(), back_inserter( metapopsnpositions)); // back_inserter add elements (push_back) to the end of metapopsnpositions https://www.cplusplus.com/reference/iterator/back_inserter/
            }

            sort(metapopsnpositions.begin(),metapopsnpositions.end());
            // Delete the repeated elements  returning an iterator to the element that should be considered its new past-the-end element
            auto last=unique(metapopsnpositions.begin(),metapopsnpositions.end());
    // See STL book (page 218) to understand why this is necessary (unique only delete the content but the vector is not resized)
            metapopsnpositions.erase(last,metapopsnpositions.end());


// This is only necessary if there are non-shared snps in data. If there are only shared (onlycommon), the shared variable is true and the procedure assumes the non-shared has been previously deleted
// At each population update data non-snp columns putting 0's in the corresponding non-common positions (that are segregating in other populations)
            if(!shared){

                uint_fast32_t p=0;
                uint_fast32_t ini{0};

                for(auto& poseachpop:snpositions){

                    uint_fast32_t acum=0; //because once inserted a position the next insertion must be shifted with respect to the non-inserted vector of snpositions.
                    for(const auto& element:metapopsnpositions){
                        if(!binary_search(poseachpop.begin(),poseachpop.end(),element)){ // insert
                            auto pos = distance(poseachpop.begin(),lower_bound(poseachpop.begin(),poseachpop.end(),element)); /// lower_bound: Returns an iterator pointing to the first element in the range [first, last) that is not less than (i.e. greater or equal to) value, or last if no such element is found

                            for(int spl=0;spl<samplesizes[p];++spl){
                                newdata[ini+spl].insert(newdata[ini+spl].begin()+pos+acum,0); // The vector is extended by inserting 0 before the element at the specified position
                            }
                            ++acum;
                        } // close if

                    }
                    ini+=samplesizes[p];
                    ++p;
                }

            } // close nonshared


            if(newdata.empty()){
                SNPSelTools::texterr(3,"Not sequences with commons",false);
                return false;
            }
            metapopsamplesize = newdata.size();
        } // close else numpops==1

        return true;

}


// Original data is copied in newdata. The number of populations correspond to the dimension of snpositions.
// The new data will store the data of all populations as if they were one.
// It considerd only common snps
// The snps positions for the whole metapop will also be updated in the metapopsnppositions vector
// Finally the samplesizes will be updated in metapopsamplesize
/// This is the version in use: specialized for shared SNPs thus assumes that the data only contains shared snps
template<typename T>
bool dataFusion(const deque< deque<bool> >&data, deque< deque<bool> >& newdata, const vector<vector<T>>& snpositions,
                vector<T>& metapopsnpositions, const vector<int>&samplesizes, int& metapopsamplesize){

        auto numpops=snpositions.size();

        if(data.empty() || numpops<1 ){
            string msg = "Empty data in dataFusion.";
            cout<<msg<<endl;
            SNPSelTools::texterr(data.size(),msg,true);
        }
       // Clean newdata if necessary
        if(!newdata.empty()){

            for(auto& seq: newdata){

                deque<bool>().swap(seq);
            }
             deque< deque<bool>>().swap(newdata);
        }

        if(!metapopsnpositions.empty()){

            vector<double>().swap(metapopsnpositions);
        }

        if(numpops==1){
            newdata=data;
            metapopsnpositions = snpositions[0];

        }
        else{

            newdata = data;

            size_t largernumsnps=snpositions[0].size();


            metapopsnpositions.reserve(numpops*largernumsnps);
// For each population copy the positions to the metatpopsnpositions vector
            for(const auto& element:snpositions){ /// note that each element here is a vector of positions
                copy(element.cbegin(),element.cend(), back_inserter( metapopsnpositions)); // back_inserter add elements (push_back) to the end of metapopsnpositions https://www.cplusplus.com/reference/iterator/back_inserter/
            }
            if(!is_sorted(metapopsnpositions.begin(),metapopsnpositions.end()))
                sort(metapopsnpositions.begin(),metapopsnpositions.end());
            // Delete the repeated elements  returning an iterator to the element that should be considered its new past-the-end element. The relative order of the elements not removed is preserved, while the elements between the returned iterator and last are left in a valid but unspecified state: https://www.cplusplus.com/reference/algorithm/unique/
            auto last=unique(metapopsnpositions.begin(),metapopsnpositions.end());
    // See STL book (page 218) to understand why this is necessary (unique only delete the content but the vector is not resized)
            metapopsnpositions.erase(last,metapopsnpositions.end());

            if(newdata.empty()){
                SNPSelTools::texterr(3,"Not sequences with commons",false);
                return false;
            }
            metapopsamplesize = newdata.size();
        } // close else numpops==1

        return true;

}


// For each population
//vector<bool> computereference(const vector<double>&snpfreqs){
//
//	vector<bool> ref(snpfreqs.size(),0);
//
//	for(int sn=0; sn< snpfreqs.size();++sn){
//
//        if (snpfreqs[sn]>0.5)
//            ref[sn]=1;
//
//	} // close for sn
//
//	return ref;
//}

// For using with the metapopulation data. The Pref population is used to compute the reference.
// This gives the right length to the reference when non-shared snps between populations are being considered (in the reference of a given population the non-shared snps will count as 0's)
// So we do not use the metapopulation snpfrequencies. The samplesize must be the Pref sample size
template<typename TBcol>
TBcol computereference(const deque < deque<bool> >&datapop,const int& Pref, const vector<int>& samplesizes){

	int numsnps=datapop[0].size(); // all data has the same length
	int count1s=0;
	int med = samplesizes[Pref]/2;
	int iniP = accumulate(samplesizes.begin(),samplesizes.begin()+Pref,0); // the sum of elements in populations previous to Pref

	vector<bool> ref(numsnps,0);

	for(int sn=0; sn<numsnps; ++sn){

		for(int sam=0; sam<samplesizes[Pref]; ++sam){
			if(datapop[iniP+sam][sn]==1) ++count1s;
		} // close for sam

		if(count1s>med) ref[sn]=1;
        count1s=0;
	} // close for sn

	return ref;
}

/// acr August 21: Compute the reference consisting in major alleles for each population
template<typename TBcol>
vector<TBcol> computereferencepops(const deque < deque<bool> >&datapop, const vector<int>& samplesizes){

	int numsnps=datapop[0].size(); // all data has the same length
	int count1s=0;
	int numpops = samplesizes.size();
	vector<vector<bool>>refs(numpops, vector<bool>(numsnps,0));
    int iniP =0; // the accumulation of samples in populations previous to each pop

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

        int med = samplesizes[pop]/2;

        for(int sn=0; sn<numsnps; ++sn){ // for each site

            for(int sam=0; sam<samplesizes[pop]; ++sam){
                if(datapop[iniP+sam][sn]==1) ++count1s;

                if(count1s>med){
                    refs[pop][sn]=1;
                    break;
                }

            } // close for sam


            count1s=0;

        } // close for sn

        iniP+=samplesizes[pop];

	} // close for pops

	return refs;
}

//// This version uses as reference the allele at higher number considering both populations jointly
//
//vector<bool> computereference(const deque < deque<bool> >&datapop){
//
//	int numsnps=datapop[0].size(); // all data has the same length
//	int count1s=0;
//	int totsample = datapop.size(); // the sum of elements in populations
//    int med = totsample/2;
//
//	vector<bool> ref(numsnps,0);
//
//	for(int sn=0; sn<numsnps; ++sn){
//
//		for(int sam=0; sam<totsample; ++sam){
//			if(datapop[sam][sn]==1) ++count1s;
//		} // close for sam
//
//		if(count1s>med) ref[sn]=1;
//        count1s=0;
//	} // close for sn
//
//	return ref;
//}

template<typename TR>
TR bufferAnalyzedSNPs(const vector<string>& Popnames, const vector<string>& rsids, const vector<vector<double>> snpositions,const vector<vector<double>> snpfreqs, string& bufferdata){
        bufferdata="\nPopulation\trsID\tPosition\tFrequency\n";

//cout<<rsids.size()<<" "<<snpositions[0].size()<<" "<<snpfreqs[1].size()<<endl;
        assert(rsids.size()==snpositions[0].size() && rsids.size()==snpfreqs[1].size());

        for(size_t pop=0; pop< Popnames.size(); ++pop){
            bufferdata+= Popnames[pop]+"\n";
            for(size_t row=0; row<rsids.size(); ++row){
                bufferdata+= "\t";
                bufferdata+= rsids[row];
                bufferdata+= "\t";
                bufferdata+= NtoS(snpositions[pop][row],10);
                bufferdata+= "\t";
                bufferdata+= NtoS(snpfreqs[pop][row]);
                bufferdata+="\n";

            }

        }
}
template<typename T=double>
void buffernvds(const vector<string>& Popnames,const int& ws, const vector<string>& rsids,const vector<double>& nvds,const vector<T> metapopsnpositions,const vector<double> metapopsnpfreqs, string& bufferdata){
        bufferdata="\nPopulation\n";

//cout<<rsids.size()<<" "<<snpositions[0].size()<<" "<<snpfreqs[1].size()<<endl;
        assert(rsids.size()==metapopsnpositions.size() && rsids.size()==metapopsnpfreqs.size() && rsids.size()==nvds.size());

        for(auto pop=0; pop< Popnames.size(); ++pop){
            bufferdata+= Popnames[pop]+"\n";
        }
        bufferdata+=NtoS(ws);
        bufferdata+="\n\trsID\tnvd\tPosition\tFrequency\n";
        for(auto row=0; row<rsids.size(); ++row){
            bufferdata+= "\t";
            bufferdata+= rsids[row];
            bufferdata+= "\t";
            bufferdata+= NtoS(nvds[row],10);
            bufferdata+= "\t";
            bufferdata+= NtoS(metapopsnpositions[row],10);
            bufferdata+= "\t";
            bufferdata+= NtoS(metapopsnpfreqs[row]);
            bufferdata+="\n";
        }
}

template<typename T, typename TR>
inline TR writeAnalyzedSNPs(const vector<T>& outnames, const vector<T>&buffers){

    assert(outnames.size() == buffers.size());
    auto file=0;
    for(auto& out: outnames){
        ofstream canalsalida(out.c_str(),ios::app);
        canalsalida<<buffers[file]<<endl;
        canalsalida.close();
        ++file;
    }

}

template<typename T=string, typename TR=void>
inline TR writenvds(const T& outname, const T& buffer){

    ofstream canalsalida(outname.c_str(),ios::app);
    canalsalida<<buffer<<endl;
    canalsalida.close();
}



/**
acr Oct 21
Write the results using the new test that compares partition variances
Oct 2023
Include the new Jeffreys test and others

*/
template<typename T, typename TR>
TR writeVartest(const string& foutname,const vector<string>& names,const vector <vector<ResultIndexes> > indexesfiles,
              vector <vector<ResultFst> > fstfiles, const vector<int>& contSLsig,const double& SL,
              double& SegSites,
              const vector<int>& numoutliers, const vector<int>& totsnps,const vector<int>& numindepsnps,const double& iniMAF,
              const vector<double>& fstmax,const vector<double>& fstmean,const vector<double>& pi0,const vector<double>& LBqval,
              vector<HapMapTable>& snpstable, const size_t& totalmaxcandidates,
              const int& numpops, const vector<vector<vector<int>>>& samplesizeP1, const vector<vector<vector<int>>>& samplesizeP2,
              const vector<vector<VarTest>>&VarTests, const vector<vector<vector<bool>>>& skip, const vector<size_t>& testtype)
    {

    assert(names.size()==fstfiles.size());
    assert(testtype.size()== numpops*indexesfiles.size()); //

 ///totalmaxcandidates is the maximum number of candidates among all runs
/// The use of totalmaxcandidates means that the accumulation is by the rank of the candidate i.e. it accumulate values for the first candidate
/// even if the first is not the same in the different files, similarly for the second candidate and so on...
// This vectors will count the significant obtained for each candidate and the averaged positions
    vector<double> significants(totalmaxcandidates,0);
    vector<double> significantsJ(totalmaxcandidates,0);
    vector<double> averagepos(totalmaxcandidates,0);
    vector<double> averagews(totalmaxcandidates,0);
    vector<size_t> numsigpos(totalmaxcandidates,0);
    vector<size_t> numsruns(totalmaxcandidates,0);
    vector<size_t> candidates(totalmaxcandidates,0);

/// Dec 21 Power and number of tests per population and candidate
    vector<vector<double>>power(numpops,vector<double>(totalmaxcandidates,0));
    size_t numruns = indexesfiles.size();
    vector<vector<size_t>> numrunspersite(numpops,vector<size_t>(totalmaxcandidates,0));
/// Dec 21 computation of the true number of tests regarding the variance analysis just count one for each candidate and run even if there are two in each population
    int totnumruns=numruns;
    for(size_t r =0; r<numruns; ++r){
            bool noPop1=true,noPop2=true;
            size_t numofcandidates =indexesfiles[r].size();
            for(size_t c=0;c<numofcandidates;++c){
                if(!skip[r][0][c]){ // count this run
                    ++numrunspersite[0][c];
                    noPop1=false;
                }
                else
                if(!skip[r][1][c]){ // count this run
                    ++numrunspersite[1][c];
                    noPop2=false;
                }
            } // close for c
            if(noPop1&&noPop2){
                --totnumruns; // discard this run
            }

    } // close for r

//cout<<numrunspersite[0][0]+numrunspersite[1][0]<<"\t"<<totnumruns<<endl;

    vector<bool>PR(numruns,0), PRB(numruns,0), potential(numruns,0); /// PR positive rate, if files are neutral this is FWER if not is the power of detecting at least 1 candidate

	//Jan 22
	vector<double>probcandidate(totalmaxcandidates,0); // the weight of each candidate as the proportion of runs where this candidate is measured.
	vector<size_t>FirstSigCandidateperRun(numruns,0);

	ofstream canalsalida(foutname.c_str(),ios::app);


    if(!indexesfiles.empty()){
/// Oct 2023 by now do not apply Bonferroni so we skip the ending headers appearing in the HacDivSel version
        if(snpstable.empty()){ // not hapmap format
            canalsalida<<"\nFile\tCandidate\tSNP\tPop\tWide\tFST\tFSTindex\tfstPval\tSig\tP1_Size\tP2_Size\tV1\tV2\tVarFtest\tPval\tJHac\tPval\tSig\tVarTestType";
        }
        else // hapmap format
            canalsalida<<"\nFile\tCandidate\trsid\tSNP\tPop\tWide\tFST\tFSTindex\tfstPval\tSig\tP1_Size\tP2_Size\tV1\tV2\tVarFtest\tPval\tJHac\tPval\tSig\tVarTestType";

        canalsalida<<endl;

        for(size_t prun=0;prun<numruns;++prun){

        /// Compute the number of candidates at each run. It is  as simply as indexesfiles[prun].size()
        ///Dec 21
            size_t numcandidates = indexesfiles[prun].size();
            double ws[2]={0,0};

            for(size_t candidate=0; candidate < numcandidates; ++candidate){

        ///
                //candidateindex =  candidate;//indexesfiles[prun][candidate].candidate-1;
                ///acr Oct 21
                bool sig1pop=false;
                int npopspersite = !skip[prun][0][candidate] + !skip[prun][1][candidate]; // for this candidate site number of populations with a minimum of minsamplesize copies of the minor allele

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

                    if(snpstable.empty()){
                        canalsalida<<names[prun]<<"\t"<<candidate+1 <<"\t"<<indexesfiles[prun][candidate].pos;
                    }
                    else{
                        canalsalida<<names[prun]<<"\t"<<candidate+1 <<"\t"<<snpstable[prun][indexesfiles[prun][candidate].pos]<<"\t"<<indexesfiles[prun][candidate].pos;
                    }

                    ws[p]=indexesfiles[prun][candidate].endW[p]-indexesfiles[prun][candidate].iniW[p];
                    canalsalida <<"\t"<<p+1 <<"\t"<<ws[p];

                    canalsalida<<"\t"<<indexesfiles[prun][candidate].fst;

                    canalsalida<<"\t"<<indexesfiles[prun][candidate].fstindex.stat;

                    canalsalida<<"\t"<<indexesfiles[prun][candidate].fstindex.pval<<"\t";

                    if(indexesfiles[prun][candidate].fstindex.pval<=SL){

                        canalsalida<<"*"<<"\t";

                    // average position only if significant
                        averagepos[candidate]+=indexesfiles[prun][candidate].pos;
                        ++numsigpos[candidate];
                    }
                    else
                        canalsalida<<"ns"<<"\t";

                    canalsalida<<samplesizeP1[prun][p][candidate]<<"\t"<<samplesizeP2[prun][p][candidate]<<"\t"<<VarTests[prun][p].v1[candidate]<<"\t"<<VarTests[prun][p].v2[candidate]<<"\t"<<VarTests[prun][p].F[candidate]<<"\t"<<VarTests[prun][p].Fpval[candidate]<<"\t"<<indexesfiles[prun][candidate].JHac.stat<<"\t"<<indexesfiles[prun][candidate].JHac.pval;
/** Dec 21
For Bonferroni we use the number of tests = number of candidates in the run = indexesfiles[prun].size()
**/
                    //int numtests = numcandidates;
                    double pBonf = min(VarTests[prun][p].Fpval[candidate]*numcandidates ,1.0);
                    double pfstBonf = min(indexesfiles[prun][candidate].fstindex.pval*numcandidates ,1.0);

                    if(!skip[prun][p][candidate]){ // for this candidate site the population had a minimum of 2 copies of the minor allele
//cout<<npopspersite<<"\t"<<ws[p]<<"\t"<<averagews[candidate]<<endl;
                        averagews[candidate]+=(ws[p]/npopspersite);
                        /// Count significants
                        if(VarTests[prun][p].Fpval[candidate]<=SL){ // Variance test
                            if(indexesfiles[prun][candidate].fstindex.pval<=SL){ //Fst test
                                if(indexesfiles[prun][candidate].JHac.pval<=SL){
                                    significants[candidate]+=0.5;// because JHac is the same for both populations and this is being computed for each population
                                    significantsJ[candidate]+=0.5;
                                    canalsalida<<"\t***"; // three test were significant

                                    //canalsalida<<"\t"<<pfstBonf<<"\t"<<pBonf;

////                                    if(!sig1pop){// count only 1 significant per population
////                                        significants[candidateindex]+=1;
////                                    /// Dec 21: count the positives just one independently of the number of candidates provided the candidate happens in enough number of runs.
////                                        if(!PR[prun]){ // only once what ever the number of candidates per run
////                                            PR[prun]=1;
////                                            if(pBonf<=SL && pfstBonf<=SL){
////                                                PRB[prun]=1;
////                                                //store the candidate where this happens
////                                                FirstSigCandidateperRun[prun]=candidateindex;
////                                            }
////
////                                        }
////                                        if(pBonf<=SL && pfstBonf<=SL){ // accumulate if there are >1 candidates
////
////                                            power[p][candidateindex]+=1.0;
////                                        }
////                                        sig1pop=true;
////                                    }
//
//                                    if(pBonf<=SL && pfstBonf<=SL){
//                                        canalsalida<<"\t*";
//                                    }
//                                    else{
//                                        canalsalida<<"\tns";
//                                    }
                                } // close J test significant
                                else{
                                    canalsalida<<"\t**";
                                    //significants[candidateindex]+=1;
                                }
                            }// close fst test
                            else if(indexesfiles[prun][candidate].JHac.pval<=SL){
                                canalsalida<<"\t*J*";
                                significants[candidate]+=0.5;
                                significantsJ[candidate]+=0.5;
                            }
                            else{

                                canalsalida<<"\t*";
                                //significants[candidateindex]+=1;
                                //canalsalida<<"\t"<<pfstBonf<<"\t"<<pBonf<<"\tns";
                            }

                        } // close if significant Fpval
                        else{

                            if(indexesfiles[prun][candidate].fstindex.pval<=SL){
                                if(indexesfiles[prun][candidate].JHac.pval<=SL){
                                    canalsalida<<"\t*J*";
                                    significants[candidate]+=0.5;
                                    significantsJ[candidate]+=0.5;
                                }
                                else
                                    canalsalida<<"\t*";

                                if(VarTests[prun][p].v1[candidate] < VarTests[prun][p].v2[candidate] )
                                    if(!potential[prun] && !PRB[prun])
                                        potential[prun]=1;
                            } // close if fstindex.paval
                            else if(indexesfiles[prun][candidate].JHac.pval<=SL){
                                canalsalida<<"\tJ*";
                                significantsJ[candidate]+=0.5;

                            }
                            else
                                canalsalida<<"\tns";
                            //canalsalida<<"\t"<<pfstBonf<<"\t"<<pBonf<<"\tns";
                        }

                    } // close if not skiprun
                    else{//skip the run for the var test of this candidate and population but not for the Fst and JHac tests
                        if(indexesfiles[prun][candidate].fstindex.pval<=SL){
                            if(indexesfiles[prun][candidate].JHac.pval<=SL){
                                canalsalida<<"\t*J*";
                                significantsJ[candidate]+=0.5;
                            }
                        }
                        else
                         if(indexesfiles[prun][candidate].JHac.pval<=SL){
                                canalsalida<<"\tJ*";
                                significantsJ[candidate]+=0.5;
                        }
                        else
                            canalsalida<<"\t---";
                        //canalsalida<<"\t"<<"---"<<"\t"<<"---"<<"\t---";
                    }

                    if(p==0){
                        switch(testtype[numpops*prun]){ // 1:Bartlett 2:Brown-Forsythe 3:F-test 4:Levene 5:LeveneHSK2 6:O'Brian
                            case 0: canalsalida<<"\t-------"; break;
                            case 1: canalsalida<<"\tBartlett"; break;
                            case 2: canalsalida<<"\tBForsythe"; break;
                            case 3: canalsalida<<"\tF-test"; break;
                            case 4: canalsalida<<"\tLevene"; break;
                            case 5: canalsalida<<"\tLeveneHSK2"; break;
                            case 6: canalsalida<<"\tO'Brian"; break;
                        }
                    }
                    else{
                        switch(testtype[numpops*prun+p]){ // 1:Bartlett 2:Brown-Forsythe 3:F-test 4:Levene 5:LeveneHSK2 6:O'Brian
                            case 0: canalsalida<<"\t-------"; break;
                            case 1: canalsalida<<"\tBartlett"; break;
                            case 2: canalsalida<<"\tBForsythe"; break;
                            case 3: canalsalida<<"\tF-test"; break;
                            case 4: canalsalida<<"\tLevene"; break;
                            case 5: canalsalida<<"\tLeveneHSK2"; break;
                            case 6: canalsalida<<"\tO'Brian"; break;
                        }
                    }
                    canalsalida<<endl;

                } // close for pops


                ///Dec 21
                if(npopspersite)
                    ++numsruns[candidate];

//cout<<probcandidate[prun]<<endl;
                candidates[candidate]=indexesfiles[prun][candidate].candidate; // this is equal to the for index candidate+1


            } // close for candidates
        } // close for pruns

        // mean number of segregating sites
        SegSites/=numruns;

        for(size_t candidate=0; candidate < significants.size(); ++candidate){

            significants[candidate]=100.0*significants[candidate]/numsruns[candidate]; // the denominator is for avoid cases with very low partition 2 size that provoke overdispersion in the variance test
            significantsJ[candidate]=100.0*significantsJ[candidate]/numruns; // here there is no partitions it just compares one population with the other
            averagews[candidate]/=numsruns[candidate];
            if(numsigpos[candidate])
                averagepos[candidate]/=numsigpos[candidate];
            else
                averagepos[candidate]=-1;
        }

//cout<<numsruns[0]<<"\t"<<totnumruns<<endl;
//assert(numsruns[0]==totnumruns);
//exit(0);
        canalsalida<<endl;
        if(numruns>1){

            canalsalida<<"SegSites\tCandidate\tProbCandidate\t#Runs\tMean_size\tMean_pos\t%Power\t%PowerJ" <<endl;

            canalsalida<<SegSites;

            /// Dec 21 Compute the power counting the non-discarded runs.
            ///Note that for each candidate the significants was count just once even if they happen in both populations, so the total power is the sum
            vector<double> totpower(totalmaxcandidates,0);
            vector<size_t>totrunspersite(totalmaxcandidates,0);
            for(size_t po=0;po<totalmaxcandidates;++po){
                totrunspersite[po]=numrunspersite[0][po]+numrunspersite[1][po];
                totpower[po]=(power[0][po]+power[1][po]) / (totrunspersite[po]);

                probcandidate[po] = (numrunspersite[0][po]+numrunspersite[1][po])/(double)totnumruns;
            }



            for(size_t candidate=0; candidate < significants.size(); ++candidate){
            // averages
                canalsalida<<"\t"<<candidates[candidate]<<"\t"<<probcandidate[candidate] <<"\t"<<totrunspersite[candidate]<<"\t"<<averagews[candidate]<<"\t"<<averagepos[candidate];
                canalsalida<<"\t"<<significants[candidate]<<"\t"<<significantsJ[candidate];
                canalsalida<<endl;
            }

            /// The maximum power is
            auto itermaxpower = max_element(begin(totpower),end(totpower));
            int candidatemax = distance(begin(totpower),itermaxpower);
            double maxpower= *itermaxpower;
//            auto p1 = max_element(begin(power[0]),end(power[0]));
//            auto p2 = max_element(begin(power[1]),end(power[1]));
//            int maxsite, maxp;
//            if(*p1>=*p2){
//                maxp=0;
//                maxsite = distance(begin(power[0]), p1);
//                maxpower =*p1;
//
//            }
//            else{
//                maxp=1;
//                maxsite = distance(begin(power[1]), p2);
//                maxpower =*p2;
//            }
//
//            maxpower/= numrunspersite[maxp][maxsite];

            double sumPR=0, sumPRB=0, sumPot=0;

            for(size_t f=0; f< numruns; ++f){
                //size_t candidatesig = FirstSigCandidateperRun[f];
                sumPR+=PR[f];
                sumPRB+=PRB[f];
                if(!PRB[f]){
                    double maxprobcandidate= *max_element(begin(probcandidate),end(probcandidate));
                    sumPot+=potential[f]*maxprobcandidate;
                }
            }

            sumPR/= (totnumruns); //
            sumPRB/= (totnumruns);
            sumPot/= (totnumruns);
            canalsalida<<"TotRuns\t%PR\t%PRBonf\t%MaxPowerBonf\tMax_Candidate_Rank\tPotential"<<endl;
            //canalsalida<<100*sumPR<<"\t"<<100*sumPRB<<"\t"<<100*maxpower<<"\t"<<candidatemax<<"\t"<<100*sumPot<<endl;
            canalsalida<<totnumruns<<"\t"<<100*sumPR<<"\t"<<100*sumPRB<<"\t"<<100*maxpower<<"\t"<<candidatemax<<"\t"<<100*sumPot<<endl;

        }
    } // close indexesfile not empty

}


/**
acr Oct 23
Write the results using the new JHac test (using Jeffreys divergence)
Indicate also the results of eos which are in the fstfiles
*/
template<typename T, typename TR>
TR writeJtest(const vector<bool>& highmeanFST,const vector<int>&totsnps,const vector<double>& fstmean, const size_t& totruns,const string& foutname,const vector<string>& names,const vector <vector<ResultIndexes> > indexesfiles,
              vector <vector<ResultFst> > fstfiles, const vector<int>& contSLsig,const double& SL,
              double& SegSites,vector<HapMapTable>& snpstable, const size_t& totalmaxcandidates,
              const int& numpops, const vector<vector<int>>& samplesizePops)
    {

    assert(names.size()==fstfiles.size());


 ///totalmaxcandidates is the maximum number of candidates among all runs
/// The use of totalmaxcandidates means that the accumulation is by the rank of the candidate i.e. it accumulate values for the first candidate
/// even if the first is not the same in the different files, similarly for the second candidate and so on...
// This vectors will count the significant obtained for each candidate and the averaged positions when there are several files analysed at once (like simulations)
    vector<double> significantsFst(totalmaxcandidates,0);
    vector<double> significantsJ(totalmaxcandidates,0);
    vector<double> FPJ(totalmaxcandidates,0); // false positives considering also the runs that had no blocks i.e. without selection, this correct a problem of low runs that could inflate false positives when many files have no blocks (and so no selection) and the program do not analyze them
    vector<double> averagepos(totalmaxcandidates,0);
    vector<double> averagews(totalmaxcandidates,0);

    vector<size_t> numsigpos(totalmaxcandidates,0);
    vector<size_t> numsruns(totalmaxcandidates,0);
    vector<size_t> candidates(totalmaxcandidates,0);

    size_t numruns = indexesfiles.size();
    vector<vector<size_t>> numrunspersite(numpops,vector<size_t>(totalmaxcandidates,0));


	ofstream canalsalida(foutname.c_str(),ios::app);

    string signfile = substring(foutname, ".", true,true); //
    signfile+="_SIG.tsv";
    ofstream canalsalidaSIG(signfile.c_str(),ios::app);

    if(!indexesfiles.empty()){

    canalsalida<<setprecision(outputprecision);
    canalsalidaSIG<<setprecision(outputprecision);
/// Oct 2023 now do not need to apply Bonferroni so we skip the ending headers appearing in the HacDivSel version
        if(snpstable.empty()){ // not hapmap format
            canalsalida<<"File\t#SNPs\tCandidate\tSNP\tWide\tIni\tEnd\tNCLASSES\tMeanFST\tFST\tFSTindex\tfstPval\tn1\tn2\tJHac\tPval\tSig";
            canalsalidaSIG<<"File\t#SNPs\tCandidate\tSNP\tWide\tIni\tEnd\tNCLASSES\tMeanFST\tFST\tFSTindex\tfstPval\tn1\tn2\tJHac\tPval\tSig";
        }
        else{ // hapmap format
            canalsalida<<"File\t#SNPs\tCandidate\trsid\tPos\tWide\tIni\tEnd\tNCLASSES\tMeanFST\tFST\tFSTindex\tfstPval\tn1\tn2\tJHac\tPval\tSig";
            canalsalidaSIG<<"File\t#SNPs\tCandidate\trsid\tPos\tWide\tIni\tEnd\tNCLASSES\tMeanFST\tFST\tFSTindex\tfstPval\tn1\tn2\tJHac\tPval\tSig";
        }
        canalsalida<<endl;
        canalsalidaSIG<<endl;
        bool sig=false;
        for(size_t prun=0;prun<numruns;++prun){

        ///April 2024
            if(highmeanFST[prun]){

                canalsalida<<"\nWARNING: THE AVERAGE FST WAS GREATER THAN 0.5 IN FILE "<<names[prun]<<"\n"<<endl;
                canalsalidaSIG<<"\nWARNING: THE AVERAGE FST WAS GREATER THAN 0.5 IN FILE "<<names[prun]<<"\n"<<endl;

            }

        /// Compute the number of candidates at each run. It is  as simply as indexesfiles[prun].size()
        ///Dec 21
            size_t numcandidates = indexesfiles[prun].size();


            for(size_t candidate=0; candidate < numcandidates; ++candidate){

        ///Feb 24
                if(indexesfiles[prun][candidate].JHac.pval<=SL)
                    sig=true;
                else
                    sig=false;
                //candidateindex =  candidate;//indexesfiles[prun][candidate].candidate-1;
                ///acr Oct 21


                if(snpstable.empty()){
                    canalsalida<<names[prun]<<"\t"<<totsnps[prun]<<"\t"<<candidate+1 <<"\t"<<indexesfiles[prun][candidate].pos;

                }
                else{
                    canalsalida<<names[prun]<<"\t"<<totsnps[prun]<<"\t"<<candidate+1 <<"\t"<<snpstable[prun][indexesfiles[prun][candidate].pos]<<"\t"<<indexesfiles[prun][candidate].pos;

                }


                canalsalida <<"\t"<<indexesfiles[prun][candidate].wide;

                /// the window size is the same in both populations
                canalsalida <<"\t"<<indexesfiles[prun][candidate].iniW[0] <<"\t"<< indexesfiles[prun][candidate].endW[0];

                canalsalida<<"\t"<<indexesfiles[prun][candidate].nclasses;

                canalsalida<<"\t"<<fstmean[prun]; ///April 2024

                canalsalida<<"\t"<<indexesfiles[prun][candidate].fst;


                canalsalida<<"\t"<<indexesfiles[prun][candidate].fstindex.stat;


                canalsalida<<"\t"<<indexesfiles[prun][candidate].fstindex.pval<<"\t";


                if(indexesfiles[prun][candidate].fstindex.pval<=SL){

                    //canalsalida<<"*"<<"\t";
                    significantsFst[candidate]+=1;

                }
                //else
                    //canalsalida<<"ns"<<"\t";

                canalsalida<<samplesizePops[prun][0]<<"\t"<<samplesizePops[prun][1]<<"\t"<<indexesfiles[prun][candidate].JHac.stat<<"\t"<<indexesfiles[prun][candidate].JHac.pval;

                 if(sig){
                        canalsalida<<"\t*";
                        significantsJ[candidate]+=1;
                        // average position only if significant J
                        averagepos[candidate]+=indexesfiles[prun][candidate].pos;
                        averagews[candidate]+=indexesfiles[prun][candidate].wide;
                        ++numsigpos[candidate];
                }
                else{
                    canalsalida<<"\tns";
                }

                canalsalida<<endl;
                ///From april 2024 esto ya no es cierto: indexesfiles[prun][candidate].candidate; does not correspond with the natural order because indexesfiles was sorted for highest JHac values
                /// Lo anterior ya no es cierto porque
                /// indexesfiles ya no se ordena por  los valores J más altos se deja el orden natural de posición (más cómodo para el usuario)
                //candidates[candidate]=indexesfiles[prun][candidate].candidate; // this is equal to the for index candidate+1
                candidates[candidate]=candidate+1;
                ///Feb 2024
                if(sig){ // write a file with only significants
                    if(snpstable.empty()){
                        canalsalidaSIG<<names[prun]<<"\t"<<totsnps[prun]<<"\t"<<candidate+1 <<"\t"<<indexesfiles[prun][candidate].pos;

                    }
                    else{
                        canalsalidaSIG<<names[prun]<<"\t"<<totsnps[prun]<<"\t"<<candidate+1 <<"\t"<<snpstable[prun][indexesfiles[prun][candidate].pos]<<"\t"<<indexesfiles[prun][candidate].pos;

                    }

                    canalsalidaSIG <<"\t"<<indexesfiles[prun][candidate].wide;

                /// the window size is the same in both populations
                    canalsalidaSIG <<"\t"<<indexesfiles[prun][candidate].iniW[0] <<"\t"<< indexesfiles[prun][candidate].endW[0];

                    canalsalidaSIG<<"\t"<<indexesfiles[prun][candidate].nclasses;

                    canalsalidaSIG<<"\t"<<fstmean[prun]; ///April 2024

                    canalsalidaSIG<<"\t"<<indexesfiles[prun][candidate].fst;


                    canalsalidaSIG<<"\t"<<indexesfiles[prun][candidate].fstindex.stat;


                    canalsalidaSIG<<"\t"<<indexesfiles[prun][candidate].fstindex.pval<<"\t";

                    canalsalidaSIG<<samplesizePops[prun][0]<<"\t"<<samplesizePops[prun][1]<<"\t"<<indexesfiles[prun][candidate].JHac.stat<<"\t"<<indexesfiles[prun][candidate].JHac.pval;


                    canalsalidaSIG<<"\t*";

                    canalsalidaSIG<<endl;
                }// close if sig

            } // close for candidates
        } // close for pruns

        // mean number of segregating sites
        SegSites/=numruns;

        for(size_t candidate=0; candidate < totalmaxcandidates; ++candidate){

            significantsFst[candidate]=100.0*significantsFst[candidate]/numruns;

            /// Enero 2024: ¿cambia el divisor numruns por totruns para incluir aquellos ficheros en los que no se contaron solo porque no tenían bloques aunque sí suf snps
            FPJ[candidate]=100.0*significantsJ[candidate]/totruns;

            significantsJ[candidate]=100.0*significantsJ[candidate]/numruns;

            if(numsigpos[candidate]){
                averagepos[candidate]/=numsigpos[candidate];
                averagews[candidate]/=numsigpos[candidate];
            }
            else
                averagepos[candidate]=-1;
        }

//cout<<numsruns[0]<<"\t"<<totnumruns<<endl;
//assert(numsruns[0]==totnumruns);
//exit(0);
        canalsalida<<endl;
        if(numruns>1){

            canalsalida<<"SegSites\tCandidate\tMean_size\tMean_pos\t%PowerJ\t%FPJJ" <<endl;

            canalsalida<<SegSites;

            for(size_t candidate=0; candidate < totalmaxcandidates; ++candidate){
            // averages
                canalsalida<<"\t"<<candidates[candidate]<<"\t"<<averagews[candidate]<<"\t"<<averagepos[candidate];
                canalsalida<<"\t"<<significantsJ[candidate]<<"\t"<<FPJ[candidate];
                canalsalida<<endl;
            }




        } // close if num runs >1
    } // close indexesfile not empty

    //canalsalida<<totruns<<endl;

}



template<typename TR>
TR writeFSTs(const string& foutname,const vector<string>& names,const vector <vector<ResultFst> > fstfiles0,
              const vector <vector<ResultFst> > fstfiles,const double& SL,
              const vector<int>& totsnps, const double& iniMAF)
    {

    assert(names.size()==fstfiles.size());

    //int inisite=0,endsite=0,previni=-1, prevend=-1,L=0,sgn=MISSINGVAL;
    //double cSL=SL, Meff=0,gamm=0, Ngvd=0, meanD=0;

	ofstream canalsalida(foutname.c_str(),ios::app);


    if(!fstfiles0.empty() || !fstfiles.empty() ){
        auto fcont = 0;
        for(auto fst: fstfiles)
            if(!fst.empty())
                ++fcont;
        double pc = (double)fcont/max(names.size(),size_t(1));
        pc*=100;
        canalsalida<<"MAF = "<<iniMAF<<endl;
        canalsalida<<"% FILES WITH K-MEANS EXTREME POSITIVE OUTLIERS = "<<NtoS(pc) <<"%\n";
    }

    for(size_t prun=0;prun<fstfiles0.size();++prun){
        canalsalida<<"\n****************************************************"<<endl;
        canalsalida<<"FILE "<<names[prun]<<"\n";
        canalsalida<<"****************************************************\n"<<endl;

        canalsalida<<"#SNPs"<<"\t"<<totsnps[prun]<<"\n";

        canalsalida<<"\n*************** Non-EOS FST OUTLIERS"<<endl;

//assert((void("Failed "+NtoS(fstfiles0[prun].size())+" vs "+NtoS((size_t)numoutliers[prun]-fstfiles[prun].size())),fstfiles0[prun].size() == (size_t)numoutliers[prun]-fstfiles[prun].size()));

        canalsalida<<fstfiles0[prun].size()<<endl;

        if(!fstfiles0[prun].empty())
            canalsalida<<"\n\t\t\tID\tPos\tFST\tpval(LK test)";

        for(int f =fstfiles0[prun].size()-1; f>=0; --f){

            canalsalida<<endl<<"\t\t\t";
            canalsalida<<setprecision(outputprecision);

            canalsalida<<fstfiles0[prun][f].name<<"\t"<<fstfiles0[prun][f].pos<<"\t"<<fstfiles0[prun][f].fst<<"\t"<<fstfiles0[prun][f].pval;



        }
        canalsalida<<endl;

        canalsalida<<"\n*************** EOS FST OUTLIERS"<<endl;

        //canalsalida<<"\n\t"<<"\t"<<"#SNPs"<<"\t"<<"#ExtremeOutliers"<<"\t"<<"Name\tSNP\tFST"<<endl;
       if(fstfiles[prun].empty()){
          canalsalida<<fstfiles[prun].size()<<endl;
          continue;
        }

        canalsalida<<fstfiles[prun].size()<<"\n";
        canalsalida<<"\n\t\t\tID\tPos\tFST\tpval(LK test)";
        for(int f =fstfiles[prun].size()-1; f>=0; --f){


            canalsalida<<endl<<"\t\t\t";
            canalsalida<<fstfiles[prun][f].name<<setprecision(outputprecision)<<"\t"<<fstfiles[prun][f].pos<<"\t"<<fstfiles[prun][f].fst<<"\t"<<fstfiles[prun][f].pval;
        }
        canalsalida<<endl;

    }


}


template <typename T>
void writeHACDistrib(const vector<T>& HAC1,const vector<T>& HAC2){

	ofstream canalsalida("HACDistrib.dat",ios::app);
	int rangemin =min(HAC1.size(),HAC2.size());
    int rangemax =max(HAC1.size(),HAC2.size());

    for(int h=0;h<rangemax;++h)
        if(HAC1.size()>=HAC2.size())
            if(h<HAC2.size())
                canalsalida<<HAC1[h]<<"\t"<<HAC2[h]<<endl;
            else
                canalsalida<<HAC1[h]<<endl;
        else
            if(h<HAC1.size())
                canalsalida<<HAC1[h]<<"\t"<<HAC2[h]<<endl;
            else
                canalsalida<<"\t\t"<<HAC2[h]<<endl;
    canalsalida<<endl;
    canalsalida.close();
}




