#include <iostream>
#include <sstream>
#include <fstream>
#include <list>
#include <map>
#include <queue>
#include <stack>
#include <cstring>
#include <cstdlib>
#include <string>
#include <algorithm>
#include "functions.h"   // contains all the rest of the methods of the program

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// SWITCHING NODES /////////////////////////////////////////////////////////////////////////////////////////////////
/* this is where everything happens
   cycles are checked if they have two switching nodes
   xml file is created with switching nodes and external switching nodes
   events file is creeated with all cycles and their sequences */
/////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

// GUS: There's still a bug in this function.

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;
  set<int> set1;
    

  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);
      direct1=it->direction;
      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 it1=adjacency_list[num1]->begin(); it1!=adjacency_list[num1]->end(); it1++){
 
      if(num2==it1->number){
	//cout<< it->direction << endl;
   

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

	}
  
     
	direct1=it1->direction;

      }

    }

  }


  if(switching.size()==2 ){
    numberOfCycles++;

    xml_file << "<Cycle " << numberOfCycles << ">"<<endl;


    file1 << "Cycle "<< numberOfCycles << endl;
    for( vector<int>::iterator it=cycles.begin(); it != cycles.end(); it++){
      file1 << correspondance[*it] << " ";
    }
    file1<< " "<<endl;
    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];
    xml_file << "<SN1>"<< correspondance[start] << "</SN1>"<<endl;
    xml_file << "<SN2>"<< correspondance[end] << "</SN2>"<<endl;

    //Find External Switching Nodes
    //Intersection of nodes of current cycle with the set of all swithcing nodes "switchnodes"
    vector<int>cycles_sorted;
    cycles_sorted=cycles;
    sort( cycles_sorted.begin(), cycles_sorted.end());
    vector<int>switch_sorted;
    switch_sorted=switchnodes;
    sort( switch_sorted.begin(), switch_sorted.end());
    set_intersection(cycles_sorted.begin(), cycles_sorted.end(), switch_sorted.begin(), switch_sorted.end(), inserter(set1, set1.end()));

    int i=0;
    for( set<int>::iterator it3=set1.begin(); it3 != set1.end() ;it3++){
      if (start != *it3 && end !=*it3){
	xml_file<< "<ESN>"<< correspondance[*it3]<<"</ESN>"<<endl;
	i++;
      }
    }
    xml_file <<"<ESN_number>"<< i <<"</ESN_number>"<<endl;
    //xml_file << endl;


    xml_file << "<all_nodes>";
    for( vector<int>::iterator it=cycles.begin(); it != cycles.end(); it++){
      xml_file << correspondance[*it] << " ";
    }
    xml_file<< "</all_nodes>"<<endl;
    //   cout << "start:" <<start << " " << end<<endl; 
    //  cout << "start:" <<correspondance[start] << " " << correspondance[end]<<endl; 
    file1 << "\nPaths: " << endl;
   



    //  Print_Path_twonodes(correspondance, predecessor, start, end);
    //get the upper and lower paths based on the switching nodes'position in path	
    Get_paths(&path1,&path2, cycles, start, end);



    //////////////////////////////////////Calculating length and coverage values for each cycle////////////////////////////////////////////





    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;
      xml_file<< "<upper>"<< expanding1 << "</upper>" <<endl; 
      file1<< ">"<< "Cycle_"<< cycles1 <<"_lower_Length:"<< length2<<"_coverage1:"<<mean1<<"_coverage2:"<< mean2 << endl;
      file1 << expanding2 << endl;
      xml_file<< "<lower>"<< expanding2 << "</lower>" <<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 lower, 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;
      xml_file<< "<upper>"<< expanding2 << "</upper>" <<endl; 

      file1<< ">"<< "Cycle_"<< cycles1 <<"_lower_Length:"<< length1<<"_coverage1:"<<mean3<<"_coverage2:"<< mean4 << endl;
      file1 << expanding1 << endl;
      xml_file<< "<lower>"<< expanding1 << "</lower>" <<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;



    }
    xml_file << "</Cycle " << numberOfCycles << ">"<<endl;
    xml_file << " "<<endl;
    path1.clear();
    path2.clear();
    count_events++; 
  }

  switching.clear();

}
    
  



