/** *************************************************************************************
 * May you live in interesting times!
 *
 *  PROGRAM: Useful template library
 *
 * 	http://acraaj.webs.uvigo.es
 *
 *  AUTHOR:
 *
 * 	Antonio Carvajal Rodriguez (acraaj@uvigo.es) April 2018
 *
 * 	Area de Genetica. Universidad de Vigo.
 *
 *
 * *************************************************************************************
 *
 * 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:
 *
 *  30/01/2019:
 *  1) Concept added: #define tstring typename
 *  2) Two new function template added toUpperCase and toLowerCase
 * ************************************************************************************* */


#pragma once


#ifndef TLIB_HPP
#define TLIB_HPP

#include <iostream>
#include <iomanip>
#include <iterator>
#include <string>
#include <fstream>
#include <sstream>
#include <stdarg.h>
#include <vector>
#include <deque>
#include <array>
#include <cmath>
#include <ctime>
#include <limits>
#include <numeric>
#include <random>
//#include <cfloat>
#include <algorithm>
#include <cstdlib>
#include <errno.h>
#include<chrono>
#include <typeinfo>
#include <ctype.h>

#define tstring typename

//#include "Traits.hpp"
const unsigned SEED=777;
const unsigned outputprecision{12};
const int MISSINGVAL = -99;
const int openfiletries{1};
const double dPI  =3.141592653589793238463;
const double raiz2 = sqrt(2.0);
const double kn = sqrt(2.0/dPI);
const double quasizero = std::numeric_limits<double>::min();
const double quasione = 0.9999999;


using namespace std;

/// Necessary declarations

template <typename T>
inline void printval(const T& v,const string& msg ,const string& name);

template <typename T>
inline void printcol(const T& v);
//
//template <typename T>
//inline double probks(const T& alam);

template <typename Type>
inline
const Type infinito()
{
    return numeric_limits<Type>::max();
}

template <typename Type>
inline
const Type toosmall()
{
    return std::numeric_limits<Type>::min();
}


template <typename T> int sign(T val) {
    return (T(0) < val) - (val < T(0));
}

// integer power because the stl pow with integer does not behaves well requires rounding e.g pow(10,2) = 99.99999 the int is 99 instead of 100
template<typename T=long int>
T pow(const T& b, const T& ne, bool error){

    error=false;
    if(b==0){

        if(ne>0)
            return 0;
        else{
            error=true;
            return toosmall<T>();
        }

    }
    else if(b==1)
        return 1;
    else if(b==-1)
        return(ne%2==0? 1: -1);
    else if(ne<0){
        error=true;
        return 0;
    }
    else{ // b>=2 n>=0
        T p=1;
        T n=ne;
        while(n>0){
            if( (p<0? -p:p) >= infinito<T>() / (b<0? -b:b)) // bas values of p and b. Manage overflow
            {
                error=true;
                return 0;
            }
            else{
                p*=b;
                --n;
            }

        } // close while
        return p;

    }

}

template <typename T= double>
T gammln(const T& xx){
	int j;
	T x,y,tmp,ser;
	static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
	-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};

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

/// Factorial and combinatory algorithms
/// NRC versions use modification of factrl NRC (p.219) algorithm which statically stores some values and calling gamml function
/// Thus if the same computations are expected to occur several times NRC versions should be used
template<typename T=long double>
inline T factorialNRC(const int& n){

    const size_t tabmax=101;
    static size_t ntop=7;
    static T a[tabmax]={1.0,1.0,2.0,6.0,24.0,120.0,720.0, 5040.0};

    if(n<0){
        printval<int>(1,"Negative factorial in routine factorialNRC","Errfact.log");
        exit(1);
    }

    if(n<=ntop) return a[n];

    if(n>(tabmax-1)) return exp(gammln<T>(n+1.0));

    size_t j;
	while(ntop<n){
		j=ntop++; // puts j equal to ntop then increases ntop
		a[ntop]=a[j]*ntop;
	}

    return a[n];
}

// This version just incorporate the table but do not call to gammln
template<typename T=long double>
inline T factorial_tab(const int& n){
    const size_t tabmax=101;
    static size_t ntop=7;
    static T a[tabmax]={1.0,1.0,2.0,6.0,24.0,120.0,720.0, 5040.0};

    if(n<0){
        printval<int>(1,"Negative factorial in routine factorialNRC","Errfact.log");
        exit(1);
    }

    size_t j;
    T an;
	if(n<tabmax){
        while(ntop<n){
            j=ntop++; // puts j equal to ntop then increases ntop
            a[ntop]=a[j]*ntop;
        }
        an = a[n];
	}
    else{
        while(ntop<(tabmax-1)){
            j=ntop++; // puts j equal to ntop then increases ntop
            a[ntop]=a[j]*ntop;
        }
        an=a[ntop];
        for(auto i=tabmax; i<=n; ++i)
            an*=i;
    }
    return an;
}

template<typename T=long double>
inline T combinatory_tab(const int& n, const int& i){

    if(n<0 || i<0){
        printval<int>(1,"Negative factorial in routine combinatory","Errfact.log");
        exit(1);
    }

    T num(1);
    T den(i);

    if(n<i) return 0;
    if(n==i || i==0 || n<2) return num;
    if(i==1 || i==n-1) return n;

    int ti = i;
    int dif =n-i;
    if(dif<i) swap(dif,ti);


    for(auto j=dif+1; j<=n; ++j)
        num*=j;
    den = factorial_tab(ti);


    return num/den;

}

template<typename T=long double>
inline T combinatoryNRC(const int& n, const int& i){

    if(n<0 || i<0){
        printval<int>(1,"Negative factorial in routine combinatory","Errfact.log");
        exit(1);
    }

    T num(1);
    T den(i);

    if(n<i) return 0;
    if(n==i || i==0 || n<2) return num;
    if(i==1 || i==n-1) return n;

    int ti = i;
    int dif =n-i;
    if(dif<i) swap(dif,ti);


    for(auto j=dif+1; j<=n; ++j)
        num*=j;
    den = factorialNRC(ti);


    return num/den;

}

template<typename T=long double>
inline T factorial(const int& n){
    if(n<0){
        printval<int>(1,"Negative factorial in routine factorial","Errfact.log");
        exit(1);
    }

    T x(1);
    for(auto i=2; i<=n; ++i)
        x*=i;
    return x;
}

template<typename T=long double>
inline T combinatory(const int& n, const int& i){

    if(n<0 || i<0){
        printval<int>(1,"Negative factorial in routine combinatory","Errfact.log");
        exit(1);
    }

    T num(1);
    T den(i);

    if(n<i) return 0;
    if(n==i || i==0 || n<2) return num;
    if(i==1 || i==n-1) return n;

    int ti = i;
    int dif =n-i;
    if(dif<i) swap(dif,ti);


    for(auto j=dif+1; j<=n; ++j)
        num*=j;
    den = factorial(ti);


    return num/den;

}

// Converts a number to string and return the string
template<typename T>
inline string NtoS(const T& number){
	ostringstream snumber;
	snumber << number;
	return snumber.str();
}
template<typename T>
inline string NtoS(const T& number, const size_t precision){
	ostringstream snumber;
	snumber<<setprecision(precision);
	snumber << number;
	return snumber.str();
}
//Search value within non-sorted  v. It returns the index if v[index]==value or -1 if value is not in v
template <typename T, typename Col=vector<T>>
inline int Tfindindx(const Col& v,const T& value){

	if(v.empty()) return -1;

    int dist=distance(v.begin(),find(v.begin(),v.end(),value));

	int res = (dist<(int)v.size())? dist: -1;

	return res;
}


//Search value within sorted  v. ** USE Tfindindx if v is non-sorted **. It returns the index if v[index]==value or -1 if value is not in v
template <typename T, typename Col=vector<T>>
inline int get_indx(const Col& v,const T& value){

	if(v.empty()) return -1;



    if(binary_search (v.begin(), v.end(), value)){

        return (lower_bound (v.begin(), v.end(), value) - v.begin() );

    }

	return -1;
}

// Given an index corresponding to a value x in a collection v1 get the index corresponding to x in another ** sorted ** collection v2
// It returns the index if v2[index]==x or -1 if value is not in v2
//** USE Tfindindxtoindx if v2 is non-sorted **.

template <typename T, typename Col=vector<T>>
inline int get_indxtoindx(const Col& v1,const int& indx, const Col& v2){

	if(v1.empty() || v2.empty()) return -1;
	if((int)v1.size()<= indx || indx<0) return -1;

    if(binary_search (v2.begin(), v2.end(), v1[indx])){

        return (lower_bound (v2.begin(), v2.end(), v1[indx]) - v2.begin() );

    }

	return -1;
}

// Given an index corresponding to a value x in a collection get the index corresponding to x in another ** non-sorted ** collection
// It returns the index if v2[index]==x or -1 if value is not in v2
//** USE get_indxtoindx if v2 is sorted **.

template <typename T, typename Col=vector<T>>
inline int Tfindindxtoindx(const Col& v1,const int& indx, const Col& v2){

	if(v1.empty() || v2.empty()) return -1;
	if(v1.size()<= indx || indx<0) return -1;

    int dist=distance(v2.begin(),find(v2.begin(),v2.end(),v1[indx]));

	int res = (dist<v2.size())? dist: -1;

	return res;
}


// the maximum value and its position
template <typename tsize=size_t, typename T, typename Tcol=vector<T>>
inline bool maxelementandpos(const Tcol& v,T& maxval,tsize &position){

	if(v.empty()) return false;

	position = distance(v.begin(),max_element(v.begin(),v.end()));
	maxval=v[position];

	return true;
}

//// The n higher elements and their positions, the shift constant is for the case when the positions do not matches the index of the vector
template <typename tsize=size_t, typename T, typename Tcol=vector<T>>
inline bool maxelementandpos(const Tcol& v,const tsize n,const tsize inishift,const tsize shift, Tcol& maxvals,vector<tsize> &positions){

	if(v.empty()) return false;

    Tcol temp= Tcol(v);
    vector<tsize> temppos; // store positions for avoiding repetition that could occur if an element is repeated but the find function always return the first encountered
    tsize position,index;

    for(auto times=0; times<n; ++times){
        tsize numiters=0;
        auto itermaxnvd = max_element(begin(temp),end(temp));
        position = distance(temp.begin(),itermaxnvd); // position gives the number of elements between begin and the max_element
        maxvals[times]=temp[position];
// the problem to solve is that if the same maximum element is repeated the find algorithm applied to the original vector will return always the same first position of that element
        auto iterini=v.begin();
        int repeats=0;
        do{
            auto iter = find(iterini,v.end(),maxvals[times]);
            index = iter-v.begin();
            if(numiters==0){
                repeats= count(temppos.begin(),temppos.end(),index); //count the number of times that index appeared before this loop
                temppos.push_back(index);
            }

            iterini=iter+1;
            ++numiters;

        }while(repeats>=numiters); // if this index is repeated get the next position
                        // firstsnp + index*shift
        positions[times]= inishift + index*shift; // the position of this element in the original sequence i.e. if shift=stepsize=50 (window 100 or 101) the element with index =0 in the set of original snps corresponds to the position 49 (the center of the window size) and the index=1 correspond to 99. If shift=1 nothing changes and the index corresponds to the real position
        temp.erase(temp.begin()+position); //
//cout<<positions[times]<<"\t"<<maxvals[times]<<endl;
    }
//printcol(v);
//exit(1);

	return true;
}


template <typename T, typename Tcol=vector<T>>
inline bool minelementandpos(const Tcol& v,T& minval,int &position){

	if(v.empty()) return false;

	position = distance(v.begin(),min_element(v.begin(),v.end()));
	minval=v[position];

	return true;
}




