Merge branch 'no-edges2'
[libdai.git] / src / weightedgraph.cpp
1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #include <algorithm>
23 #include <cassert>
24 #include <dai/weightedgraph.h>
25 #include <dai/util.h>
26
27
28 namespace dai {
29
30
31 using namespace std;
32
33
34 ostream & operator << (ostream & os, const DEdgeVec & rt) {
35 os << "[";
36 for( size_t n = 0; n < rt.size(); n++ )
37 os << rt[n] << ", ";
38 os << "]";
39 return os;
40 }
41
42
43 /// Calculate rooted tree from a tree T and a root
44 DEdgeVec GrowRootedTree( const Graph & T, size_t Root ) {
45 DEdgeVec result;
46 if( T.size() == 0 )
47 return result;
48 else {
49 // Make a copy
50 Graph Gr = T;
51
52 // Nodes in the tree
53 set<size_t> treeV;
54
55 // Start with the root
56 treeV.insert( Root );
57
58 // Keep adding edges until done
59 while( !(Gr.empty()) )
60 for( Graph::iterator e = Gr.begin(); e != Gr.end(); ) {
61 bool e1_in_treeV = treeV.count( e->n1 );
62 bool e2_in_treeV = treeV.count( e->n2 );
63 assert( !(e1_in_treeV && e2_in_treeV) );
64 if( e1_in_treeV ) {
65 // Add directed edge, pointing away from the root
66 result.push_back( DEdge( e->n1, e->n2 ) );
67 treeV.insert( e->n2 );
68 // Erase the edge
69 Gr.erase( e++ );
70 } else if( e2_in_treeV ) {
71 result.push_back( DEdge( e->n2, e->n1 ) );
72 treeV.insert( e->n1 );
73 // Erase the edge
74 Gr.erase( e++ );
75 } else
76 e++;
77 }
78
79 return result;
80 }
81 }
82
83
84 UEdgeVec RandomDRegularGraph( size_t N, size_t d ) {
85 // Algorithm 1 in "Generating random regular graphs quickly"
86 // by A. Steger and N.C. Wormald
87 //
88 // Draws a random graph with size N and uniform degree d
89 // from an almost uniform probability distribution over these graphs
90 // (which becomes uniform in the limit that d is small and N goes
91 // to infinity).
92
93 assert( (N * d) % 2 == 0 );
94
95 bool ready = false;
96 UEdgeVec G;
97
98 size_t tries = 0;
99 while( !ready ) {
100 tries++;
101
102 // Start with N*d points {0,1,...,N*d-1} (N*d even) in N groups.
103 // Put U = {0,1,...,N*d-1}. (U denotes the set of unpaired points.)
104 vector<size_t> U;
105 U.reserve( N * d );
106 for( size_t i = 0; i < N * d; i++ )
107 U.push_back( i );
108
109 // Repeat the following until no suitable pair can be found: Choose
110 // two random points i and j in U, and if they are suitable, pair
111 // i with j and delete i and j from U.
112 G.clear();
113 bool finished = false;
114 while( !finished ) {
115 random_shuffle( U.begin(), U.end() );
116 size_t i1, i2;
117 bool suit_pair_found = false;
118 for( i1 = 0; i1 < U.size()-1 && !suit_pair_found; i1++ )
119 for( i2 = i1+1; i2 < U.size() && !suit_pair_found; i2++ )
120 if( (U[i1] / d) != (U[i2] / d) ) {
121 // they are suitable
122 suit_pair_found = true;
123 G.push_back( UEdge( U[i1] / d, U[i2] / d ) );
124 U.erase( U.begin() + i2 ); // first remove largest
125 U.erase( U.begin() + i1 ); // remove smallest
126 }
127 if( !suit_pair_found || U.empty() )
128 finished = true;
129 }
130
131 if( U.empty() ) {
132 // G is a graph with edge from vertex r to vertex s if and only if
133 // there is a pair containing points in the r'th and s'th groups.
134 // If G is d-regular, output, otherwise return to Step 1.
135
136 vector<size_t> degrees;
137 degrees.resize( N, 0 );
138 for( UEdgeVec::const_iterator e = G.begin(); e != G.end(); e++ ) {
139 degrees[e->n1]++;
140 degrees[e->n2]++;
141 }
142 ready = true;
143 for( size_t n = 0; n < N; n++ )
144 if( degrees[n] != d ) {
145 ready = false;
146 break;
147 }
148 } else
149 ready = false;
150 }
151
152 return G;
153 }
154
155
156 } // end of namespace dai