#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cstring>
#include <cmath>
#include <stdlib.h>
#include <algorithm>
#include <sstream>
#include <iterator>
#include <map>
#include <time.h>
#include "Array"
#include "Core"
#include "Dense"
using namespace Eigen;
using namespace std;
/*


Contains the source code for haploPS alpha version. 
//csvline_populate is used to parse input with a delimiter.
//
required library: Eigen. 
If you want to compile the source code, please install Eigen. The compilation command using g++ is as follows:
 g++ -L ~/library_path/Eigen/ -Wall HaploPS_version2.cpp 

Any question or command, please email: xuanyao.liu@nus.edu.sg



*/

//function 1, parsing###################################################################################3
//void _csvline_populate(vector<string> &record, const string& line, char delimiter);

//////function parsing delimited files//////////////////

//////function parsing delimited files//////////////////
void _csvline_populate(vector<string> &record, const string& line, char delimiter)
{
    int linepos=0;
    int inquotes=false;
    char c;
    int linemax=line.length();
    string curstring;
    record.clear();
       
    while(line[linepos]!=0 && linepos < linemax)
    {
       
        c = line[linepos];
       
        if (!inquotes && curstring.length()==0 && c=='"')
        {
            //beginquotechar
            inquotes=true;
        }
        else if (inquotes && c=='"')
        {
            //quotechar
            if ( (linepos+1 <linemax) && (line[linepos+1]=='"') ) 
            {
                //encountered 2 double quotes in a row (resolves to 1 double quote)
                curstring.push_back(c);
                linepos++;
            }
            else
            {
                //endquotechar
                inquotes=false; 
            }
        }
        else if (!inquotes && c==delimiter)
        {
            //end of field
            record.push_back( curstring );
            curstring="";
        }
        else if (!inquotes && (c=='\r' || c=='\n') )
        {
            record.push_back( curstring );
            return;
        }
        else
        {
            curstring.push_back(c);
        }
        linepos++;
    }
    record.push_back( curstring );
    return;
}

////////////////end of parse program//////////////////////

////////////////end of parse program//////////////////////
////////////////function2 maxcount program/////////////////////////
struct by_second
{
    template <typename Pair>
    bool operator()(const Pair& a, const Pair& b)
    {
        return a.second < b.second;
    }
};


