dag.h

/**************************************************************************
 * Copyright (c) 2004-2011 T. M. Murali                                   *
 * Copyright (c) 2004 Greg Grothaus                                       *
 *                                                                        *
 * 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 _DAG_H
#define _DAG_H

#include<iostream>
#include<map>
#include<set>
#include<vector>
#include<iterator>
#include<set>
#include<string>
#include<sstream>
#include "setTemplates.C"
#include <queue>

using namespace std;

// forward declaration.
template< typename NodeType > class DAG;


template< typename NodeInfo > class DAGNode
{
  friend class DAG< DAGNode < NodeInfo > >;

private:
  NodeInfo _nodeInfo;

protected:
  set< DAGNode< NodeInfo > *> parents[3],ancestors[3];  //0-all, 1-is_a, 2-part_of      
  set< DAGNode< NodeInfo > *> children[3], descendants[3];      //0-all, 1-is_a, 2-part_of      

public:
  DAGNode()
      : _nodeInfo(), parents(), ancestors(), children(), descendants()
    {}
  DAGNode(const NodeInfo &info)
      : _nodeInfo(info), parents(), ancestors(), children(), descendants()
    {}

  DAGNode(const DAGNode& rhs)
      : _nodeInfo(rhs._nodeInfo), parents(), ancestors(), children(), descendants()
    {}


  virtual ~DAGNode()
    {}

  virtual void addChild(DAGNode< NodeInfo > *child, unsigned int type = 0)
    {
#ifdef DEBUG
//        cout << "DAGNode::addChild(): node is " << _nodeInfo.getId()
//             << ". Child is " << child->_nodeInfo.getId()
//             << " Edge type is " << type
//             << endl;
#endif
      // check for loops.
      if (isParent(child, type))
        cerr << "DAG::addChild(): adding child " << child->_nodeInfo.getId()
             << " for node " << _nodeInfo.getId() << " but child is already a parent!"
             << endl;
      children[type].insert(child);
      // also insert child in type 0 (which stores all children).
      children[0].insert(child);
      
    }
  
  virtual void addParent(DAGNode< NodeInfo > *parent, unsigned int type = 0)
    {
#ifdef DEBUG
//       cout << "DAGNode::addParent(): node is " << _nodeInfo.getId()
//            << ". Parent is " << parent->_nodeInfo.getId() << endl;
#endif      
      // check for loops.
      if (isChild(parent, type))
        cerr << "DAG::addParent(): adding parent " << parent->_nodeInfo.getId()
             << " for node " << _nodeInfo.getId() << " but parent is already a child!"
             << endl;
      parents[type].insert(parent);
      // also insert parent in type 0 (which stores all parents).
      parents[0].insert(parent);
    }
  
  virtual set< DAGNode< NodeInfo > *> getParents(unsigned int type = 0)
    {
      return(parents[type]);
    }

  virtual unsigned int getParentEdgeType(DAGNode< NodeInfo > *parent) 
    {
      for (unsigned int i = 1; i < 3; i++)
        if (isParent(parent, i))
          return(i);
      // i do not know any type for this parent.
      cerr << "DAG::getParentEdgeType(): No edge type known for node " << _nodeInfo.getId()
           << " and parent " << parent->_nodeInfo.getId() << endl;
      return(0);
    }
  
//  virtual
  set< DAGNode< NodeInfo > *> getChildren(unsigned int type = 0)
    {
      return(children[type]);
    }
  
  virtual set< DAGNode< NodeInfo > *> getAncestors(int type = 0)
    {
      
      //check if already computed/no parents
      if(ancestors[type].size()>0)
        return ancestors[type];
      ancestors[type]= parents[type];
      for(typename set< DAGNode< NodeInfo > *>::iterator i = parents[type].begin();
          i != parents[type].end(); i++)
        ancestors[type]= setunion((*i)->getAncestors(type), ancestors[type]);
  
      return ancestors[type];
    }

  virtual set< DAGNode< NodeInfo > *> getDescendants(int type = 0)
    {
      //check if already computed or no descendants.
      if(descendants[type].size()>0)
        return descendants[type];
      descendants[type] = children[type];
      for(typename set< DAGNode< NodeInfo > *>::iterator i = children[type].begin();
          i != children[type].end(); i++)
        descendants[type]= setunion((*i)->getDescendants(type), descendants[type]);
  
      return descendants[type];
    }


  virtual NodeInfo getNodeInfo() const
    {
      return(_nodeInfo);
    }
  
  
  virtual NodeInfo &getNodeInfo() 
    {
      return(_nodeInfo);
    }
  
  
  virtual bool isChild(DAGNode< NodeInfo > *node, unsigned int type = 0)
    {
      return(children[type].end() != children[type].find(node));
      
    }
  
  virtual bool isParent(DAGNode< NodeInfo > *node, unsigned int type = 0) 
    {
      return(parents[type].end() != parents[type].find(node));
      
    }
  
};

  


template< typename NodeType > class DAG
{

protected:

  vector< NodeType* > _nodes;
  // map telling me in which index of _nodes each node is stored.
  map< const NodeType*, unsigned int > _nodesToIndices;

protected:

  // copy rhs into *this.
  void _copy(const DAG< NodeType > &rhs);

  // copy the nodes of rhs into *this.
  void _copyNodes(const DAG< NodeType > &rhs);
  
  // methods to get the index from a NodeType*.
  unsigned int _getNodeIndex(const NodeType *node) const
    {
      return(_nodesToIndices.find(node)->second);
    }
  
  
  
public:
  DAG()
      : _nodes(), _nodesToIndices()
    {}

  DAG(const DAG< NodeType > &rhs)
    {
      _copy(rhs);
    }
  
  
  DAG(istream *in);

  virtual ~DAG()
    {
      clear();
    }


  const DAG< NodeType > &operator=(const DAG &rhs)
    {
      if (this != &rhs)
        {
          clear();
          _copy(rhs);
        }
      return(*this);
    }
  
  void clear()
    {
      for (unsigned int i = 0; i < _nodes.size(); i++)
        delete _nodes[i];
      _nodes.clear();
      _nodesToIndices.clear();
    }
  

  virtual void addNode(NodeType *node)
    {
#ifdef DEBUG
//       cout << "DAG::addNode(): node is " << node->_nodeInfo.getId()
//            << "." << endl;
#endif      
      _nodes.push_back(node);
      _nodesToIndices[node] = _nodes.size() - 1;
    }
  
//   virtual void addNode(const NodeType)
//     {
//       _nodes.push_back(node);
//     }
  

  void computeDifference(DAG< NodeType> &other, DAG< NodeType> &difference);

  
  virtual void computePower(unsigned int exponent, DAG< NodeType> &powerDAG);
  
  virtual void computeLeaves(vector< NodeType* > &leaves, unsigned int type = 0);

  virtual void computeRoots(vector< NodeType* > &roots, unsigned int type = 0);
  
  
  virtual void computeTopologicalSort(vector< NodeType* > &sortedFunctions,
                                      unsigned int type = 0);
  

  virtual void computeTransitiveClosure(DAG< NodeType> &closedDAG);

  virtual void computeTransitiveReduction(DAG< NodeType> &reducedDAG);

  

//  virtual void read(string file);

//  virtual void read(istream *in);
  
//  virtual set<NodeType*> getAncestors(NodeType* obj, int type=0);

//  virtual set<NodeType*> getParents(NodeType* obj, int type=0);

//  virtual set<NodeType*> getDescendants(NodeType* obj, int type=0);

//  virtual set<NodeType*> getChildren(NodeType* obj, int type=0);
  
//   NodeType* getFunctionById(int id) const;
  
//   NodeType* getFunctionById(string id) const
//     {
//       unsigned int intID;
//       stringstream ss;
//       ss << id;
//       ss >> intID;
//       return(getFunctionById(intID));
//     }
  
  unsigned int getNumEdges() const;
  
  unsigned int getNumNodes() const;
  
  void print(string file) const;
  
  void print(ostream &ostr) const;
  

};


// template< typename NodeType >
// set< NodeType *> DAG< NodeType >::getAncestors(NodeType * obj, int type)
// {
//   return(obj->getAncestors(type));
// }

// template< typename NodeType >
// set< NodeType* > DAG< NodeType >::getDescendants(NodeType* obj, int type)
// {
//   return(obj->getDescendants(type));
// }

// template< typename NodeType >
// set< NodeType* > DAG< NodeType >::getParents(NodeType* obj, int type)
// {
//   return obj->getParents(type);
// }

// template< typename NodeType >
// set< NodeType* > DAG< NodeType >::getChildren(NodeType* obj, int type)
// {
//   return obj->getChildren(type);
// }

template< typename NodeType >
void DAG< NodeType >::_copyNodes(const DAG< NodeType > &rhs)
{
  // add nodes from rhs to *this.
  for (unsigned int i = 0; i < rhs._nodes.size(); i++)
    {
      addNode(new NodeType(*rhs._nodes[i]));
//       _nodes.push_back(new NodeType(*rhs._nodes[i]));
//       _nodesToIndices[_nodes.back()] =  i;
    }
}

template< typename NodeType >
void DAG< NodeType >::_copy(const DAG< NodeType > &rhs)
{
  _copyNodes(rhs);
  unsigned int parentIndex, edgeType;
  NodeType *node, *parent;
  typename set< NodeType *>::iterator pitr;
  // loop over all edges of rhs and add them to *this.
  for (unsigned int i = 0; i < rhs._nodes.size(); i++)
    {
      node = _nodes[i];
      set< NodeType *> parents = rhs._nodes[i]->getParents();
      for (pitr = parents.begin(); pitr != parents.end(); pitr++)
        {
          parentIndex = rhs._getNodeIndex(*pitr);
          parent = _nodes[parentIndex];
          edgeType = rhs._nodes[i]->getParentEdgeType(*pitr);
          // add edge to *this.
          node->addParent(parent, edgeType);
          parent->addChild(node, edgeType);
        }
    }
}

template< typename NodeType >
unsigned int DAG< NodeType >::getNumEdges() const
{
  // loop over the nodes.
  typename vector< NodeType* >::const_iterator nitr;
  unsigned int numEdges = 0;
  for (nitr = _nodes.begin(); nitr != _nodes.end(); nitr++)
    {
      set< NodeType *> parents = (*nitr)->getParents();
      numEdges += parents.size();
    }
  return(numEdges);
}

template< typename NodeType >
unsigned int DAG< NodeType >::getNumNodes() const
{
  return(_nodes.size());
}

template< typename NodeType >
void DAG< NodeType >::print(string file) const
{
  ofstream fstr(file.c_str());
  print(fstr);
  fstr.close();
}

template< typename NodeType >
void DAG< NodeType >::print(ostream &ostr) const
{
  // loop over the nodes.
  typename vector< NodeType* >::const_iterator nitr;
  typename set< NodeType *>::const_iterator sitr;
  for (nitr = _nodes.begin(); nitr != _nodes.end(); nitr++)
    {
      set< NodeType *> parents = (*nitr)->getParents();
#ifdef DEBUG
//       cout << "Node " << (*nitr)->_nodeInfo.getId() << " has "
//            << parents.size() << " parents." << endl;
#endif
      for (sitr = parents.begin(); sitr != parents.end(); sitr++)
        ostr << (*nitr)->_nodeInfo.getId()
             << "\t" << (*sitr)->_nodeInfo.getId()
             << "\t" << (*nitr)->getParentEdgeType(*sitr)
             << endl;
    }
  
}

template< typename NodeType >
void DAG< NodeType >::computeDifference(DAG< NodeType> &other,
                                        DAG< NodeType> &difference)
{
  // let us add all the nodes in *this to difference.
  difference._copyNodes(*this);
//   for (unsigned int i = 0; i < _nodes.size(); i++)
//     {
//       difference._nodes.push_back(new NodeType(*_nodes[i]));
//       difference._nodesToIndices[difference._nodes.back()] =  i;
//     }
  

  typename set< NodeType * >::iterator pitr;
  NodeType *nodeInOther, *nodeInDifference, *parentInOther, *parentInDifference;
  unsigned int parentIndex, parentEdgeType;
  // loop over all edges in *this.
  for (unsigned int i = 0; i < _nodes.size(); i++)
    {
      nodeInOther = other._nodes[i];
      nodeInDifference = difference._nodes[i];
      set< NodeType *> parents = _nodes[i]->getParents();
      for (pitr = parents.begin(); pitr != parents.end(); pitr++)
        {
          parentIndex = _getNodeIndex(*pitr);
          parentInOther = other._nodes[parentIndex];
          parentEdgeType = _nodes[i]->getParentEdgeType(*pitr);
          if (nodeInOther->isParent(parentInOther))
            // this edge exists in other, so i ignore it.
            continue;
          // add edge between node and other to difference.
          
          // locate parent in difference.
          parentInDifference = difference._nodes[parentIndex];
          nodeInDifference->addParent(parentInDifference, parentEdgeType);
          parentInDifference->addChild(nodeInDifference, parentEdgeType);
        }
    }
  
}

template< typename NodeType >
void DAG< NodeType >::computeLeaves(vector< NodeType* > &leaves, unsigned int type)
{
  for (unsigned int i = 0; i < _nodes.size(); i++)
    {
      if (0 == _nodes[i]->getChildren(type).size())
        // _nodes[i] has no leaves.
        leaves.push_back(_nodes[i]);
    }
}

template< typename NodeType >
void DAG< NodeType >::computePower(unsigned int exponent, DAG< NodeType> &powerDAG)
{
  if (0 == exponent)
    return;
  if (1 == exponent)
    {
      powerDAG = *this;
      return;
    }
  DAG< NodeType > grandparentDAG;
  if (2 < exponent)
    // make a recursive call.
    {
      computePower(exponent - 1, grandparentDAG);
    }
  else
    grandparentDAG = *this;

  // invariant: at this stage, grandparentDAG = (*this)^{exponent - 1}.
  
  // now join *this and grandparentDAG to form powerDAG. 

  // first, add all nodes in *this to powerDAG.
  powerDAG._copyNodes(*this);
//   for (unsigned int i = 0; i < _nodes.size(); i++)
//     powerDAG._nodes.push_back(new NodeType(*_nodes[i]));

  // next, for each node in *this, find its parents in *this, find the
  // parents of the parents in grandparentDAG, and in powerDAG,
  // connect the node to each of its grandparents.
  NodeType *nodeInDAG, *parentInGrandparentDAG, *grandparentInPowerDAG;
  typename set< NodeType *>::iterator pitr, gitr;
  unsigned int grandparentIndex;
  
  for (unsigned int i = 0; i < _nodes.size(); i++)
    {
//      nodeInDAG = _nodes[i];
      set< NodeType *> parents = _nodes[i]->getParents();
      for (pitr = parents.begin(); pitr != parents.end(); pitr++)
        {
          // find the parent in grandparentDAG.
          parentInGrandparentDAG = grandparentDAG._nodes[_getNodeIndex(*pitr)];
          // find the parents of parentInGrandparentDAG.
          set< NodeType *> grandparents = parentInGrandparentDAG->getParents();
          for (gitr = grandparents.begin(); gitr != grandparents.end(); gitr++)
            {
              // in powerDAG, find the counterpart of *gitr. i use the
              // fact that node indices are identical in powerDAG and
              // grandparentDAG.
              grandparentIndex = grandparentDAG._getNodeIndex(*gitr);
              grandparentInPowerDAG = powerDAG._nodes[grandparentIndex];
              powerDAG._nodes[i]->addParent(grandparentInPowerDAG);
              grandparentInPowerDAG->addChild(powerDAG._nodes[i]);
            }
        }
    }
}


template< typename NodeType >
void DAG< NodeType >::computeRoots(vector< NodeType* > &roots, unsigned int type)
{
  for (unsigned int i = 0; i < _nodes.size(); i++)
    {
      if (0 == _nodes[i]->getParents(type).size())
        // _nodes[i] has no parents.
        roots.push_back(_nodes[i]);
    }
}

template< typename NodeType >
void DAG< NodeType >::computeTopologicalSort(vector< NodeType* > &sortedNodes,
                                             unsigned int type)
{
  map< NodeType *, unsigned int > numVisitedParents;
  set< NodeType * > roots;
  deque< NodeType * > nodeQueue;
  for (unsigned int i = 0; i < _nodes.size(); i++)
    {
      // initialise
      numVisitedParents[_nodes[i]] = 0;
      if (0 == _nodes[i]->getParents(type).size())
        // _nodes[i] has no parents.
        roots.insert(_nodes[i]);
//        nodeQueue.push_back(_nodes[i]);
    }
  NodeType *nextNode;
  // pop elements from the queue until it is empty.
  while (!nodeQueue.empty() || !roots.empty())
    {
      if (nodeQueue.empty())
        // insert a node from from roots and remove it from roots
        // so that it is not inserted. this way the nodes in each
        // component of the DAG appear consecutively in the sorted
        // order.
        {
          nodeQueue.push_back(*(roots.begin()));
          roots.erase(roots.begin());
        }
      
      nextNode = nodeQueue.front();
      nodeQueue.pop_front();
      sortedNodes.push_back(nextNode);
      // process children of nextNode.
      set< NodeType *> children = nextNode->getChildren(type);
      typename set< NodeType* >::const_iterator itr;
      
      for (itr = children.begin();
           itr != children.end(); itr++)
        {
          numVisitedParents[*itr]++;
          // push child into queue if all parents are visited.
          if (numVisitedParents[*itr] == (*itr)->getParents().size())
            nodeQueue.push_back(*itr);
        }
    }
}

template< typename NodeType >
void DAG< NodeType >::computeTransitiveClosure(DAG< NodeType> &closedDAG)
{
  cout << "DAG has " << getNumNodes() << " nodes and " << getNumEdges()
       << " edges." << endl;
  // i must operate on the nodes in closedDAG, so first make closedDAG
  // a copy of *this.
  closedDAG = *this;
  cout << "Copy of DAG has " << closedDAG.getNumNodes() << " nodes and "
       << closedDAG.getNumEdges() << " edges." << endl;
//   // add each node in *this to closedDAG.
//   for (unsigned int i = 0; i < _nodes.size(); i++)
//     closedDAG.addNode(_nodes[i]);

  unsigned int numClosureEdges = 0;
  typename set< NodeType *>::iterator aitr;
  // loop over each node in *this.
  for (unsigned int i = 0; i < closedDAG._nodes.size(); i++)
    {
      cout << "DAG::computeTransitiveClosure(): Processing " << closedDAG._nodes[i]->_nodeInfo.getId() << endl;
      // find ancestors of each node.
      set< NodeType *> ancestors = closedDAG._nodes[i]->getAncestors();
      for (aitr = ancestors.begin(); aitr != ancestors.end(); aitr++)
        // add link between ancestor and node to closedDAG.
        {
          cout << "\tChild " << closedDAG._nodes[i]->_nodeInfo.getId()
               << ", parent " << (*aitr)->_nodeInfo.getId() << endl;
          closedDAG._nodes[i]->addParent(*aitr);
          (*aitr)->addChild(closedDAG._nodes[i]);
          numClosureEdges++;
        }
    }
  cout << "DAG< NodeType >::computeTransitiveClosure(): added "
       << numClosureEdges << " edges to the transitive closure." << endl;
  cout << "Transitive closure of DAG has " << closedDAG.getNumNodes() << " nodes and "
       << closedDAG.getNumEdges() << " edges." << endl;
}



template< typename NodeType >
void DAG< NodeType >::computeTransitiveReduction(DAG< NodeType> &reducedDAG)
{
  cout << "DAG has " << getNumNodes() << " nodes and "
       << getNumEdges() << " edges." << endl;
  DAG< NodeType > closedDAG, squareClosedDAG;
  // compute transitive closure of *this.
//  computeTransitiveClosure(closedDAG);
  closedDAG = *this;
  cout << "Copy of DAG has "
       << closedDAG.getNumNodes() << " nodes and "
       << closedDAG.getNumEdges() << " edges." << endl;
  
  // compute square of transitive closure.
  closedDAG.computePower(2, squareClosedDAG);
  cout << "Square of transitive closure of DAG has "
       << squareClosedDAG.getNumNodes() << " nodes and "
       << squareClosedDAG.getNumEdges() << " edges." << endl;
  
  // reducedDAG consists of edges in transitive closure that are not in the square.
  closedDAG.computeDifference(squareClosedDAG, reducedDAG);
  cout << "Transitive reduction of DAG has "
       << reducedDAG.getNumNodes() << " nodes and "
       << reducedDAG.getNumEdges() << " edges." << endl;
}


#endif // _DAG_H