Fixed regression in TFactor::partSum
[libdai.git] / include / dai / mr.h
1 /* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
2 Radboud University Nijmegen, The Netherlands /
3 Max Planck Institute for Biological Cybernetics, Germany
4
5 This file is part of libDAI.
6
7 libDAI is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 libDAI is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with libDAI; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22
23 /// \file
24 /// \brief Defines class MR
25 /// \todo Improve documentation
26
27
28 #ifndef __defined_libdai_mr_h
29 #define __defined_libdai_mr_h
30
31
32 #include <vector>
33 #include <string>
34 #include <dai/factorgraph.h>
35 #include <dai/daialg.h>
36 #include <dai/enum.h>
37 #include <dai/properties.h>
38 #include <dai/exceptions.h>
39 #include <boost/dynamic_bitset.hpp>
40
41
42 namespace dai {
43
44
45 /// Approximate inference algorithm by Montanari and Rizzo
46 class MR : public DAIAlgFG {
47 private:
48 bool supported; // is the underlying factor graph supported?
49
50 std::vector<size_t> con; // con[i] = connectivity of spin i
51 std::vector<std::vector<size_t> > nb; // nb[i] are the neighbours of spin i
52 std::vector<std::vector<double> > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
53 std::vector<double> theta; // theta[i] is the local field on spin i
54 std::vector<std::vector<double> > M; // M[i][_j] is M^{(i)}_j
55 std::vector<std::vector<size_t> > kindex; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
56 std::vector<std::vector<std::vector<double> > > cors;
57
58 static const size_t kmax = 31;
59 typedef boost::dynamic_bitset<> sub_nb;
60
61 size_t N;
62
63 std::vector<double> Mag;
64
65 double _maxdiff;
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 DAI_ENUM(UpdateType,FULL,LINEAR)
73
74 /// Enumeration of different ways of initializing the cavity correlations
75 DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT)
76
77 /// Verbosity
78 size_t verbose;
79
80 /// Tolerance
81 double tol;
82
83 /// Update equations
84 UpdateType updates;
85
86 /// How to initialize the cavity correlations
87 InitType inits;
88 } props;
89
90 /// Name of this inference method
91 static const char *Name;
92
93 public:
94 /// Default constructor
95 MR() : DAIAlgFG(), supported(), con(), nb(), tJ(), theta(), M(), kindex(), cors(), N(), Mag(), _maxdiff(), _iters(), props() {}
96
97 /// Copy constructor
98 MR( const MR &x ) : DAIAlgFG(x), supported(x.supported), con(x.con), nb(x.nb), tJ(x.tJ), theta(x.theta), M(x.M), kindex(x.kindex), cors(x.cors), N(x.N), Mag(x.Mag), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
99
100 /// Assignment operator
101 MR& operator=( const MR &x ) {
102 if( this != &x ) {
103 DAIAlgFG::operator=(x);
104 supported = x.supported;
105 con = x.con;
106 nb = x.nb;
107 tJ = x.tJ;
108 theta = x.theta;
109 M = x.M;
110 kindex = x.kindex;
111 cors = x.cors;
112 N = x.N;
113 Mag = x.Mag;
114 _maxdiff = x._maxdiff;
115 _iters = x._iters;
116 props = x.props;
117 }
118 return *this;
119 }
120
121 /// Construct from FactorGraph fg and PropertySet opts
122 MR( const FactorGraph &fg, const PropertySet &opts );
123
124
125 /// @name General InfAlg interface
126 //@{
127 virtual MR* clone() const { return new MR(*this); }
128 virtual MR* create() const { return new MR(); }
129 virtual std::string identify() const;
130 virtual Factor belief( const Var &n ) const;
131 virtual Factor belief( const VarSet &/*ns*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
132 virtual std::vector<Factor> beliefs() const;
133 virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
134 virtual void init() {}
135 virtual void init( const VarSet &/*ns*/ ) { DAI_THROW(NOT_IMPLEMENTED); }
136 virtual double run();
137 virtual double maxDiff() const { return _maxdiff; }
138 virtual size_t Iterations() const { return _iters; }
139 //@}
140
141
142 /// @name Additional interface specific for MR
143 //@{
144 //@}
145
146 private:
147 void init(size_t Nin, double *_w, double *_th);
148 void makekindex();
149 void init_cor();
150 double init_cor_resp();
151 void solvemcav();
152 void solveM();
153
154 double _tJ(size_t i, sub_nb A);
155
156 double Omega(size_t i, size_t _j, size_t _l);
157 double T(size_t i, sub_nb A);
158 double T(size_t i, size_t _j);
159 double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
160 double Gamma(size_t i, size_t _l1, size_t _l2);
161
162 double appM(size_t i, sub_nb A);
163 void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
164
165 double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
166
167 void setProperties( const PropertySet &opts );
168 PropertySet getProperties() const;
169 std::string printProperties() const;
170 };
171
172
173 } // end of namespace dai
174
175
176 #endif