// open a file and left it opened in the stream datafile
template <typename T=bool>
T openfile(const string& path, const string& filename0,ifstream& datafile, bool boolalfa){

    /// Open file

    int tries=0;

	string filename=path+filename0;
	if(datafile.is_open())
        datafile.close();

	// check input file name
	do{
		datafile.clear();
		datafile.open(filename.c_str());
		if(!datafile.fail())
		{
			tries=0;
			break; // the open file operation succeed
		}
		else{
            cout<<"\nCould not open file: "<<filename<<endl;
			switch(tries){
				case 0:	transform(filename.begin(),filename.end(), filename.begin(), (int(*)(int)) tolower);
cout<<"\nTRYING FILE NAME "<<filename<<endl;
						tries++;
						break;
				case 1: transform(filename.begin(),filename.end(), filename.begin(), (int(*)(int)) toupper);
cout<<"\nTRYING FILE NAME "<<filename<<endl;
						tries++;
						break;
				default: tries++;
			}

		}

	}while(tries<openfiletries);

	if(tries>=openfiletries)
	{
		string msg = "\nERROR: can not find input file after "+ NtoS(tries+1)+" tries\n";
		cout<<msg<<endl;
		return false;

	}

    if(boolalfa) // set the boolalpha flag (for reading bool values as their textual representations true or false instead of integer 1 or 0
        datafile >> boolalpha; // what is changed is the str datafile

    datafile.clear(); // set the internal error state flags
    return true;
}

// Data is changed by sorting or nth_element

template <typename Tel=double,typename Tcol=vector<double>>
inline Tel odd_median(Tcol& seq)
{
    size_t n = seq.size() / 2;

    nth_element(seq.begin(), seq.begin()+n, seq.end());
    return seq[n];
}

template <typename Tel=double,typename Tcol=vector<double>>
inline Tel median(Tcol& seq)
{
    if(seq.size()%2!=0) return odd_median<Tel>(seq);
    else{

        size_t n=seq.size()/2;
        sort(seq.begin(),seq.end());
        return((seq[n-1]+seq[n])/2);
    }

}

// Data are left constant

template <typename Tel=double,typename Tcol=vector<double>>
inline Tel c_odd_median(const Tcol& seq)
{
    size_t n = seq.size() / 2;
    Tcol temp = seq;
    nth_element(temp.begin(), temp.begin()+n, temp.end());
    return temp[n];
}

// Data are left constant
template <typename Tel=double,typename Tcol=vector<double>>
inline Tel c_median(const Tcol& seq)
{
    if(seq.size()%2!=0) return c_odd_median<Tel>(seq);
    else{
        Tcol temp = seq;
        size_t n=seq.size()/2;
        sort(temp.begin(),temp.end());
        return((temp[n-1]+temp[n])/2);

    }

}

template <typename Tel,typename Tcol>
inline void quartiles(const Tcol & seq, Tel& Q1, Tel& Q2, Tel& Q3){
    Tcol temp = seq;
    auto const pQ1 = temp.size() / 4;
    auto const pQ2 = temp.size() / 2;
    auto const pQ3 = pQ1 + pQ2;
    std::nth_element(temp.begin(),          temp.begin() + pQ1, temp.end());
    std::nth_element(temp.begin() + pQ1 + 1, temp.begin() + pQ2, temp.end());
    std::nth_element(temp.begin() + pQ2 + 1, temp.begin() + pQ3, temp.end());

    Q1 = temp[pQ1];
    Q2 = temp[pQ2];
    Q3 = temp[pQ3];
}

template <typename Tel,typename Tcol, typename iT>
inline Tel quartil(const Tcol & seq, const iT& q){
    Tcol temp= seq;
    auto pQ=iT();
    switch(q){
        case 1: pQ = temp.size() / 4; break;
        case 2: pQ = temp.size() / 2; break;
        case 3: pQ = temp.size() / 4 + temp.size() / 2; break;
        default: exit(4);
    }

    std::nth_element(temp.begin(), temp.begin() + pQ, temp.end());

    return(temp[pQ]);

}

// IQR: interquartile range. The seq is partially sorted as done by nth_element
template <typename Tel,typename Tcol>
inline void iqr(Tcol & seq, Tel& Q1, Tel& Q3, Tel& IQR, int& posQ1, int& posQ3){

    auto const pQ1 = seq.size() / 4;
    auto const pQ2 = seq.size() / 2;
    auto const pQ3 = pQ1 + pQ2;
    std::nth_element(seq.begin(),          seq.begin() + pQ1, seq.end());
    std::nth_element(seq.begin() + pQ1 + 1, seq.begin() + pQ3, seq.end());

    Q1 = seq[pQ1];
    Q3 = seq[pQ3];

    /// The interquartile range is
    IQR = Q3-Q1;

    posQ1= pQ1;
    posQ3 = pQ3;
}

// IQR: just return interquartile range without position information. The seq is partially sorted as done by nth_element
template <typename Tel=double,typename Tcol=vector<double>>
inline Tel iqr(Tcol & seq){

    auto const pQ1 = seq.size() / 4;
    auto const pQ2 = seq.size() / 2;
    auto const pQ3 = pQ1 + pQ2;
    std::nth_element(seq.begin(),          seq.begin() + pQ1, seq.end());
    std::nth_element(seq.begin() + pQ1 + 1, seq.begin() + pQ3, seq.end());

    Tel Q1 = seq[pQ1];
    Tel Q3 = seq[pQ3];

    Tel IQR = Q3-Q1;
    /// The interquartile range is
    return IQR;
}



///// K-means algorithm
//template <typename T=double, typename Tcoll= vector<T>, typename Tcollofcoll = vector<Tcoll> >
//inline void kmeans(const Tcoll& objects, Tcollofcoll& classes,const double& maxd,  const size_t k)
//{
//    if(k<=1) {printval<int>(k,"K equal or less than 1 in kmeans","kmeans.log");return;}
//    if(k>=objects.size()) {printval<int>(k,"K equal or greater than objects in kmeans","kmeans.log");return;}
//
//    classes.resize(k,Tcoll(1));
//
//    classes[0][0] = objects[0];
//    classes[k-1][0] = (objects.back());
//
//    for(auto clss=1; clss< k-1; ++clss){
//        classes[clss].push_back(objects[clss]);
//    }
//
//    Tcoll centroids(k);
//    for(auto clss=0; clss< k; ++clss)
//        centroids[clss]=classes[clss][0];
//
//    bool end = false;
//    int clase = -1;
//    double dist =maxd;
//    do{
//        for(auto obj: objects){ // for each object
//            for(auto ctr=0; ctr< k; ++ctr){ // compute the distance between the object and the centroids
//                double d = abs(obj-centroids[ctr]);
//                if( d< dist){
//                    dist = d;
//                    clase = ctr;
//                }
//            } // close for centroids
//            classes[clase].push_back(obj);
//            dist=maxd;
//        } // close for objects
//
//        // recompute centroids
//        end =true;
//        for(auto cl=0; cl<classes.size(); ++cl){
//            double z=0;
//            for(auto x:classes[cl]){
//                z+=x;
//            }
//            z/=classes[cl].size();
//            if(centroids[cl]!= z){
//                centroids[cl]=z;
//                end=false;
//            }
//
//        }
//        if(!end){ // clean classes for the next iteration
//            for(auto cl=0; cl<classes.size(); ++cl){
//                classes[cl].resize(1);
//            }
//            classes[0][0]= objects[0];
//            classes[k-1][0]= objects.back();
//            for(auto clss=1; clss< classes.size()-1; ++clss){
//                classes[clss]=(objects[clss]);
//            }
//        }
//
//
//    }while(!end);
//
//}


// rIQR: right extreme of the interquartile range. Q3+3*IQR
//The seq is partially sorted as done by nth_element

template <typename Tel,typename Tcol>
inline void riqr(Tcol & seq, Tel& rIQR){

    auto const pQ1 = seq.size() / 4;
    auto const pQ3 = pQ1 + seq.size() / 2;

    std::nth_element(seq.begin(),          seq.begin() + pQ1, seq.end());
    std::nth_element(seq.begin() + pQ1 + 1, seq.begin() + pQ3, seq.end());

    rIQR = seq[pQ3]+ 3*(seq[pQ3]-seq[pQ1]);
}

//MAD: http://en.wikipedia.org/wiki/Median_absolute_deviation
// The expectation of scale=1.4826 times the MAD for large samples of normally distributed Xi is approximately equal to
//the population standard deviation. sigma = scale*MAD
// Compute Median Absolute Deviation. Values in data are sorted or partially sorted at output
// It has assymptotically 37% of efficiency at the normal (has higher variance (Var_optimal_estimator/Var_MAD = 0.37)
//than the optimal estimator)
// However MAD is very robust to outliers it has a breakdown of 50% (maximum possible breakdown)
template <typename Tel,typename Tcol=vector<Tel>>
inline Tel MAD(const Tcol& cdata, const double& scale)
{

    size_t n;
    bool odd = true;
    Tel mad(0);
    Tcol data = cdata;

    if(data.size()%2!=0) mad = odd_median<Tel>(data);
    else{
        odd=false;
        n=data.size()/2;
        sort(data.begin(),data.end());
        mad = (data[n-1]+data[n])/2;
    }

    for(auto& element:data)
        element= abs(element-mad);

    if(odd)
        mad = odd_median<Tel>(data);
    else{

        sort(data.begin(),data.end());
        mad = (data[n-1]+data[n])/2;

    }

    return mad*scale;


}
// Robust estimator of scale: Croux & Rousseeuw 92
// This is the naive algorithm
template <typename Tel=double,typename Tcol=vector<Tel>>
inline double Qn(const Tcol& data)
{
    auto n = data.size(); // n
    auto dim=n*(n-1); // n*(n-1)
    dim/=2; // n*(n-1)/2;
    auto k = (n/2 +1)*(n/2);
    k/=2;
    Tcol q(dim);
    int c=0;
    for(auto i=0; i< data.size()-1; ++i)
        for(auto j=i+1; j< data.size(); ++j){
            q[c] = abs(data[i] - data[j]);
            ++c;
        }

    nth_element(q.begin(), q.begin()+k, q.end());
    double qn = q[k];

    // correction factor
    double dn(0);
    switch(n){
    case 1: // error
        exit(13);
    case 2:
        dn = 0.399;
        break;
    case 3:
        dn = 0.994;
        break;
    case 4:
        dn = 0.512;
        break;
    case 5:
        dn = 0.844;
        break;
    case 6:
        dn = 0.611;
        break;
    case 7:
        dn = 0.857;
        break;
    case 8:
        dn = 0.669;
        break;
    case 9:
        dn = 0.872;
        break;
    default:
        if(n%2) // odd
            dn= n/(n+1.4);
        else
            dn= n / (n+3.8);

    }

    qn *= (2.2219*dn);
    return qn;
}

// Mean (Ave) Absolute Deviation
template <typename T>
inline void AAD(const T& data,double& ave, double& ad){

    auto n=data.size();
	double ep=0.0,s=0.0;

    for(const auto& element:data)
        s+= element;
    ave=s/n;
	ad=0.0;

	for(const auto& element:data)
        ad+= abs(element - ave);

    ad/=n;

}
/**
Anova 1 way statistic if between group variances are homogeneous then it follows F with K-1 and sum(ni) -K degrees of freedom
See Brown and Forsythe 1974 or Ramsey&Ramsey 2007
*/

template <typename T=double,typename Tcol=vector<T> ,typename Tcolcol=vector<Tcol> >
inline T Anova1wayStat(const int& K,const int& N, const Tcolcol& z,const  Tcol& ni){

    vector<T>avezi(K,0.0);
    T accumaccum =0.0;
//    for(int i=0; i<K; ++i){
//        avezi[i] = accumulate(z[i].begin(),z[i].end(),0.0);
//        accumaccum+=avezi[i];
//
//        avezi[i]/=ni[i];
//
//    }

    for(int i=0; i<K; ++i){
        for(int j=0; j< ni[i]; ++j){
            avezi[i]+= z[i][j];
        }
        accumaccum+=avezi[i];
        avezi[i]/=ni[i];
    }

// from From https://stackoverflow.com/questions/58681103/how-to-sum-up-a-vector-of-vector-int-in-c-without-loops
//    double accumaccum = std::accumulate(
//        z.begin(), z.end(),                       // iterators for the outer vector
//        0.0,                                            // initial value for summation - 0
//        [](T init, const vector<T>& intvec){ // binaryOp that sums a single vector<int>
//            return accumulate(
//                intvec.begin(), intvec.end(), // iterators for the inner vector
//                init);                        // current sum
//                                              // use the default binaryop here
//        }
//    );

    T aveave = accumaccum/N;
//cout<<accumaccum<<" \t"<<aveave<<endl;
    ///Anova statistic W see Brown and Forsythe 74 or Wang et al 17 (Table 1)
    T numW = 0,S2=0;
    for(int i=0;i<K;++i){
        S2 = avezi[i]-aveave;
        S2*=S2;
        numW+= ni[i]*S2;
    }

    double denW=0;
    for(int i=0;i<K;++i){
        for(int j=0; j<ni[i];++j){
            S2= z[i][j]-avezi[i];
            S2*=S2;
            denW+= S2;
        }
    }

    T W=((N-K)*numW)/(denW*(K-1));

    return W;

}

