#include <iostream>
#include <sstream>
#include <fstream>
#include <list>
#include <map>
#include <queue>
#include <stack>
#include <cstring>
#include <cstdlib>
#include <string>
#include <algorithm>
#include <set>

#define MAX 1024
#define MAX_NB_NODES 5000000

using namespace std;

 typedef struct node{
    int number;
   char *direction;
  }NODE;


  ofstream file1; // for all the cycles
 ofstream myfile; // for the cycles in fasta format
 ofstream tabfile; // for the cycles in tab format file
ofstream xml_file; // for the cycles in xml format
  list<NODE> **adjacency_list;
  list<NODE> **adjacency_list_copy;
 map<string,int> labels; //hold labels of nodes
 map<int, int> coverage1;
 map<int, int> coverage2;
map<int, string> ingoing;
map<int,string> outgoing;
vector<int>cyclesnew;
map<int,string>allcycles; 
vector<int>cycles;
vector<int>cycle_reverse;
vector<int>switching;
vector<int>path1;
vector<int>path2;
map<int,string>sequences;
vector<int>switchnodes; // NEW switching nodes from all the graph
vector<int>switch_RF; // NEW switching nodes from all the graph of type RF
vector<int>switch_FR; // NEW switching nodes from all the graph of type FR
int cycles1=0;
int count_events=0;
int snps=0;
int k_value;
string expanding1;
 string expanding2;
/*Added this for cycle detection using Pilu's code*/
stack<int> marked;
stack<int> point;
bool* mark;
bool** adjacency;
int numberOfCycles;

/////////////////////////////////////////////////////// FUNCTION PROTOTYPES /////////////////////////////////////////////////
int read_nodes(char **correspondance);
char * Trim_char(char *input);
string *read_graph(list<NODE> **adjacency_list, int *nb_nodes, char *file_name);
int connected_components(int *E, list<int> **adjacency_list, int nb_nodes);
void find_nodes(vector<int>path, int start, int end, list<NODE>**adjacency_list, string *correspondance);
void allPaths(vector<int>path, int start, int end, list<NODE>**adjacency_list, string *correspondance);
int valid_path(char *direction1, char *direction2);
void Print_Path_twonodes(string * correspondance, int* predecessor, int start, int end);
void  Get_paths(vector<int>*path1,vector<int>*path2, vector<int>cycles, int start,int end);
int Get_Length(vector<int>path1, map<int,string>&sequence, string *correspondance);
void Switching_Nodes(vector<int>cycles, list<NODE>**adjacency_list, string *correspondance, int numofcycles);
void Print_Path(string * correspondance, int* predecessor, int i, int it);
void DFS_Visit(int i, list<NODE> **adjacency_list, int nb_nodes, string* color, int* predecessor,int time, string *correspondance, int* discovered, int* finished);
void DFS(list<NODE> **adjacency_list, int nb_nodes, string * correspondance);

void Get_sequences(string *correspondance, char *filename, map<int,string>&sequences, map<int,int>&coverage1,map<int,int>&coverage2 );



/////////////////////////////////////////////////////// FUNCTION DEFINITIONS  ////////////////////////////////////////////////
int read_nodes(char **correspondance)
{
  return 4;
}




///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// VALID PATH  /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


