Adopted contributions by Christian.
[libdai.git] / mr.h
1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #ifndef __MR_H__
23 #define __MR_H__
24
25
26 #include <vector>
27 #include "factorgraph.h"
28 #include "daialg.h"
29 #include "enum.h"
30
31
32 namespace dai {
33
34
35 using namespace std;
36
37
38 class sub_nb;
39
40
41 class MR : public DAIAlgFG {
42 private:
43 bool supported; // is the underlying factor graph supported?
44
45 vector<size_t> con; // con[i] = connectivity of spin i
46 vector<vector<size_t> > nb; // nb[i] are the neighbours of spin i
47 vector<vector<double> > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
48 vector<double> theta; // theta[i] is the local field on spin i
49 vector<vector<double> > M; // M[i][_j] is M^{(i)}_j
50 vector<vector<size_t> > kindex; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
51 vector<vector<vector<double> > > cors;
52
53 static const size_t kmax = 31;
54
55 size_t N;
56
57 vector<double> Mag;
58
59 public:
60 ENUM2(UpdateType,FULL,LINEAR)
61 ENUM3(InitType,RESPPROP,CLAMPING,EXACT)
62
63 UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
64 InitType Inits() const { return GetPropertyAs<InitType>("inits"); }
65
66 MR( const FactorGraph & fg, const Properties &opts );
67 void init(size_t _N, double *_w, double *_th);
68 void makekindex();
69 void read_files();
70 double init_cor();
71 double init_cor_resp();
72 void solvemcav();
73 void solveM();
74 double run();
75 Factor belief( const Var &n ) const;
76 Factor belief( const VarSet &ns ) const { assert( 0 == 1 ); }
77 vector<Factor> beliefs() const;
78 Complex logZ() const { return NAN; }
79 void init() { assert( checkProperties() ); }
80 static const char *Name;
81 string identify() const;
82 double _tJ(size_t i, sub_nb A);
83
84 double Omega(size_t i, size_t _j, size_t _l);
85 double T(size_t i, sub_nb A);
86 double T(size_t i, size_t _j);
87 double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
88 double Gamma(size_t i, size_t _l1, size_t _l2);
89
90 double appM(size_t i, sub_nb A);
91 void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
92
93 double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
94 MR* clone() const { assert( 0 == 1 ); }
95
96 bool checkProperties();
97
98 };
99
100
101 // represents a subset of nb[i] as a binary number
102 // the elements of a subset should be thought of as indices in nb[i]
103 class sub_nb {
104 private:
105 size_t s;
106 size_t bits;
107
108 public:
109 // construct full subset containing nr_elmt elements
110 sub_nb(size_t nr_elmt) {
111 #ifdef DEBUG
112 assert( nr_elmt < sizeof(size_t) / sizeof(char) * 8 );
113 #endif
114 bits = nr_elmt;
115 s = (1U << bits) - 1;
116 }
117
118 // copy constructor
119 sub_nb( const sub_nb & x ) : s(x.s), bits(x.bits) {}
120
121 // assignment operator
122 sub_nb & operator=( const sub_nb &x ) {
123 if( this != &x ) {
124 s = x.s;
125 bits = x.bits;
126 }
127 return *this;
128 }
129
130 // returns number of elements
131 size_t size() {
132 size_t size = 0;
133 for(size_t bit = 0; bit < bits; bit++)
134 if( s & (1U << bit) )
135 size++;
136 return size;
137 }
138
139 // increases s by one (for enumeration in lexicographical order)
140 sub_nb operator++() {
141 s++;
142 if( s >= (1U << bits) )
143 s = 0;
144 return *this;
145 }
146
147 // return i'th element of this subset
148 size_t operator[](size_t i) {
149 size_t bit;
150 for(bit = 0; bit < bits; bit++ )
151 if( s & (1U << bit) ) {
152 if( i == 0 )
153 break;
154 else
155 i--;
156 }
157 #ifdef DEBUG
158 assert( bit < bits );
159 #endif
160 return bit;
161 }
162
163 // add index _j to this subset
164 sub_nb &operator +=(size_t _j) {
165 s |= (1U << _j);
166 return *this;
167 }
168
169 // return copy with index _j
170 sub_nb operator+(size_t _j) {
171 sub_nb x = *this;
172 x += _j;
173 return x;
174 }
175
176 // delete index _j from this subset
177 sub_nb &operator -=(size_t _j) {
178 s &= ~(1U << _j);
179 return *this;
180 }
181
182 // return copy without index _j
183 sub_nb operator-(size_t _j) {
184 sub_nb x = *this;
185 x -= _j;
186 return x;
187 }
188
189 // empty this subset
190 sub_nb & clear() {
191 s = 0;
192 return *this;
193 }
194
195 // returns true if subset is empty
196 bool empty() { return (s == 0); }
197
198 // return 1 if _j is contained, 0 otherwise ("is element of")
199 size_t operator&(size_t _j) { return s & (1U << _j); }
200
201 friend std::ostream& operator<< (std::ostream& os, const sub_nb x) {
202 if( x.bits == 0 )
203 os << "empty";
204 else {
205 for(size_t bit = x.bits; bit > 0; bit-- )
206 if( x.s & (1U << (bit-1)) )
207 os << "1";
208 else
209 os << "0";
210 }
211 return os;
212 }
213 };
214
215
216 }
217
218
219 #endif