template <typename T=double,typename Tcol=vector<T> ,typename Tcolcol=vector<Tcol> >
inline T Anova1wayStat(const int& K,const int& N, const Tcolcol& z,const  Tcol& ni, const Tcol& avei){

//    T accumaccum =0;
//    for(int i=0; i<K; ++i){
//        for(int j=0; j< ni[i]; ++j){
//            accumaccum+=z[i][j];
//
//        }
//
//    }
// from From https://stackoverflow.com/questions/58681103/how-to-sum-up-a-vector-of-vector-int-in-c-without-loops
    double accumaccum = std::accumulate(
        z.begin(), z.end(),                       // iterators for the outer vector
        0.0,                                            // initial value for summation - 0
        [](T init, const vector<T>& intvec){ // binaryOp that sums a single vector<int>
            return accumulate(
                intvec.begin(), intvec.end(), // iterators for the inner vector
                init);                        // current sum
                                              // use the default binaryop here
        }
    );

    T aveave = accumaccum/N;
//cout<<accumaccum<<" \t"<<aveave<<endl;
    ///Anova statistic W see Brown and Forsythe 74 or Wang et al 17 (Table 1)
    T numW = 0,S2=0;
    for(int i=0;i<K;++i){
        S2 = avei[i]-aveave;
        S2*=S2;
        numW+= ni[i]*S2;
    }

    double denW=0;
    for(int i=0;i<K;++i){
        for(int j=0; j<ni[i];++j){
            S2= z[i][j]-avei[i];
            S2*=S2;
            denW+= S2;
        }
    }

    T W=((N-K)*numW)/(denW*(K-1));

    return W;

}


/*
// simple example shrinkage estimator for variance: http://en.wikipedia.org/wiki/Mean_squared_error#Variance

        double hdn = HACdistri1.size();
        var1*= (hdn-1)/(hdn + 1 +(hdn -1)*kurt1/hdn);
        hdn = HACdistri2.size();
        var2*= (hdn-1)/(hdn + 1 +(hdn -1)*kurt2/hdn);

other robust estimators of sigma are MAD see Ripley04

Robust estimator for the location is hodges lehmann estimator:http://en.wikipedia.org/wiki/Hodges%E2%80%93Lehmann_estimator:

In the simplest case, the "Hodges–Lehmann" statistic estimates the location parameter for a univariate population.
Its computation can be described quickly.
For a dataset with n measurements, its Cartesian product with itself has n(n + 1)/2 pairs, including the pair of each item and itself.
For each such pair, the mean is computed;
finally, the median of these n(n + 1)/2 averages is defined to be the Hodges–Lehmann estimator of location.
It is a robust statistic that has a breakdown point of 0.29, which means that the statistic remains bounded even if nearly 30 percent of
the data have been contaminated. This robustness is an important advantage over the sample mean, which has a zero breakdown point,
being proportional to any single observation and so liable to being misled by even one outlier.
The sample median is even more robust, having a breakdown point of 0.50.
The Hodges–Lehmann estimator is much better than the sample mean when estimating mixtures of normal distributions, also.

For symmetric distributions, the Hodges–Lehmann statistic has greater efficiency (lower variance)
than does the sample median. For the normal distribution, the Hodges-Lehmann statistic is nearly as efficient as the sample mean.
For the Cauchy distribution (Student t-distribution with one degree of freedom), the Hodges-Lehmann is infinitely more efficient
than the sample mean, which is not a consistent estimator of the median.

For non-symmetric populations, the Hodges-Lehmann statistic estimates the population's "pseudo-median" a location parameter that is
closely related to the median. The difference between the median and pseudo-median is relatively small, and so this distinction is
neglected in elementary discussions.

*/
// Hodges–Lehmann estimator of the location parameter of a distribution
template <typename T=double,typename Tcol=vector<T>>
inline double HLe(const Tcol& data){
    size_t n= data.size();
    size_t dim = n*(n+1)/2;
    size_t k=0;
    Tcol temp(dim,0);
    for(auto i=0; i<n; ++i)
        for(auto j=i; j<n; ++j){
            temp[k] = (data[i]+data[j])/2;
            ++k;
        }
    return median<T>(temp);
}

// Compute cuasivariance of a STL container using NRC two-step algorithm
template <typename T, typename Tcol>
inline void avesamplevar(const Tcol& data,T& ave, T& var){
	auto n=data.size();
	T ep=0.0,s=0.0;
	ave=0;
	var=0;
    if(n==0) return;
    for(auto& element:data)
        s+= element;
    ave=s/n;
    if(n==1) return;
	for(const auto& element:data){

		s=element-ave;
		ep+= s;
		var+= s*s;
	}

	var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
}

// Given an already computed mean value it computes cuasivariance of a STL container using second step of NRC two-step algorithm
template <typename T=double, typename Tcol=vector<T>>
inline void samplevar(const Tcol& data,const T& ave, T& var){
	auto n=data.size();
	T ep=0.0,s=0.0;
	var=0;
    if(n==0) return;
    if(n==1) return;

	for(const auto& element:data){
		s=element-ave;
		ep+= s;
		var+= s*s;
	}

	var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
}


// Compute cuasivariance of a STL container using NRC two-step algorithm
// use only a partial range, from inir to endr
template <typename T, typename Tcol>
inline void avesamplevarpart(const Tcol& data,T& ave, T& var, const int& inir, const int& endr){
	auto n=max(0,endr-inir);
	T ep=0.0,s=0.0;
	ave=0;
	var=0;
    if(n==0) return;
    for(int i = inir; i <endr; ++i){
            s+= data[i];
    }
    ave=s/n;
    if(n==1) return;
	for(int i = inir; i <endr; ++i){

		s=data[i]-ave;
		ep+= s;
		var+= s*s;
	}

	var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
}

// Compute cuasivariance and kurtosis excess of a STL container using NRC two-step algorithm p. 618
/// skewness and kurtosis are computed in the b1 and b2 forms as defined in Joanes 98, and Joanes and Gill 98
template <typename T, typename Tcol>
inline void avesamplevar(const Tcol& data,T& ave, T& var, T& skew, T& kurt){
	auto n=data.size();
	T ep=0.0,s=0.0,p=0.0;
	ave=0;
	var=0;
    if(n==0) return;
    for(const auto& element:data)
        s+= element;
    ave=s/n;
	skew=0.0;
	kurt=0.0;
    if(n==1) return;
	for(const auto& element:data){

		s=element-ave;
		ep+= s;
		var+= (p=s*s);
		skew+=(p*=s);
		kurt+=(p*=s);
	}

	var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
	if(var!=0){
        skew/= (n*var*sqrt(var)); // this skewness is b1=m3/s³  in Joanes 98 (var = s²)
        kurt=kurt/(n*var*var)-3.0; // kurtosis b2 = m4/s⁴ - 3
	}
	//if(var<0.0){var=0.0;SNPSelTools::texterr("negative var",false);}
}

// Compute cuasivariance and kurtosis excess of a STL container using NRC two-step algorithm p. 618
/// if Pearson=true skewness and kurtosis are computed as defined in Pearson 1931 (see Joanes 1998 or Mandalos 2010)
/// if Pearson =false then n the b1 and b2 forms are computed as defined in Joanes 98, and Joanes and Gill 98
template <typename T, typename Tcol>
inline void avesamplevar(const bool Pearson, const Tcol& data,T& ave, T& var, T& skew, T& kurt){
	auto n=data.size();
	T ep=0.0,s=0.0,p=0.0;
	ave=0;
	var=0;
    if(n==0) return;
    for(const auto& element:data)
        s+= element;
    ave=s/n;
	skew=0.0;
	kurt=0.0;
    if(n==1) return;
	for(const auto& element:data){

		s=element-ave;
		ep+= s;
		var+= (p=s*s);
		skew+=(p*=s);
		kurt+=(p*=s);
	}

	var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
	if(var!=0){
        skew/= (n*var*sqrt(var)); // this skewness is b1=m3/s³  in Joanes 98 (var = s²)
        kurt=kurt/(n*var*var)-3.0; // this is kurtosis b2 = m4/s⁴ - 3

        if(Pearson){ // see for example p. 184 Joanes 98
            double p = 3.0/2.0;
            double g= (n/(n-1.0));

            skew*=pow(g,p);

            kurt = (kurt+3)*pow(g,2.0) -3.0;

        }
	}
	//if(var<0.0){var=0.0;SNPSelTools::texterr("negative var",false);}
}

// Compute cuasivariance and kurtosis excess of a STL container using NRC two-step algorithm p. 618

// use only a partial range, from inir to endr
template <typename T, typename Tcol>
inline void avesamplevar(const Tcol& data,T& ave, T& var, T& skew, T& kurt, const int& inir, const int& endr){
	auto n=max(0,endr-inir);
	T ep=0.0,s=0.0,p=0.0;
	ave=0;
	var=0;
    if(n==0) return;
    for(int i = inir; i <endr; ++i){
            s+= data[i];
    }

    ave=s/n;

	skew=0.0;
	kurt=0.0;
    if(n<=1) return;
	for(int i = inir; i <endr; ++i){

        s=data[i]-ave;
        ep+= s;
        var+= (p=s*s);
        skew+=(p*=s);
        kurt+=(p*=s);

	}

    var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
	if(var!=0){
        skew/= (n*var*sqrt(var));
        kurt=kurt/(n*var*var)-3.0;
	}
	//if(var<0.0){var=0.0;SNPSelTools::texterr("negative var",false);}
}

// use only a partial range, from inir to endr, of the positive data
template <typename T, typename Tcol>
inline void positavesamplevar(const Tcol& data,T& ave, T& var, T& skew, T& kurt, const int& inir, const int& endr){
	auto n=0,counter=0;
	T ep=0.0,s=0.0,p=0.0;
	ave=0;
	var=0;
    if(n==0) return;
    for(int i = inir; i <endr; ++i)
        if(data[i]>0){
            s+= data[i];
            ++counter;
        }
    n=max(counter,1);
    ave=s/n;
	var=0.0;
	skew=0.0;
	kurt=0.0;
    if(n==1) return;
	for(int i = inir; i <endr; ++i){
        if(data[i]>0){
            s=data[i]-ave;
            ep+= s;
            var+= (p=s*s);
            skew+=(p*=s);
            kurt+=(p*=s);
        }
	}

    var=(var-ep*ep/n)/(n-1);
	if(var<0.0){var=0.0;}
	if(var!=0){
        skew/= (n*var*sqrt(var));
        kurt=kurt/(n*var*var)-3.0;
	}
	//if(var<0.0){var=0.0;SNPSelTools::texterr("negative var",false);}
}


// Compute variance and kurtosis excess of a STL container using NRC two-step algorithm p. 618
template <typename T, typename Tcol>
inline void avevar(const Tcol& data,T& ave, T& var, T& skew, T& kurt){
	auto n=data.size();
	T ep=0.0,s=0.0,p=0.0;

    for(const auto& element:data)
        s+= element;
    ave=s/n;
	var=0.0;
	skew=0.0;
	kurt=0.0;

	for(const auto& element:data){

		s=element-ave;
		ep+= s;
		var+= (p=s*s);
		skew+=(p*=s);
		kurt+=(p*=s);
	}

	var=(var-ep*ep/n)/(n);
	if(var<0.0){var=0.0;}
	if(var!=0){
        skew/= (n*var*sqrt(var));
        kurt=kurt/(n*var*var)-3.0;
	}
	//if(var<0.0){var=0.0;SNPSelTools::texterr("negative var",false);}
}


