treedex.h

/**************************************************************************
 * Copyright (c) 2002-2011 T. M. Murali                                   *
 * Copyright (c) 2007-2011 Naveed Massjouni                               *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm 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.                                    *
 *                                                                        *
 * Biorithm 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      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef TREEDECOMPOSITION_H
#define TREEDECOMPOSITION_H

#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <set>
#include <sstream>
#include <stack>
#include "graph.h"
#include "treewidth.h"
using namespace  std;

typedef MyNode Node;
typedef MyEdge Edge;
typedef MyNode::EdgeIterator EdgeIterator;
typedef MyGraph Graph;
typedef MyGraph::NodeIterator NodeIterator;
typedef string NodeId;
typedef string BagId;

class TreeDecomposition {

private:
  map< BagId, set<NodeId> > bagMap; // maps bags to set of nodes they contain
  map< BagId, set<NodeId> > blackNodesMap; 
  Graph _tree; // underlying graph representing the tree structure
  const Graph _origGraph; // original graph to be decomposed
  vector<BagId> bagOrder; // the order the bags are added
  vector<NodeId> nodeOrder; // the order the nodes are processed
  int bagCount;
  int nodeCount;

  vector<NodeId> getNodes(Graph graph) {
    vector<NodeId> nodes;
    NodeIterator nodesItr = graph.nodes();
    while (nodesItr.hasNext()) {
      Node node = nodesItr.next();
      nodes.push_back(node.getId());
    }
    return nodes;
  }

  vector<Edge> getEdges(Graph& graph) {
    vector<Edge> edges;
    Graph::MyEdgeIterator edgeItr = graph.edges();
    while (edgeItr.hasNext()) {
      Edge edge = edgeItr.next();
      edges.push_back(edge);
    }
    return edges;
  }

  vector<NodeId> getNeighbors(Graph& graph, NodeId nodeId) {
    vector<NodeId> neighbors;
    EdgeIterator edgeItr = graph.incidentEdges(nodeId);
    while (edgeItr.hasNext()) {
      const Edge& edge = *edgeItr.next();
      Node* nodePtr = edge.opposite(nodeId);
      neighbors.push_back(nodePtr->getId());
    }
    return neighbors;
  }

  set<NodeId> getNeighborSet(Graph& graph, NodeId nodeId) {
    cout << endl << "getNeighborSet called: ";
    set<NodeId> result;
    vector<NodeId> neighbors = getNeighbors(graph, nodeId);
    for (int i = 0; i < neighbors.size(); i++) {
      result.insert(neighbors[i]);
    }
    cout << result << endl;
    return result;
  }

  // Connect the neighbors of a node to each other.
  void connectNeighbors(Graph& graph, NodeId nodeId) {
    vector<NodeId> neighbors = getNeighbors(graph, nodeId);
    for (int i = 0; i < neighbors.size(); i++) {
      for (int j = i+1; j < neighbors.size(); j++) {
        graph.addEdge(neighbors[i], neighbors[j]);
      }
    }
  }

  // Determine if graph is a tree.
  bool isTree(Graph& graph) {
    if (graph.numNodes() != graph.numEdges() + 1)
      return false;
    Treewidth treewidth(graph);
    int prevSize = 0;
    while (treewidth.numNodes() != prevSize) {
      prevSize = treewidth.numNodes();
      treewidth.runRule(Treewidth::TWIG_RULE);
    }
    return treewidth.numNodes() == 1;
  }

  // Determine if graph is a clique.
  bool isClique(Graph& graph) {
    unsigned int numNodes = graph.numNodes();
    unsigned int numEdges = graph.numEdges();
    if (numNodes*(numNodes - 1)/2 == numEdges)
      return true;
    return(false);
  }

  // Find a bag in bagOrder with the minimum index
  // that contains every node in neighborSet.
  BagId findBagWithNeighbors(set<NodeId> neighborSet, int minIndex) {
    BagId bagId;
    set<BagId> bagSet;
    vector<BagId> intersection(neighborSet.size());
    vector<BagId>::iterator endOfIntersection;
    for (int i = minIndex; i < bagOrder.size(); i++) {
      bagId = bagOrder[i];
      set<BagId> bagSet = bagMap[bagId];
      endOfIntersection = set_intersection(
        neighborSet.begin(), neighborSet.end(),
        bagSet.begin(), bagSet.end(), intersection.begin());
      int intersectionSize = endOfIntersection - intersection.begin();
      if (intersectionSize == neighborSet.size())
        return bagId;
    }
    return "";
  }

  BagId addBag(const set<NodeId>& nodeSet) {
      ostringstream oss;
      oss << "bag" << (1000000 + bagCount);
      BagId bagId = oss.str();
      _tree.addNode(bagId);
      bagMap[bagId] = nodeSet;
      bagCount++;
      updateBlackNodes(bagId, nodeSet);
      return bagId;
  }

  void updateBlackNodes(BagId bagId, const set<NodeId>& nodeSet) {
    Graph induced;
    _origGraph.computeInducedSubgraph(nodeSet, induced);

    NodeIterator itr = induced.nodes();
    while (itr.hasNext()) {
      Node node = itr.next();
      if (node.degree() == 0)
        blackNodesMap[bagId].insert(node.getId());
    }
  }

public:
  TreeDecomposition()
    : bagCount(0), nodeCount(0) {}

  TreeDecomposition(Graph graph)
    : bagCount(0), nodeCount(0), _origGraph(graph) {}

  void printBags() {
    for (map< BagId, set<NodeId> >::iterator itr = bagMap.begin();
         itr != bagMap.end(); itr++) {
      //cout << itr->first << "|" << itr->second << endl;
      BagId bagId = itr->first;
      cout << bagId << "|";
      set<NodeId> nodes = itr->second;
      for (set<NodeId>::iterator nodeItr = nodes.begin();
           nodeItr != nodes.end(); nodeItr++) {
        NodeId nodeId = *nodeItr;
        cout << nodeId;
        set<NodeId> blackNodes = blackNodesMap[bagId];
        if (blackNodes.find(nodeId) != blackNodes.end())
          cout << "*"; // is black node, so mark with *
        cout << ' ';
      }
      cout << endl;
    }
  }

  int getMaxBagSize() {
    int max = 0;
    map< BagId, set<NodeId> >::iterator itr = bagMap.begin();
    while (itr != bagMap.end()) {
      if (itr->second.size() > max)
        max = itr->second.size();
      itr++;
    }
    return max;
  }

  BagId getMaxBagId() {
    BagId maxBagId = "undefined";
    int max = 0;
    map< BagId, set<NodeId> >::iterator itr = bagMap.begin();
    while (itr != bagMap.end()) {
      if (itr->second.size() > max) {
        max = itr->second.size();
        maxBagId = itr->first;
      }
      itr++;
    }
    return maxBagId;
  }

  set<NodeId> getBlackNodesForBagId(BagId bagId) {
    return blackNodesMap[bagId];
  }

  class TreeWidthNodeSelector
  {
  protected:
    // the graph the node selector will operate on. it is a reference
    // since the TreeDecomposition class will also operate on the same instance.
    MyGraph &_graph;
  public:
    TreeWidthNodeSelector(MyGraph & g)
        : _graph(g)
      {}

    virtual ~TreeWidthNodeSelector()
      {}

    virtual void initialise() = 0;

    virtual MyNodeId getNextNode() = 0;
    
    virtual void update(MyNodeId deletedNodeId, set< MyNodeId > &neighbours)=0;
    
  };

  class TreeWidthNodeSelectorMinDegree : public TreeWidthNodeSelector
  {
  protected:
    // the current minimum degree.
    unsigned int _minDegree;
    
    // store nodes in order of degree. for each possible degree (the
    // index of the vector, store a list of nodes with that degree.
    vector< list< MyNodeId > > _degreeLists;
    // for each node, store a pointer to the list where it is stored.
    // this map will be useful for efficiently moving nodes from one
    // list to another.
    map< MyNodeId, list< MyNodeId >::iterator > _nodePositions;
  public:
    TreeWidthNodeSelectorMinDegree(MyGraph & g)
        : TreeWidthNodeSelector(g), _minDegree(0),
          _degreeLists(_graph.numNodes() + 1),
          _nodePositions()
      {}

    virtual ~TreeWidthNodeSelectorMinDegree()
      {}

    //virtual void initialise();
    virtual void initialise() {
      _graph._createDegreeLists(_minDegree, _degreeLists, _nodePositions);
    }
    
    virtual MyNodeId getNextNode() {
      static unsigned int numNodes = 0;
      // remove the node with smallest degree. get the element at the
      // front of the right list and also pop it.
      string nodeId = _degreeLists[_minDegree].front();
      _degreeLists[_minDegree].pop_front();
      cout << ++numNodes << " Returning node " << nodeId << " with degree "
           << _minDegree << ". ";
      return(nodeId);
    }


    virtual void  update(MyNodeId deletedNodeId, set<NodeId> &neighbours) {
      // decrement the degree of all its neighbours.
      _graph._decrementDegrees(_degreeLists, _nodePositions, neighbours);
      // update minDegree. it either decreases by 1, stays the same,
      // or increases.
      if ((0 != _minDegree) && (0 != _degreeLists[_minDegree - 1].size()))
        _minDegree--;
      else if (0 != _degreeLists[_minDegree].size()) {} // do nothing.
      else while ((_minDegree <= _graph.numNodes()) &&
                  (0 == _degreeLists[++_minDegree].size()));
    }

  };

  struct MyNodeFillInRatioPQNode
  {
    MyNT _fillInRatio;
    unsigned int _fillIn;
    MyNodeId _id;
  };

  // function object for comparing two node PQ nodes. i want nodes with
  // smaller weight at top, so the comparison is not < but >
  // (priority_queue::top() returns the max element.)
  struct compareFillInValues : public
    binary_function<
      const MyNodeFillInRatioPQNode&, const MyNodeFillInRatioPQNode&, bool> 
  {
    bool operator()(
        const MyNodeFillInRatioPQNode& __x,
        const MyNodeFillInRatioPQNode& __y) const
    {
      return(__x._fillIn > __y._fillIn);
    }
  };

  class TreeWidthNodeSelectorMinFillIn : public TreeWidthNodeSelector
  {
  protected:
    typedef priority_queue< MyNodeFillInRatioPQNode,
                            vector< MyNodeFillInRatioPQNode >,
                            compareFillInValues >
    MyFillInPQ;
    MyFillInPQ _npq;

    vector<NodeId> getNeighbors(Graph& graph, NodeId nodeId) {
      vector<NodeId> neighbors;
      EdgeIterator edgeItr = graph.incidentEdges(nodeId);
      while (edgeItr.hasNext()) {
        const Edge& edge = *edgeItr.next();
        Node* nodePtr = edge.opposite(nodeId);
        neighbors.push_back(nodePtr->getId());
      }
      return neighbors;
    }

    // Compute the number of pairs of neighbors of a node that are not connected
    unsigned int computeFillIn(Graph& graph, NodeId nodeId) {
      vector<NodeId> neighbors = getNeighbors(graph, nodeId);
            unsigned int fillIn = 0;
      for (int i = 0; i < neighbors.size(); i++) {
        for (int j = i+1; j < neighbors.size(); j++) {
                      if (!graph.isEdge(neighbors[i], neighbors[j]))
                        fillIn++;
                        }
      }
            return(fillIn);
    }

  public:
    TreeWidthNodeSelectorMinFillIn(MyGraph & g)
        : TreeWidthNodeSelector(g), _npq() {}

    virtual ~TreeWidthNodeSelectorMinFillIn() {}

    virtual void initialise() {
      MyGraph::NodeIterator nitr = _graph.nodes();
      MyNodeFillInRatioPQNode npqNode;
      MyNT quotient;
      unsigned int degree;
      while (nitr.hasNext())
        {
          MyNode &node = nitr.next();
          npqNode._id = node.getId();
          npqNode._fillIn = computeFillIn(_graph, npqNode._id);
          degree = node.degree();
          quotient = degree*(degree - 1)/2.0;
          if (0 == quotient)
            npqNode._fillInRatio = 0;
          else
            npqNode._fillInRatio = npqNode._fillIn/quotient;
          _npq.push(npqNode);
        }
    }
    
    virtual MyNodeId getNextNode() {
      static unsigned int numNodes = 0;
      MyNodeFillInRatioPQNode pqNode;
      do
        {
          pqNode = _npq.top();
          _npq.pop();
        }
      // the pq may have incorrect and old info for nodes i have already
      // processed, so i keep popping nodes until i find a node in _graph.
      while (!_graph.isNode(pqNode._id));
      cout << ++numNodes << " Returning node " << pqNode._id
           << " with fill in " << pqNode._fillIn
           << " and fill in ratio " << pqNode._fillInRatio << ". ";

      return(pqNode._id);
    }

    virtual void  update(MyNodeId deletedNodeId, set<NodeId> &neighbours) {
      set<NodeId>::const_iterator nitr;
      MyNodeFillInRatioPQNode npqNode;
      unsigned int degree;
      MyNT quotient;
      for (nitr = neighbours.begin(); nitr != neighbours.end(); nitr++)
        {
          npqNode._id = *nitr;
          npqNode._fillIn = computeFillIn(_graph, npqNode._id);
          degree = _graph._getNode(npqNode._id).degree();
          quotient = degree*(degree - 1)/2.0;
          if (0 == quotient)
            npqNode._fillInRatio = 0;
          else
            npqNode._fillInRatio = npqNode._fillIn/quotient;
          _npq.push(npqNode);
        }
    }

  };

  // function object for comparing two edge PQ nodes for computing
  // minimum spanning trees. i want nodes with smaller weight at top, so
  // the comparison is not < but > (priority_queue::top() returns the
  // max element.) 
  struct compareFillInRatioValues
    : public binary_function< const MyNodeFillInRatioPQNode&,
                              const MyNodeFillInRatioPQNode&, bool> 
  {
    bool operator()(const MyNodeFillInRatioPQNode& __x,
                    const MyNodeFillInRatioPQNode& __y) const
    {
      return(__x._fillInRatio > __y._fillInRatio);
    }
  };

  class TreeWidthNodeSelectorMinFillInRatio
    : public TreeWidthNodeSelectorMinFillIn
  {
  protected:
    // the only change between TreeWidthNodeSelectorMinFillInRatio and
    // TreeWidthNodeSelectorMinFillIn is the compare function in the pq.
    typedef priority_queue<
      MyNodeFillInRatioPQNode,
      vector< MyNodeFillInRatioPQNode >,
      compareFillInRatioValues > MyFillInRatioPQ;

    MyFillInRatioPQ _npq;
  public:
    TreeWidthNodeSelectorMinFillInRatio(MyGraph & g)
        : TreeWidthNodeSelectorMinFillIn(g), _npq()
      {}

    virtual ~TreeWidthNodeSelectorMinFillInRatio()
      {}
  };

  // making pq tree efficient (avoid the "if (!graph.isNode(nodeId))"
  // check). store MyNodeFillInPQNodes in a two sets, one sorted by _id
  // and the other by _fillIn. pop from the _fillIn map. pop same _id
  // from the _id map. create MyNodeFillInPQNodes for all the
  // neighbours. for each neighbour, find the entry in the _id map,
  // update the entry, find the entry in the _fillIn map, delete the
  // entry, and insert the new entry.

  // Recursively find tree decomposition.
  void decompose() {
    Graph graph = _origGraph;

    unsigned int maxBagSize  = 0;
    MyNodeId id, nodeId;
    
    //TreeWidthNodeSelectorMinFillInRatio twns(graph);
    //TreeWidthNodeSelectorMinFillIn twns(graph);
    TreeWidthNodeSelectorMinDegree twns(graph);
    twns.initialise();
    
    while (!isTree(graph) && !isClique(graph)) {
      nodeId = twns.getNextNode();
      nodeOrder.push_back(nodeId);

#ifdef DEBUG
      // we may have already processed this node.
      if (!graph.isNode(nodeId)) {
        cerr << "Node id " << nodeId
          << " is not a graph node but the node selector returned it."
          << endl;
        continue;
      }
#endif      
      set<NodeId> bagNodes = getNeighborSet(graph, nodeId);
      connectNeighbors(graph, nodeId);
      graph.deleteNode(nodeId);
      
//     cout << ++nodeCount << " Deleted node " << nodeId << " with "
//          << bagNodes.size() << " neighbors" << endl;

      twns.update(nodeId, bagNodes);

      bagNodes.insert(nodeId);
      BagId newBagId = addBag(bagNodes);
      cout << "adding bag: " << bagNodes << endl;
      bagOrder.push_back(newBagId);
      if (bagNodes.size() > maxBagSize)
        maxBagSize = bagNodes.size();
      cout << " Graph size is (" << graph.numNodes()
           << ", " << graph.numEdges() << "). Treewidth is "
           << maxBagSize << "." << endl;
    }

    cout << "Base Case (tree/clique): " << endl;
    //graph.printEdges(cout);
    set<NodeId> nodeSet;
    graph.getNodeSet(nodeSet);
    BagId newBagId = addBag(nodeSet);
    bagOrder.push_back(newBagId);

    // Note that bagOrder should have 1 more element than nodeOrder.
    // Process the bags in the reverse order that they were added.
    // For each bag b_i, there is an associated node n_i,
    // except for the last bag that was added.
    // For a given bag b_i, consider only bags b_j: j > i.
    // Add an edge from a bag b1 in b_i to a bag b2 in b_j
    // s.t. (b1 - n1) is a subset of b2.
    cout << "Going to add edges" << endl;
    for (int i = nodeOrder.size()-1; i >= 0; i--) {
      cout << ".";
      set<NodeId> naybears = bagMap[bagOrder[i]];
      naybears.erase(nodeOrder[i]);
      BagId bagToConnect = findBagWithNeighbors(naybears, i+1);
      _tree.addEdge(bagOrder[i], bagToConnect);
    }
    cout << endl;

  }

  void computeInducedSubtree(const NodeId nodeId, TreeDecomposition& subtree) {
    set<NodeId> bag;
    bag.insert(nodeId);
    set<BagId> bagsWithNode;

    NodeIterator bagItr = _tree.nodes();
    while (bagItr.hasNext()) {
      BagId bagId = bagItr.next().getId();
      set<NodeId> nodeSet = bagMap[bagId];
      if (nodeSet.find(nodeId) != nodeSet.end()) {
        cout << "found bag: " << bagId << endl;
        bagsWithNode.insert(bagId);
        subtree.bagMap[bagId] = bag;
      }
    }

    _tree.computeInducedSubgraph(bagsWithNode, subtree._tree);
  }

};

#endif