Changed license from GPL v2+ to FreeBSD (aka BSD 2-clause) license
[libdai.git] / include / dai / weightedgraph.h
index 9423059..edb1422 100644 (file)
@@ -1,22 +1,15 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
-    Radboud University Nijmegen, The Netherlands
-    
-    This file is part of libDAI.
-
-    libDAI 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 2 of the License, or
-    (at your option) any later version.
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
 
-    libDAI 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 libDAI; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+/** \file
+ *  \brief Defines some utility functions for (weighted) undirected graphs, trees and rooted trees.
+ *  \todo Improve general support for graphs and trees (in particular, a good tree implementation is needed).
+ */
 
 
 #ifndef __defined_libdai_weightedgraph_h
 #include <map>
 #include <iostream>
 #include <set>
-#include <cassert>
 #include <limits>
+#include <climits>   // Work-around for bug in boost graph library
+#include <dai/util.h>
+#include <dai/exceptions.h>
+#include <dai/graph.h>
 
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/prim_minimum_spanning_tree.hpp>
-
+#include <boost/graph/kruskal_min_spanning_tree.hpp>
 
 
 namespace dai {
 
 
-/// Directed edge
+/// Represents a directed edge
 class DEdge {
     public:
-        size_t  n1, n2;
-    
-        DEdge() {}
-        DEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
-        bool operator==( const DEdge &x ) const {
-            return ((n1 == x.n1) && (n2 == x.n2));
-        }
-        bool operator!=( const DEdge &x ) const {
-            return !(*this == x);
-        }
+        /// First node index (source of edge)
+        size_t first;
+
+        /// Second node index (target of edge)
+        size_t second;
+
+        /// Default constructor
+        DEdge() : first(0), second(0) {}
+
+        /// Constructs a directed edge pointing from \a m1 to \a m2
+        DEdge( size_t m1, size_t m2 ) : first(m1), second(m2) {}
+
+        /// Tests for equality
+        bool operator==( const DEdge &x ) const { return ((first == x.first) && (second == x.second)); }
+
+        /// Smaller-than operator (performs lexicographical comparison)
         bool operator<( const DEdge &x ) const {
-            return( (n1 < x.n1) || ((n1 == x.n1) && (n2 < x.n2)) );
+            return( (first < x.first) || ((first == x.first) && (second < x.second)) );
         }
+
+        /// Writes a directed edge to an output stream
         friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
-            os << "(" << e.n1 << "," << e.n2 << ")";
+            os << "(" << e.first << "->" << e.second << ")";
             return os;
         }
 };
 
 
-/// Undirected edge
+/// Represents an undirected edge
 class UEdge {
     public:
-        size_t  n1, n2;
-    
-        UEdge() {}
-        UEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
-        UEdge( const DEdge & e ) : n1(e.n1), n2(e.n2) {}
+        /// First node index
+        size_t first;
+
+        /// Second node index
+        size_t second;
+
+        /// Default constructor
+        UEdge() : first(0), second(0) {}
+
+        /// Constructs an undirected edge between \a m1 and \a m2
+        UEdge( size_t m1, size_t m2 ) : first(m1), second(m2) {}
+
+        /// Construct from DEdge
+        UEdge( const DEdge &e ) : first(e.first), second(e.second) {}
+
+        /// Tests for inequality (disregarding the ordering of the nodes)
         bool operator==( const UEdge &x ) {
-            return ((n1 == x.n1) && (n2 == x.n2)) || ((n1 == x.n2) && (n2 == x.n1));
+            return ((first == x.first) && (second == x.second)) || ((first == x.second) && (second == x.first));
         }
+
+        /// Smaller-than operator
         bool operator<( const UEdge &x ) const {
-            size_t s = n1, l = n2;
-            if( s > l )
-                std::swap( s, l );
-            size_t xs = x.n1, xl = x.n2;
-            if( xs > xl )
-                std::swap( xs, xl );
+            size_t s = std::min( first, second );
+            size_t l = std::max( first, second );
+            size_t xs = std::min( x.first, x.second );
+            size_t xl = std::max( x.first, x.second );
             return( (s < xs) || ((s == xs) && (l < xl)) );
         }
+
+        /// Writes an undirected edge to an output stream
         friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
-            if( e.n1 < e.n2 )
-                os << "{" << e.n1 << "," << e.n2 << "}";
+            if( e.first < e.second )
+                os << "{" << e.first << "--" << e.second << "}";
             else
-                os << "{" << e.n2 << "," << e.n1 << "}";
+                os << "{" << e.second << "--" << e.first << "}";
             return os;
         }
 };
 
 