/**
// moments: Compute the mean and the 2nd moment .... until rth moment about the mean. Recall that the 1st moment about the mean is 0
Returns a collection of r elements with type Tcol with first element the mean and the remaining if any moments until the rth (included)

m(r) = sum(xi - ave)^r / n

thus m(2) returns the moments with moments[0]=mean and moments[1] the second moment
m(1) returns the mean in moment[0]

*/
template <typename T=double, typename Tcol=vector<T>>
inline void moments(const int& r, const Tcol& data,Tcol& moments){


	auto n=data.size();

	ToZeroResize((size_t)r, moments); // set the size and initiallize to 0
	T s=0.0;

    for(const auto& element:data)
        s+= element;

    moments[0]=s/n; // the mean

    Tcol dd = data;

	// deviate from the mean
	if(r>1){
        for(size_t d=0; d<n; ++d){

            dd[d]-=moments[0];

            // for each desired moment
            for(size_t m=2; m<=r; ++m){
                moments[m-1]+= pow(dd[d],(double)m);

            }
        }

        // Get the moments
        for(auto& m:moments){
            m/=n;
        }

    }

}

// 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 at each sample
// the common elements at lower freq than maf are stored in lowcommons
// if used with maf=0 then copy all common elements and left lowcommons empty
template <typename T= double, typename Tcol= vector<T>, typename Tcolofcol = vector<Tcol> >
inline void copycommon(const Tcolofcol& data,const Tcolofcol& freqs, Tcol& commons, Tcol& lowcommons, const double& maf){

    bool consider=true, considerlow=false;

    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: data)
        if(minel>sizepos.size())
            minel =sizepos.size();

    if(!commons.empty())

         Tcol().swap(commons);

    if(!lowcommons.empty())
         Tcol().swap(lowcommons);

    int numel=0,element2=-1;

// Search elements present in all collections

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

         for(size_t p2=1;p2<data.size();++p2){ // compare with elements in other collections

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

                break;
             }
             // if the element is common but not reaches maf frequency
             if( (freqs[p2][element2] < (maf)) || ( (1.0-freqs[p2][element2]) < (maf)) )
                considerlow=true;
         }

         if(consider){
            // there is more than maf frequency in population 0
            if( (freqs[0][numel] >= (maf)) && ( (1.0-freqs[0][numel]) >= (maf)) ){ // more than ncopies copy in the sample

                if(considerlow){ // but only ncopies copy in the next collection p2

                    if(!binary_search(lowcommons.begin(),lowcommons.end(),element)){
                        lowcommons.push_back(element); // add to the list
                        //cout<<"pop 2"<<endl;
                    }
                }
                else{ // ok is common with more than ncopies  at each population
                    if(!binary_search(commons.begin(),commons.end(),element))
                        commons.push_back(element);
                }

            }
            else{ // the element is at frequency less than maf in population 0
                if(!binary_search(lowcommons.begin(),lowcommons.end(),element)){
                    lowcommons.push_back(element);
                    //cout<<"pop 1: "<<freqs[0][numel]<<endl;
                }
            }

        } // close if consider

        consider=true;
        considerlow=false;
        ++numel;

        if( (commons.size()+lowcommons.size()) ==minel) break; // cannot exists more common snps
    } // close for elements
    sort(commons.begin(), commons.end());
    sort(lowcommons.begin(), lowcommons.end());
}


/***********************************************************************************
AUTHOR: A Carvajal-Rodriguez
DATE:	081016
FUNCTION: correlxy(container<T> x, container<T>y)
INPUT: two STL containers of size n
PERFORM: computes unbiased estimation of correlation following biometry 15.3 p.559-560. Performs linear or inverse regression estimating a and b which are returned in the pair<rT,rT>
OUTPUT: the standard or unbiased estimate rho of correlation coefficient
************************************************************************************/
template <typename rT=double,typename colT=double, typename Col=vector<colT>>
rT correlxy(const colT& x0, const colT&y,const bool& unbiased,const bool& inverse, pair<rT,rT>& regression){

	int n;
	rT sumx{0}, sumy{0},sumxy=0, sumx2=0,sumy2=0, avex,avey;
	rT rho,denom;
	rT x2desv, b,a;

	if(x0.size()!=y.size()){
		throw invalid_argument("Sample size different in regcorrelxy");
	}
	if(x0.empty()){
        throw invalid_argument("Empty vector in regcorrelxy");
	}

	n=x0.size();
    int nm_n{n}; // nonmissing_n
    colT x{x0};
    if(inverse){
        for(auto& el:x){
            el = (el>0)? (1.0/el): MISSINGVAL;
            if(el<=MISSINGVAL) throw invalid_argument("division by 0 in inverse regression\n");
        }
    }

//for(int ix=0; ix<n;++ix)
//    cout<< x0[ix]<<"\t"<<x[ix]<<endl;
//exit(0);

	for(int i=0;i <n;++i){
		if(y[i]!=MISSINGVAL){
            sumx+=x[i];
            sumy+=y[i];
            sumxy+= x[i]*y[i];
            sumx2+=x[i]*x[i];
            sumy2+=y[i]*y[i];
		}
		else --nm_n;
	}

	avex= sumx/nm_n;
	avey=sumy/nm_n;

	rho = sumxy - nm_n*avex*avey;

	x2desv = (sumx2-nm_n*avex*avex);

	// regression
	b = rho / x2desv ;
    a = avey - b*avex;
    regression = make_pair(a,b);

	denom = x2desv*(sumy2-nm_n*avey*avey);

	if(denom<=0){
		throw invalid_argument("Zero or Negative denominator in regcorrelxy");
    }

	denom=sqrt(denom);

    // correlation
	rho/=denom;

    if(unbiased){
        if(nm_n>4)
            rho=rho*(1 + (1-rho*rho)/(2*(nm_n-4)));
        else
            throw invalid_argument("Too low sample size in regcorrelxy");
    }

	return rho;
}

// deletelist must be sorted from minor to major
// delete the data elements at positions as indicated by the values stored in deletelist
// delete first the higher positions so that index update is not necessary i.e. AFTER DELETING 7TH THE 4TH IS STILL THE 4TH
/// OJO que data.size() si que va cambiando
template <typename T= double, typename Tcol= vector<T>, typename Tcol2= vector<int>>
inline void del_elements(Tcol& data,const Tcol2& deletelist){
    if(deletelist.empty()) return;
    for(int de=(int)deletelist.size()-1; de>=0; --de)
        if( (deletelist[de]>=0) && deletelist[de]< (int)data.size()){
            data.erase(data.begin()+deletelist[de]);
    //cout<<"Deleted "<<deletelist[de]<<endl;
        }

}
/// https://www.techiedelight.com/erase-elements-vector-cpp/
template <typename T= int, typename Tcol= vector<T>>
void erase(Tcol &v, int key)
{
    for (auto it = v.begin(); it != v.end();)
    {
        if (*it == key) {
            it = v.erase(it); // note if we erase the iterator is invalidated and we restore it as l-value by the return of erase. Also because we erasee we do not increase it.
        }
        else {
            ++it;
        }
    }
}


// Convert bits in integer
/**
Usage:
deque<bool> test(100,0); test[0]=1;test[3]=1;test[1]=0;
int sum = bitstoint<deque<bool>>(test);
cout<<sum<<endl;
exit(0);

**/
template <typename Tcol= vector<bool>>
inline int bitstoint(Tcol data){
    return accumulate(data.rbegin(), data.rend(), 0, [](int ini, int bit) { return( (ini << 1) + bit); });
}

template <typename T=int>
void texterr(const T& val,const string filename,const string error_txt,bool finish)
{
/*
	Author: ACR
	Writes in ERROR.txt the error message 'error_txt'
*/
	if(finish){
        printval<int>(val,error_txt,filename);
        exit(val);
	}
	else{

        printval<int>(val,error_txt,filename);
	}
}

// store all words of a given line in the vector words
// return assumes RVO (p. 23  c++SL book)
template <typename Tcol= vector<bool>>
vector<string> get_words(const string& line)
{
	vector<string> words;

	string first_w ("");

	for(size_t i = 0; i < line.length(); i++)
	{
		if( line[i] !=' ' && line[i] !='\t' && line[i] !='\n' && line[i] !='\r')
		{
			first_w += line.substr(i, 1);

			if( line[i+1] ==' ' || line[i+1] =='\t' || line[i+1] =='\n'|| line[i+1] =='\r'
					|| line[i+1] =='\0')
			{
				words.push_back(first_w);

				first_w = "";
			}

		}
	}

	return words;
}

template <typename T= bool>
vector<string> get_words(const string& line, const char& separator)
{
	vector<string> words;

	string first_w ("");

	for(size_t i = 0; i < line.length(); i++)
	{
		if( line[i] !=separator && line[i] !='\n' && line[i] !='\r')
		{
			first_w += line.substr(i, 1);

			if( line[i+1] ==separator || line[i+1] =='\n'|| line[i+1] =='\r'
					|| line[i+1] =='\0')
			{
				words.push_back(first_w);

				first_w = "";
			}

		}
	}

	return words;
}

template <typename Tcol= vector<bool>>
void portablemkdir(const string& dir){

	#ifdef __unix
        mkdir(dir.c_str(),S_IRWXU);
	#elif __APPLE__
		mkdir(dir.c_str(),S_IRWXU);
	#elif macintosh
		mkdir(dir.c_str(),S_IRWXU);
	#elif __linux__
		mkdir(dir.c_str(),S_IRWXU);
	#elif linux
		mkdir(dir.c_str(),S_IRWXU);
    #elif _WIN32
        mkdir(dir.c_str());
    #else
        cout<<"error detecting OS"<<endl;
    #endif


}
template <typename Tcol= vector<bool>>
bool checkWordisInt(string line){
	bool negative=false;
	for(size_t pos = 0; pos < line.length(); pos++){
		while(isspace(line[pos])) ++pos; // ignore initial spaces
		if(isdigit(line[pos]) || (line[pos]=='-' && !negative)) //acr09, allowing negative numbers
			negative=true; // if negative - should appear just 1 time
		else
			return(false);
	}
	return true;

}

/// check is natural number
template <typename Tcol= vector<bool>>
bool isNatNum(string line){

	for(size_t pos = 0; pos < line.length(); pos++){
		while(isspace(line[pos])) ++pos; // ignore initial spaces
		if(!isdigit(line[pos]) )
			return(false);
	}
	return true;

}

/*
*	NRC++ routines

*/
template <typename T= double>
void gcf(T& gammcf, const T& a, const T& x, T& gln){
	const int ITMAX=5000;
	const double EPS=numeric_limits<double>::epsilon();
	//const double EPS=numeric_limits<double>::epsilon() -> DBL_EPSILON;
	const double FPMIN=numeric_limits<double>::min()/EPS;
	int i;
	double an,b,c,d,del,h;

	gln=gammln(a);
	b=x+1.0-a;
	c=1.0/FPMIN;
	d=1.0/b;
	h=d;
	for(i=1;i<=ITMAX;++i){
		an= -i*(i-a);
		b+=2.0;
		d=an*d+b;
		if(fabs(d)<FPMIN) d=FPMIN;
		c=b+an/c;
		if(fabs(c)<FPMIN) c=FPMIN;
		d=1.0/d;
		del=d*c;
		h*=del;
		if(fabs(del-1.0) <= EPS) break;
	}

	if(i > ITMAX){
		string msg("a too large, ITMAX too small in routine gcf");
		cout<<msg<<" x = "<<x<<endl;
		texterr(x,"TLib_ERROR.log", msg,true);
	}
	gammcf=exp(-x+a*log(x)-gln)*h; // Put factors in front
}
template <typename T= double>
void gser(T& gamser, const T& a, const T& x, T& gln){
	const int ITMAX=5000;
	const double EPS= std::numeric_limits<double>::epsilon();
	//const double EPS=numeric_limits<double>::epsilon() -> DBL_EPSILON;
	int n;
	double sum,del,ap;

	gln=gammln(a);
	if(x<=0.0){
		if(x<0.0) texterr(8,"TLib_ERROR.log","x less than 0 in gser",true);
		gamser=0.0;
		return;
	}
	else{
		ap=a;
		del=sum=1.0/a;
		for(n=0;n<ITMAX;++n){
			++ap;
			del*=x/ap;
			sum+=del;
			if(fabs(del)<fabs(sum)*EPS){
				gamser=sum*exp(-x+a*log(x)-gln);
				return;
			}
		}
		string msg("a too large, ITMAX too small in routine gser");
		cout<<msg<<" x = "<<x<<endl;
		texterr(9,"TLib_ERROR.log",msg,true);

		return;
	}
}