template <typename Fwd>
typename std::map<typename std::iterator_traits<Fwd>::value_type, int>::value_type
most_frequent_element(Fwd begin, Fwd end)
{
    std::map<typename std::iterator_traits<Fwd>::value_type, int> count;

    for (Fwd it = begin; it != end; ++it)
        ++count[*it];

    return *std::max_element(count.begin(), count.end(), by_second());
}
////////////////////////////end of maxcount////////////////////////////////////////
////////////////////////main program///////////////////////////////////////////////////
int main(int argc,char *argv[]){
            clock_t t1,t2;
            t1=clock();
            string inputfile;//input genotype file
            string legendfile;
            float freq;
            string outputfile;
            assert(argc ==9);
            
            inputfile=argv[2]; //genotype file
            legendfile=argv[4]; //legend file
            freq=atof(argv[6]); //search haplotype frequency
            outputfile=argv[8];
            
            ofstream outf;
            outf.open(outputfile.c_str());
            outf<<"chr"<<"\t"<<"startpos"<<"\t"<<"endpos"<<"\t"<<"genetic_dist"<<"\t"<<"physical_dist"<<"\t"<<"n.snp"<<endl;


            int rowcount =0;
            stringstream ss;
            int snptotal;
            int posflag; 
            int extra;
            int extraleft=0;
            int extraright=0;
            int numchr;
            double thres=0.98;//threshold for similarity index
            vector<string> row;
            string line;
            string header;//header of table
            int i;
            int j;
            //int m;
//////read in haplotype file///////////////
                  ifstream in(inputfile.c_str());
                  if (in.fail())  { cout << "Input genotype file not found "<<inputfile <<endl; return 0; }
///////read in first line & number of columns///////////
                  getline(in,header);//remove header
                  _csvline_populate(row, header, '\t'); //count number of SNPs if tab delimited
                  if(row.size()==1) _csvline_populate(row, header, ' '); //count number of SNPs if space delimited
                  snptotal=row.size();
////read every line of the haplotype file, & number of rows///////////
	     		  while(getline(in, line)  && !in.eof() )/////////need to break when exceeding window size
			                     {  
			        
			                                      rowcount++;
                                                    }
                      in.close();
                      numchr=rowcount+1; //number of chromosomes


///////////////////////////////read in haplotype file into 2d dynamic arrays//////////////////////////////////


                  int **hapData=0;
                  int numSample=rowcount+1;
                  int numMarker=snptotal;

                  hapData=new int *[numSample];
                  for(i=0;i<numSample;i++)
                      hapData[i]=new int[numMarker];


                      FILE *infp;


                      if((infp=fopen(inputfile.c_str(),"r"))==NULL){
                              printf("cannot open input file%s","hap");
                              exit(1);
                      }


                 for(i=0;i<numSample;i++)
                     {

                              for(int k=0;k<numMarker;k++)
                              {
                                   fscanf(infp,"%d",&hapData[i][k]);

                              }
                     }

                fclose(infp);

 
                     posflag=0;
                     while(posflag<snptotal-1){
//extend the haplotype to left and right at extra=10*i////////////////////////////////////////////////////////////////
                     extra=10;
                     float crithap=floor(freq*numchr);
                     double simindex;
                     int simcount=numchr;                     
                     while(simcount>crithap){  
//make extras into string & find the highest occurrence and the rownumber & canonical haplotype
                     string haplo;
                     string canonicalhap;
                     int canonicalhaprow;
                     vector<string> hapblock;
                     hapblock.clear();
                     int haplo_start=max(posflag-extra,0);//rows in c++
                     int haplo_end=min(posflag+extra,snptotal-1);
                  
                     int ncolmatrix=haplo_end-haplo_start+1; 
                     MatrixXi hapblock_matrix(numchr,ncolmatrix);
////read every line of the haplotype file, & number of rows///////////
                     rowcount=0;
	     		      while(rowcount<numSample)
			                     { 
		                              for(i=haplo_start;i<=haplo_end;i++){
                                                                        haplo+=hapData[rowcount][i];
                                                                        hapblock_matrix(rowcount,i-haplo_start)=hapData[rowcount][i];
                                                                          }
                                              hapblock.push_back(haplo);
                                  
                                              haplo.clear();
                                              rowcount++;
                                                    }//end of while
                      std::pair<string, int> x = most_frequent_element(hapblock.begin(), hapblock.end());
                      canonicalhap=x.first;
                      canonicalhaprow=std::find(hapblock.begin(),hapblock.end(),canonicalhap)-hapblock.begin();
                      
//calculate distances, and get frequency
                      simcount=0;
                      for(i=0;i<numchr;++i){
                                          simindex=0;
                                          for(j=0;j<hapblock_matrix.cols();++j){
                                                                                simindex+=(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))*(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j));
                                                                                }      
                                                simindex=1-simindex/hapblock_matrix.cols();
                                                if(simindex>thres) {simcount++;}
                                                }
                                              extra=extra+10;  
                      }//end of left right extra