//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    SWITCH DIRECTIONS  ///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int switch_directions(char *dir1, char *dir2)
{
  if(dir1[0]== '0') return 0;

  if(dir1[1]!=dir2[0])
    return 1;

  return 0;

}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    ENUMERATE ALL CYCLES  ///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



////////////////////////////////////////////printArray - not used //////////////////////////////////////////////////////////////////
void printArray(bool *array, int size)
{
  for (int i=0;i<size;i++)
    cout<<array[i]<<" ";
  cout<<endl;
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    PRINT CIRCUIT ///////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

int tot2 = 0;

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

    }

  if(cyclesnew.size()!=2){

    cout << "\n" << "Larger than 2: " << ++tot2 << "\n";
    for (int i = 0; i < (int)cyclesnew.size(); i++)
      cout << correspondance[cyclesnew[i]] << " ";
    cout << "\n";

    //reverse cycle
    for(vector<int>::iterator i=cyclesnew.end()-1;i!=cyclesnew.begin()-1;i--){

      cycles.push_back(*i);
    }

    Switching_Nodes(cycles, adjacency_list, correspondance,numberOfCycles);
    file1<<"\n"<< endl;
  }
  //}
  cyclesnew.clear();
  cycles.clear();

  //ELEMENTS BACK TO STACK POINT AFTER CYCLE IS PRINTED
  while (!stack2.empty())  
    {
      int u=stack2.top();
      stack2.pop();
      stack1.push(u);
    }
  //cout<<endl;

}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////   backtrack  ///////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

int tot = 0;

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

  //  strcpy(dir2,"00");


  //cout << " num at begin : " << num  <<endl;
  //bool g;
  //bool increased=false;
  cout << "bactrack s " << correspondance[s]<< " v is "<< correspondance[v]<< endl;
  bool f = false;
  bool removed = false;
  point.push(v);
  mark[v] = true;
  marked.push(v);
  // num=stack_SN(point,correspondance);
  int w;

  for (list<NODE>::iterator it=adjacency_list_copy[v]->begin(); it!=adjacency_list_copy[v]->end(); removed ? it : it++)
  { 
    removed = false;

    w = it->number;
    //string dir1 (it->second);
    dir2 = it->direction;
    
    // GUS: I'm removing the switching node count for the moment, till
    //we are sure that the standard algorithm is correct. Will add it
    //again later.
    //int value=switch_directions(dir1,dir2);
    
    //if (value==1) {increased=true;}
    
    //cout << value << num << dir1 << dir2 << endl;
    
    if(w < s)
    {
      //erase w from the adjacency list of v
      cout << "w is smaller than s: " <<  correspondance[s] << " w " << correspondance[w] << endl;
      
      // GUS: We should not increment the iterator in the outer
      // for-loop if an element was removed!  When you use it =
      // list.erase(it) the iterator already points to the next
      // element in the list. By doing this we may wrongly skip some
      // elements of adjacency_list[v].
      removed = true;
      
      it = adjacency_list_copy[v]->erase(it);
      
      // GUS: Tarjan's algorithm works for directed graphs, we cannot
      // remove the other edge, because for the algorithm they are
      // different.

      //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);
      //   }
    }
    
    //else if (w == s && (num+value) !=3) {
    else if (w == s)
    {
      //int SN=stack_SN(marked,correspondance);
      cout << "Aqui!!" << tot++ << "\n";
      printCircuit(point, correspondance);
      f = true;
    }
    else if (!mark[w]) 
    {
      //if((num+value) != 3){
      bool g = backtrack(s, w, mark, marked, point, nb_nodes, correspondance, num+value, dir2);
      f = f || g;  
      //}
    }
    //num=num-1; 
  }
  
  if (f) 
  {
    while (marked.top() != v) 
    {
      int u = marked.top();
      //cout <<" pop  marked : " << u << " " << correspondance[u] << endl;
      marked.pop();
      mark[u] = false;
    }
    marked.pop();
    //cout <<" pop  marked : " << k << " " << correspondance[k] << endl;
    mark[v] = false;
  }
  //cout << "pop point :  " << point.top() << " " << correspondance[point.top()] <<endl;
  point.pop();

  return f;
}



