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