-typedef std::vector<UEdge>  UEdgeVec;
-typedef std::vector<DEdge>  DEdgeVec;
-std::ostream & operator << (std::ostream & os, const DEdgeVec & rt);
+/// Represents an undirected graph, implemented as a std::set of undirected edges
+class GraphEL : public std::set<UEdge> {
+    public:
+        /// Default constructor
+        GraphEL() {}
+
+        /// Construct from range of objects that can be cast to UEdge
+        template <class InputIterator>
+        GraphEL( InputIterator begin, InputIterator end ) {
+            insert( begin, end );
+        }
+
+        /// Construct from GraphAL
+        GraphEL( const GraphAL& G ) {
+            for( size_t n1 = 0; n1 < G.nrNodes(); n1++ )
+                foreach( const Neighbor n2, G.nb(n1) )
+                    if( n1 < n2 )
+                        insert( UEdge( n1, n2 ) );
+        }
+};
+
+
+/// Represents an undirected weighted graph, with weights of type \a T, implemented as a std::map mapping undirected edges to weights
 template<class T> class WeightedGraph : public std::map<UEdge, T> {};
-typedef std::set<UEdge>     Graph;
 
 
-/// Use Prim's algorithm to construct a minimal spanning tree from the weighted graph Graph
-/// The weights should be positive. Use implementation in Boost Graph Library.
-template<typename T> DEdgeVec MinSpanningTreePrims( const WeightedGraph<T> & Graph ) {
-    DEdgeVec result;
-    if( Graph.size() > 0 ) {
+/// Represents a rooted tree, implemented as a vector of directed edges
+/** By convention, the edges are stored such that they point away from 
+ *  the root and such that edges nearer to the root come before edges
+ *  farther away from the root.
+ */
+class RootedTree : public std::vector<DEdge> {
+    public:
+        /// Default constructor
+        RootedTree() {}
+
+        /// Constructs a rooted tree from a tree and a root
+        /** \pre T has no cycles and contains node \a Root
+         */
+        RootedTree( const GraphEL &T, size_t Root );
+};
+
+
+/// Constructs a minimum spanning tree from the (non-negatively) weighted graph \a G.
+/** \param G Weighted graph that should have non-negative weights.
+ *  \param usePrim If true, use Prim's algorithm (complexity O(E log(V))), otherwise, use Kruskal's algorithm (complexity O(E log(E))).
+ *  \note Uses implementation from Boost Graph Library.
+ *  \note The vertices of \a G must be in the range [0,N) where N is the number of vertices of \a G.
+ */
+template<typename T> RootedTree MinSpanningTree( const WeightedGraph<T> &G, bool usePrim ) {
+    RootedTree result;
+    if( G.size() > 0 ) {
         using namespace boost;
         using namespace std;
-        typedef adjacency_list< vecS, vecS, undirectedS, property<vertex_distance_t, int>, property<edge_weight_t, double> > boostGraph;
-        typedef pair<size_t, size_t> E;
+        typedef adjacency_list< vecS, vecS, undirectedS, no_property, property<edge_weight_t, double> > boostGraph;
 
         set<size_t> nodes;
-        vector<E> edges;
+        vector<UEdge> edges;
         vector<double> weights;
-        edges.reserve( Graph.size() );
-        weights.reserve( Graph.size() );
-        for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
+        edges.reserve( G.size() );
+        weights.reserve( G.size() );
+        for( typename WeightedGraph<T>::const_iterator e = G.begin(); e != G.end(); e++ ) {
             weights.push_back( e->second );
-            edges.push_back( E( e->first.n1, e->first.n2 ) );
-            nodes.insert( e->first.n1 );
-            nodes.insert( e->first.n2 );
+            edges.push_back( e->first );
+            nodes.insert( e->first.first );
+            nodes.insert( e->first.second );
         }
 
+        size_t N = nodes.size();
+        for( set<size_t>::const_iterator it = nodes.begin(); it != nodes.end(); it++ )
+            if( *it >= N )
+                DAI_THROWE(RUNTIME_ERROR,"Vertices must be in range [0..N) where N is the number of vertices.");
+
         boostGraph g( edges.begin(), edges.end(), weights.begin(), nodes.size() );
-        vector< graph_traits< boostGraph >::vertex_descriptor > p( num_vertices(g) );
-        prim_minimum_spanning_tree( g, &(p[0]) );
-
-        // Store tree edges in result
-        result.reserve( nodes.size() - 1 );
-        size_t root = 0;
-        for( size_t i = 0; i != p.size(); i++ )
-            if( p[i] != i )
-                result.push_back( DEdge( p[i], i ) );
-            else
-                root = i;
-
-        // We have to store the minimum spanning tree in the right
-        // order, such that for all (i1, j1), (i2, j2) in result,
-        // if j1 == i2 then (i1, j1) comes before (i2, j2) in result.
-        // We do this by reordering the contents of result, effectively
-        // growing the tree starting at the root. At each step, 
-        // result[0..N-1] are the edges already added to the tree,
-        // whereas the other elements of result still have to be added.
-        // The elements of nodes are the vertices that still have to
-        // be added to the tree.
-
-        // Start with the root
-        nodes.erase( root );
-        size_t N = 0;
-
-        // Iteratively add edges and nodes to the growing tree
-        while( N != result.size() ) {
-            for( size_t e = N; e != result.size(); e++ ) {
-                bool e1_in_tree = !nodes.count( result[e].n1 );
-                if( e1_in_tree ) {
-                    nodes.erase( result[e].n2 );
-                    swap( result[N], result[e] );
-                    N++;
-                    break;
-                }
+        size_t root = *(nodes.begin());
+        GraphEL tree;
+        if( usePrim ) {
+            // Prim's algorithm
+            vector< graph_traits< boostGraph >::vertex_descriptor > p(N);
+            prim_minimum_spanning_tree( g, &(p[0]) );
+
+            // Store tree edges in result
+            for( size_t i = 0; i != p.size(); i++ ) {
+                if( p[i] != i )
+                    tree.insert( UEdge( p[i], i ) );
+            }
+        } else {
+            // Kruskal's algorithm
+            vector< graph_traits< boostGraph >::edge_descriptor > t;
+            t.reserve(  N - 1 );
+            kruskal_minimum_spanning_tree( g, std::back_inserter(t) );
+
+            // Store tree edges in result
+            for( size_t i = 0; i != t.size(); i++ ) {
+                size_t v1 = source( t[i], g );
+                size_t v2 = target( t[i], g );
+                if( v1 != v2 )
+                    tree.insert( UEdge( v1, v2 ) );
             }
         }
-    }
 
+        // Direct edges in order to obtain a rooted tree
+        result = RootedTree( tree, root );
+    }
     return result;
 }
 
 
-/// Use Prim's algorithm to construct a maximal spanning tree from the weighted graph Graph
-/// The weights should be positive. Use implementation in Boost Graph Library.
-template<typename T> DEdgeVec MaxSpanningTreePrims( const WeightedGraph<T> & Graph ) {
-    T maxweight = Graph.begin()->second;
-    for( typename WeightedGraph<T>::const_iterator it = Graph.begin(); it != Graph.end(); it++ )
-        if( it->second > maxweight )
-            maxweight = it->second;
-    // make a copy of the graph
-    WeightedGraph<T> gr( Graph );
-    // invoke MinSpanningTreePrims with negative weights
-    // (which have to be shifted to satisfy positivity criterion)
-    for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
-        it->second = maxweight - it->second;
-    return MinSpanningTreePrims( gr );
-}
-
-
-/// Calculate rooted tree from a tree T and a root
-DEdgeVec GrowRootedTree( const Graph & T, size_t Root );
-
-
-UEdgeVec RandomDRegularGraph( size_t N, size_t d );
-
-
-/// Use Dijkstra's algorithm to solve the shortest path problem for the weighted graph Graph and source vertex index s
-template<typename T>
-std::map<size_t,size_t> DijkstraShortestPaths( const WeightedGraph<T> & Graph, size_t s ) {
-/*
- * from wikipedia 
- *
-In the following algorithm, u := Extract_Min(Q) searches for the vertex u in the vertex set Q that has the least d[u] value. That vertex is removed from the set Q and returned to the user.
-
-  function Dijkstra(G, w, s)
-     for each vertex v in V[G]                        // Initializations
-           d[v] := infinity                                 // Unknown distance function from s to v
-           previous[v] := undefined
-     d[s] := 0                                        // Distance from s to s
-     S := empty set                                   // Set of all visited vertices
-     Q := V[G]                                        // Set of all unvisited vertices
-     while Q is not an empty set                      // The algorithm itself
-           u := Extract_Min(Q)                        // Remove best vertex from priority queue
-           S := S union {u}                           // Mark it 'visited'
-           for each edge (u,v) outgoing from u
-                  if d[u] + w(u,v) < d[v]             // Relax (u,v)
-                        d[v] := d[u] + w(u,v) 
-                        previous[v] := u
-
-To find the shortest path from s to t:
-
-u := t
-while defined previous[u]
-      insert u to the beginning of S
-      u := previous[u]
-       
-*/
-
-    // Calculate set of nodes in G
-    std::set<size_t> nodes;
-    for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
-        nodes.insert( e->n1.n1 );
-        nodes.insert( e->n1.n2 );
+/// Constructs a minimum spanning tree from the (non-negatively) weighted graph \a G.
+/** \param G Weighted graph that should have non-negative weights.
+ *  \param usePrim If true, use Prim's algorithm (complexity O(E log(V))), otherwise, use Kruskal's algorithm (complexity O(E log(E))).
+ *  \note Uses implementation from Boost Graph Library.
+ *  \note The vertices of \a G must be in the range [0,N) where N is the number of vertices of \a G.
+ */
+template<typename T> RootedTree MaxSpanningTree( const WeightedGraph<T> &G, bool usePrim ) {
+    if( G.size() == 0 )
+        return RootedTree();
+    else {
+        T maxweight = G.begin()->second;
+        for( typename WeightedGraph<T>::const_iterator it = G.begin(); it != G.end(); it++ )
+            if( it->second > maxweight )
+                maxweight = it->second;
+        // make a copy of the graph
+        WeightedGraph<T> gr( G );
+        // invoke MinSpanningTree with negative weights
+        // (which have to be shifted to satisfy positivity criterion)
+        for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
+            it->second = maxweight - it->second;
+        return MinSpanningTree( gr, usePrim );
     }
-    if( !nodes.count( s ) )
-        return std::map<size_t, size_t>();
-
-    std::map<size_t, double> d;
-    std::map<size_t, size_t> previous;
-    for( std::set<size_t>::const_iterator n = nodes.begin(); n != nodes.end(); n++ )
-        d[*n] = std::numeric_limits<double>::infinity();
-
-    d[s] = 0;
-    std::set<size_t> S;
-    std::set<size_t> Q = nodes;
-
-    while( Q.size() ) {
-        double least = d[*Q.begin()];
-        std::set<size_t>::iterator u_least = Q.begin();
-        for( std::set<size_t>::iterator _u = Q.begin(); _u != Q.end(); _u++ )
-            if( d[*_u] < least ) {
-                u_least = _u;
-                least = d[*_u];
-            }
-        size_t u = *u_least;
-        Q.erase( u_least );
-            
-        S.insert( u );
-        for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
-            size_t v;
-            if( e->n1.n1 == u )
-                v = e->n1.n2;
-            else if( e->n1.n2 == u )
-                v = e->n1.n1;
-            else
-                continue;
-            if( d[u] + e->n2 < d[v] ) {
-                d[v] = d[u] + e->n2;
-                previous[v] = u;
-            }
-        }
-    }
-
-    return previous;
 }