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