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

#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
list<NODE> **adjacency_list;
list<NODE> **adjacency_list_copy;
map<string,int> labels;
map<int, int> coverage1;
map<int, int> coverage2;
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;
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;
}


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

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// 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 /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// TODO: replace by an array of lenght 256 : return revcomp[b]
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) { // FIXME: remove this
    // 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;
    
      //untabbed.erase(std::remove(untabbed.begin(),untabbed.end(),'\t'),untabbed.end());
      //      strcpy(mynode->direction,"asa");
      //length=str.copy(buffer,6,5);
      mynode->direction=new char[2];
      char *trimmed=Trim_char(str3);
      strcpy(mynode->direction,trimmed);
      //cout <<mynode->number << " "<< mynode->direction << "\t"<< endl;
      /*add edge*/
      adjacency_list[n1]->push_back(*mynode);
       
      // NODE *mynode2=new NODE;
      mynode->number=n1;
       
      // strcpy(mynode2->direction,mynode->direction);
      //  adjacency_list[n2]->push_back(*mynode);

      //      adjacency_list-[n1].number->push_back(n2);
      //adjacency_list[n2]->push_back(n1);
 
     
    }



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



  /*fill correspondance array*/
  string *correspondance = new string[labels.size()];
  // direction=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*/
  /*I think that the number of nodes was not correct*/
  (*nb_nodes)=k-1;

  return correspondance;




}