template <typename T= double>
double gammp(const T& a, const T& x){
	double gamser,gammcf,gln;

	if(x<0.0||a<=0.0) texterr(10,"TLib_ERROR.log","Invalid arguments in gammp",true);
	if(x<a+1.0){ // use the series representation
		gser(gamser,a,x,gln);
		return gamser;
	}
	else{ // use the continued fraction representation
		gcf(gammcf,a,x,gln);
		return 1.0-gammcf; // and take its complement
	}
}

template <typename T= double>
double gammq(const T& a, const T& x){
	double gamser,gammcf,gln;
	if(x<0.0||a<=0.0){cout<<"x: "<<x<<" a: "<<a<<endl; texterr(11,"TLib_ERROR.log","Invalid arguments in gammq",true);}
	if(x<a+1.0){ // use the series representation
		gser(gamser,a,x,gln);
		return 1.0 - gamser;
	}
	else{ // use the continued fraction representation
		gcf(gammcf,a,x,gln);
		return gammcf; // and take its complement
	}
}

/*
*	From NRC++ pg 232. Incomplete beta function betai and the auxiliar betcf

*/
template <typename T= double>
double betacf(const T& a, const T& b, const T& x){

	const int MAXIT=100;
	const double EPS= std::numeric_limits<double>::epsilon();//Returns the machine epsilon, that is, the difference between 1.0 and the next value representable by the floating-point type T
	const double FPMIN= std::numeric_limits<double>::min()/EPS; // number close to the smaller representable positive floating point number
	int m,m2;
	double aa,c,d,del,h,qab,qam,qap;

	qab=a+b;
	qap=a+1.0;
	qam=a-1.0;
	c=1.0;
	d=1.0-qab*x/qap;
	if(fabs(d) < FPMIN) d = FPMIN;
	d=1.0/d;
	h=d;
	for(m=1;m<=MAXIT;++m){
		m2=2*m;
		aa=m*(b-m)*x/((qam+m2)*(a+m2));
		d=1.0+aa*d;
		if(fabs(d) < FPMIN) d = FPMIN;
		c=1.0+aa/c;
		if(fabs(c) < FPMIN) c = FPMIN;
		d=1.0/d;
		h*= d*c;
		aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
		d=1.0+aa*d;
		if(fabs(d) < FPMIN) d = FPMIN;
		c=1.0+aa/c;
		if(fabs(c) < FPMIN) c = FPMIN;
		d=1.0/d;
		del=d*c;
		h*=del;
		if(fabs(del-1.0) <= EPS) break;
	}// close for

	if(m>MAXIT){
        string msg ="a "+NtoS(a)+" or b "+NtoS(b) +" too big or MAXIT "+NtoS(MAXIT)+" too small in betafc";
        cerr<<msg<<endl;
        texterr(1,"TLib_ERROR.log",msg,false);
	}
	return h;
}
template <typename T= double>
double betai(const T& a, const T& b, const T& x){

	double bt;
	if(x<0.0 || x>1.0)  texterr(12,"TLib_ERROR.log","Bad x in routine betai",true);
	if(x==0.0 || x==1.0) bt=0.0;
	else
		bt = exp(gammln(a+b) - gammln(a) - gammln(b)+ a*log(x)+ b*log(1.0-x));
	if( x < (a+1.0)/(a+b+2.0))
		return bt* betacf(a,b,x)/a;
	else return 1.0 - bt*betacf(b,a,1.0-x)/b;
}

template <typename T= double>
double erffc(const T& x){
	return x < 0.0 ? 1.0 + gammp(0.5,x*x): gammq(0.5,x*x);

}
// Two tailed  p-value of 0 is 0.5 NRC++ pag 641 x is Z as rsqrt(n) in 14.5.2
template <typename T= double>
double getpval_Z(const T& x, const bool twotail){
    double p=0.5*erffc(x/raiz2);
    if(x<0) p=1.0-p;
    if(twotail) p*=2.0;
    return p;
}


/*
*	From NRC++ pg 621 in function ttest. Given a test t value with dg degrees of freedom
*	get the 2-tailed probability P which  is betai( dg/2, 1/2, dg/(dg + t^2) ) pg. 233

*/
template <typename T= double>
double  getpval_t(const T& t,const int & df, const bool twotail){

    double p = betai(0.5*df, 0.5, df/(df+t*t));

    if(!twotail) p/=2.0;
	return p;

}

// Get the probability density of t distribution. Useful for some outliers tests

template <typename T= double>
double t_density(const T& x, const int& gl){

    bool par = (gl%2)==0;
    double Gt=0.0, df = gl;
    double num =df-1, den=sqrt(df);
    int i=3,j=2;

    while(df-j > 1){
        num*=(df-i);
        den*=(df-j);
        i+=2;
        j+=2;
    }

    if(par){
        den*=2;

    }
    else{
        den*= dPI;
    }

    Gt = num/ den;

    double e = (df+1.0) / 2.0;

    e = pow( (1.0 + x*x / df), e );
    e= 1.0/e;

    return Gt*e;

}

/*
*	Modified from NRC++ pg 624. Two-sided ftest
Link about two vs one-sided tests: https://stats.stackexchange.com/questions/267192/doubling-or-halving-p-values-for-one-vs-two-tailed-tests

*/
template <typename T= double>
void getpval_ftest(const T & fval, const T& df1, const T& df2, T& prob){

    prob = 2.0*betai(0.5*df2,0.5*df1,df2/(df2+df1*fval));
    if(prob>1.0) prob=2.0-prob;
}

/// Function version
template <typename T= double>
T getpval_ftest(const T & fval, const T& df1, const T& df2){

    T prob = 2.0*betai(0.5*df2,0.5*df1,df2/(df2+df1*fval));
    if(prob>1.0) prob=2.0-prob;

    return prob;
}

/// August 2021 version
/*
*	From NRC++ pg 624. Ftest

Link about two vs one-sided tests: https://stats.stackexchange.com/questions/267192/doubling-or-halving-p-values-for-one-vs-two-tailed-tests
*/
template <typename T= double>
void getpval_ftest(const bool twotail,const T & fval, const T& df1, const T& df2, T& prob){

    // one-tail
    prob = betai(0.5*df2,0.5*df1,df2/(df2+df1*fval));

    if(twotail){
        if(prob>0.5) prob=1.0-prob;
        prob*=2.0;
    }
}

/// Function version
template <typename T= double>
T getpval_ftest(const bool twotail,const T & fval, const T& df1, const T& df2){

    // one-tail
    T prob = betai(0.5*df2,0.5*df1,df2/(df2+df1*fval));

    if(twotail){
        if(prob>0.5) prob=1.0-prob;
        prob*=2.0;
    }
    return prob;
}


// One tail F-test to check is VA is larger than VB.
//df1 should be n1 - 1 with n1 the range of the data where var1 was obtained. Similarly df2.
template <typename T= double>
void Ftest(const T & Va, const T & Vb,const T& df1, const T& df2, T& fval, T& prob){

    if(Vb==0.0){fval = MISSINGVAL;prob=0;return;}
    fval = Va/Vb;
    if(fval<=1){prob=1;return;}
    prob = betai(0.5*df2,0.5*df1,df2/(df2+df1*fval));
}

// Given a chissquare value and the degrees of freedom return the probability P
//ACR July 21 change df to type T2=int
template <typename T= double, typename T2=int>
double getChiPval(const T& chsk,const T2 & df){

	if(chsk==0.0) return 1.0;
	return gammq(0.5*df,0.5*chsk);
}

// DAgostino K2 ommnibus test is goodness of fit measure of departure from normality http://en.wikipedia.org/wiki/D%27Agostino%27s_K-squared_test
/// As defined in DAgostino et al 90 it uses the Pearson definitions for Skewness and Kurtosis (equations 4 and 5 in that article)

//asinh(x) = ln(x+sqrt(x²+1) in C++: double asinh (double x); http://www.cplusplus.com/reference/cmath/asinh/
template <typename T= double>
double DAK2test(const T& skew, const T& kurt, const int& in){

	if(in<=7){ // cannot be computed
        string msg="DAK2test test cannot be performed because sample size n = "+ NtoS(in)+" is less than 8";
        cout<<msg<<endl;
        texterr(2,"TLib_ERROR.log",msg,false);
        return 0; // significant departure is assumed
	}
	double n = in;
	double G2g1 = (36*(n-7)*(n*n + 2*n -5))/ ( (n-2)*(n+5)*(n+7)*(n+9) );
	double G1g2 = ((6*(n*n - 5*n +2))/((n+7)*(n+9))) *sqrt((6*(n+3)*(n+5))/(n*(n-2)*(n-3)));
	double w2=sqrt(2*G2g1 +4)-1;
	double w = sqrt(w2);
	double delta = 1 / sqrt(log(w));
	double a2 = 2/(w2-1);
	double m2g1 = 6*(n-2)/((n+1)*(n+3));
	double x = skew/( sqrt(a2)*sqrt(m2g1) );

	/// Skewness test
	double Z1=delta*asinh(x); //Z1 = delta* log( skew/( sqrt(a2)*sqrt(m2g1) ) + sqrt((skew*skew)/(a2*m2g1) +1) );

	getpval_Z(Z1,true); // true means two tail

	double m1g2 = 6/(n+1);
	double m2g2 = (24*n*(n-2)*(n-3)) / ((n+1)*(n+1)*(n+3)*(n+5));

	double A = 6 + (8/G1g2)*(2/G1g2 +sqrt(1+4/(G1g2*G1g2)));

	double Z2part1 = sqrt((9*A)/2);
	double Z2part2 = ( 1-2/(9*A));
	double Z2part3num = (1-2/A);
	double Z2part3denom = (1+((kurt-m1g2)/sqrt(m2g2))*sqrt(2/(A-4)));
    double Z2part3 = pow(Z2part3num/Z2part3denom, 1.0/3.0);

    /// Kurtosis test
    double Z2 = Z2part1*(Z2part2- Z2part3);


    getpval_Z(Z2,false); // false means one tail

    /// Omnibus test
    double K2 = Z1*Z1 + Z2*Z2;

//if(Z2part3denom>0)
//
//else{
//	string msg="Z2 could not be computed in DAK2test. Value into cubic root was "+ NtoS(Z2part3denom)+". Only Skewness test is performed.";
//	cout<<msg<<endl;
//	texterr(2,"TLib_ERROR.log",msg,false);
//	Z2=0;
//}



	return getChiPval(K2,2);

}

// DAgostino K2 test is goodness of fit measure of departure from normality http://en.wikipedia.org/wiki/D%27Agostino%27s_K-squared_test
/// As defined in DAgostino et al 90 it uses the Pearson definitions for Skewness and Kurtosis (equations 4 and 5 in that article)
/// THIS VERSION besides the ommnibus K2 test also performs the skewness and kurtosis test separately

