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