treewidth.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/>.      *
 *                                                                        *
 **************************************************************************/

/*
  Treewidth is a class that reduces a graph by a series of "safe" rules.
  It also generates a minimum value for the treewidth of the graph.
  The rules that have been implemented so far are:
  * Islet Rule
  * Twig Rule
  * Series Rule
  * Triangle Rule
  * Simplicial Rule
  * Almost Simplicial Rule

  Example usage:
    Treewidth treewidth("graph.in");
    TODO: fill this in
    cout << "low = " << treewidth.getLow() << endl;
*/

#ifndef TREEWIDTH_H
#define TREEWIDTH_H

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

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

class Treewidth {

  private:

  int low;
  Graph graph;

  public:

  enum RuleType {ISLET_RULE = 1, TWIG_RULE = 2, SERIES_RULE = 4,
    TRIANGLE_RULE = 8, SIMPLICIAL_RULE = 16, ALMOST_SIMPLICIAL_RULE = 32,
    ALL_RULES_MASK = 0xffff};

  // Constructors
  Treewidth(string graphFile) : graph(graphFile, "asdf"), low(0) {}
  Treewidth(Graph graphIn) : graph(graphIn), low(0) {}

  void getGraph(Graph &g) const {
    g = graph;
  }

  // Returns a copy of the graph
  Graph getGraphCopy() {
    Graph g = graph;
    return g;
  }

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

  int getDegree(NodeId nodeId) {
    return graph.degree(nodeId);
  }

  int getLow() {
    return low;
  }

  int numNodes() {
    return graph.numNodes();
  }

  int numEdges() {
    return graph.numEdges();
  }

  void printEdges(ostream& ostr) {
    graph.printEdges(ostr);
  }

  void runRule(RuleType rule) {
    stack<string> nodesToProcess;
    set<string> deletedNodes;
    NodeIterator nodesItr = graph.nodes();
    while (nodesItr.hasNext()) {
      Node node = nodesItr.next();
      nodesToProcess.push(node.getId());
    }
    while (!nodesToProcess.empty()) {
      NodeId nodeId = nodesToProcess.top();
      nodesToProcess.pop();
      // Skip this node if we have already deleted it.
      if (deletedNodes.find(nodeId) != deletedNodes.end())
        continue;
      //cout << "looking at node: " << nodeId << endl;
      int degree = graph.degree(nodeId);
      switch(rule) {
        case ISLET_RULE:
          if (degree == 0) {
            //cout << "deleting it by islet rule" << endl;
            graph.deleteNode(nodeId);
            deletedNodes.insert(nodeId);
          }
          break;
        case TWIG_RULE:
          if (degree == 1) {
            pushNeighbors(nodesToProcess, nodeId);
            //cout << "deleting it by twig rule" << endl;
            graph.deleteNode(nodeId);
            deletedNodes.insert(nodeId);
            updateLow(degree);
          }
          break;
        case SERIES_RULE:
          if (degree == 2) {
            pushNeighbors(nodesToProcess, nodeId);
            connectNeighbors(nodeId);
            //cout << "deleting it by series rule" << endl;
            graph.deleteNode(nodeId);
            deletedNodes.insert(nodeId);
            updateLow(degree);
          }
          break;
        case TRIANGLE_RULE:
          if (degree == 3) {
            vector<NodeId> neighbors = getNeighbors(nodeId);
            if (graph.isEdge(neighbors[0], neighbors[1])
                || graph.isEdge(neighbors[1], neighbors[2])
                || graph.isEdge(neighbors[0], neighbors[2])) {
              pushNeighbors(nodesToProcess, nodeId);
              connectNeighbors(nodeId);
              //cout << "\tdeleting it by triangle rule" << endl;
              graph.deleteNode(nodeId);
              deletedNodes.insert(nodeId);
              updateLow(degree);
            }
          }
          break;
        case SIMPLICIAL_RULE:
          if (isSimplicial(nodeId)) {
            pushNeighbors(nodesToProcess, nodeId);
            //cout << "deleting it by simplicial rule" << endl;
            graph.deleteNode(nodeId);
            deletedNodes.insert(nodeId);
            updateLow(degree);
          }
          break;
        case ALMOST_SIMPLICIAL_RULE:
          if (isAlmostSimplicial(nodeId)) {
            pushNeighbors(nodesToProcess, nodeId);
            connectNeighbors(nodeId);
            //cout << "deleting it by almost simp. rule" << endl;
            graph.deleteNode(nodeId);
            deletedNodes.insert(nodeId);
            updateLow(degree);
          }
          break;
      }
    }
  }

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

