/* ***************************************************************************
 *
 *                              KisSplice
 *      de-novo calling alternative splicing events from RNA-seq data.
 *
 * ***************************************************************************
 *
 * Copyright INRIA 
 *  contributors :  Vincent Lacroix
 *                  Pierre Peterlongo
 *                  Gustavo Sacomoto
 *                  Vincent Miele
 *                  Alice Julien-Laferriere
 *                  David Parsons
 *
 * pierre.peterlongo@inria.fr
 * vincent.lacroix@univ-lyon1.fr
 *
 * This software is a computer program whose purpose is to detect alternative
 * splicing events from RNA-seq data.
 *
 * This software is governed by the CeCILL license under French law and
 * abiding by the rules of distribution of free software. You can  use,
 * modify and/ or redistribute the software under the terms of the CeCILL
 * license as circulated by CEA, CNRS and INRIA at the following URL
 * "http://www.cecill.info".

 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.

 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the
 * same conditions as regards security.
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL license and that you accept its terms.
 */
 
 
 
 
 
// ===========================================================================
//                               Include Libraries
// ===========================================================================
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <map>
#include <list>
#include <utility>
#include <cassert>




// ===========================================================================
//                             Include Project Files
// ===========================================================================
#include "NGraph.h"
#include "debug.h"
#include "CycleCompression.h"
#include "Utils.h"




// ===========================================================================
//                             Declare Used Namespaces
// ===========================================================================
using namespace std;




// ===========================================================================
//                                   Constant
// ===========================================================================
#define MAX 1024



//############################################################################
//                                                                           #
//                                Class NGraph                               #
//                                                                           #
//############################################################################

// ===========================================================================
//                                  Constructors
// ===========================================================================
NGraph::NGraph( int kValue, int outputtedSnps )
{
  _kValue = kValue;
  _nbOutput = outputtedSnps;

}

/*!
  TODO : check : Create a NGraph for the BCC passed in edges ? all_edges ?
*/
NGraph::NGraph( CGraph& cgraph, vector<char *>& seqs,
    vector<LabelledCEdge>& all_edges, vector<CEdge>& edges )
{
  _nbOutput = 0;
  for ( int i = 0 ; i < (int)edges.size() ; i++ )
  {
    int u = edges[i].getFirst();
    int v = edges[i].getSecond();

    
    //~ printf( "contructing %s %s (0x%x, 0x%x)\n", seqs[u], seqs[v], seqs[u], seqs[v] );
    //~ getchar();
    insert_node(std::to_string(u), seqs[u]);
    insert_node(std::to_string(v), seqs[v]);

    insert_bidirected_edges( all_edges, edges[i] );
    edges[i].swap_ends();
    insert_bidirected_edges( all_edges, edges[i] );
  }
  _kValue = cgraph.k_value;
  expand_parallel_edges();
}


// ===========================================================================
//                                 Public Methods
// ===========================================================================
void NGraph::insert_empty_node( string u )
{
  _strToNode[u] = _nodeToStr.size(); // _nodeToStr.size() corresponds to the current number of nodes in the graph
  _nodeToStr.push_back(u);
  
  _nodes.push_back( NNode() );
}

// Insert the new edge u->v (label).
// If the edge u->v is already present with a different label, concatenate <label> to the existing one.
// If its the same label do nothing.
void NGraph::insert_edge( string u, string v, string label )
{
  // If the node u doesn't exist, create it
  if ( _strToNode.find(u) == _strToNode.end() )
  {
    insert_empty_node(u);
  }
  
  // If the node v doesn't exist, create it
  if ( _strToNode.find(v) == _strToNode.end() )
  {
    insert_empty_node(v);
  }

  
  int u_idx = _strToNode[u];
  int v_idx = _strToNode[v];
  
  list<NEdge>::iterator it;
  for ( it = _nodes[u_idx].getMAdjList().begin() ; it != _nodes[u_idx].getMAdjList().end() ; it++ )
  {
    if ( it->get_node() == v_idx )
    {
      int offset;
      
      for ( offset = 0 ; offset+2 <= (int)it->get_labels().size() ; offset += 2 )
      {
        if ( it->get_labels().substr(offset, 2) == label )
        {
          break;
        }
      }
      
      if ( offset == (int)it->get_labels().size() )
      {
        it->add_label( label );
      }
      
      break;
    }
  }
  
  // The node 'v' is not in the adjList of _nodes[u_idx]
  if ( it == getAdjList(u_idx).end() )
  {
    NEdge new_edge( v_idx, label );
    _nodes[u_idx].add_to_adj_list( new_edge );
  }
}