/**
 * 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;
  
}



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


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

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// ALL PATHS between 2 nodes //////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void find_nodes(vector<int>path, int start, int end, list<NODE>**adjacency_list, string *correspondance){
  if(start==end){
    for( vector<int>::iterator it1=path.begin(); it1 != path.end(); it1++)
      cout << *it1 << " " ;
    return;
  }

  for (list<NODE>::iterator it=adjacency_list[start]->begin(); it!=adjacency_list[start]->end(); it++){
    path.push_back(it->number);
    find_nodes(path, it->number, end, adjacency_list, correspondance);
  
    path.pop_back();
    return;
  }

  return;
}



void allPaths(vector<int>path, int start, int end, list<NODE>**adjacency_list, string *correspondance){
  cout << start << end << endl;
  path.push_back(start);
  find_nodes (path,start,end, adjacency_list, correspondance) ;
}




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

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// 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;
    }
}

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






///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////GET 2 PATHS //////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
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 ////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


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){
    //    file1 << direction << endl;
    //   file1 << "seq1:"<< sequence1 << endl;
    //    file1 << "seq2:" << seq2<<endl;
    //  file1 <<"seq1_common:" << sequence1.substr(sequence1.length()-k+1)<< endl;
    //  file1 << "seq2_common:"<< seq2.substr(0,k-1)<< endl;
    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 /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


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;

}




///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// SWITCHING NODES /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


void Switching_Nodes(vector<int>cycles, list<NODE>**adjacency_list, string *correspondance, int cycles1){
  vector<int>::iterator iterator=cycles.begin();
  vector<int>::iterator p=cycles.begin()+1;
  //list<NODE>::iterator it;
  int vertex1=0;
  int vertex2=0;  
  char *direct1=new char[2];
  strcpy(direct1,"00");
  char *directfirst=new char[2];
  strcpy(directfirst,"00");
  char *directlast=new char[2];
  strcpy(directlast,"00");
  int valid;  
  int one=cycles.front();
  int total1=0;
  double mean1=0;
  int total2=0;
  double mean2=0;
  int total3=0;
  int mean3=0;
  int total4=0;
  int mean4=0;

    

  int lastnum=cycles[cycles.size()-1];
  //file1 << "last " <<correspondance[lastnum]<<endl;
  int firstnum=cycles[0];

  //file1 << correspondance[firstnum]<<endl;
  for (list<NODE>::iterator it=adjacency_list[lastnum]->begin(); it!=adjacency_list[lastnum]->end(); it++){
    if(firstnum==it->number){
      strcpy(direct1,it->direction);
      strcpy(directfirst,it->direction);

      //       file1<< "first direction :" << direct1 <<endl;
    }
  }
  
  for(int i=0;i+1 < cycles.size();i++){
    int num1=cycles[i];
    int num2=cycles[i+1];
    //cout <<"num1: " <<correspondance[num1]<< " " <<"num2: "<< correspondance[num2] << endl;
    for (list<NODE>::iterator it=adjacency_list[num1]->begin(); it!=adjacency_list[num1]->end(); it++){
 
      if(num2==it->number){
	//cout<< it->direction << endl;
   

	if( (valid=valid_path(direct1, it->direction))==0){
	  // file1 << direct1 << " "<< it->direction << endl;
	  //cout<< direct1 << " "<< it->direction << endl;
	  switching.push_back(num1);

	}
  
     
	strcpy(direct1,it->direction);

      }

    }

  }

  //if( (valid=valid_path(directlast, directfirst))==0){
  //    cout << directlast << " "<< directfirst << endl;
  // switching.push_back(1);
  //}

  file1 << "No of switching nodes:"<< switching.size() << endl;
  file1 << "Number of nodes in cycle:"<< cycles.size() <<endl;
  // cout << "The switching nodes are: ";
  // for( vector<int>::iterator it1=switching.begin(); it1 != switching.end(); it1++)
  //   cout <<correspondance[ *it1] << " " ;

  if(switching.size()==2 ){
    file1 << "The switching nodes are: ";
    for( vector<int>::iterator it1=switching.begin(); it1 != switching.end(); it1++)
      file1 <<correspondance[ *it1] << " " ;

    int start=switching[0];
    int end=switching[1];
    //   cout << "start:" <<start << " " << end<<endl; 
    //  cout << "start:" <<correspondance[start] << " " << correspondance[end]<<endl; 
    file1 << "\nPaths: " << endl;
   
    //  Print_Path_twonodes(correspondance, predecessor, start, end);
    Get_paths(&path1,&path2, cycles, start, end);

    int  length1=Get_Length(path1,sequences, correspondance, adjacency_list,expanding1);
    // cout << "Next path" << endl;
    int length2=Get_Length(path2,sequences, correspondance,adjacency_list,expanding2);
    if (length1>=length2){
      file1 <<"Upper path:"<<endl;
      file1 <<"Node\tCoverage_1\tCoverage_2"<< endl;
      for( vector<int>::iterator it1=path1.begin(); it1 != path1.end(); it1++){
    
	file1 <<correspondance[*it1] << "\t " ;
	int num=atoi(correspondance[*it1].c_str());
        file1<< coverage1[num] <<"\t";
        total3=total3+ coverage1[num];
	file1<< coverage2[num]<<"\t";
	total4=total4+coverage2[num];
	file1 << " " << endl;
      }
      file1 << " " << endl;
      if(length1-k_value-1==0){
	mean3=0;
	mean4=0;
      }
      else{
	mean3=total3/(length1-k_value-1);
	mean4=total4/(length1-k_value-1);
      }
      file1 << "Average coverage of upper path in condition 1:"<< mean3 <<endl;
      file1 << "Average coverage of upper path in condition 2:"<< mean4 <<endl;
      file1 << " " << endl;

      //cout <<" "<< endl;
      file1<<"Lower path: "<<endl;
      ;
      file1 <<"Node\tCoverage_1\tCoverage_2"<< endl;

      for( vector<int>::iterator it1=path2.begin(); it1 != path2.end(); it1++){
	file1 <<correspondance[*it1] << "\t " ;
	int num=atoi(correspondance[*it1].c_str());
        file1<< coverage1[num] <<"\t";
        total1=total1+coverage1[num];
	file1<< coverage2[num]<<"\t";
	total2=total2+coverage2[num];

	file1 << " " << endl;

      }
      file1 << " " << endl;

      if(length2-k_value-1==0){
	mean1=0;
	mean2=0;
      }
      else{
	mean1=total1/(length2-k_value-1);
	mean2=total2/(length2-k_value-1);
      }
      file1 << "Average coverage of lower path in condition 1:"<< mean1 <<endl;
      file1 << "Average coverage of lower path in condition 2:"<< mean2 <<endl;
      file1 << " " << endl;
      file1 <<"Upper Length:" << length1 << endl;
      file1 <<"Lower Length:" << length2 << endl;
      file1<< ">"<< "Cycle_"<< cycles1 <<"_upper_Length:"<< length1<<"_coverage1:"<<mean3<<"_coverage2:"<< mean4 << endl;
      file1 << expanding1 << endl;
      file1<< ">"<< "Cycle_"<< cycles1 <<"_lower_Length:"<< length2<<"_coverage1:"<<mean1<<"_coverage2:"<< mean2 << endl;
      file1 << expanding2 << endl;
      file1 <<" "<< endl;
      tabfile<<cycles1<<"\t"<<expanding1<<"\t"<<expanding2<<"\t"<<length1<<"\t"<<length2<<"\t"<<mean3<<"\t"<<mean4<<"\t"<<mean1<<"\t"<<mean2<<"\t"<<endl;


      myfile << ">"<< "Cycle_"<< cycles1 <<"_upper_Length:"<< length1<<"_coverage1:"<<mean3<<"_coverage2:"<< mean4 << endl;
      myfile << expanding1 << endl;
      myfile<< ">"<< "Cycle_"<< cycles1 <<"_lower_Length:"<< length2 <<"_coverage1:"<<mean1<<"_coverage2:"<< mean2 << endl;
      myfile << expanding2 << endl;
      //  myfile <<" "<< endl;

      if(length1==length2) snps++;

    }






    //else largest length is upper, reverse them
    else{
      file1 <<"Upper path:"<< endl;
      file1 <<"Node\tCoverage_1\tCoverage_2"<< endl;
      for( vector<int>::iterator it1=path2.begin(); it1 != path2.end(); it1++){
	file1 <<correspondance[*it1] << "\t " ;
	int num=atoi(correspondance[*it1].c_str());
        file1<< coverage1[num] <<"\t";
        total1=total1+coverage1[num];
	file1<< coverage2[num]<<"\t";
	total2=total2+coverage2[num];

	file1 << " " << endl;
      }
      file1 << " " << endl;

      if(length2-k_value-1==0){
	mean1=0;
	mean2=0;
      }
      else{
	mean1=total1/(length2-k_value-1);
	mean2=total2/(length2-k_value-1);
      }
      file1 << "Average coverage of upper path in condition 1:"<< mean1 <<endl;
      file1 << "Average coverage of upper path in condition 2:"<< mean2 <<endl;
      file1 << " " << endl;



      //cout <<" "<< endl;
      file1<<"Lower path: "<<endl;

      for( vector<int>::iterator it1=path1.begin(); it1 != path1.end(); it1++){
	file1 <<correspondance[*it1] << "\t " ;
	int num=atoi(correspondance[*it1].c_str());
        file1<< coverage1[num] <<"\t";
        total3=total3+ coverage1[num];
	file1<< coverage2[num]<<"\t";
	total4=total4+coverage2[num];
	file1 << " " << endl;
      }

      file1 << " " << endl;
      if(length1-k_value-1==0){
	mean3=0;
	mean4=0;
      }
      else{
	mean3=total3/(length1-k_value-1);
	mean4=total4/(length1-k_value-1);
      }
      file1 << "Average coverage of lower path in condition 1:"<< mean3 <<endl;
      file1 << "Average coverage of lower path in condition 2:"<< mean4<<endl;
      file1 << " " << endl;


      file1 <<"Upper Length:" << length2 << endl;
      file1 <<"Lower Length:" << length1 << endl;
      file1<< ">"<< "Cycle_"<< cycles1 <<"_upper_Length:"<< length2<<"_coverage1:"<<mean1<<"_coverage2:"<< mean2 << endl;
      file1 << expanding2 << endl;
      file1<< ">"<< "Cycle_"<< cycles1 <<"_lower_Length:"<< length1<<"_coverage1:"<<mean3<<"_coverage2:"<< mean4 << endl;
      file1 << expanding1 << endl;
      file1 <<" "<< endl;
      tabfile<<cycles1<<"\t"<<expanding2<<"\t"<<expanding1<<"\t"<<length2<<"\t"<<length1<<"\t"<<mean1<<"\t"<<mean2<<"\t"<<mean3<<"\t"<<mean4<<"\t"<<endl;


      myfile << ">"<< "Cycle_"<< cycles1 <<"_upper_Length:"<< length2<<"_coverage1:"<<mean1<<"_coverage2:"<< mean2 << endl;
      myfile << expanding2 << endl;
      myfile<< ">"<< "Cycle_"<< cycles1 <<"_lower_Length:"<< length1<<"_coverage1:"<<mean3<<"_coverage2:"<< mean4 << endl;
      myfile << expanding1 << endl;

    }

    path1.clear();
    path2.clear();
    count_events++; 
  }

  switching.clear();
  // Print_path(vertex1, vertex2, adjacency_list,correspondance);
}
    
  





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

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////// PRINT PATH //////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////

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);
  }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// DFS_VISIT /////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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)  ){
      //list<int>*cycleslist;
      //cycleslist->push_front(it->number);
      Print_Path(correspondance, predecessor,it->number,i);
      // Switching_Nodes(cycleslist, adjacency_list);
      //cout<<endl;
      //cout <<  correspondance[it->number]<<" "<< correspondance[i];
      file1<<"\n"<<endl;

      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);
      //reverse(cycles.begin(),cycles.end());
      // Switching_Nodes(cycles, adjacency_list, correspondance,cycles1);
      cycles1++;
      cycles.clear();
      // cycle_reverse.clear();

      //cout << cycles<< " cycle"<< endl;
    }

  }

  color[i]="black";
  finished[i]=time++;
  //cout<< "finished:"<< correspondance[i] << endl;
  //cout << cycles<<endl;
}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    DFS  ///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
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];
//}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    ENUMERATE ALL CYCLES  ///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
void printArray(bool *array, int size)
{
  for (int i=0;i<size;i++)
    cout<<array[i]<<" ";
  cout<<endl;
}


void printCircuit(stack<int> stack1, string *correspondance) {
  
  stack<int> stack2;  

  while (!stack1.empty())
    {
      int u=stack1.top();
      stack1.pop();
      // cout<<u<<" ";
      stack2.push(u);
      cyclesnew.push_back(u);
      //cout << correspondance[u] << " ";
    }

  if(cyclesnew.size()!=2){
    numberOfCycles++;
    for(vector<int>::iterator i=cyclesnew.end()-1;i!=cyclesnew.begin()-1;i--){
      //cout << correspondance[*i]<< " ";
      //sum=sum+correspondance[*i];

      cycles.push_back(*i);
    }
    //cout << " " << endl;
  
    file1 << "Cycle:" << numberOfCycles <<endl ;
    for(vector<int>::iterator i=cyclesnew.begin();i!=cyclesnew.end();i++){
      file1 << correspondance[*i]<< " ";
    }
    file1 << endl;
    Switching_Nodes(cycles, adjacency_list, correspondance,numberOfCycles);
    file1<<"\n"<< endl;
  }

  cyclesnew.clear();
  cycles.clear();
  while (!stack2.empty())  
    {
      int u=stack2.top();
      stack2.pop();
      stack1.push(u);
    }
  //cout<<endl;

}

bool backtrack(int s, int v, bool *mark,stack<int>& marked, stack<int>& point, int nb_nodes, string *correspondance) {

  bool f = false;
  point.push(v);
  mark[v] = true;
  marked.push(v);
  int w;
  list<NODE>::iterator it1;

  /*This section could be optimised using adjacency lists instead of adjacency matrix, as in the original Tarjan publication*/

  for (list<NODE>::iterator it=adjacency_list_copy[v]->begin(); it!=adjacency_list_copy[v]->end(); it++)
    {
      w=it->number;
      
      if(w < s)
	{
	  //erase v from the adjacency list of w
	  for (list<NODE>::iterator it2=adjacency_list_copy[w]->begin(); it2!=adjacency_list_copy[w]->end(); it2++)
	    {
	      if (v==it2->number)
		it2=adjacency_list_copy[w]->erase(it2);
	    }
	  //erase w from the adjacency list of v
	  it=adjacency_list_copy[v]->erase(it);
	}
      else if (w == s) 
	{
	  printCircuit(point, correspondance);
	  f = true;
	} 
      else if (!mark[w]) {
	bool g = backtrack(s, w,mark, marked, point, nb_nodes,correspondance);
	f = f || g;
      }
    }
 

  //for (int w = 0; w < nb_nodes; w++) {

  //  if (adjacency[v][w]) {
  //   if (w < s) {
  //	adjacency[v][w] = false;
  //	adjacency[w][v]=false;
  //   } else if (w == s) {
  //	printCircuit(point, correspondance);


  //	f = true;
  //   } else if (!mark[w]) {
  //	bool g = backtrack(s, w,mark, marked, point, nb_nodes,correspondance, adjacency_list);
  //	f = f || g;
  //    }
  //  }
  //  }

  
  if (f) 
    {
      while (marked.top() != v) 
	{
	  int u = marked.top();
	  marked.pop();
	  mark[u] = false;
	}
      marked.pop();
      mark[v] = false;
		    
    }

  point.pop();
  return f;
}