//iterates
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//extend haplotype to the right at 10 snps
                  extra=extra-20;
                  int positiveflag=1;
                  while(positiveflag>0){
                 
                  string haplo;
                  string canonicalhap;
                  int canonicalhaprow;
                  vector<string> hapblock;
                  hapblock.clear();
                  int haplo_start=max(posflag-extra,0);//rows in c++
                  int haplo_end=min(posflag+extra+10*positiveflag,snptotal-1);
                  if(haplo_end-haplo_start+1>1){
                           int ncolmatrix=haplo_end-haplo_start+1; 
                           MatrixXi hapblock_matrix(numchr,ncolmatrix);
                
                           rowcount=0;
	     		   while(rowcount<numSample)
			                     { 
		                              for(i=haplo_start;i<=haplo_end;i++){
                                                                        haplo+=hapData[rowcount][i];
                                                                        hapblock_matrix(rowcount,i-haplo_start)=hapData[rowcount][i];
                                                                          }
                                              hapblock.push_back(haplo);
                                  
                                              haplo.clear();
                                              rowcount++;
                                                    }//end of while
                    
                std::pair<string, int> x = most_frequent_element(hapblock.begin(), hapblock.end());
                canonicalhap=x.first;
                canonicalhaprow=std::find(hapblock.begin(),hapblock.end(),canonicalhap)-hapblock.begin();
                      
                      
//calculate distances, and get frequency
                      simcount=0;
                      for(i=0;i<numchr;++i){
                                          simindex=0;
                                          for(j=0;j<hapblock_matrix.cols();++j){
                                                                                simindex+=(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))*(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j));
                                                                                }      
                                                simindex=1-simindex/hapblock_matrix.cols();
                                                if(simindex>thres) {simcount++;}
                                                }
                       if(simcount>crithap)   positiveflag++;
                       if(simcount<=crithap || haplo_end==snptotal-1) {
                                            extraright=extra+10*(positiveflag-1);
                                            positiveflag=-1;                                                            
                                            }
             }//end of if nsnp>!
             if(haplo_end-haplo_start+1==1){
                                                       positiveflag=-1;
                                                       extraright=extra;
                                                       }                     
                                             }//while positive flag

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
//extra to the left at 10 snps

                  int negativeflag=1;
                  while(negativeflag>0){
                  string haplo;
                  string canonicalhap;
                  int canonicalhaprow;
                  vector<string> hapblock;
                  hapblock.clear();
                  int haplo_start=max(posflag-extra-10*negativeflag,0);//rows in c++
                  int haplo_end=min(posflag+extraright,snptotal-1);
                  if(haplo_end-haplo_start+1>1){
                            int ncolmatrix=haplo_end-haplo_start+1; 
                            MatrixXi hapblock_matrix(numchr,ncolmatrix);
                
////read every line of the genotype file, & number of rows///////////
                  rowcount=0;
	    	  while(rowcount<numSample)
			                     { 
	        	                          for(i=haplo_start;i<=haplo_end;i++){
                                                                        haplo+=hapData[rowcount][i];
                                                                        hapblock_matrix(rowcount,i-haplo_start)=hapData[rowcount][i];
                                                                          }
                                                  hapblock.push_back(haplo);
                                  
                                                  haplo.clear();
                                                  rowcount++;
                                              }//end of while
                    
                      std::pair<string, int> x = most_frequent_element(hapblock.begin(), hapblock.end());
                      canonicalhap=x.first;
                      canonicalhaprow=std::find(hapblock.begin(),hapblock.end(),canonicalhap)-hapblock.begin();
                      
                      
//calculate distances, and get frequency
                      simcount=0;
                      for(i=0;i<numchr;++i){
                                          simindex=0;
                                          for(j=0;j<hapblock_matrix.cols();++j){
                                                                                simindex+=(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))*(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j));
                                                                                }      
                                                simindex=1-simindex/hapblock_matrix.cols();
                                                if(simindex>thres) {simcount++;}
                                                }
                       if(simcount>crithap)   negativeflag++;
                       if(simcount<=crithap || haplo_start==0) {
                                            extraleft=extra+10*(negativeflag-1);
                                            negativeflag=-1;                                                            
                                            }
                    }//end of nsnp
                        if(haplo_end-haplo_start+1==1){
                                                       negativeflag=-1;
                                                       extraleft=extra;
                                                       }                     
                                                                                                            
                                             
                                             }//while negative flag
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//extra to the right at 1 snp
                      
                       positiveflag=1;
                       while(positiveflag>0){
                   string haplo;
                  string canonicalhap;
                  int canonicalhaprow;
                  vector<string> hapblock;
                  hapblock.clear();
                  int haplo_start=max(posflag-extraleft,0);//rows in c++
                  int haplo_end=min(posflag+extraright+positiveflag,snptotal-1);
                  if(haplo_end-haplo_start+1>1){
                  int ncolmatrix=haplo_end-haplo_start+1; 
                  MatrixXi hapblock_matrix(numchr,ncolmatrix);
                
////read every line of the genotype file, & number of rows///////////
                  rowcount=0;
	     		  while(rowcount<numSample)
			                     { 
		                          for(i=haplo_start;i<=haplo_end;i++){
                                                                        haplo+=hapData[rowcount][i];
                                                                        hapblock_matrix(rowcount,i-haplo_start)=hapData[rowcount][i];
                                                                          }
                                   hapblock.push_back(haplo);
                                  
                                   haplo.clear();
                                  //matrix of hapblock
                                  /**/ 
                                  rowcount++;
                                                    }//end of while
                      std::pair<string, int> x = most_frequent_element(hapblock.begin(), hapblock.end());
                      canonicalhap=x.first;
                      canonicalhaprow=std::find(hapblock.begin(),hapblock.end(),canonicalhap)-hapblock.begin();
                      //cout<<canonicalhap<<" happens at position "<<canonicalhaprow<<endl;
                      
                      
//calculate distances, and get frequency
                      simcount=0;
                      for(i=0;i<numchr;++i){
                                          simindex=0;
                                          for(j=0;j<hapblock_matrix.cols();++j){
                                                                                //cout<<hapblock_matrix(i,j);
                                                                                simindex+=(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))*(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j));
                                                                                //(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))
                                                                                }      
                                                                                //cout<<endl;
                                                simindex=1-simindex/hapblock_matrix.cols();
                                                //cout<<simindex<<endl;
                                                if(simindex>thres) {simcount++;}
                                                }
                       if(simcount>crithap)   positiveflag++;
                       if(simcount<=crithap || haplo_end==snptotal-1) {
                                            extraright=extraright+(positiveflag-1);
                                            positiveflag=-1;                                                            
                                            }
             }//end of if nsnp>!
             if(haplo_end-haplo_start+1==1){
                                                       positiveflag=-1;
                                                       extraright=extraright;
                                                       }                     
                                                             

                                             
                                             }//while positive flag