void NGraph::insert_node( string v, string seq )
{
  // If the node doesn't exist yet, create it and insert it
  if ( _strToNode.find(v) == _strToNode.end() )
  {
    insert_empty_node( v );
  }
  
  _nodes[_strToNode[v]].setSequence( seq );
}



void NGraph::insert_bidirected_edges( vector<LabelledCEdge>& all_edges, CEdge edge)
{
  char small[3];
  char large[3];
  
  small[0] = (char)0;
  small[1] = (char)0;
  small[2] = '\0';  // Anything strictly smaller than (int)'F'
  
  large[0] = (char)127;
  large[1] = (char)0;
  large[2] = '\0';  // Anything strictly greater than (int)'R'

  
  vector<LabelledCEdge>::iterator low, upper, it;
  low   = lower_bound( all_edges.begin(), all_edges.end(), LabelledCEdge( edge, small ) );
  upper = upper_bound( all_edges.begin(), all_edges.end(), LabelledCEdge( edge, large ) );

  for ( it = low ; it != upper ; it++ )
  {
    int u = it->getFirst();
    int v = it->getSecond();
    insert_edge(std::to_string(u), std::to_string(v), it->label);
  }
}




/*!
\brief Modifies the structure of the adjacency list to reflect the parallel
 edges of the graph.
*/
void NGraph::expand_parallel_edges( void )
{
  // For each node in the graph...
  for ( int v = 0 ; v < (int)_nodes.size() ; v++ )
  {
    list<NEdge> to_add;
    list<NEdge>::iterator it;
    // For each neighbour of the current node
    for ( it = _nodes[v].getMAdjList().begin() ; it != _nodes[v].getMAdjList().end() ; it++ )
    {
      for ( int offset = 2 ; offset+2 <= (int)it->get_labels().size() ; offset += 2 )
      {
        to_add.push_back( NEdge( it->get_node(), it->get_labels().substr(offset, 2) ) );
      }
      
      // Keep just the first label
      if ( (int)it->get_labels().size() > 2 )
      {
        it->del_all_but_first();
      }
    }

    // Add the parallel edges
    if ( to_add.size() != 0 )
    {
      _nodes[v].getMAdjList().splice( _nodes[v].getMAdjList().begin(), to_add );
    }
  }
}
/*!
 * \brief Searchs and compress bubbles into NGraph
 * \param *compressed_bubbles number of compressed bubbles
 * \param snp_log_file
 * \param bccid: index of the bcc being searched and compressed.
 * \param output_context
 */
void NGraph::compress_all_bubbles( int *compressed_bubbles, FILE * snp_log_file, const int bccid, const bool output_context )
{
  bool *removed = new bool[2 * (int)_nodes.size() + MAX];
  for ( int i = 0; i < (2*(int)_nodes.size() + MAX); i++ )
  {
    removed[i] = false;
  }

  // Some vertices may be added to graph 'g', but we need to traverse
  // just the original ones.
  int initial_size = (int)_nodes.size();

  // Tries to close a bubble from every node with out_degree = 2
  for (int u = initial_size - 1; u >= 0 ; u--)
  {
    if ( not removed[u] )
    {
      DEBUG(cerr << "compress_all_bubbles, u = " << u << "\n";)

      // Looks for a bubble with opening node u
      find_and_compress_bubble(  u, 'F', removed, snp_log_file, bccid, output_context);
      
      // Looks for a bubble with opening node R(u)
      find_and_compress_bubble(  u, 'R', removed, snp_log_file, bccid, output_context);
    }
  }

  (*compressed_bubbles) = _nodes.size() - initial_size;
  *this = remove_marked_nodes( *this, removed, false );
  
  delete [] removed;
}



// Removes the arcs of (open -> u), (open -> l); adds the arc (open ->
// new_idx) and (R(new_idx) -> R(open))
/*!
 * \brief Fix_bubble_neighborhood
 * \param &new_node
 * \param new_idx
 * \param new_dir
 * \param open
 * \param u
 * \param l
 * \param *removed
 */
