$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2003-2014 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#include "qub_findloops.h"

/*

*/

#include <math.h>

#include <iostream>
#include <list>
#include <map>
#include <vector>

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/prim_minimum_spanning_tree.hpp>

using namespace std;
using namespace boost;

bool findloops_findpath(map< size_t, list<size_t> >& edges, vector<int>& path, size_t from, size_t to, set<size_t>& visited)
{
  visited.insert(from);
  for (list<size_t>::iterator ei = edges[from].begin(); ei != edges[from].end(); ++ei) {
    path.push_back( (int) *ei );
    if ( *ei == to )
      return true;
    if ( (visited.find(*ei) == visited.end()) && findloops_findpath(edges, path, *ei, to, visited) )
      return true;
    path.pop_back();
  }
  visited.erase(from);
  return false;
}

bool findloops_findpath(map< size_t, list<size_t> >& edges, vector<int>& path, size_t from, size_t to)
{
  path.push_back((int) from);
  set<size_t> visited;
  return findloops_findpath(edges, path, from, to, visited);
}



/*
Finds the fundamental cycles in an undirected graph.  Returns the number of cycles.

  Nvertex: a.k.a. Nstate
  Nedge: actually twice the number of edges, counting (a, b) and (b, a)
  in_edges: connected vertex pairs, alternating from0, to0, from1, to1; include both (a, b) and (b, a)
  out_loopsizes: on exit, the number of states in each fundamental cycle; allocate at least Nedge
  out_loopstates: on exit, the states in each cycle i are listed starting at out_loopstates[i*Nvertex]

  See David Colquhoun, Kathryn Dowsland, Marco Beato and Andrew Plested. (2004). How to impose microscopic reversibility in complex reaction mechanisms. (with appendix by Kathryn. A. Dowsland and Frank G. Ball)  Biophysical Journal, 86, 3510–3518.

*/

extern "C" QUBFAST_API int    qub_findloops(int Nvertex, int Nedge, int *in_edges, int *out_loopsizes, int *out_loopstates)
{
  typedef adjacency_list < vecS, vecS, undirectedS,
    property<vertex_distance_t, int>, property < edge_weight_t, int > > Graph;
  typedef pair<int, int> E;
  
  vector<E> edges;
  vector<int> weights;
  for (int i=0; i<Nedge; ++i) {
    edges.push_back( E( in_edges[2*i], in_edges[2*i+1] ) );
    weights.push_back(1);
  }
  
  Graph g(Nvertex);
  property_map<Graph, edge_weight_t>::type weightmap = get(edge_weight, g);
  for (size_t jj = 0; jj < edges.size(); ++jj) {
    graph_traits<Graph>::edge_descriptor e; bool inserted;
    tie(e, inserted) = add_edge(edges[jj].first, edges[jj].second, g);
    weightmap[e] = weights[jj];
  }
  
  vector< graph_traits< Graph >::vertex_descriptor > pred_of_v( num_vertices(g) );
  prim_minimum_spanning_tree(g, &pred_of_v[0]);
  
  map<size_t, list<size_t> > keptEdges;
  for (size_t i=0; i<pred_of_v.size(); ++i) {
    if ( i == pred_of_v[i] ) // no predecessor: unconnected, or tree root
      continue;
    
    keptEdges[i].push_back( pred_of_v[i] );
    keptEdges[pred_of_v[i]].push_back( i );

    // for (size_t j = edges.size()-1; j >= 0; --j) {
    size_t j = edges.size();
    if ( j-- ) {
      do {
	if ( (edges[j].first == ((int) pred_of_v[i]) && edges[j].second == ((int) i))
	     || (edges[j].first == ((int) i) && edges[j].second == ((int) pred_of_v[i])) ) {
	  edges.erase( edges.begin() + j );
	  break;
	}
      } while ( j-- );
    }
  }
  // now edges has the erased ones (that complete the loops)
  int Nloop = 0;
  
  for (vector<E>::iterator ei=edges.begin(); ei != edges.end(); ++ei) {
    std::vector<int> path;
    if ( findloops_findpath(keptEdges, path, ei->first, ei->second) ) {
      if ( path.size() > 2 ) {
	out_loopsizes[Nloop] = (int) path.size();
	for (size_t j=0; j<path.size(); ++j)
	  out_loopstates[Nloop*Nvertex + j] = path[j];
	++Nloop;
      }
    } else {
      cerr << "missing loop!" << endl;
    }
  }
  
  return Nloop;
}