////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////
//extra to the left at 1 snp
                       negativeflag=1;
                       while(negativeflag>0){
                   string haplo;
                  string canonicalhap;
                  int canonicalhaprow;
                  vector<string> hapblock;
                  hapblock.clear();
                  int haplo_start=max(posflag-extraleft-negativeflag,0);//rows in c++
                  int haplo_end=min(posflag+extraright,snptotal-1);
            if(haplo_end-haplo_start+1>1){
                  int ncolmatrix=haplo_end-haplo_start+1; 
                  MatrixXi hapblock_matrix(numchr,ncolmatrix);
                  //cout<<snptotal<<endl;
                  //cout<<haplo_start<<" "<<haplo_end<<" "<<rowcount+1<<" "<<ncolmatrix<<endl; 
                
////read every line of the genotype file, & number of rows///////////
                  rowcount=0;
	     	  while(rowcount<numSample)
			                     { 
		                          for(i=haplo_start;i<=haplo_end;i++){
                                                                        haplo+=hapData[rowcount][i];
                                                                        hapblock_matrix(rowcount,i-haplo_start)=hapData[rowcount][i];
                                                                          }
                                   hapblock.push_back(haplo);
                                  
                                   haplo.clear();
                                  //matrix of hapblock
                                  /**/ 
                                  rowcount++;
                                                    }//end of while 
                      std::pair<string, int> x = most_frequent_element(hapblock.begin(), hapblock.end());
                      canonicalhap=x.first;
                      canonicalhaprow=std::find(hapblock.begin(),hapblock.end(),canonicalhap)-hapblock.begin();
                      //cout<<canonicalhap<<" happens at position "<<canonicalhaprow<<endl;
                      
                      
//calculate distances, and get frequency
                      simcount=0;
                      for(i=0;i<numchr;++i){
                                          simindex=0;
                                          for(j=0;j<hapblock_matrix.cols();++j){
                                                                                //cout<<hapblock_matrix(i,j);
                                                                                simindex+=(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))*(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j));
                                                                                //(hapblock_matrix(i,j)-hapblock_matrix(canonicalhaprow,j))
                                                                                }      
                                                                                //cout<<endl;
                                                simindex=1-simindex/hapblock_matrix.cols();
                                                //cout<<simindex<<endl;
                                                if(simindex>thres) {simcount++;}
                                                }
                       if(simcount>crithap)   negativeflag++;
                       if(simcount<=crithap || haplo_start==0) {
                                            extraleft=extraleft+(negativeflag-1);
                                            negativeflag=-1;                                                            
                                            }
                    }//end of nsnp
                        if(haplo_end-haplo_start+1==1){
                                                       negativeflag=-1;
                                                       extraleft=extraleft;
                                                       }                     
                                                                                                            
                                                           
                                             
                                             }//while negative flag
////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////

int hapfinal_start=max(0,posflag-extraleft);
int hapfinal_end=min(max(min(snptotal, posflag + extraright), posflag + 1), snptotal);

rowcount=0;
int startpos=0;
int endpos=0;
float genetic_start=0;
float genetic_end=0;

///read in legend fiile
                ifstream in3(legendfile.c_str());
                if (in3.fail())  { cout << "Input legend file not found: "<<legendfile <<endl; 
                                   return 0; 
                                   }
                getline(in3,header);//remove header
                     while(getline(in3, line)  && in3.good() )/////////need to break when exceeding window size
                     {  
  
                      if(rowcount==hapfinal_start){
                                             _csvline_populate(row, line, '\t');
                                              if(row.size()==1)   _csvline_populate(row, line, ' ');  
                                               startpos=atoi(row[3].c_str());
                                               genetic_start=atof(row[2].c_str());                                           
                                                   }//if start
                      if(rowcount==hapfinal_end){
                                             _csvline_populate(row, line, '\t');
                                              if(row.size()==1)   _csvline_populate(row, line, ' ');  
                                               endpos=atoi(row[3].c_str());
                                               genetic_end=atof(row[2].c_str());                                           
                                                   }//if start
                                                       
                      
                       rowcount++;
                                  
                           }//end of while read
                        in3.close();//close infile
    outf <<row[0]<<"\t"<<startpos<<"\t"<<endpos<<"\t"<<genetic_end-genetic_start<<"\t"<<endpos-startpos;
    outf<<"\t"<<hapfinal_end-hapfinal_start+1<<endl;

    if(hapfinal_end+1<snptotal-1) posflag=hapfinal_end+2;
    else posflag=snptotal*10;
}
////////////////////////////free memory/////////////////////////


for(i=0;i<numSample;i++)
delete [] hapData[i];
delete [] hapData;



t2=clock();
return EXIT_SUCCESS;
}