void NGraph::fix_bubble_neighborhood( NNode& new_node, int new_idx,
    char new_dir, idx_dir open, idx_dir u, idx_dir l, bool *removed )
{
  int   v     = open.first;
  char  v_dir = open.second;

  // Find the arcs correspond to open -> u and open -> l, with the
  // direction 'open.second'.
  list<NEdge>::iterator it;
  list<NEdge>::iterator to_erase[2];
  int idx = 0;

  for ( it = _nodes[v].getMAdjList().begin() ; it != _nodes[v].getMAdjList().end() ; it++ )
  {
    if (!removed[it->get_node()] && it->get_labels()[0] == v_dir)
    {
      idx_dir curr = make_pair(it->get_node(), it->get_labels()[1]);

      if ( is_same_node(curr, u) || is_same_node(curr, l) )
      {
        to_erase[idx++] = it;
      }

      if ( idx == 2 )
      {
        break;
      }
    }
  }

  _nodes[v].getMAdjList().erase(to_erase[0]);
  _nodes[v].getMAdjList().erase(to_erase[1]);

  // Adds the arc (v -> new_idx), with label (v_dir + new_dir)
  _nodes[v].add_to_adj_list( NEdge( new_idx, v_dir, new_dir ) );

  // Adds the reverse arc (new_idx -> v), with label (R(new_dir) + R(v_dir))
  new_node.add_to_adj_list( NEdge( v, reverse_dir(new_dir), reverse_dir(v_dir) ) );
}

// Looks for a bubble from 'u' in the direction 'dir'. If found,
// compress the bubble, adding a new node. Necessary with the arc (u ->
// new_node) with label 'dir + dir'; and (new_node -> u) with label
// 'R(dir) + R(dir)'.
/*!
 * \brief Find and compress bubble if it is a SNPs
 * \param u the node tried
 * \param dir the direction to go where
 * \param removed the nodes to be removed
 * \param snp_log_file indicator
 * \param bccid index of the bcc being compressed
 */
void NGraph::find_and_compress_bubble( int u, char dir, bool *removed,
				       FILE * snp_log_file, const int bccid, const bool output_context)
{
  // I'm assuming that the switching node that opens the bubble should
  // have out degree exactly 2.
  if ( get_out_degree(u, dir) != 2 )
  {
    return;
  }

  idx_dir open = make_pair(u, dir);

  // Find the two lips
  idx_dir lips[2];
  int n_lips = 0;
  list<NEdge>::const_iterator it;
  for ( it = getAdjList(u).begin() ; it != getAdjList(u).end() ; it++ )
  {
    if ( not removed[it->get_node()] && it->get_labels()[0] == dir )
    {
      lips[n_lips++] = make_pair(it->get_node(), it->get_labels()[1]);
    }
  }

  // Upper and lower lips should have the same length
  if ( _nodes[lips[0].first].getSequence().size() != _nodes[lips[1].first].getSequence().size() )
  {
    return;
  }

  // Upper and lower lips cannot be the same node
  if (lips[0].first == lips[1].first)
  {
    return;
  }

  // Check if it's possible to close the bubble.
  idx_dir close = close_bubble( open, lips[0], lips[1]);
  DEBUG(cerr << "find_and_compress_bubble, close = " << close.first << "\n";)

  if (close.first != -1)
  {
    DEBUG(fprintf(stderr, "find_and_compress_bubble, open = (%d, %c), l[0] = (%d, %c), l[1] = (%d,%c), close = (%d, %c)\n", open.first, open.second, lips[0].first, lips[0].second, lips[1].first, lips[1].second, close.first, close.second);)

    NNode new_node;
    int new_idx = (int)_nodes.size();

    // Fix neighborhoods
    fix_bubble_neighborhood( new_node, new_idx, dir, open, lips[0], lips[1], removed );
    fix_bubble_neighborhood( new_node, new_idx, reverse_dir(dir), reverse_dir(close),
         reverse_dir(lips[0]), reverse_dir(lips[1]), removed );

    // Remove the old nodes
    removed[lips[0].first] = true;
    removed[lips[1].first] = true;

    // Add the new node to _nodes and _nodeSequences
    string  new_seq = merge_nodes( dir, lips[0], lips[1] );
    new_node.setSequence( new_seq );
    _nodes.push_back( new_node );

    // Output the bubble in the log file
    output_bubble( dir, lips[0], lips[1], open, close, snp_log_file, bccid, output_context);

    // Incrementing the bubble count
    incrementNbOutput();
  }
}

/*!
 *\brief Merges two nodes creating a new one composed with the merge sequences
 *\brief  and print the corresponding lower and upper path into
 * the snp file.
 * \param dir
 * \param u
 * \param l
 * \param snp_log_file
 * \param bccid
 */
string NGraph::merge_nodes( char dir, idx_dir u, idx_dir l )
{
  string new_seq;
  string u_seq = _nodes[u.first].getSequence();
  string l_seq = _nodes[l.first].getSequence();

  if (u.second != dir) u_seq = reverse_complement(u_seq);
  if (l.second != dir) l_seq = reverse_complement(l_seq);

  new_seq = merge_sequences(u_seq, l_seq);

  return new_seq;
}