void enumerateAllCircuits(int nb_nodes, string *correspondance) {
  cout << "Cycles found:\n"<<endl;
  mark = new bool[nb_nodes];
  for (int i = 0; i < nb_nodes; i++) {
    mark[i] = false;
  }
  for (int s = 0; s < nb_nodes; s++) {
    backtrack(s, s, mark, marked, point,nb_nodes,correspondance);
    while (!marked.empty()) {
      int u = marked.top();
      marked.pop();
      mark[u] = false;
    }
  }


}

void printAdjacencyMatrix(int nb_nodes) {
  for (int i=0;i<nb_nodes;i++)
    {
      for (int j=0;j<nb_nodes;j++)
	cout<<adjacency[i][j]<<" ";
      cout<<endl;
    }
}

void printAdjacencyList(int nb_nodes,list<NODE> **adjacency_list, string * correspondance) {
  for (int i=0;i<nb_nodes;i++)
    {
      for (list<NODE>::iterator it=adjacency_list[i]->begin(); it!=adjacency_list[i]->end(); it++){
	cout<<it->number<<" ";
	//cout<< correspondance[it->number] << " ";
      }
      cout<<endl;
    }
}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    GET 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);
      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);
      //cout << sequence <<" "<<endl;
      // cout << "NEW"<<endl;

      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);

      // cout << cov2 <<endl;

      // if((labels.find(str1)) !=  (labels.end()))
      //	{
      //	  int  n1=it->second;
      //  cout << n1 << " n1" << endl;
      // cout << it->first << " "<<endl; ;
      //sequences[it]=new string();
      
      // coverage1[n1]=new int[2];
      // coverage1[n1]=new int;
      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));
      //mymap.insert ( pair<char,int>('a',100) );

      //  cout << n1 << " " << sequences[n1]<<" " << endl;
      // strcpy(sequences[n1],sequence);
      // for(map<string,int>::iterator iter = labels.begin();iter!=labels.end();iter++)
      //correspondance[(*iter).second] = (*iter).first;
    
    }

  file_op.close();

 

}


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

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

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// MAIN FUNCTION /////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char ** argv)
{
 

  cout<< "Finding cycles"<<endl;

  //list<NODE> **sarah;

  /*allocation of necessary structures*/

  adjacency_list = new list<NODE> *[MAX_NB_NODES];
  adjacency_list_copy = new list<NODE> *[MAX_NB_NODES];
  for (int i=0;i<MAX_NB_NODES;i++){
    adjacency_list[i]=new list<NODE>;
    adjacency_list_copy[i]=new list<NODE>;
  }



  if (argc != 7) { // Check the value of argc. If not enough parameters have been passed, inform user and exit.
    std::cout << "Usage is ./paths graph_edges_file  graph_nodes_file  k_value outputfile_1.fa  outputfile_cycles.txt output_tab_file.txt \n";
    std::cin.get();
    exit(0);
  }
  else{

    /*read graph*/
    int nb_nodes=0;
    string *direction;
    string *correspondance = read_graph(adjacency_list,&nb_nodes,argv[1]); //reads graph file and creates adjaceny list
    // string *correspondance1 = read_graph(adjacency_list_copy,&nb_nodes,argv[1]);
    //adjacency_list_copy->merge(adjacency_list);
    cout<<"the number of nodes is "<<nb_nodes<<endl;
    string value=argv[3];  // k value
    char *filefasta=argv[4];  //outpu file1 -fasta formatted results
    char *fileresults=argv[5]; //output file 2 -txt output
    char *filetabs=argv[6];
    myfile.open (filefasta);
    file1.open(fileresults);
    tabfile.open(filetabs);
    k_value=atoi(value.c_str());


    tabfile<< "Cycle\tUpperpath\tLowerPath\tLength_Upper\tLength_Lower\tUpper_Coverage_1\tUpper_Coverage_2\tLow_Coverage_1\tLow_Coverage_2\tPrediction"<<endl;

    /*vector indicating who is in what class*/
    int *E = new int[nb_nodes];
    int *newE= new int[nb_nodes];
    int *cyc=new int[nb_nodes];
  
    //get the string of each node and the coverages and save it to sequences map<int,string>
    Get_sequences(correspondance, argv[2],sequences, coverage1, coverage2);
  
  
    cout << "Finding events"<< endl;
    //for
    //file1 << sequences[12] << endl;
    file1<<"============================================================================"<<endl;
    file1<< "Summary of resutls can be found at the end of this file." << endl;
    file1<<"============================================================================"<<endl;
    file1<< "The cycles in the graph are: " << endl;
    file1<<"============================================================================"<<endl;
  
    /*New method to find cycles*/
    /*First store the graph as an adjacency matrix*/
    //cout << "create adjacency matrix "<<endl;
    //  adjacency = new bool*[nb_nodes];
    // for(int i = 0; i < nb_nodes; ++i)
    //   {
    //    adjacency[i] = new bool[nb_nodes];
    //    for (int j=0;j<nb_nodes;j++)
    //	adjacency[i][j]=0;
    //  }
    // adjacency=new bool[nb_nodes][nb_nodes];
  
    // printAdjacencyMatrix(nb_nodes);
    //for (int i=0;i<nb_nodes;i++)
    //  for (list<NODE>::iterator it=adjacency_list[i]->begin(); it!=adjacency_list[i]->end(); it++)
    //    adjacency[i][it->number]=1;
  
    /*Then query it to get all circuits*/
    //printAdjacencyMatrix(nb_nodes);

    // for (int i=0;i<MAX_NB_NODES;i++){
    //copy (adjacency_list[i]->begin(), adjacency_list[i]->end(),     // source
    //         adjacency_list_copy[i]->begin());
    //}
    adjacency_list_copy=adjacency_list;
    cout<<endl;
    cout << "entering enumerateAllcircuits" << endl;
    //printAdjacencyList(nb_nodes,adjacency_list,correspondance);
    enumerateAllCircuits(nb_nodes,correspondance);
    cout<<"Number of cycles is "<< numberOfCycles <<endl;
    //exit(0);
    cyclesnew.clear();

    // Find cycles with depth first search
    //DFS(adjacency_list,nb_nodes,correspondance);
 
    file1<<"\n\n============================================================================"<<endl;
    file1<<"Summary of results:"<<endl;
    file1 << "No of cycles:" << numberOfCycles <<endl;
    file1 << "No of cycles investigated for events(having exactly two switching nodes): " << count_events << endl;
    file1 << "No of cycles having equal lower and upper paths, identified as SNPs: " << snps <<endl;
    file1 << "No of cycles having unequal upper and lower paths: " << count_events-snps << endl;
    file1<<"============================================================================"<<endl;

    cout<<"============================================================================"<<endl;
    cout<< "Summary of resutls" << endl;
    cout <<"============================================================================"<<endl;

    cout << "No of cycles:" << numberOfCycles <<endl;
    cout << "No of cycles investigated for events(having exactly two switching nodes): " << count_events << endl;
    cout << "No of cycles having equal lower and upper paths, identified as SNPs: " << snps <<endl;
    cout << "No of cycles having unequal upper and lower paths: " << count_events-snps << endl;
    cout<<"============================================================================"<<endl;

    myfile.close();
    file1.close();
    tabfile.close();


    // cout<<"The number of connected components is "<<nb_connected_components<<endl;





    //for (int j=0;j<nb_connected_components;j++)
    //  {// cout<<j<<"\t";
    //   for (int i=0;i<nb_nodes;i++)
    //	if (E[i]==j && newE[i]==j){
    //	  cout<<correspondance[i]<<"\t";
    //}
    //     cout<<endl;
    //   }
    // cout<<endl;

  }


}