weightedgraph.h

Go to the documentation of this file.

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


#ifndef __defined_libdai_weightedgraph_h
#define __defined_libdai_weightedgraph_h


#include <vector>
#include <map>
#include <iostream>
#include <set>
#include <limits>
#include <climits>   // Work-around for bug in boost graph library

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

#include <dai/util.h>
#include <dai/exceptions.h>
#include <dai/graph.h>


namespace dai {


class DEdge {
    public:
        size_t first;

        size_t second;

        DEdge() : first(0), second(0) {}

        DEdge( size_t m1, size_t m2 ) : first(m1), second(m2) {}

        bool operator==( const DEdge &x ) const { return ((first == x.first) && (second == x.second)); }

        bool operator<( const DEdge &x ) const {
            return( (first < x.first) || ((first == x.first) && (second < x.second)) );
        }

        friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
            os << "(" << e.first << "->" << e.second << ")";
            return os;
        }

        std::string toString() const {
            std::stringstream ss;
            ss << *this;
            return ss.str();
        }
};


class UEdge {
    public:
        size_t first;

        size_t second;

        UEdge() : first(0), second(0) {}

        UEdge( size_t m1, size_t m2 ) : first(m1), second(m2) {}

        UEdge( const DEdge &e ) : first(e.first), second(e.second) {}

        bool operator==( const UEdge &x ) {
            return ((first == x.first) && (second == x.second)) || ((first == x.second) && (second == x.first));
        }

        bool operator<( const UEdge &x ) const {
            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)) );
        }

        friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
            if( e.first < e.second )
                os << "{" << e.first << "--" << e.second << "}";
            else
                os << "{" << e.second << "--" << e.first << "}";
            return os;
        }

        std::string toString() const {
            std::stringstream ss;
            ss << *this;
            return ss.str();
        }
};


class GraphEL : public std::set<UEdge> {
    public:
        GraphEL() {}

        template <class InputIterator>
        GraphEL( InputIterator begin, InputIterator end ) {
            insert( begin, end );
        }

        GraphEL( const GraphAL& G ) {
            for( size_t n1 = 0; n1 < G.nrNodes(); n1++ )
                bforeach( const Neighbor n2, G.nb(n1) )
                    if( n1 < n2 )
                        insert( UEdge( n1, n2 ) );
        }
};


template<class T> class WeightedGraph : public std::map<UEdge, T> {};



class RootedTree : public std::vector<DEdge> {
    public:
        RootedTree() {}


        RootedTree( const GraphEL &T, size_t Root );
};



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, no_property, property<edge_weight_t, double> > boostGraph;

        set<size_t> nodes;
        vector<UEdge> edges;
        vector<double> weights;
        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->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() );
        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;
}



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 );
    }
}


} // end of namespace dai


#endif