void NGraph::output_bubble( char dir, idx_dir u, idx_dir l, idx_dir open, idx_dir close,
			    FILE * snp_log_file, const int bccid, const bool output_context)
{
  string u_seq = _nodes[u.first].getSequence();
  string l_seq = _nodes[l.first].getSequence();
  
  if (u.second != dir) u_seq = reverse_complement(u_seq);
  if (l.second != dir) l_seq = reverse_complement(l_seq);
  int numCycle = getNbOutput();
  int contextL = 0;
  int contextR = 0;

  string open_seq  = open.second  == 'F' ? _nodes[open.first].getSequence()  : reverse_complement(_nodes[open.first].getSequence());
  string close_seq = close.second == 'F' ? _nodes[close.first].getSequence() : reverse_complement(_nodes[close.first].getSequence());

    // This is a bit misleading. The bubble is always output in the
    // forward direction. So when dir (the direction in which the
    // bubble was found) is 'R' we need to revert everything to
    // preserve the inner nodes in forward direction 'F'.
    if (dir == 'R')
    {
      swap(open_seq, close_seq);
      open_seq  = reverse_complement(open_seq);
      close_seq = reverse_complement(close_seq);
    }
    int k = getKValue();
    if (output_context)
    {
      contextL = open_seq.size() - k;
      contextR = close_seq.size() - k;
      u_seq = open_seq.substr(0, open_seq.size() - (k-1) ) + u_seq + close_seq.substr(k-1);
      l_seq = open_seq.substr(0, open_seq.size() - (k-1) ) + l_seq + close_seq.substr(k-1);
    }
    else
    {
      u_seq = open_seq.substr(open_seq.size()-k, 1 ) + u_seq + close_seq.substr( k-1, 1);
      l_seq = open_seq.substr(open_seq.size()-k, 1 ) + l_seq + close_seq.substr( k-1, 1);
    }

    output_sequences( u_seq, l_seq,  snp_log_file, bccid, numCycle, contextL , contextR, k );
}


// From 'curr', check if there is only one out-neighbor (in both
// directions), different from 'prev'. If it's true, return 'true' and
// put this node in 'next'.
/*!
 * \brief Look for an out-neighbor of the current node different of prev.
 * Return TRUE if one is found
 * \param prev the previous neighbor
 * \param curr the current node& direction being studied
 * \param &next the outgoing possible neighbor if found
 */
bool NGraph::follow_path(idx_dir prev, idx_dir curr, idx_dir &next) const
{
  if ( get_out_degree( curr.first, curr.second ) != 1
    || get_out_degree(  curr.first, reverse_dir( curr.second ) ) != 1 )
  {
    return false;
  }

  list<NEdge>::const_iterator it;
  for ( it = getAdjList(curr.first).begin() ; it != getAdjList(curr.first).end() ; it++ )
  {
    if ( it->get_labels()[0] == curr.second && it->get_node() != prev.first )
    {
      next = make_pair( it->get_node(), it->get_labels()[1] );
      return true;
    }
  }

  return false;
}

/*!
 * \brief Looks for the index and direction of the node closing the current
 * bubble defined by open, upper and lower
 */
idx_dir NGraph::close_bubble( idx_dir open, idx_dir upper, idx_dir lower) const
{
  idx_dir close = make_pair(-1, '-'), u_close, l_close;

  if (follow_path( open, upper, u_close) && follow_path( open, lower, l_close))
  {
    if (u_close.first == l_close.first && u_close.second == l_close.second)
    {
      close = u_close;
    }
  }
  return close;
}
// ============== Path compression methods
 // Compress all linear path and output the graph to the files.
/*!
 * \brief  Compress all linear paths of the NGraph object
 */
void NGraph::compress_all_paths( void )
{
  bool* visited = new bool[2 * _nodes.size() + MAX];
  bool* removed = new bool[2 * _nodes.size() + MAX];

  for ( int i = 0; i < (2 * (int)_nodes.size() + MAX); i++ )
  {
    visited[i] = removed[i] = false;
  }
  
  list<path_t> paths_list;
  for (int i = 0; i < (int)_nodes.size(); i++)
  {
    if (!visited[i])
    {
      DEBUG(cerr << "Explore node: " << i << "\n";)

      path_t path = find_maximal_linear_path( i, 'F');
      paths_list.push_back(path);

      // Mark all visited nodes. It may be just a single node.
      path_t::iterator it;
      for (it = path.begin(); it != path.end(); it++)
      {
        visited[it->first] = true;
      }
    }
  }
  delete [] visited;
  
  list<path_t>::iterator it;
  for (it = paths_list.begin(); it != paths_list.end(); it++)
  {
    compress_path( *it, removed);
  }
  DEBUG(cerr << "compress_all_path, out\n";)

  *this = remove_marked_nodes( *this, removed, false);
  
  delete [] removed;
}