int valid_path(char *direction1, char *direction2){
  //  cout<< "To check:" << direction1 <<" " << direction2<<endl;
  if(strcmp(direction1,"RR")==0 || (strcmp(direction1,"FR")==0))
      if(strcmp(direction2,"RR")==0 || (strcmp(direction2,"RF")==0))
	return 1;
      else{ //cout << direction1 << direction2<<endl;
 return 0;
      }

 if(strcmp(direction1,"FF")==0 || (strcmp(direction1,"RF")==0))
      if(strcmp(direction2,"FF")==0 || (strcmp(direction2,"FR")==0))
	return 1;

      else{// cout << direction1 << direction2<<endl;
 return 0;
      }
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// DFS_VISIT -not used /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void DFS_Visit(int i, list<NODE> **adjacency_list, int nb_nodes, string* color, int* predecessor,int time, string *correspondance, int* discovered, int* finished){
color[i]="gray";
time++;

discovered[i]=time;
//cout << "discovered: "<< correspondance[i]<<endl;

for (list<NODE>::iterator it=adjacency_list[i]->begin(); it!=adjacency_list[i]->end(); it++){
if(color[it->number].compare("white")==0){
predecessor[it->number]=i;
 DFS_Visit(it->number,adjacency_list, nb_nodes,color,predecessor,time, correspondance, discovered, finished);
}

if(color[it->number].compare("gray")==0 && (predecessor[i]!=it->number)  ){

 Print_Path(correspondance, predecessor,it->number,i);


 file1 << "Cycle:" << cycles1 <<endl ;
 for(vector<int>::iterator i=cycles.begin();i!=cycles.end();i++){
   file1 << correspondance[*i]<< " ";
 }
 file1 <<" "<<endl;

 
 Switching_Nodes(cycles, adjacency_list, correspondance,cycles1);

cycles1++;
 cycles.clear();

}

}

color[i]="black";
finished[i]=time++;

}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    DFS - not used  ///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
void DFS(list<NODE> **adjacency_list, int nb_nodes, string * correspondance ){
string* color=new string[nb_nodes];
int* predecessor=new int[nb_nodes];
int *discovered=new int[nb_nodes];
int *finished=new int[nb_nodes];
for (int i=0;i<nb_nodes;i++) {
color[i]="white";
predecessor[i]=-1;
discovered[i]=-1;
finished[i]=-1;
}
int time=0;
int cycles;
for (int i=0;i<nb_nodes;i++) {
if(color[i].compare("white")==0)
  DFS_Visit(i, adjacency_list, nb_nodes,color,predecessor,time,correspondance, discovered, finished);

}


}


//void print_correspondance(int i, string *correspondance){
//cout << correspondance[i];
//}