//asinh(x) = ln(x+sqrt(x²+1) in C++: double asinh (double x); http://www.cplusplus.com/reference/cmath/asinh/
template <typename T= double>
T DAK2test(const bool verbose, const T& skew, const T& kurt, const int& in, const bool Stwotail,const bool Ktwotail, double&Spval,double& Kpval){

	if(in<=7){ // cannot be computed
        string msg="DAK2test test cannot be performed because sample size n = "+ NtoS(in)+" is less than 8";
        if(verbose){
            cout<<msg<<endl;
            texterr(2,"TLib_ERROR.log",msg,false);
        }

        Spval=0;
        Kpval=0;
        return 0; // significant departure is assumed
	}
	double n = in;
	double G2g1 = (36*(n-7)*(n*n + 2*n -5))/ ( (n-2)*(n+5)*(n+7)*(n+9) );
	double G1g2 = ((6*(n*n - 5*n +2))/((n+7)*(n+9))) *sqrt((6*(n+3)*(n+5))/(n*(n-2)*(n-3)));
	double w2=sqrt(2*G2g1 +4)-1;
	double w = sqrt(w2);
	double delta = 1 / sqrt(log(w));
	double a2 = 2/(w2-1);
	double m2g1 = 6*(n-2)/((n+1)*(n+3));
	double x = skew/( sqrt(a2)*sqrt(m2g1) );

	/// Skewness test
	double Z1=delta*asinh(x); //Z1 = delta* log( skew/( sqrt(a2)*sqrt(m2g1) ) + sqrt((skew*skew)/(a2*m2g1) +1) );

	Spval=getpval_Z(Z1,Stwotail); // Stwotail=true means two tail

	double m1g2 = 6/(n+1);
	double m2g2 = (24*n*(n-2)*(n-3)) / ((n+1)*(n+1)*(n+3)*(n+5));

	double A = 6 + (8/G1g2)*(2/G1g2 +sqrt(1+4/(G1g2*G1g2)));

	double Z2part1 = sqrt((9*A)/2);
	double Z2part2 = ( 1-2/(9*A));
	double Z2part3num = (1-2/A);
	double Z2part3denom = (1+((kurt-m1g2)/sqrt(m2g2))*sqrt(2/(A-4)));
    double Z2part3 = pow(Z2part3num/Z2part3denom, 1.0/3.0);


double Z2=0;
if(Z2part3denom<=0){
	string msg="Z2 could not be computed in DAK2test. Value into cubic root was "+ NtoS(Z2part3denom)+". Only Skewness test is performed.";
	if(verbose){
        cout<<msg<<endl;
        texterr(2,"TLib_ERROR.log",msg,false);
    }
	Kpval=1;
}
else{
    /// Kurtosis test
    Z2 = Z2part1*(Z2part2- Z2part3);
    Kpval=getpval_Z(Z2,Ktwotail); // false means one tail
}

    /// Omnibus test
    double K2 = Z1*Z1 + Z2*Z2;

	return getChiPval(K2,2);

}

// Modification of the DAgostino K2 test to test only skewness
template <typename T= double>
double DAK1test(const T& skew, const int& in){

	double G2g1,w2,delta,m2g1,a2,Z1,K1;
	double n = in;

	G2g1 = (36*(n-7)*(n*n + 2*n -5))/ ( (n-2)*(n+5)*(n+7)*(n+9) );

	w2 =sqrt(2*G2g1 +4)-1;

	delta = 1 / sqrt(log(sqrt(w2)));

	a2 = 2/(w2-1);

	m2g1 = 6*(n-2)/((n+1)*(n+3));

	Z1 = delta* log( skew/( sqrt(a2)*sqrt(m2g1) ) + sqrt((skew*skew)/(a2*m2g1) +1) );


	return getpval_Z(Z1,true); // see D'Agostino et al 90

}


/// ******************* VARIANCE HOMOGENEITY TESTS


/// Adaptive location estimators (Keselman eta al 08; Othman et al 12)
/// Hogg measures and Hinge estimators. Z must be sorted

template <typename T= double>
inline void ULBCDE(const T& alpha, const vector<T>& Z, T& U, T& L, T& B, T& C, T& D, T& E){

    if(!is_sorted(Z.begin(),Z.end(),std::greater<int>())){
        string msg("Measure ULBCDE requires sorted vector in decreasing order");
		cout<<msg<<" Alpha = "<<alpha<<endl;
		texterr(7,"TLib_ERROR.log", msg,true);
    }

    if(alpha>0.5){
        string msg("Measure ULBCDE requires alpha not greater than 0.5");
		cout<<msg<<" Alpha = "<<alpha<<endl;
		texterr(7,"TLib_ERROR.log", msg,true);
    }

    int n =(int)Z.size();
    int nAlpha = n*alpha;

    if(nAlpha<1){
        string msg("Measure ULBCDE requires N*alpha greater than 0");
		cout<<msg<<"n, Alpha = "<<n<<", " <<alpha<<endl;
		texterr(7,"TLib_ERROR.log", msg,true);
    }


    U=0,L=0;
    /// alpha*(1 + nB + nC +nD +nE + 1) = 1 and nB=nE, nC=nD with nC=2nB => 2+6nB=1/alpha => nB= (1-2*alpha)/(6*alpha)
    /// accumulate the values
    int nB=  n*(1.0-2*alpha)/(6*alpha);
    int nC= 2*n*(1.0-2*alpha)/(6*alpha);

    /// accumulate the values
    int accB = nAlpha +nB;
    int accC = accB +nC;
    int accD = accC +nC;
    int accE = accD +nB;

    for(int i=0;i< n; ++i){

        if(i<nAlpha){
            U+=Z[i]; // Largest nalpha observations
            L+= Z[n-i-1]; // smallest nalpha observations
        }
        else if(i<accB){
            B+=Z[i];
        }
        else if(i<accC){
            C+=Z[i];
        }
        else if(i<accD){
            D+=Z[i];
        }
        else if(i<accE){
            E+=Z[i];
        }

    }

    U/=nAlpha;
    L/=nAlpha;
    if(nB>0){
        B/=nB;
        E/=nB;
    }
    if(nC>0){
        C/=nC;
        D/=nC;
    }


}


template <typename T= double>
inline void HQ1(const T& alpha,const T&X1,const T&Xn, const T& median, T& alpha1, T& alpha2){

    ;

}

/// X1 should be the highest value of the sample Xn is the lowest
template <typename T= double>
inline void HSK2(const T& alpha,const T&X1,const T&Xn, const T& median, T& alpha1, T& alpha2){

    T U = X1-median;
    T L = median - Xn;
    alpha1= alpha*(U/(U+L));
    alpha2=alpha-alpha1;

}
/// Z must be sorted
template <typename T= double>
inline T trimmedMean(const T& alpha1,const T& alpha2, const int& n,const vector<T>& Z){
    int g1= alpha1*n;
    int g2= alpha2*n;
    int h = n-g1-g2;
    int fin = n-g2;
    T mt =0;
    for(int i=g1; i<fin; ++i){
        mt+= Z[i];
    }
    return mt/h;
}