/* \brief Compress the current path removing the inter-nodes (non branching nodes)
 * and updates the NGraph object
 * \param path: the current path
 * \param removed: nodes to be removed
 */
void NGraph::compress_path( list< pair<int,char> > path, bool removed[] )
{
  // Nothing to compress
  if (path.size() == 1) return;

  DEBUG(fprintf(stderr, "compress_path\n");)
  DEBUG(print_path(path);)

  // Mark nodes to remove
  list<pair<int,char> >::iterator it;
  for ( it = path.begin(); it != path.end(); it++ )
  {
    removed[it->first] = true;
  }

  // Fix the neighborhoods
  int first = path.back().first, last = path.front().first;
  char first_dir = path.back().second, last_dir = path.front().second;
  
  NNode new_node;
  int new_node_idx = (int)_nodes.size();

  DEBUG(cerr << "new_node_idx = " << new_node_idx << "\n";)

  // TODO: Fix the neighbors when 'path' is a cycle, currently we
  // don't add the self-loop
  fix_neighborhood( new_node, new_node_idx, first, first_dir, removed);
  fix_neighborhood( new_node, new_node_idx, last, last_dir, removed);

  // Fix the sequences
  string new_seq = get_path_sequence( path);

  DEBUG(cerr << "compress_path, new_node_nb.size = " << new_node.getAdjList().size() << "\n";)

  new_node.setSequence( new_seq );
  _nodes.push_back( new_node );
}


/*!
\brief Search the maximal linear path that includes the vertex 'start'
considering the direction 'dir'and returns it.
\param start the vertex to be included
\param dir the direction
*/
path_t NGraph::find_maximal_linear_path( int start, char dir ) const
{
  path_t path;

  path.push_back(make_pair(start, dir));

  // Explore left side (forward) of the path:
  // o <- o <- start
  int left = start;
  char left_dir = dir;
  while ((left = next_node( left, left_dir, &left_dir)) != -1 && left != start)
    path.push_front(make_pair(left, left_dir));

  // It's a cycle. We don't need to explore the other side.
  if (left == start)
    return path;

  // Explore right side (reverse) of the path:
  // start <- o <- o
  // It's equivalent to go forward on the reverse of dir.
  int right = start;
  char right_dir = reverse_dir(dir);
  while ((right = next_node( right, right_dir, &right_dir)) != -1)
    path.push_back(make_pair(right, reverse_dir(right_dir)));

  return path;
}


/*!
 \brief Checks if the path given by 'start' and direction 'dir' can be
 extended to left. Returns the next node or -1.
 \param start the starting node
 \param dir the direction to go to
 \param next_dir the future direction to go to
 */
int NGraph::next_node( int start, char dir, char *next_dir ) const
{
  int next = -1;

  if ( get_out_degree(start,dir) == 1 && get_in_degree(start,dir) == 1 )
  {
    DEBUG( fprintf( stderr, "out_degree = %d\n", get_out_degree(start,dir) ); )
    
    list<NEdge>::const_iterator it;
    
    for ( it = getAdjList(start).begin() ; it != getAdjList(start).end() ; it++ )
    {
      if ( it->get_labels()[0] == dir )
      {
        next = it->get_node();
        *next_dir = it->get_labels()[1];
        break;
      }
    }

    DEBUG(cerr << "next = " << _nodeToStr[next] << "\n";)
    
    if ( get_in_degree(next, *next_dir) == 1 && get_out_degree(next, *next_dir) == 1 )
    {
      return next;
    }
  }
  
  return -1;
}


/*!
 * \brief Concatenates the sequences of the path respecting the orientations
\param path the path
 */
string NGraph::get_path_sequence( path_t path )
{
  if (path.size() == 1)
  {
    string curr_node_seq = _nodes[path.front().first].getSequence();

    if (path.front().second == 'R')
    {
      curr_node_seq = reverse_complement(curr_node_seq);
    }

    return curr_node_seq;
  }

  // path: curr <- previous ...
  int curr = path.front().first;
  char dir = path.front().second;
  string previous;
  string curr_node;

  curr_node = ( dir == 'R' ? reverse_complement(_nodes[curr].getSequence()) : _nodes[curr].getSequence() );
  path.pop_front();
  previous = get_path_sequence( path );

  // Removes the prefix of size (k-1) from the current node
  curr_node.erase(0, _kValue - 1);

  previous += curr_node;

  return previous;
}