  // Push the neighbors of a node onto the given stack.
  void pushNeighbors(stack<NodeId>& s, NodeId nodeId) {
    vector<NodeId> neighbors = getNeighbors(nodeId);
    for (int i = 0; i < neighbors.size(); i++) {
      s.push(neighbors[i]);
    }
  }

  void updateLow(int degree) {
    low = degree > low ? degree : low;
  }

  // A node is simplicial if its neighbors induce a clique.
  bool isSimplicial(NodeId nodeId) {
    if (graph.degree(nodeId) <= 1)
      return true;
    vector<NodeId> neighbors = getNeighbors(nodeId);
    int size = neighbors.size();
    for (int i = 0; i < size; i++) {
      for (int j = i+1; j < size; j++) {
        if (!graph.isEdge(neighbors[i], neighbors[j])) {
          return false;
        }
      }
    }
    return true;
  }

  // A node v almost simplicial if there is a neighbor w (lonely node) of v
  // such that all other neighbors of v form a clique.
  bool isAlmostSimplicial(NodeId nodeId) {
    if (graph.degree(nodeId) <= 2)
      return true;
    map< NodeId, set<NodeId> > neigborsMap;
    vector<NodeId> neighbors = getNeighbors(nodeId);
    int numNeighbors = neighbors.size();
    // Map a list of the neighbors' neighbors to each neighbor.
    // Only consider neighbors of the original nodeId that was passed in.
    for (int i = 0; i < numNeighbors ; i++) {
      set<NodeId> neighborsNeighbors;
      neighborsNeighbors.insert(neighbors[i]);
      for (int j = 0; j < numNeighbors; j++) {
        if (i == j) continue;
        if (graph.isEdge(neighbors[i], neighbors[j])) {
          neighborsNeighbors.insert(neighbors[j]);
        }
      }
      neigborsMap[neighbors[i]] = neighborsNeighbors;
    }
    bool foundLonely = false;
    NodeId lonelyNodeId = "";
    for (int i = 0; i < numNeighbors ; i++) {
      if (neigborsMap[neighbors[i]].size() < numNeighbors-1) {
        if (foundLonely && neighbors[i] != lonelyNodeId) {
          return false; // can have only 1 lonely node
        } else {
          foundLonely = true;
          lonelyNodeId = neighbors[i];
        }
      } else if (neigborsMap[neighbors[i]].size() == numNeighbors-1) {
        for (int j = 0; j < numNeighbors; j++) {
          // If node is missing from list of neighbors
          if (neigborsMap[neighbors[i]].find(neighbors[j]) ==
              neigborsMap[neighbors[i]].end()) {
            if (foundLonely && neighbors[j] != lonelyNodeId) {
              return false; // can have only 1 lonely node
            } else {
              foundLonely = true;
              lonelyNodeId = neighbors[j];
            }
          }
        }
      }
    }
    return true;
  }

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

  void core(unsigned int k) {
    graph.core(k, graph, false, &cout);
  }

  // Returns a clique containing nodeId or an empty vector
  vector<NodeId> getClique(NodeId nodeId) {
    vector<NodeId> nodes;
    if (isSimplicial(nodeId)) {
      nodes = getNeighbors(nodeId);
      nodes.push_back(nodeId);
      return nodes;
    }
    //cout << "Did not find clique for " << nodeId << endl;
    return nodes;
  }

  void reduce() {
    unsigned int prevSize = 0;
    while (numNodes() != prevSize) {
      prevSize = numNodes();

      cout << "running islet rule...          ";
      runRule(ISLET_RULE);
      cout << " num nodes: " << numNodes()
           << "\tlow: " << getLow() << endl;

      cout << "running twig rule...           ";
      runRule(TWIG_RULE);
      cout << " num nodes: " << numNodes()
           << "\tlow: " << getLow() << endl;
      if (numNodes() != prevSize) continue;

      cout << "running series rule...         ";
      runRule(SERIES_RULE);
      cout << " num nodes: " << numNodes()
           << "\tlow: " << getLow() << endl;
      if (numNodes() != prevSize) continue;

      cout << "running triangleRule...        ";
      runRule(TRIANGLE_RULE);
      cout << " num nodes: " << numNodes()
           << "\tlow: " << getLow() << endl;
      if (numNodes() != prevSize) continue;

      cout << "running simplicialRule...      ";
      runRule(SIMPLICIAL_RULE);
      cout << " num nodes: " << numNodes()
           << "\tlow: " << getLow() << endl;
      if (numNodes() != prevSize) continue;

      cout << "running almostSimplicialRule...";
      runRule(ALMOST_SIMPLICIAL_RULE);
      cout << " num nodes: " << numNodes()
           << "\tlow: " << getLow() << endl;

    }
  }

  int lowerBound() {
    reduce();
    return low;
  }

};

#endif