/// One-side Levene test (e.g. to test if v1>v2) based upon group assymetric trimmed means.
/**
This is the cHSK2 method in Keselman et al 2008. These authors recommend alpha = 0.2
The rows of the vector Z become sorted from major to minor.
**/
template <typename T= double>
T LeveneHSK2(const T& alpha, const int& K, const int& N, const vector<T> &ni,vector<vector<T>>& Z, vector<T>& median, T& F){

    if(K!=(int)median.size()){
		string msg("K and vector size are different in routine Levene");
		cout<<msg<<" v size = "<<median.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine Levene");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    vector<vector<T>> Zs=Z;

    double alpha1=0,alpha2=0;
    for(int i=0; i<K;++i){
        sort(Z[i].begin(),Z[i].end(),std::greater<int>()); // sort from major to minor

        HSK2(alpha,Z[i][0], Z[i][ni[i]-1], median[i], alpha1, alpha2);
        double mt = trimmedMean(alpha1,alpha2, ni[i],Z[i]);

        for(int j=0;j<ni[i];++j){
            Zs[i][j]= fabs(Z[i][j]-mt);
        }
    }

    F=Anova1wayStat(K,N, Zs,ni);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}


/// Bartlett test
template <typename T= double>
T Bartlett(const int& K,const int& N, const vector<T> &ni, const vector<T>& v, T& F){

    if(K!=(int)v.size()){
		string msg("K and vector size are different in routine Bartlett");
		cout<<msg<<" v size = "<<v.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine Bartlett");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }

	T S2=0;
    // Compute pooled sample variance
	for(int i=0; i<K; ++i){
        S2+=(ni[i]-1)*v[i];
	}

    int NmK=N-K;
	S2/=NmK;

	//Compute Bartlett statistic
    T B=0, D=0;
    for(int i=0; i<K; ++i){
        B+=(ni[i]-1)* log(v[i]);
        D+= 1.0/(ni[i]-1);
	}

	B= NmK*log(S2) - B;

	D-= 1.0/(NmK);

	D*= 1.0/(3.0*(K-1));
	D+=1;

    F=B/D;

	return getChiPval(B,K-1); // s

}

/// One-side Levene test (e.g. to test if v1>v2)
template <typename T= double>
T Levene(const int& K, const int& N, const vector<T> &ni,const vector<vector<T>>& Z, const vector<T>& ave, T& F){

    if(K!=(int)ave.size()){
		string msg("K and vector size are different in routine Levene");
		cout<<msg<<" v size = "<<ave.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine Levene");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    vector<vector<T>> Zs=Z;

    for(int i=0; i<K;++i){
        for(int j=0;j<ni[i];++j){
            Zs[i][j]= fabs(Zs[i][j]-ave[i]);
        }
    }

    F=Anova1wayStat(K,N, Zs,ni);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}

/// A variation in the Levene for low and assymmetric sample sizes
template <typename T= double>
T Levene2(const int& K, const int& N, const vector<T> &ni,const vector<vector<T>>& Z, const vector<T>& ave, T& F){

    if(K!=(int)ave.size()){
		string msg("K and vector size are different in routine Levene");
		cout<<msg<<" v size = "<<ave.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine Levene");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    vector<vector<T>> Zs=Z;

    for(int i=0; i<K;++i){
        for(int j=0;j<ni[i];++j){
            Zs[i][j]= fabs(Zs[i][j]-ave[i]);
        }
    }

    F=Anova1wayStat(K,N, Zs,ni,ave);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}


/// One-side BForsythe test (e.g. to test if v1>v2)
template <typename T= double>
T BForsythe(const int& K, const int& N, const vector<T> &ni, const vector<vector<T>>& Z, const vector<T>& median, T& F){

    if(K!=(int)median.size()){
		string msg("K and median vector size are different in routine BForsythe");
		cout<<msg<<" v size = "<<median.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine BForsythe");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }

    vector<vector<T>> Zs=Z;

    for(int i=0; i<K;++i){
        for(int j=0;j<ni[i];++j){
            Zs[i][j]= fabs(Zs[i][j]-median[i]);
        }
    }

    F=Anova1wayStat(K,N, Zs,ni);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}

// especial version with original avei
template <typename T= double>
T BForsythe2(const int& K, const int& N, const vector<T> &ni, const vector<vector<T>>& Z, const vector<T>& median, T& F, const vector<T>& avei){

    if(K!=(int)median.size()){
		string msg("K and median vector size are different in routine BForsythe");
		cout<<msg<<" v size = "<<median.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine BForsythe");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }

    vector<vector<T>> Zs=Z;

    for(int i=0; i<K;++i){
        for(int j=0;j<ni[i];++j){
            Zs[i][j]= fabs(Zs[i][j]-median[i]);
        }
    }

    F=Anova1wayStat(K,N, Zs,ni,avei);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}


/// One-side OBrien test (e.g. to test if v1>v2)
template <typename T= double>
T OBrien(const int& K, const int& N, const vector<T> &ni,const vector<vector<T>>& Z, const vector<T>& avei, const vector<T>& v, T& F){

    if(K!=(int)v.size()){
		string msg("K and vector size are different in routine OBrien");
		cout<<msg<<" v size = "<<v.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine OBrien");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }

     if(int n = *min_element(begin(ni),end(ni)) < 3 ){
     		string msg("OBrien test cannot be calculated for n<3");
            cout<<msg<<" Some n is = "<<n<<endl;
            texterr(7,"TLib_ERROR.log", msg,true);
     }
    vector<vector<T>> Zs=Z;

    for(int i=0; i<K;++i){
        T ari = (ni[i]-1.5)*ni[i];
        T vri=v[i]*(ni[i]-1)/2;
        ///
        T denomri= (ni[i]-1)*(ni[i]-2);

        for(int j=0;j<ni[i];++j){
            T S=Zs[i][j]-avei[i];
            S*=S;
            Zs[i][j]= (ari*S - vri)/denomri;
        } // close for j
    } // close for i

    F=Anova1wayStat(K,N, Zs,ni);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}

template <typename T= double>
T OBrien2(const int& K, const int& N, const vector<T> &ni,const vector<vector<T>>& Z, const vector<T>& avei, const vector<T>& v, T& F){

    if(K!=(int)v.size()){
		string msg("K and vector size are different in routine OBrien");
		cout<<msg<<" v size = "<<v.size()<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }
    if(N<=K){
		string msg("N <= K in routine OBrien");
		cout<<msg<<" N = "<<N<<" K = "<<K <<endl;
		texterr(7,"TLib_ERROR.log", msg,true);

    }

     if(int n = *min_element(begin(ni),end(ni)) < 3 ){
     		string msg("OBrien test cannot be calculated for n<3");
            cout<<msg<<" Some n is = "<<n<<endl;
            texterr(7,"TLib_ERROR.log", msg,true);
     }
    vector<vector<T>> Zs=Z;

    for(int i=0; i<K;++i){
        T ari = (ni[i]-1.5)*ni[i];
        T vri=v[i]*(ni[i]-1)/2;
        ///
        T denomri= (ni[i]-1)*(ni[i]-2);

        for(int j=0;j<ni[i];++j){
            T S=Zs[i][j]-avei[i];
            S*=S;
            Zs[i][j]= (ari*S - vri)/denomri;
        } // close for j
    } // close for i

    F=Anova1wayStat(K,N, Zs,ni,avei);
    T df1=K-1, df2 = N-K;
    // one-side F-test
    return getpval_ftest(false,F, df1, df2);
}


/// ******************* END OF VARIANCE HOMOGENEITY TESTS



// kolmogorov Smirnov probability function
template <typename T= double>
inline double probks(const T& alam){

	const double EPS1= 1.0e-6, EPS2= 1.0e-16;

	double a2, fac=2.0, sum=0.0, term, termbf=0.0;

	a2= -2.0*alam*alam;
	for(int j=1;j<=300; ++j){
		term= fac*exp(a2*j*j);
		sum+=term;
		if(fabs(term) <= EPS1*termbf || fabs(term)<= EPS2*sum) return sum;
		fac=-fac; // Alternating signs in sum
		termbf=fabs(term);
	}
	return 1.0; // Get here only by failing to converge
}

/*
* Function KStwo for Kolmogorov-Smirnov in p. 630 NRC++.
* In this version the vectors inside data are sorted.
*
*/
template <typename Tcol= vector<bool>>
void KStwo(vector< vector<double> >& data,vector<double>& d, vector<double> &prob){

    int n = data.size();
	int numpairs = n * (n-1) / 2;
	int np=0;

    d.resize(numpairs);
    prob.resize(numpairs);


	size_t j1=0, j2=0;
	int d1,d2;
	double dt, en1, en2, en, fn1=0.0, fn2=0.0;

    // sort each vector in data
    for(auto& popdata: data)
        sort(popdata.begin(),popdata.end());

    for(size_t p1 = 0; p1 < data.size(); ++p1)
        for(size_t p2 = p1+1; p2<data.size(); ++p2){

            size_t n1 = data[p1].size();
            size_t n2 = data[p2].size();

            en1=n1; // store sizes as double
            en2=n2;

            d[np]=0.0;

            while(j1<n1 && j2<n2){
                if( (d1=data[p1][j1]) <= (d2=data[p2][j2]) ) fn1 = j1++ / en1;
                if(d2 <= d1) fn2 = j2++/en2;
                if( (dt=fabs(fn2-fn1)) > d[np]) d[np] = dt;
            }

            en = sqrt(en1*en2/(en1+en2));
            prob[np] = probks( (en+0.12+0.11/en)*d[np]); // compute significance
            ++np;
        }

}



///
/*
the cumulative density function (CDF) of a standard normal (Gaussian) random variable.
It is closely related to the error function
*/
//double phi(double x)
//{
//    // constants
//    double a1 =  0.254829592;
//    double a2 = -0.284496736;
//    double a3 =  1.421413741;
//    double a4 = -1.453152027;
//    double a5 =  1.061405429;
//    double p  =  0.3275911;
//
//    // Save the sign of x
//    int sign = 1;
//    if (x < 0)
//        sign = -1;
//    x = fabs(x)/sqrt(2.0);
//
//    // A&S formula 7.1.26
//    double t = 1.0/(1.0 + p*x);
//    double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
//
//    return 0.5*(1.0 + sign*y);
//}





/** What follows:

Code modified from the book C++ Templates - The Complete Guide
by David Vandevoorde and Nicolai M. Josuttis, Addison-Wesley, 2002

Code modified from the book C++ Standard Library
by Nicolai M. Josuttis, Addison-Wesley, 2012 2nd Ed.

Usage:
deque<bool> test(100,0);
printcol(test);
vector<double> freqs(10, 0.1);
printcol(freqs);


**/


/*
printElements will print in stdout or file a single basic type value or collection depending on the Policy. By default is the collection in the stdout.

*/


class baseTPrintPolicy{
public:
    template<typename T, typename To = std::ostream>
    static void print(const T& value, To& out){
        out<<value<<std::endl;
        out<< std::flush;
    }

//    template<typename T, typename To = std::ostream>
//    static void print(const T& value, const string& msg, To& out){
//        if(value==T(0))
//            out<<msg<<std::endl;
//        else
//        out<<msg<<" "<<value<<std::endl;
//        out<< std::flush;
//    }
    /// 28 abril 2026
    template<typename T, typename To = std::ostream>
    static void print(const T& value, const string& msg, To& out)
    {
        if (value == T(0))
            out << msg << '\n';
        else
            out << msg << " " << value << '\n';
    }



};
/* from: http://stackoverflow.com/questions/9407367/determine-if-a-type-is-an-stl-container-at-compile-time
CHECK IF A TYPE IS AN STL CONTAINER OR NOT

template <typename T>
struct is_cont {
  static const bool value = false;
};

template <typename T,typename Alloc>
struct is_cont<std::vector<T,Alloc> > {
  static const bool value = true;
};

// usage

if (is_cont<X>::value) { ... }


*/

    template <typename T>
    struct element_is_coll { // by default value is false
        static const bool value = false;
    };

    template <typename T,typename Alloc>
    struct element_is_coll<std::vector<T,Alloc> > { // this specialization when the type of the struct is a vector the value is true
        static const bool value = true;
    };

    template <typename T,typename Alloc>
    struct element_is_coll<std::deque<T,Alloc> > { // this specialization when the type of the struct is a deque the value is true
        static const bool value = true;
    };

template <typename typeelement>
class collTPrintPolicy{
public:

    template<typename T, typename To = std::ostream>
    static void print(const T& seq, To& out){


        if( typeid(bool) == typeid(seq[0]) )
            out << noboolalpha;

        for(const auto& elem : seq){
            out<<elem<<" ";
        }
        out<<std::endl;

    }

    /// Oct 2023
    template<typename T, typename To = std::ostream>
    static void print(const T& seq, To& out,const bool absval){


        if( typeid(bool) == typeid(seq[0]) )
            out << noboolalpha;

        if(absval){
            for(const auto& elem : seq){
                out<<fabs(elem)<<" ";
            }
        }
        else
            for(const auto& elem : seq){
                out<<elem<<" ";
            }
        out<<std::endl;

    }

};



// specialization of collTPrintPolicy
template <>
class collTPrintPolicy<deque<bool>>{
public:

    template<typename T, typename To = std::ostream>
    static void print(const T& seq, To& out){



        out << noboolalpha;

        for(const auto& elem : seq){
            out<<elem<<" ";
        }
        out<<std::endl;

    }

};

template <typename typeelement,typename Alloc>
class collTPrintPolicy<deque<typeelement,Alloc> >{
public:

    template<typename T, typename To = std::ostream>
    static void print(const T& seq, To& out){

        //if (element_is_coll<typeelement>::value){
            if( typeid(bool) == typeid(seq[0][0]) )
                out << noboolalpha;
            for(const auto& elem : seq){

                //print(elem, out);
                for(const auto& eleminelemen : elem){
                    out<<eleminelemen<<" ";
                }
                out<<std::endl;
            }
        //}
    }

};


// Print in stdout. Default policy is for collections
template <typename T, typename Policy=collTPrintPolicy<T>>
inline void printElements(const T& X){

    Policy::print(X,std::cout);

}

// Oct 2023
template <typename T, typename Policy=collTPrintPolicy<T>>
inline void printElements(const T& X, const bool absval){

    Policy::print(X,std::cout,absval);

}

// Print in std::ofstream by default stdout
template <typename T, typename Policy=collTPrintPolicy<T>, typename To = std::ostream>
inline void printElements(const T& X, To& out){

    Policy::print(X,out);

}

// Overload if string is the second argument: the compiler Code::Blocks was not able to resolve without adding the 3rd argument
template <typename T, typename Policy=collTPrintPolicy<T>>
inline void printElements(const T& X, const string& name, const bool append=false){
    ofstream fout;
    if(append) fout.open(name, ofstream::app);
    else  fout.open(name, ofstream::out);
    Policy::print(X,fout);

    fout.close();
}

//template <typename T, typename Policy=collTPrintPolicy<T>>
//inline void printElements(const T& X, const string& msg, const string& name, const bool append=false){
//    ofstream fout;
//    if(append) fout.open(name, ofstream::app);
//    else  fout.open(name, ofstream::out);
//    Policy::print(X,msg,fout);
//    fout<< std::flush;
//    fout.close();
//}

///28 abril 2026
template <typename T, typename Policy=collTPrintPolicy<T>>
inline void printElements(const T& X, const string& msg, const string& name, const bool append=false)
{
    ofstream fout;

    if (append)
        fout.open(name, ios::out | ios::app);
    else
        fout.open(name, ios::out | ios::trunc);

    if (!fout) {
        cerr << "ERROR: could not open log file: " << name << endl;
        return;
    }

    Policy::print(X, msg, fout);
}


//Simpler sintaxis

template <typename T>
inline void printcol(const T& v){
    printElements(v);
}


//Oct 2023
template <typename T>
inline void printcol(const T& v, const bool absval){
    printElements(v,absval);
}


template <typename T>
inline void printcol(const deque<T>& v){
    printElements<deque<T>>(v);
}

template <typename T>
inline void printcol(const T& v, const string& name){
    printElements(v, name);
}

template <typename T>
void printval(const T& v){
    printElements<T,baseTPrintPolicy>(v);
}

template <typename T>
inline void printval(const T& v, const string& name){

    printElements<T,baseTPrintPolicy>(v, name);
}

template <typename T>
inline void printval(const T& v,const string& msg ,const string& name){

    printElements<T,baseTPrintPolicy>(v, msg, name,false); /// May 2024 if true then it appends and seems to be a hang out sometimes, do not know why
}


/* resamplwithout: sample numelements integers (not reals) without replacement from a range [min-max]
Example of use:
vector<int>v;
resamplwithout(1, 11, 3,v,default_random_engine(seed) ); // generate a vector with 3 integers between 1 and 10
printcol(v);
exit(0);
*/
// WARNING: If this function is used inside a for the engine g should be initialized incorporating the index in order to avoid the same seed in different for iterations
//e.g. default_random_engine(std::chrono::system_clock::now().time_since_epoch().count()+i+j) if the engine is inside a double for

template <typename T=size_t, typename T2=size_t>
inline bool resamplwithout(const T& min, const T& max, const T2& numelements, vector<T>& sample, default_random_engine&& g){
    if(max<=min) return false; // dar aviso del error o resolver
    if(numelements<1) return false; // dar aviso del error o resolver
    auto n = max-min;
    if(n<numelements){
        string msg="Cannot sample  "+ NtoS(numelements)+" elements from "+NtoS(n);
        cout<<msg<<endl;
        texterr(13,"TLib_ERROR.log",msg,false);
    }

    T c=min;
    // Define the vector of possible values
    vector<T> universe(n);
    // Charge the vector
    for_each(universe.begin(),universe.end(), [&](T& element){element=c;++c;});
    // Shuffle
    shuffle ( universe.begin(), universe.end(), g );
    // Clean vector if necessary
    if(!sample.empty()){

        vector<T>().swap(sample);
    }
    // Get the desired values
    copy(universe.cbegin(),universe.cbegin()+numelements, back_inserter( sample));

    return true;
}
//return the number of different elements in a collection excluding missing val if any
// no elements in v are rearranged
// WARNING: NOTE THAT after unique, the elements between the returned iterator and last are left in a valid but unspecified state
/// it is for this that v should be passed as constant or by value
template <typename T=int, typename Tcol=vector<T>>
int count_distinct_els(Tcol v, T missingval)
{
    sort(v.begin(), v.end()); // Average case O(n log n), worst case O(n^2) (usually implemented as quicksort.
  // To guarantee worst case O(n log n) replace with make_heap, then sort_heap.

  // Unique will take a sorted range, and move things around to get duplicated
   // items to the back and returns an iterator to the end of the unique section of the range
    auto unique_end = unique(v.begin(), v.end()); // Again n comparisons

    return (distance(v.begin(), unique_end) - count(v.begin(), unique_end,missingval)); // Constant time for random access iterators (like vector's)
}

// Same as before but now the position of the first repeated element is returned within the variable firstrepeatedpos
template <typename T=int, typename Tcol=vector<T>>
int count_distinct_els(Tcol v, T missingval, T& firstrepeatedpos)
{
    //Tcol vt= v;
    sort(v.begin(), v.end()); // Average case O(n log n), worst case O(n^2) (usually implemented as quicksort.
  // To guarantee worst case O(n log n) replace with make_heap, then sort_heap.

  // adjacent find just detect the first pair of repeated elements and return an iterator to the repeated position
      firstrepeatedpos = *adjacent_find(v.begin(), v.end());


  // Unique will take a sorted range, and move things around to get duplicated
   // items to the back and returns an iterator to the end of the unique section of the range
    auto unique_end = unique(v.begin(), v.end()); // Again n comparisons
    int pos = distance(v.begin(), unique_end);

    return (pos - count(v.begin(), unique_end,missingval)); // Constant time for random access iterators (like vector's)
}

/// Go over a collection of collections and check values for a user defined column looking for repetitions. The row with the repeated item is deleted
///Same as before but now the position of the first repeated element is returned within the variable firstrepeatedpos
//template <typename T=int, typename Tcol=vector<T>, typename Tcolofcol=vector<Tcol> >
//int filter_rep_SNPs(Tcolofcol& tv, T missingval, T& firstrepeatedpos)
//{
//
//    for(auto& e: cv)
//
//    sort(v.begin(), v.end()); // Average case O(n log n), worst case O(n^2) (usually implemented as quicksort.
//  // To guarantee worst case O(n log n) replace with make_heap, then sort_heap.
//
//  // adjacent find just detect the first pair of repeated elements and return an iterator to the repeated position
//      firstrepeatedpos = *adjacent_find(v.begin(), v.end());
//
//
//  // Unique will take a sorted range, and move things around to get duplicated
//   // items to the back and returns an iterator to the end of the unique section of the range
//    auto unique_end = unique(v.begin(), v.end()); // Again n comparisons
//    int pos = distance(v.begin(), unique_end);
//
//    return (pos - count(v.begin(), unique_end,missingval)); // Constant time for random access iterators (like vector's)
//}

// return the number of different elements in a collection of collections excluding missing val if any
// store a sorted list of these elements
template <typename T=int, typename Tcol=vector<T>, typename Tcolofcol=vector<Tcol> >
int do_unique_list(const Tcolofcol& cv, T missingval, Tcol& uniquelist)
{
    Tcol v, templist;

    for(auto& e: cv){
        v = e;
        sort(v.begin(), v.end());
        auto unique_end = unique(v.begin(), v.end()); //
        for (auto it = begin(v); it!=unique_end; ++it){ // for each unique element
            if(*it!= missingval){
                templist.push_back(*it);
            }
        }

        v.clear();
    } // close for cv

    sort(templist.begin(), templist.end());
    auto unique_end = unique(templist.begin(), templist.end());

    if(!uniquelist.empty()) uniquelist.clear();
    uniquelist.reserve(templist.size());

    for (auto it = begin(templist); it!=unique_end; ++it){ // for each unique element

        uniquelist.push_back(*it);

    }

  // Unique will take a sorted range, and move things around to get duplicated
   // items to the back and returns an iterator to the end of the unique section of the range

    //return (distance(templist.begin(), unique_end) ); // Constant time for random access iterators (like vector's)
    return(uniquelist.size());
}

// Same as before but now the position of the first repeated element is returned within the variable firstrepeatedpos
template <typename T=int, typename Tcol=vector<T>, typename Tcolofcol=vector<Tcol> >
int do_unique_list(const Tcolofcol& cv, T missingval, Tcol& uniquelist, T& firstrepeatedpos )
{
    Tcol v, templist;

    for(auto& e: cv){
        v = e;
        sort(v.begin(), v.end());
        auto unique_end = unique(v.begin(), v.end()); //
        firstrepeatedpos=*unique_end;
        for (auto it = begin(v); it!=unique_end; ++it){ // for each unique element
            if(*it!= missingval){
                templist.push_back(*it);
            }
        }

        v.clear();
    } // close for cv

    sort(templist.begin(), templist.end());
    auto unique_end = unique(templist.begin(), templist.end());

    if(!uniquelist.empty()) uniquelist.clear();
    uniquelist.reserve(templist.size());

    for (auto it = begin(templist); it!=unique_end; ++it){ // for each unique element

        uniquelist.push_back(*it);

    }

  // Unique will take a sorted range, and move things around to get duplicated
   // items to the back and returns an iterator to the end of the unique section of the range

    //return (distance(templist.begin(), unique_end) ); // Constant time for random access iterators (like vector's)
    return(uniquelist.size());
}


// the vector alelids will contain the identifiers of the different alleles except the missing if any
// the vector freqs will contain the absolute number of copies of each element except the missing value
// the freqs positions corresponds to the sorted elements i.e. the first element in freq correspond to the minor element  and the freqs.back() corresponds to the highest-valued element in cv
// it stores in  n the sum of total elements including the missing value, so the difference between n - accumulate(freqs) will indicate if the missing value was present and the number of copies that were present
// it returns the number of distinct elements minus 1 if the missing val is present

// WARNING: NOTE THAT after unique, the elements between the returned iterator and last are left in a valid but unspecified state
// it it because of this that unique is not applied to cv but to a copy
template <typename T=int, typename Tcol=vector<T>, typename Tfreq>
int count_distinct_els_andfreqs(const Tcol& cv, T missingval,Tcol& alid, vector<Tfreq>& freqs, T&n)
{

    Tcol v=cv;
    n=v.size();
    sort(v.begin(), v.end());

  // Unique will take a sorted range, and move things around to get duplicated
   // items to the back and returns an iterator to the end of the unique section of the range
    auto unique_end = unique(v.begin(), v.end()); // Again n comparisons
    auto numids = distance(v.begin(), unique_end) - count(v.begin(), unique_end,missingval);
    if(!alid.empty()) alid.clear();
    alid.reserve(numids);

    if(!freqs.empty()) freqs.clear();

    for (auto it = begin(v); it!=unique_end; ++it){ // for each unique element
        if(*it!= missingval){
            alid.push_back(*it);
            freqs.push_back(count(begin(cv),end(cv), *it)); // count the number of occurrences
        }

    }

    return (numids); // Constant time for random access iterators (like vector's)
}

// returns the value occurring more times
template <typename T=int, typename Tcol=vector<T>>
T max_occurrence_el(const Tcol& v)
{

    int cont(0), maxc{0};
    T maxval;
    for(auto& val:v){

        if(cont< (maxc=count(v.begin(), v.end(),val)) ){
            cont=maxc;
            maxval=val;
        }
    }

    return maxval;
}

template <typename T=int>
inline void trim_left(string& s, const string& delimiters =" \f\n\r\t\v"){

    s.erase(0, s.find_first_not_of(delimiters) );

}

/// from stackoverflow: https://stackoverflow.com/questions/216823/whats-the-best-way-to-trim-stdstring?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa
// trim from start (in place)
static inline void ltrim(std::string &s) {
    s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](unsigned char ch) {
        return !std::isspace(ch);
    }));
}

// trim from end (in place)
static inline void rtrim(std::string &s) {
    s.erase(std::find_if(s.rbegin(), s.rend(), [](unsigned char ch) {
        return !std::isspace(ch);
    }).base(), s.end());
}

// trim from both ends (in place)
static inline void trim(std::string &s) {
    ltrim(s);
    rtrim(s);
}


//remove spaces in the string (in place)
// Code from: second upvoted in https://stackoverflow.com/questions/83439/remove-spaces-from-stdstring-in-c
static inline void clean_spaces(string& str){

    string::iterator end_pos = std::remove(str.begin(), str.end(), ' ');
    str.erase(end_pos, str.end());
}

// extract the part of the word at the left or right of the tag
/*
Usage
bool left=true;
cout<< substring(ofilenametxt, ".", left)<<endl;
cout<<substring(ofilenamexls, ".", !left)<<endl;

*/
/// Modified in Feb 2024 to add the last argument that allows to detect the tag from the beginning (last=false) or for the end (last=true)
template<typename T=string>
inline T substring(const string& word, const string& tag, const bool left=true, const bool last=false)
{
    string sub=word;
    size_t pos;

    if(last){
        pos= word.rfind(tag);
    }
    else{
        pos = word.find(tag);
    }
    if(pos!=string::npos){ // the tag is present
        if(left)
            sub=word.substr(0, pos);
        else
            sub=word.substr(pos);
    }

    return sub;
}


template <tstring T=string>
T ToUpperCase(const T& texto){
    T temp = texto;
    transform(temp.begin(),temp.end(), temp.begin(), (int(*)(int)) toupper); // transform to upper case
    return temp;
}

template <tstring T=string>
T ToLowerCase(const T& texto){
    T temp = texto;
    transform(temp.begin(),temp.end(), temp.begin(), (int(*)(int)) tolower); // transform to lower case
    return temp;
}

template <typename Tcol>
inline void ToZero(Tcol& col){
    fill(col.begin(), col.end(), 0);
}

// differences between assing and resize (resize only assing the new elements not the old ones): https://stackoverflow.com/questions/19673311/stl-vector-resize-and-assign
template <typename Tcol>
inline void ToZeroResize(const size_t newsize, Tcol& col){
    col.assign(newsize, 0);
}

/*
La interpolación lineal se basa en la idea de que, para interpolar un valor faltante en una serie de datos, se puede utilizar
la recta que pasa por los dos puntos más cercanos al valor faltante.
La fórmula para interpolar un valor y_i en una serie de datos y se basa en los dos puntos más cercanos
 x_{i-1} y x_{i+1} en la serie de datos x y es la siguiente:
y_i = y_{i-1} + (y_{i+1} - y_{i-1}) * (x_i - x_{i-1}) / (x_{i+1} - x_{i-1})
Donde i es el índice del valor faltante y x_i es el valor correspondiente en la serie de datos x.


#include <iostream>
#include <vector>

std::vector<int> interpolar(std::vector<int> x, std::vector<int> y) {
    std::vector<int> y_interp;

    // Iterar sobre los valores de y
    for (int i = 0; i < x.size(); i++) {
        // Si hay un valor faltante en y, interpolarlo
        if (i >= y.size()) {
            int y_i = y.back() + (y.back() - y[y.size()-2]) * (x[i] - x[y.size()-2]) / (x[y.size()] - x[y.size()-2]);
            y_interp.push_back(y_i);
        } else {
            y_interp.push_back(y[i]);
        }
    }

    return y_interp;
}

int main() {
    std::vector<int> x = {1, 2, 3, 4, 5, 6};
    std::vector<int> y = {1, 2, 3, 4};

    std::vector<int> y_interp = interpolar(x, y);

    for (int i = 0; i < y_interp.size(); i++) {
        std::cout << y_interp[i] << " ";
    }
    std::cout << std::endl;

    return 0;
}

*/
//se asume y,size() <x.size() y ambos ordenados de menor a  mayor
template <typename T, typename Tcol>
std::vector<T> interpolarY(Tcol x, Tcol y) {

    std::vector<int> y_interp;

    // Iterar sobre los valores de y
    for (int i = 0; i < x.size(); ++i) {
        // Si hay un valor faltante en y, interpolarlo
        if (i >= y.size()) {
            int y_i = y.back() + (y.back() - y[y.size()-2]) * (x[i] - x[y.size()-2]) / (x[y.size()] - x[y.size()-2]);
            y_interp.push_back(y_i);
        } else {
            y_interp.push_back(y[i]);
        }
    }

    return y_interp;
}


template <typename T, typename Tcol>
T QQplot(Tcol& x, Tcol &y){
    sort(x.begin(),x.end());
    sort(y.begin(),y.end());
}



#endif