// Changes all the incoming arcs of 'node' to point to 'new_idx', if
// necessary flip the orientations. Returns all the the outgoing arcs
// of 'node' in 'new_node_nb'.
/*!
 * \brief update the neighborhood after the compressions of the paths
 * \param new_node a new node
 * \param new_idx index of the new node
 * \param node_id neighbor node
 * \param dir direction of the neighbor node labeled
 * \param removed the removed nodes
 */
void NGraph::fix_neighborhood( NNode& new_node, int new_idx, int node_id,
    char dir, bool removed[] )
{
  list<NEdge>::iterator it;
  for ( it = _nodes[node_id].getMAdjList().begin() ; it != _nodes[node_id].getAdjList().end() ; it++ )
  {
    if ( not removed[it->get_node()] )
    {
      DEBUG(cerr << "Fix neighborhood, node = "<< node_id << "\n";)

      // Fixing outgoing arc: 'node -> (it->node)'
      NEdge new_edge = *it;
      
      // If 'node' is used in reverse direction, we need to flip the
      // directions of the out-going arcs.
      if ( dir == 'R' )
      {
        new_edge.revert_dir( 0 );
      }

      new_node.getMAdjList().push_back( new_edge );

      // Fixing incoming arcs: 'v -> node' (there can be more than one)
      int v = it->get_node();
      list<NEdge>::iterator iit;
      for ( iit = _nodes[v].getMAdjList().begin() ; iit != _nodes[v].getMAdjList().end() ; iit++ )
      {
        if ( iit->get_node() == node_id )
        {
          DEBUG(cerr << "iit->node = " << iit->get_node() << "\n";)

          iit->set_node( new_idx );
          // If 'node' is used in reverse direction, we need to flip the
          // directions of the incoming arcs
          if ( dir == 'R' )
          {
            iit->revert_dir(1);
          }
        }
      }
    }
  }
}



// ================ Writing methods
void NGraph::print_graph_edges( FILE *stream, string (*node_label)(int), bool *filter, bool value ) const
{
  for (int u = 0; u < (int)_nodes.size(); u++ )
  {
    if ( filter == NULL || filter[u] == value )  
    {
      list<NEdge>::const_iterator it;
      for ( it = getAdjList(u).begin() ; it != getAdjList(u).end() ; it++ )
      {
        if ( filter == NULL || filter[it->get_node()] == value )
        {  
          string u_lbl = (node_label == NULL ? _nodeToStr[u] : node_label(u));
          string v_lbl = (node_label == NULL ? _nodeToStr[it->get_node()] : node_label(it->get_node()));
         
          fprintf( stream, "%s\t%s\t%s\n", u_lbl.c_str(), v_lbl.c_str(), it->get_labels().c_str() );
        }
      }
    }
  }
}

void NGraph::print_graph_nodes( FILE *stream, string (*node_label)(int), bool *filter, bool value) const
{
  for (int u = 0; u < (int)_nodes.size(); u++)
  {
    if (filter == NULL || filter[u] == value)
    {
      string u_lbl = ( node_label == NULL ? _nodeToStr[u] : node_label(u) );
      string seq   = _nodes[u].getSequence();
      string seq_r = reverse_complement( _nodes[u].getSequence() );

      fprintf( stream, "%s\t%s\t%s\n", u_lbl.c_str(), seq.c_str(), seq_r.c_str() );
    }
  }
}

//new code
void NGraph::print_graph_edges_new(  int *lines_written, FILE *stream1, FILE *stream, string (*node_label)(int), bool *filter, bool value ) const
{
  for (int u = 0; u < (int)_nodes.size(); u++ )
  {
    if ( filter == NULL || filter[u] == value )  
    {
      list<NEdge>::const_iterator it;
      for ( it = getAdjList(u).begin() ; it != getAdjList(u).end() ; it++ )
      {
        if ( filter == NULL || filter[it->get_node()] == value )
        {  
          string u_lbl = (node_label == NULL ? _nodeToStr[u] : node_label(u));
          string v_lbl = (node_label == NULL ? _nodeToStr[it->get_node()] : node_label(it->get_node()));
         
          fprintf( stream, "%s\t%s\t%s\n", u_lbl.c_str(), v_lbl.c_str(), it->get_labels().c_str() );
	  *lines_written = *lines_written + 1;
        }
      }
    }
  }
  fprintf( stream1, "%d\n",*lines_written );
}

