#include <iostream>
#include <algorithm>
#include <iterator>
int main(int argc, char **argv)
{
const char *infilename = argv[1];
vector<Sequence::Fasta> data;
Sequence::Alignment::GetData(data,infilename);
Sequence::Alignment::validForPolyAnalysis(data.begin(),data.end()) )
{
cout << "Original position order:\n";
copy(SNPs.pbegin(),SNPs.pend(),
ostream_iterator<double>(cout," "));
cout << endl;
random_shuffle(SNPs.pbegin(),SNPs.pend());
cout << "Permuted positions:\n";
copy(SNPs.pbegin(),SNPs.pend(),
ostream_iterator<double>(cout," "));
cout <<'\n'<< endl;
cout << "Original data:\n";
copy(SNPs.begin(),SNPs.end(),ostream_iterator<string>(cout,"\n"));
cout <<'\n'<< endl;
random_shuffle(SNPs.begin(),SNPs.end());
cout << "Permuted data:\n";
copy(SNPs.begin(),SNPs.end(),ostream_iterator<string>(cout,"\n"));
cout << endl;
double pi = 0.;
Sequence::PolySites::const_site_iterator itr = SNPs.sbegin();
while( itr < SNPs.send() )
{
unsigned A = count(itr->second.begin(),itr->second.end(),'A');
unsigned G = count(itr->second.begin(),itr->second.end(),'G');
unsigned C = count(itr->second.begin(),itr->second.end(),'C');
unsigned T = count(itr->second.begin(),itr->second.end(),'T');
unsigned N = count(itr->second.begin(),itr->second.end(),'N');
double SH = 1.;
unsigned n = itr->second.length() - N;
SH -= (A>0) ? (double(A)/double(n))*(double(A-1)/double(n-1)) : 0.;
SH -= (G>0) ? (double(G)/double(n))*(double(G-1)/double(n-1)) : 0.;
SH -= (C>0) ? (double(C)/double(n))*(double(C-1)/double(n-1)) : 0.;
SH -= (T>0) ? (double(T)/double(n))*(double(T-1)/double(n-1)) : 0.;
pi += SH;
++itr;
}
cout << "Pi (for the entire region) = "<< pi << endl;
}
}