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