void NGraph::print_graph_nodes_new(  int *lines_written, FILE *stream1, FILE *stream, string (*node_label)(int), bool *filter, bool value) const
{
  for (int u = 0; u < (int)_nodes.size(); u++)
  {
    if (filter == NULL || filter[u] == value)
    {
      string u_lbl = ( node_label == NULL ? _nodeToStr[u] : node_label(u) );
      string seq   = _nodes[u].getSequence();
      string seq_r = reverse_complement( _nodes[u].getSequence() );

      fprintf( stream, "%s\t%s\t%s\n", u_lbl.c_str(), seq.c_str(), seq_r.c_str() );
      *lines_written = *lines_written + 1;
    }
  }
  fprintf( stream1, "%d\n",*lines_written );
}

//end new code

list<NEdge>::iterator NGraph::eraseAdjList( int node_id, int w )
{
  list<NEdge>::iterator it;
  for ( it = _nodes[node_id].getMAdjList().begin() ; it != _nodes[node_id].getMAdjList().end() ; it++ )
  {
    if ( it->get_node() == w )
    {
      return _nodes[node_id].getMAdjList().erase( it );
    }
  }
  
  // Should never be reached
  assert( false );
  return it;
}


// ===========================================================================
//                               Non inline accessors
// ===========================================================================
int NGraph::get_in_degree( int node, char dir ) const
{
  return get_out_degree( node, reverse_dir(dir) );
}

int NGraph::get_out_degree( int node, char dir ) const
{
	int degree = 0;

	list<NEdge>::const_iterator it;
  
	for ( it = getAdjList(node).begin() ; it != getAdjList(node).end() ; it++ )
  {
		if ( it->get_labels()[0] == dir )
    {
			degree++;
    }
  }
      
	return degree;
}


// ===========================================================================
//                          Utilities graph's functions
// ===========================================================================
// TODO: Move to graph_utils

void read_graph_edges(NGraph& g, FILE* edge_file)
{
  char *buffer = new char[100 * MAX], *u = new char[MAX], *v = new char[MAX], *label = new char[MAX];

  while (fgets(buffer, 100 * MAX, edge_file) != NULL)
  {
    char *p;

    // out going node
    p = strtok(buffer, "\t\n");
    strcpy(u, p);

    // in coming node
    p = strtok(NULL, "\t\n");
    strcpy(v, p);

    // edge label
    p = strtok(NULL, "\t\n");
    strcpy(label, p);

    // Coverage values, ignore them for now
    //while (p != NULL)
    //  p = strtok(NULL, "\t\n");

    g.insert_edge(u, v, label);
  }
}

// Sequence size is limited to 100 * MAX
void read_graph_nodes(NGraph& g, FILE* node_file)
{
  char *buffer = new char[100 * MAX], *str = new char[MAX], *seq = new char[100*MAX];
  fprintf(stdout, "NGraph : read_graph_nodes \n");

  while (fgets(buffer, 100 * MAX, node_file) != NULL)
  {
    char *p;
    fprintf(stdout, "NGraph : buffer nodes: %s\n",buffer);

    if (strlen(buffer) == 100 * MAX)
    {  
      p = strtok(buffer, "\t\n");
      fprintf(stdout, "ERROR: node %s with sequence larger than %d!", p, 100 * MAX);
      exit(0);
    }

    // Node label
    p = strtok(buffer, "\t\n");
    strcpy(str, p);

    // Node seq
    p = strtok(NULL, "\t\n");
    strcpy(seq, p);

    g.insert_node(str, seq); 
  }
}
/*!
 * \brief Returns a NGraph like g without some nodes (marked nodes)
 *
 * \param g the original NGraph object
 * \param filter pointer to a boolean array, contains one value per nodes
 * \param value the removal value (T/F). If for a node i of the _nodes vector attribute, filter[i] equals to value
 * then the node i will not be into the remaining graph.
 */
