c8bfc80b00a91d0dbe284c823b81cafda23ba007
[libdai.git] / include / dai / mr.h
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, 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 class MR, which implements loop corrections as proposed by Montanari and Rizzo
11
12
13 #ifndef __defined_libdai_mr_h
14 #define __defined_libdai_mr_h
15
16
17 #include <vector>
18 #include <string>
19 #include <dai/factorgraph.h>
20 #include <dai/daialg.h>
21 #include <dai/enum.h>
22 #include <dai/properties.h>
23 #include <dai/exceptions.h>
24 #include <dai/graph.h>
25 #include <boost/dynamic_bitset.hpp>
26
27
28 namespace dai {
29
30
31 /// Approximate inference algorithm by Montanari and Rizzo [\ref MoR05]
32 /** \author Bastian Wemmenhove wrote the original implementation before it was merged into libDAI
33 */
34 class MR : public DAIAlgFG {
35 private:
36 /// Is the underlying factor graph supported?
37 bool supported;
38
39 /// The interaction graph (Markov graph)
40 GraphAL G;
41
42 /// tJ[i][_j] is the hyperbolic tangent of the interaction between spin \a i and its neighbour G.nb(i,_j)
43 std::vector<std::vector<Real> > tJ;
44 /// theta[i] is the local field on spin \a i
45 std::vector<Real> theta;
46
47 /// M[i][_j] is \f$ M^{(i)}_j \f$
48 std::vector<std::vector<Real> > M;
49 /// Cavity correlations
50 std::vector<std::vector<std::vector<Real> > > cors;
51
52 /// Type used for managing a subset of neighbors
53 typedef boost::dynamic_bitset<> sub_nb;
54
55 /// Magnetizations
56 std::vector<Real> Mag;
57
58 /// Maximum difference encountered so far
59 Real _maxdiff;
60
61 /// Number of iterations needed
62 size_t _iters;
63
64 public:
65 /// Parameters for MR
66 struct Properties {
67 /// Enumeration of different types of update equations
68 /** The possible update equations are:
69 * - FULL full updates, slow but accurate
70 * - LINEAR linearized updates, faster but less accurate
71 */
72 DAI_ENUM(UpdateType,FULL,LINEAR);
73
74 /// Enumeration of different ways of initializing the cavity correlations
75 /** The possible cavity initializations are:
76 * - RESPPROP using response propagation ("linear response")
77 * - CLAMPING using clamping and BP
78 * - EXACT using JunctionTree
79 */
80 DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
81
82 /// Verbosity (amount of output sent to stderr)
83 size_t verbose;
84
85 /// Tolerance for convergence test
86 Real tol;
87
88 /// Update equations
89 UpdateType updates;
90
91 /// How to initialize the cavity correlations
92 InitType inits;
93 } props;
94
95 public:
96 /// Default constructor
97 MR() : DAIAlgFG(), supported(), G(), tJ(), theta(), M(), cors(), Mag(), _maxdiff(), _iters(), props() {}
98
99 /// Construct from FactorGraph \a fg and PropertySet \a opts
100 /** \param fg Factor graph.
101 * \param opts Parameters @see Properties
102 * \note This implementation only deals with binary variables and pairwise interactions.
103 * \throw NOT_IMPLEMENTED if \a fg has factors depending on three or more variables or has variables with more than two possible states.
104 */
105 MR( const FactorGraph &fg, const PropertySet &opts );
106
107
108 /// \name General InfAlg interface
109 //@{
110 virtual MR* clone() const { return new MR(*this); }
111 virtual MR* construct( const FactorGraph &fg, const PropertySet &opts ) const { return new MR( fg, opts ); }
112 virtual std::string name() const { return "MR"; }
113 virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
114 virtual Factor belief( const VarSet &/*vs*/ ) const;
115 virtual Factor beliefV( size_t i ) const;
116 virtual std::vector<Factor> beliefs() const;
117 virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
118 virtual void init() {}
119 virtual void init( const VarSet &/*ns*/ ) { DAI_THROW(NOT_IMPLEMENTED); }
120 virtual Real run();
121 virtual Real maxDiff() const { return _maxdiff; }
122 virtual size_t Iterations() const { return _iters; }
123 virtual void setProperties( const PropertySet &opts );
124 virtual PropertySet getProperties() const;
125 virtual std::string printProperties() const;
126 //@}
127
128 private:
129 /// Initialize cors
130 Real calcCavityCorrelations();
131
132 /// Iterate update equations for cavity fields
133 void propagateCavityFields();
134
135 /// Calculate magnetizations
136 void calcMagnetizations();
137
138 /// Calculate the product of all tJ[i][_j] for _j in A
139 /** \param i variable index
140 * \param A subset of neighbors of variable \a i
141 */
142 Real _tJ(size_t i, sub_nb A);
143
144 /// Calculate \f$ \Omega^{(i)}_{j,l} \f$ as defined in [\ref MoR05] eqn. (2.15)
145 Real Omega(size_t i, size_t _j, size_t _l);
146
147 /// Calculate \f$ T^{(i)}_A \f$ as defined in [\ref MoR05] eqn. (2.17) with \f$ A = \{l_1,l_2,\dots\} \f$
148 /** \param i variable index
149 * \param A subset of neighbors of variable \a i
150 */
151 Real T(size_t i, sub_nb A);
152
153 /// Calculates \f$ T^{(i)}_j \f$ where \a j is the \a _j 'th neighbor of \a i
154 Real T(size_t i, size_t _j);
155
156 /// Calculates \f$ \Gamma^{(i)}_{j,l_1l_2} \f$ as defined in [\ref MoR05] eqn. (2.16)
157 Real Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
158
159 /// Calculates \f$ \Gamma^{(i)}_{l_1l_2} \f$ as defined in [\ref MoK07] on page 1141
160 Real Gamma(size_t i, size_t _l1, size_t _l2);
161
162 /// Approximates moments of variables in \a A
163 /** Calculate the moment of variables in \a A from M and cors, neglecting higher order cumulants,
164 * defined as the sum over all partitions of A into subsets of cardinality two at most of the
165 * product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
166 * entries of the partitions.
167 *
168 * \param i variable index
169 * \param A subset of neighbors of variable \a i
170 */
171 Real appM(size_t i, sub_nb A);
172
173 /// Calculate sum over all even/odd subsets B of \a A of _tJ(j,B) appM(j,B)
174 /** \param j variable index
175 * \param A subset of neighbors of variable \a j
176 * \param sum_even on return, will contain the sum over all even subsets
177 * \param sum_odd on return, will contain the sum over all odd subsets
178 */
179 void sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd);
180 };
181
182
183 } // end of namespace dai
184
185
186 #endif