Merged Generalized Loop Correction code kindly provided by Siamak Ravanbakhsh
[libdai.git] / include / dai / glc.h
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 /// \file
10 /// \brief Defines classes GLC and Cobweb, which implement the "Generalized Loop Correction method"
11 /// \todo Fix the init of GLC
12
13
14 #ifndef __defined_libdai_glc_h
15 #define __defined_libdai_glc_h
16
17
18 #include <algorithm>
19 #include <set>
20 #include <iostream>
21 #include <string>
22 #include <dai/util.h>
23 #include <dai/daialg.h>
24 #include <dai/enum.h>
25 #include <dai/cobwebgraph.h>
26 #include <dai/properties.h>
27 #include <dai/exceptions.h>
28
29
30 namespace dai {
31
32
33 /// Represents a region for the GLC algorithm.
34 /** This region contains all the proper factors from the factor graph. It also contains a factor for incoming messages and the
35 * cavity distribution over this region. The outgoing messages also receive a factor mainly so that we can easily calculate the marginal
36 * of the region over the outgoing message variables.
37 *
38 * \author Siamak Ravanbakhsh
39 */
40 class Cobweb {
41 public:
42 /// Pointer to the approximate inference algorithm and factor graph over region R (provided in parameters of GLC)
43 DAIAlgFG *cav;
44
45 /// Global to local factor index conversion.
46 /** Factor(I) in original graph corresponds to cav->factor(_g2l[I]) in this region
47 * the j'th incoming message to this region M(R,j) corresponds to cav->factor(_g2l[-j-1])
48 * the j'th cavity factor (for pairwise case),_cavitydists[R][j] corresponds to cav->factor(_g2l[-m-j-1])
49 * when m = M(R).size()
50 * note that for full cavity the _cavitydists[R] has a single factor and for uniform cavity it is empty
51 * the j'th outgoing message from this region outM(R,j) corresponds to cav->factor(_g2l[-c-m-j-1])
52 * when m = M(R).size() and c = _cavitydists[R].size()
53 */
54 std::map<int, size_t> _g2l;
55
56 /// Default constructor
57 Cobweb(){}
58
59 /// Construct from global-to-local factor index conversion structure
60 Cobweb( const std::map<int, size_t>& g2l ): cav(0), _g2l(g2l) {}
61
62 /// Default destructor
63 virtual ~Cobweb(){ delete cav; }
64
65 /// Sets the factor graph and inference algorithm for this region
66 void setInfAlg( DAIAlgFG* alg ) {
67 cav = alg;
68 }
69
70 /// Returns the factor of this region globaly indexed by \a I (mainly used to access messages to and from this region)
71 Factor factor( size_t I ) {
72 return cav->factor( _g2l[I] );
73 }
74
75 /// Initializes the inference algorithm over this region
76 void initialize() {
77 cav->init();
78 }
79
80 /// Calculates the marginal over \a ns after temporarily removing the factors indexed by \a rmgind
81 Factor marginal( const VarSet& ns, const std::vector<size_t>& rmgind ) {
82 // Set the factors in rmgind to one
83 VarSet vs;
84 std::vector<size_t> rmlinds; // local index
85 rmlinds.reserve( rmgind.size() );
86 for( std::vector<size_t>::const_iterator it = rmgind.begin(); it != rmgind.end(); it++ ) {
87 vs |= cav->factor( _g2l[*it] ).vars();
88 rmlinds.push_back( _g2l[*it] );
89 }
90 cav->makeRegionCavity( rmlinds, true );
91 // Initialize if necessary
92 cav->init( vs );
93 cav->run();
94 Factor result;
95 result = cav->belief( ns );
96 cav->restoreFactors();
97 return result;
98 }
99
100 /// Calculates the marginal over the variables in outgoing message indexed by \a outmsgind after temporarily removing the factors indexed by \a rmgind
101 /** This function is used to calculate outgoing messages from this region
102 */
103 Factor marginal( const int& outmsgind, const std::vector<size_t>& rmgind ) {
104 // set the factors in rmgind to one
105 VarSet vs;
106 std::vector<size_t> rmlinds; // local index
107 rmlinds.reserve( rmgind.size() );
108 for( std::vector<size_t>::const_iterator it = rmgind.begin(); it != rmgind.end(); it++ ) {
109 vs |= cav->factor( _g2l[*it] ).vars();
110 rmlinds.push_back( _g2l[*it] );
111 }
112 cav->makeRegionCavity( rmlinds, true );
113 // Initialize if necessary
114 cav->init( vs );
115 cav->run();
116 Factor result;
117 result = cav->beliefF( _g2l[outmsgind] );
118 cav->restoreFactors();
119 return result;
120 }
121
122 /// Sets (or multiplies, if \a multiply == \c true) the factor indexed by \a gind by factor \a msg
123 void updateFactor( int gind, const Factor& msg, bool multiply=false ) {
124 if( !multiply )
125 cav->setFactor( _g2l[gind], msg, false );
126 else
127 cav->setFactor( _g2l[gind], (msg*(cav->factor( _g2l[gind] ))).normalized(), false );
128 }
129
130 /// Runs inference and returns belief of variable \a v
131 Factor belief( const Var v ) {
132 cav->run();
133 return cav->belief( v );
134 }
135
136 /// Runs inference and returns belief of variables \a vs
137 Factor belief( const VarSet& vs ) {
138 cav->run();
139 return cav->belief( vs );
140 }
141 };
142
143
144 /// Approximate inference algorithm "Generalized Loop Correction" [\ref RYG12]
145 /** Inherits from a CobwebGraph which in turn inherits from a FactorGraph.
146 * Although a CobwebGraph represents the graph with cobweb clusters (for design reasons?),
147 * the actual vector of cobwebs is contained in GLC class. All the other functionalities
148 * needed to represent the cobwebgraph is in its own class. The cobweb provides the function
149 * of marginalization using the provided inference method. This can be beneficial when using
150 * uniform cavities (or a subset of pairwise cavities, but this is not implemented yet).
151 * However, in general the performance does not improve using Junction Tree instead of exact marginalization...
152 *
153 * \author Siamak Ravanbakhsh
154 */
155 class GLC : public DAIAlgCG {
156 private:
157 /// Stores for each variable the approximate cavity distribution
158 std::vector<Cobweb> _CWs;
159
160 /// Gives the starting global index for the factors of outgoing messages of each region
161 std::vector<int> _outOffset;
162
163 /// Cavity distributions
164 /** For each region it contains:
165 * - Nothing if using UniCAV
166 * - One Factor if using FullCAV
167 * - All pairwise caivity factors if using PairCAV or Pair2CAV
168 */
169 std::vector<std::vector<Factor> > _cavitydists;
170
171 /// Single variable beliefs
172 std::vector<Factor> _beliefs;
173
174 /// Beliefs over factors
175 std::map<VarSet, Factor> _factorBeliefs;
176
177 /// Maximum difference encountered so far
178 Real _maxdiff;
179
180 /// Number of iterations needed
181 size_t _iters;
182
183 public:
184 /// Parameters for GLC
185 struct Properties {
186 /// Enumeration of possible ways to initialize the cavities
187 /** The following initialization methods are defined:
188 * - FULL calculates the marginal using calcMarginal()
189 * - PAIR calculates only second order interactions using calcPairBeliefs() with \a accurate == \c false
190 * - PAIR2 calculates only second order interactions using calcPairBeliefs() with \a accurate == \c true
191 * - UNIFORM uses a uniform distribution
192 */
193 DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM);
194
195 /// Enumeration of different possible types of region (\ominus r in the paper)
196 /** The following initialization methods are defined:
197 * - SINGLE using single variable
198 * - FACTOR partitioning regions of that form corresponding to GLC (in the paper)
199 * - LOOP partitioning regions of that form corresponding to GLC (in the paper)
200 * - DELTA partitioning regions of that form corresponding to GLC (in the paper)
201 * - OVFACTOR, OVLOOP, OVDELTA the corresponding form of region in the overlapping form corresponding to GLC+ (in the paper)
202 * For example, using OVFACTOR each factor becomes a region (\ominus r), while using FACTOR the regions contain only
203 * non-overlapping factors and if some variables can't be covered by any remaining factors single variable regions
204 * are constructed for them.
205 * With this convention OVDELTA here corresponds to DELTA using CVM (HAK).
206 */
207 DAI_ENUM(RegionType,SINGLE,FACTOR,OVFACTOR,LOOP,OVLOOP,DELTA,OVDELTA)
208
209 /// Enumeration of different update schedules
210 /** The following update schedules are defined:
211 * - SEQFIX sequential fixed schedule
212 * - SEQRND sequential random schedule
213 */
214 DAI_ENUM(UpdateType,SEQFIX,SEQRND);
215
216 /// Verbosity (amount of output sent to stderr)
217 size_t verbose;
218
219 // Maximum length of the loops to use in construction of loopy regions
220 size_t loopdepth;
221
222 /// Maximum number of iterations
223 size_t maxiter;
224
225 /// The algorithm will quit and report the current result at reaching this runtime
226 Real maxtime;
227
228 /// Tolerance for convergence test
229 Real tol;
230
231 /// Complete or partial reinitialization of the messages in calculating the cavity distribution?
232 bool reinit;
233
234 /// Whether or not to include auxiliary factors in a region.
235 /** The auxilary factors are the ones that only depend on variables on the boundary of a region (\oplus r - \ominus r).
236 * This sometimes produces orders of magnitude improvement (e.g., in grids), but it is not discussed in the paper.
237 */
238 bool auxfactors;
239
240 /// How to initialize the cavities
241 CavityType cavity;
242
243 // Type of region
244 RegionType rgntype;
245
246 /// What update schedule to use
247 UpdateType updates;
248
249 /** Indicates what happens if a subset of variables in the boundary of a region (\ominus r_p) is shared by
250 * some neighbors such that one (\ominus r_{p,q1}) is a subset of another (\ominus r_{p,q2}).
251 * - ALL all such neighbors are included in the the updates.
252 * - TOP only (\ominus r_{p,q2}) is included unless (\ominus r_{p,q2} = \ominus r_{p,q1}) in which case both are included
253 * - CLOSEST (default): similar to TOP but in case of a tie the the region r_q with largest r_q \cap r_p is considered
254 * \note Not important in perfomance!
255 */
256 CobwebGraph::NeighborType neighbors;
257
258 /// Name of the algorithm used to initialize the cavity distributions
259 std::string cavainame;
260
261 /// Parameters for the algorithm used to initialize the cavity distributions
262 PropertySet cavaiopts;
263
264 // Name of the algorithm used for marginalization over a region (EXACT is fine)
265 std::string inainame;
266
267 // Parameters for algorithm used for marginalization over a region
268 PropertySet inaiopts;
269 } props;
270
271
272 public:
273 /// Default constructor
274 GLC() : DAIAlgCG(), _CWs(),_outOffset(), _cavitydists(), _beliefs(), _factorBeliefs(), _maxdiff(), _iters(), props() {}
275
276 /// Construct from FactorGraph \a fg and PropertySet \a opts
277 /** \param fg Factor graph.
278 * \param opts Parameters @see Properties
279 */
280 GLC( const FactorGraph &fg, const PropertySet &opts );
281
282 /// \name General InfAlg interface
283 //@{
284 virtual GLC* clone() const { return new GLC(*this); }
285 virtual GLC* construct( const FactorGraph &fg, const PropertySet &opts ) const { return new GLC( fg, opts ); }
286 virtual std::string name() const { return "GLC"; }
287 virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
288 virtual Factor belief( const VarSet &vs ) const;
289 virtual Factor beliefV( size_t i ) const { return _beliefs[i]; }
290 virtual std::vector<Factor> beliefs() const;
291 virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
292 virtual void init(){ initCWs(); }
293 virtual void init( const VarSet &/*ns*/ ) { init(); }
294 virtual Real run();
295 virtual Real maxDiff() const { return _maxdiff; }
296 virtual size_t Iterations() const { return _iters; }
297 virtual void setMaxIter( size_t maxiter ) { props.maxiter = maxiter; }
298 virtual void setProperties( const PropertySet &opts );
299 virtual PropertySet getProperties() const;
300 virtual std::string printProperties() const;
301 //@}
302
303 /// \name Additional interface specific for GLC
304 //@{
305 /// Returns constant reference to outer region \a alpha
306 const Cobweb& CW( size_t alpha ) const {
307 DAI_DEBASSERT( alpha < nrCWs() );
308 return _CWs[alpha];
309 }
310
311 /// Returns reference to outer region \a alpha
312 Cobweb& CW( size_t alpha ) {
313 DAI_DEBASSERT( alpha < nrCWs() );
314 return _CWs[alpha];
315 }
316
317 /// Approximates the cavity distribution of variable \a i, using the inference algorithm \a name with parameters \a opts
318 Real CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts );
319
320 /// Approximates all cavity distributions using inference algorithm \a name with parameters \a opts
321 Real InitCavityDists( const std::string &name, const PropertySet &opts );
322
323 /// Updates the belief of the Markov blanket of region \a R using its intersection with \a _R2'th neighbor (partitioning case).
324 void NewPancake( size_t R, size_t _R2 );
325
326 /// Updates the belief of the Markov blanket of region \a R using all its overlapping neighbors (non-partitioning case).
327 /** \note The implementation uses a more compact form than the message passing over neighborhood region-graph explained in the paper.
328 */
329 void OVNewPancake( size_t R );
330
331 /// Calculates the belief of variable \a i: if \a isFinal == \c true, it returns the geometric average given by all regions that include \a i
332 void CalcBelief( size_t i, bool isFinal = false );
333
334 /// Calculates the belief of factor \a I
335 void CalcFactorBelief( size_t I );
336
337 /// Designates the variables in each region of the graph based on the rgntype
338 std::vector<SmallSet<size_t> > calcRegions();
339
340 /// Initializes the cobweb regions and all the messages between regions.
341 void initCWs();
342
343 /// Initializes all the cobweb regions.
344 /** Defines their factors using the cavity distribution (as factor(s))
345 * and all the anticipated messages to be exchanged with other regions (treating a message as a factor),
346 * and also it builds a mapping from cobweb factors to factors of the original factor graph and messages.
347 */
348 void setCWs( const std::string &name, const PropertySet &opts );
349
350 /// Helper function for partitioning case to build the regions by finding loops of maximum given size.
351 /** If a variable is contained in a loop it is not considered for inclusion in any other loop.
352 */
353 void findLoopClusters( SmallSet<size_t>& remaining, std::set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars );
354
355 /// Helper function for general case to build the regions by finding loops of maximum given size.
356 /** All loops of max-size are found and variables that are not covered are contained in smaller loops or in factor-sized or single-variable regions.
357 */
358 void findOVLoopClusters( SmallSet<size_t>& remaining, std::set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars );
359 //@}
360 };
361
362
363 } // end of namespace dai
364
365
366 #endif