//new code
void read_graph_edges_new(NGraph& g, FILE* info_file,  FILE* contents_file, FILE* data_file, int *required_sequence, int *file_index)
{

  char *buffer = new char[100 * MAX], *u = new char[MAX], *v = new char[MAX], *label = new char[MAX];

  int written_lines_before, written_lines_for_record, written_lines;
  // have to scan the  
  fprintf(stdout, "NGraph : NEW read_graph_edges read required sequence %d file_index = %d\n",*required_sequence,*file_index);
  //normally would skip directly to the record, but for the time being use the formatted file
  for (int i = 0; i < *required_sequence+*file_index-2 ; i++) fgets(buffer, 100 * MAX, contents_file);
  fscanf(contents_file, "%d \n",&written_lines_before );
  fscanf(contents_file, "%d \n",&written_lines_for_record );
  fprintf(stdout, "\n written_lines_before = %d written_lines_for_record = %d\n",written_lines_before,written_lines_for_record);
  // if a graph edge  was not written to disk   
  if ( written_lines_before == written_lines_for_record) 
    { fprintf(stderr, "The required record %d does not exist in the bcc graph !\n", *required_sequence); exit(0); }
  
  // where to search the record
  written_lines =  written_lines_for_record - written_lines_before;
  //normally would skip directly to the record, but for the time being use the formatted file
  for (int i = 0; i < written_lines_before ; i++) fgets(buffer, 100 * MAX, data_file);
  
  for (int i = 0; i < written_lines ; i++) {
    char *p;
    fgets(buffer, 100 * MAX, data_file);
    fprintf(stdout, "NGraph : buffer edges: %s\n",buffer);
    // out going node
    p = strtok(buffer, "\t\n");
    strcpy(u, p);

    // in coming node
    p = strtok(NULL, "\t\n");
    strcpy(v, p);

    // edge label
    p = strtok(NULL, "\t\n");
    strcpy(label, p);

    // Coverage values, ignore them for now
    //while (p != NULL)
    //  p = strtok(NULL, "\t\n");

    g.insert_edge(u, v, label);
  }
    
}

// Sequence size is limited to 100 * MAX
void read_graph_nodes_new(NGraph& g, FILE* info_file,  FILE* contents_file, FILE* data_file, int *required_sequence, int *file_index)
{
  char *buffer = new char[100 * MAX], *str = new char[MAX], *seq = new char[100*MAX];

  int written_lines_before, written_lines_for_record, written_lines;
  // have to scan the  
  fprintf(stdout, "NGraph : NEW read_graph_nodes read required sequence %d file_index = %d\n",*required_sequence,*file_index);
  //normally would skip directly to the record, but for the time being use the formatted file
  for (int i = 0; i < *required_sequence+*file_index-2 ; i++) fgets(buffer, 100 * MAX, contents_file);
  fscanf(contents_file, "%d \n",&written_lines_before );
  fscanf(contents_file, "%d \n",&written_lines_for_record );
  fprintf(stdout, "\n written_lines_before = %d written_lines_for_record = %d\n",written_lines_before,written_lines_for_record);
  // if a graph node  was not written to disk   
  if ( written_lines_before == written_lines_for_record) 
    { fprintf(stderr, "The required record %d does not exist in the bcc graph !\n", *required_sequence); exit(0); }
  
  // where to search the record
  written_lines =  written_lines_for_record - written_lines_before;
  //normally would skip directly to the record, but for the time being use the formatted file
  for (int i = 0; i < written_lines_before ; i++) { fgets(buffer, 100 * MAX, data_file); }

  for (int i = 0; i < written_lines ; i++) {
    char *p;
    fgets(buffer, 100 * MAX, data_file);
    fprintf(stdout, "NGraph : buffer nodes: %s\n",buffer);
    
    if (strlen(buffer) == 100 * MAX)
      {  
	p = strtok(buffer, "\t\n");
	fprintf(stdout, "ERROR: node %s with sequence larger than %d!", p, 100 * MAX);
	exit(0);
      }
    
    // Node label
    p = strtok(buffer, "\t\n");
    strcpy(str, p);
    
    // Node seq
    p = strtok(NULL, "\t\n");
    strcpy(seq, p);
    
    g.insert_node(str, seq); 
  }

}
// end new code
NGraph remove_marked_nodes( NGraph& g, bool* filter, bool value )
{
  if ( filter == NULL )
  {
    return g;
  }

  NGraph new_g( g.getKValue(), g.getNbOutput() );

  for ( int u = 0 ; u < (int)g.getNbNodes() ; u++ )
  {
    if ( filter[u] == value )
    {
      string u_sequence = g.getSequence(u);
      new_g.insert_node(std::to_string(u), u_sequence);
      if ( filter[u] == value )  
      {
        list<NEdge>::const_iterator it;
        for ( it = g.getAdjList(u).begin() ; it != g.getAdjList(u).end() ; it++ )
        {
          if ( filter[it->get_node()] == value )
          {
            new_g.insert_edge(std::to_string(u), std::to_string(it->get_node()), it->get_labels());
          }
        }
      }
    }
  }
  // For the new_g, need to re-expand the possible parallel edges
  new_g.expand_parallel_edges();
  return new_g;
}

// Only use in DEBUG
void print_path(list<pair<int,char> > path)
{
  if (path.size() != 1)
  {
    fprintf (stderr, "Path: ");
    list<pair<int,char> >::iterator it;
    for (it = path.begin(); it != path.end(); it++)
      fprintf(stderr, "(%d, %c) ", it->first, it->second);
    fprintf(stderr, "\n");
  }
}