///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    enumerate all circuits
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void enumerateAllCircuits(int nb_nodes, string *correspondance,vector<int>switchnodes) {
  cout << "Enumerate, nb_nodes=" << nb_nodes << "\n";

  mark = new bool[nb_nodes];
  bool f=true;
  for (int i = 0; i < nb_nodes; i++) {
    mark[i] = false;
  }  
  int num=0; 
  char *dir1=new char[3];
  
  dir1[0]='0';
  dir1[1]='0';
  dir1[2]='\n';
  //for( vector<int>::iterator it=switchnodes.begin(); it != switchnodes.end(); it++){
  //mark[*it] = false;
  //}

  // GUS: We need to actually permute the order to have the SN at the
  // beginning of the adjacency list, it's not correct to just call
  // them first. Remember Tarjan's needs a total order, the problem
  // would be in the comparion "w < s". I'm switching back to the
  // original version for the moment.

  //  for( vector<int>::iterator it=switchnodes.begin(); it != switchnodes.end(); it++){
  // cout << "Enumerate, SN:" << correspondance[*it] << "\n";
  
  //  num=0;
    for (int s = 0; s < nb_nodes; s++) {

    backtrack(s, s, mark, marked, point,nb_nodes,correspondance,num,dir1);
    //backtrack(*it,*it , mark, marked, point,nb_nodes,correspondance,num, dir1);

    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++)
    {
      cout << correspondance[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;
    }
}


///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////    check if a node is a switching node     ////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



int check_if_switching_node(map<int, string> ingoing,map <int,string>outgoing, string *correspondance){


  for( map<int,string>::iterator it1=ingoing.begin(); it1 != ingoing.end(); it1++){
    for( map<int,string>::iterator it2=outgoing.begin(); it2 != outgoing.end(); it2++){
      if(it1->first!=it2->first){
	//incoming FF outgoing R type FR
	string dir1 (it1->second);
	string dir2 (it2->second);

	if(strcmp(dir1.c_str(), "FF")==0){
	  if(strcmp(dir2.c_str(), "RR")==0 || strcmp(dir2.c_str(), "RF")==0){
	    //cout << "switching node found!!! "<<"incoming : " << it1->first << " " << dir1 << " outgoing " << it2->first << " " << dir2 << endl; 
	    return 2;
	  }
	}

	//incoming RR outgoing F, type RF
	if(strcmp(dir1.c_str(), "RR")==0){
	  if(strcmp(dir2.c_str(), "FF")==0 || strcmp(dir2.c_str(), "FR")==0){
	    //cout << "switching node found!!! "<<"incoming : " << it1->first << " " << dir1 << " outgoing " << it2->first << " " << dir2 << endl; 
	    return 1;
	  }
	}

	//incoming FR outgoing F, type RF
	if(strcmp(dir1.c_str(), "FR")==0){
	  if(strcmp(dir2.c_str(), "FF")==0 || strcmp(dir2.c_str(), "FR")==0){
	    //cout << "switching node found!!! "<<"incoming : " << it1->first << " " << dir1 << " outgoing " << it2->first << " " << dir2 << endl; 
	    return 1;
	  }
	}

	//incoming RF outgoing R tpye FR


	if(strcmp(dir1.c_str(), "RF")==0){
	  if(strcmp(dir2.c_str(), "RF")==0 || strcmp(dir2.c_str(), "RR")==0){
	    //cout << "switching node found!!! "<<"incoming : " << it1->first << " " << dir1 << " outgoing " << it2->first << " " << dir2 << endl; 
	    return 2;
	  }
	}

      }

    }
  }

  return 0;
}






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

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

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char ** argv)
{
 
  cout<< "Finding cycles"<<endl;



  /*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 != 8) { // 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 xml_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];
    char *file_xml=argv[7];
    myfile.open (filefasta);
    file1.open(fileresults);
    tabfile.open(filetabs);
    xml_file.open(file_xml);
    k_value=atoi(value.c_str());
    char *dir1=new char[2];
    char *dir2=new char[2];


    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);
  
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////Find switching nodes in the graph
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    for(int i=0; i< nb_nodes; i++){
      //cout << correspondance[i] << endl;

      for (list<NODE>::iterator it=adjacency_list[i]->begin(); it!=adjacency_list[i]->end(); it++){
	//strcpy(dir1, it->direction);
	dir1=it->direction;

	int n1= atoi(correspondance[it->number].c_str());

	ingoing[n1]=new char[2];
	outgoing[n1]=new char[2];
        

	outgoing[n1]=it->direction;



	//set inverse direction for each node in the graph

	if(strcmp(dir1, "RR")==0)
	  strcpy(dir2,"FF");
	else if(strcmp(dir1, "FF")==0)
	  strcpy(dir2,"RR");
	else if(strcmp(dir1, "FR")==0)
	  strcpy(dir2,"FR");
	else if(strcmp(dir1, "RF")==0)
	  strcpy(dir2,"RF");

	//cout << "direction 2 "<< dir2 << endl;
	ingoing[n1]=dir2;

	outgoing.insert(pair<int, string>(n1,dir1)); //insert outgoing edge to map outgoing
	ingoing.insert(pair<int, string>(n1,dir2)); //insert ingoing edge to map ingoing
	//ingoing.insert(pair<int, string>(correspondance[it->number], it->direction)); 

      }

      int graphnode=atoi(correspondance[i].c_str());

      //check if current node is switching by comparing its ingoing edges against all its outgoing edges for a switch in direction
      int sn=check_if_switching_node(ingoing,outgoing,correspondance);


      if(sn==1)
	switch_RF.push_back(i);   //reverse-forward switching nodes
      if(sn==2)
	switch_FR.push_back(i);  // forward-reverse switching nodes
      if (sn==1 || sn==2){
	//all switching nodes
	switchnodes.push_back(i);  // ALL switching nodes are put in this vector
	//cout << "SN: " << correspondance[i] << "\n";
      }
      ingoing.clear();
      outgoing.clear();
    }
    int numberswitch=switchnodes.size();

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


    cout << "Number of total nodes: "<< nb_nodes<<endl;
    cout << "Number of switching nodes: "<<numberswitch<<endl;

    cout<<"============================================================================"<<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;
  




    //copy adjacency_list because we delete nodes in Tarjan's algorithm and we want to keep the original for finding the paths later on
    adjacency_list_copy=adjacency_list;
    cout<<endl;
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////Run Tarjan's algorithm
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    cout << "Running Tarjan's algorithm" << endl;
    xml_file << "<Cycles> "<<endl;
    printAdjacencyList(nb_nodes,adjacency_list,correspondance);
    enumerateAllCircuits(nb_nodes,correspondance,switchnodes);
    cout<<"Number of cycles is "<< numberOfCycles <<endl;

    cyclesnew.clear();

    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////Output
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    file1<<"\n\n============================================================================"<<endl;
    file1<<"Summary of results:"<<endl;
    //file1 << "No of cycles:" << numberOfCycles <<endl;
    file1 << "No of cycles with 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 with 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;
    xml_file << "</Cycles>"<<endl;
    myfile.close();
    file1.close();
    tabfile.close();


  }


}