//////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// REMOVE \t  from string/////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
char * Trim_char(char *input){
  int i=0;
  int length=strlen(input);
  //cout<< length;
  char *trimmed=new char[length];
  int j=0;
  for(i=0;i<length;i++){
    if(input[i]!='\t'){
    trimmed[j]=input[i];
    //  cout << trimmed[i];
j++;
    }
    
  }
  trimmed[j]='\0';
 
  // int length2= strlen(trimmed);
  // cout << trimmed << length2 <<endl;
  return trimmed;
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// REVERSE COMPLEMENT /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

static char complement(char b)
{
switch(b)
{
case 'A': return 'T';
case 'T': return 'A';
case 'G': return 'C';
case 'C': return 'G';

case 'a': return 't';
case 't': return 'a';
case 'g': return 'c';
case 'c': return 'g';
 }
return '?';
 
}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// REVERSE COMPLEMENT /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

string  revcomp(string seq){

string s(seq.begin(),seq.end());

    // iterate through all of the characters
    string::iterator pos;
    for (pos = s.begin(); pos != s.end(); ++pos) {
      // cout << *pos;
    }
    // cout << endl;

    
    reverse (s.begin(), s.end());
    //  cout << "reverse:       " << s << endl;
 
    for(pos=s.begin();pos!=s.end();++pos){

      *pos=complement(*pos);
    }

 return s;
}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// READ GRAPH /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*reads a graph and returns the number of nodes*/
string *read_graph(list<NODE> **adjacency_list, int *nb_nodes, char *file_name)
{
NODE *mynode=new NODE;
  char *str = new char[2*MAX];
  char *str1 = new char[MAX];
  char *str2 = new char[MAX];
 char *sequence = new char[MAX];
char *str3=new char[3];
 string untabbed;
 char *strdirect=new char[2];
  char *tmp_str;
  int k=0; //nb_node iterator
  int n1,n2=0;

 
map<string, int> directions;
 
  std::ifstream file_op(file_name,ios::in);
  if (!file_op.is_open())
    {
      cerr<<"unable to open file, did you specify one ?"<<endl;
      cerr<<"usage is ./cyclepaths input_file1 input_file2"<<endl;
      cerr<<endl;
      exit(0);
    }



  while(!file_op.eof())
    {
      file_op.getline(str,2 * MAX);
      int i=0;
      /*get first string*/
      while (str[i]!='\t')// || str[i]==' ')
	{
	  str1[i]=str[i];
	  i++;
	}
      str1[i]='\0';
      i++;
      /*get_second string*/
      int j=0;
      while (str[i]!='\t')
	{
	  str2[j]=str[i];
	  i++;
	  j++;
	}
      str2[j]='\0';

 /*get third string RF, RR, FF. FR*/
int w=0;
while (str[i]!='\0')
  { 
	  str3[w]=str[i];
	  i++;
	  w++;
	}
      str3[w]='\0';




    
      map<string,int>::iterator it;
      it = labels.find(str1);
      // cout << str1 << "\t"<< str2 <<"\t";
      //  cout<< str3<<endl;
  /*process string 1*/
      if((labels.find(str1)) ==  (labels.end()))
	{
	  /*the node has not been seen before, we add it*/
	  tmp_str=new char[MAX];
	  tmp_str=strcpy(tmp_str,str1);
	  labels[tmp_str]=k;
	  n1=k;
	  k++;
	}

      else{
	  n1=it->second;
      }
      /*process string 2*/
      it=labels.find(str2);
      if(it ==  (labels.end()))
	{
	  /*the node has not been seen before, we add it*/
	  tmp_str=new char[MAX];
	  tmp_str=strcpy(tmp_str,str2);
	  labels[tmp_str]=k;
	  n2=k;
	  k++;
	}

      else{
	n2=it->second;
      }


      
      mynode->number=n2;
    

       mynode->direction=new char[2];
       char *trimmed=Trim_char(str3);
       strcpy(mynode->direction,trimmed);

      /*add edge*/
       adjacency_list[n1]->push_back(*mynode);
       

       mynode->number=n1;
            
    }



      map<string,int>::iterator it;
      it = labels.find(str1);



  /*fill correspondance array*/
  string *correspondance = new string[labels.size()];

  for(map<string,int>::iterator iter = labels.begin();iter!=labels.end();iter++)
    correspondance[(*iter).second] = (*iter).first;

 
  file_op.close();

  //(*nb_nodes)=k;

  /*VL220410*/

  (*nb_nodes)=k-1;

  return correspondance;




}




//////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// connected components- not used//////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/**
 * This method gives the number of connected components of a graph. 
 * It  is based on a breadth-first search.
  * @param E array of integers indicating which connected component each node belongs to. This array is initialised to 0 and is filled during the bfs.
  * @param adjacency_list graph stored as adjacency lists
  * @param nb_nodes number of nodes of the graph
  * @return nb_connected_components number of connected components of the graph
 */

int connected_components(int *E, list<int> **adjacency_list, int nb_nodes)
{
  queue<int> Q;
  bool* visited=new bool[nb_nodes];
bool* finished= new bool[nb_nodes];

  for (int i=0;i<nb_nodes;i++)  
    visited[i]=0;

  int c=0;

  for (int i=0;i<nb_nodes;i++)
    if (!visited[i])
      {
	/*bfs*/
	Q.push(i);
	visited[i]=true;
	while(!Q.empty())
	  {
	    int node_to_process=Q.front();
	    Q.pop();
	    /*process node, i.e. determine its component*/
	   E[node_to_process]=c;
	   finished[node_to_process]=true;
	    /*add unvisited neighbors to the queue*/	    
	    for (list<int>::iterator it=adjacency_list[node_to_process]->begin(); it!=adjacency_list[node_to_process]->end(); it++){
	     
		if (!visited[*it])
		  {
		    visited[*it]=true;
		    Q.push(*it);
		
		  }
	      }
	  }
	/*end of bfs*/
	/*next iteration of this loop will correspond to a new connected component*/
	c++;
      }

  return c;
  
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////GET 2 PATHS //////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
/* Find the upper and lower path, given the start switching node and the end switching nodes. Returns path1 and path2, upper and lower path,  */

void  Get_paths(vector<int>*path1, vector<int>*path2, vector<int>cycles, int start,int end){
  bool found=0;  
  vector<int>path3;
  vector<int>path4;
  vector<int>::iterator it=cycles.begin();

 while(*it!=start && it!=cycles.end()){
      path3.push_back(*it);
      it++;
    }
  if(*it==start){
    path3.push_back(*it);
    path2->push_back(*it);
    it++;
  }
 
  while(*it!=end && it!=cycles.end()){
  path2->push_back(*it);
  //cout << *it<< " ";
  it++;
  }
  if(*it==end){
    path4.push_back(*it);
    path2->push_back(*it);  
    it++;
}
  while(it!=cycles.end()){
    path4.push_back(*it);
  
      it++;
  }

  

   for( vector<int>::iterator it=path3.end()-1; it != path3.begin()-1; it--)
     path1->push_back(*it);

 for( vector<int>::iterator it=path4.end()-1; it != path4.begin()-1; it--){

     path1->push_back(*it);
 }


}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// Merge Sequence ////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* Merges the sequence of the nodes in the path -used by next function, Get_Length*/

void Merge_Sequence(string& seq1, string sequence1, string seq2, string direction,int k){

 string add;
  string reverse;
  string reverse2;
  if (strcmp(direction.c_str(), "FF")==0){
 
    //if(sequence1.substr(sequence1.length()-k+1)!=seq2.substr(0,k-1)) cout << "NOT EQUAL"<<endl;
 add=seq2.substr(k-1);
//file1 << " toadd  " <<  add << endl;
   seq1.append(add);
//   file1 <<  "seq1: " << seq1<< endl;
  }

 if (strcmp(direction.c_str(), "FR")==0){
    // file1  << direction << endl;
 reverse=revcomp(seq2);
 //file1 << "seq1:"<< sequence1<< endl;
 //file1 << "seq2:" << seq2<<endl;
 //file1 << "seq2reverse:" << reverse <<endl;
 //file1 <<"seq1_common:" << sequence1.substr(sequence1.length()-k+1)<< endl;
 //file1 << "seq2_common:"<< reverse.substr(0,k-1)<< endl;
 //if(sequence1.substr(sequence1.length()-k+1)!=reverse.substr(0,k-1)) cout << "NOT EQUAL"<<endl;
 add=reverse.substr(k-1);
 //file1<< " toadd  " <<  add << endl;
   seq1.append(add);
 //  file1 << "seq1: " << seq1<< endl;
  }
 if (strcmp(direction.c_str(), "RF")==0){
 // file1 << direction << endl;
 add=seq2.substr(k-1);
 reverse=revcomp(sequence1);
 //file1 << "seq1:"<< sequence1 <<" " << reverse << endl;
 //file1 << "seq2:" << seq2 << endl;
 //file1 <<"seq1_common:" << reverse.substr(reverse.length()-k+1)<< endl;
 //file1 << "seq2_common:"<< seq2.substr(0,k-1)<< endl;
 //if(reverse.substr(reverse.length()-k+1)!=seq2.substr(0,k-1)) cout << "NOT EQUAL"<<endl;
 //file1 << " toadd  " <<  add << endl;
   seq1.append(add);
  // file1 << "seq1: " << seq1<< endl;


  }
 if (strcmp(direction.c_str(), "RR")==0){
 //  file1 <<  direction << endl;
 reverse=revcomp(seq2);

 add=reverse.substr(k-1);
 reverse2=revcomp(sequence1);
//  file1 << "seq1:"<< sequence1 << " "<< reverse2 << endl;
// file1 << "seq2:" << seq2<< " " << reverse << endl;
//file1 <<"seq1_common:" << reverse2.substr(reverse2.length()-k+1)<< endl;
// file1 << "seq2_common:"<< reverse.substr(0,k-1)<< endl;
 //if(reverse2.substr(reverse2.length()-k+1)!=reverse.substr(0,k-1)) cout << "NOT EQUAL"<<endl;
 //file1 << " toadd  " <<  add << endl;
   seq1.append(add);
 //  file1 << "seq1: " << seq1<< endl;
  }



}



///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// Get Length /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* finds the length of a path of a cycle */

int Get_Length(vector<int> path,  map<int, string>&sequences, string *correspondance,  list<NODE>**adjacency_list, string& expanding)
{

   int length=0;
  int current=0;  
   map<int, char*>::iterator iter;
  map<string,int>::iterator it1;
for( vector<int>::iterator it=path.begin(); it != path.end(); it++){
  
   
  int nu1=atoi(correspondance[*it].c_str());

 }

char *direct1=new char[2];
  strcpy(direct1,"00");
  int first_time=0;  
  int one=cycles.front();
  //  cycles.push_back(one);
for(int i=0;i+1 < path.size();i++){
  int num1=path[i];
  int num2=path[i+1];
  // cout <<"num1: " <<correspondance[num1]<< " " <<"num2: "<< correspondance[num2] << endl;
int n1=atoi(correspondance[num1].c_str());
int n2=atoi(correspondance[num2].c_str());
// cout << sequences[n1] << endl;
 
 for (list<NODE>::iterator it=adjacency_list[num1]->begin(); it!=adjacency_list[num1]->end(); it++){
 
  if(num2==it->number){
    //   cout<< it->direction << endl;
   
    direct1=it->direction;
    if(first_time==0){
      if (strcmp((it->direction), "RF")==0 || strcmp((it->direction), "RR")==0 ){
 	expanding=revcomp(sequences[n1]);}
    else {expanding=sequences[n1];}
     // file1 << expanding <<endl;
 first_time++;
    }
    // cout << "dsd"<< first_time << endl;
    Merge_Sequence(expanding,sequences[n1], sequences[n2], it->direction,  k_value);
   
  
  }    
 

 
 }

 }
length=expanding.length();
//cout << "Path  sequence: " << expanding << endl;
 

//cout << " with length "<< length << endl;
  return length;

}



///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////// PRINT PATH - not used //////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////

void Print_Path(string * correspondance, int* predecessor, int i, int it){
//cout << "Print_path"<<endl;
  if (i==it){
    //  cout <<  correspondance[i]<< " ";
 cycles.push_back(i);
  }
else if (predecessor[it]==-1)
  { //cout << correspondance[it]<< " ";
cycles.push_back(it);
  // cout << "direction: " <<  direction[it]<<" ";
}
else {

  Print_Path(correspondance, predecessor, i, predecessor[it]);

  //cout << correspondance[it] <<" ";
cycles.push_back(it);
}
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    GET SEQUENCES -/////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* gets the sequences of each node from the graph file containing the nodes and their sequences */

void Get_sequences(string *correspondance, char *filename, map<int,string>& sequences, map<int,int>& coverage1,map<int,int>& coverage2)
{
  char *str = new char[10*MAX];
  char *str1 = new char[10*MAX];
  char *str2 = new char[10*MAX];
  char *sequence = new char[10*MAX];
  char *number1=new char[10*MAX];
  char *number2=new char[10*MAX];
 
  std::ifstream file_op(filename,ios::in);

  if (!file_op.is_open())
  {
    cerr<<"unable to open file, did you specify one ?"<<endl;
    cerr<<"usage is ./paths graph_edges_400k25.txt graph_nodes_400k25.txt 25 sequences_25.fa events_new25.txt events_tab_25.txt"<<endl;
    cerr<<endl;
    exit(0);
  }

  while(!file_op.eof())
  {
    file_op.getline(str,10 * MAX);
    if (strlen(str) != 0)
    {  
        int i=0;
	/*get first string label*/
	while (str[i]!='\t')// || str[i]==' ')
	{
	  str1[i]=str[i];
	  i++;
	}
	str1[i]='\0';
	i++;
	/*get_second string sequence*/
	int j=0;
	while (str[i]!='\t')
	{
	  sequence[j]=str[i];
	  i++;
	  j++;
	}
	i++;
	sequence[j]='\0';
	int lengthofseq=strlen(sequence);

	map<string,int>::iterator it;
	it = labels.find(str1);

	int n1 =atoi(str1); 
    
	// third string-ignore
	while (str[i]!='\t')
	{
	  i++;
	}
	//get fourth string-coverage1
	i++;
	int w=0;
	while (str[i]!='\t')
	{
	  number1[w]=str[i];
	  i++;
	  w++;
	}
	number1[w]='\0';
	// cout << number1 <<endl;
	int cov1=atoi(number1);
      
	// cout << cov1 <<endl;

	//get fifth string-coverage2
	int z=0;
	while (str[i]!='\0')
	{
	  number2[z]=str[i];
	  i++;
	  z++;
	}
	number2[z]='\0';
	// cout << number1 <<endl;
	int cov2=atoi(number2);

	coverage1.insert(pair<int, int>(n1, cov1));
	coverage2.insert(pair<int, int>(n1, cov2));
	//	  n1=it->second;
	sequences[n1]=new char[lengthofseq];
	// cout << "coverage1 " << cov1;
	// cout << "coverage2 " << cov2;
	sequences[n1]=sequence;
	sequences.insert(pair<int, string>(n1,sequence));
    }
  }

  file_op.close();

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////   END OF FILE functions.h /////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////