Misc small changes:
[libdai.git] / include / dai / 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 __defined_libdai_mr_h
23 #define __defined_libdai_mr_h
24
25
26 #include <vector>
27 #include <string>
28 #include <dai/factorgraph.h>
29 #include <dai/daialg.h>
30 #include <dai/enum.h>
31
32
33 namespace dai {
34
35
36 class sub_nb;
37
38
39 class MR : public DAIAlgFG {
40 private:
41 bool supported; // is the underlying factor graph supported?
42
43 std::vector<size_t> con; // con[i] = connectivity of spin i
44 std::vector<std::vector<size_t> > nb; // nb[i] are the neighbours of spin i
45 std::vector<std::vector<double> > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
46 std::vector<double> theta; // theta[i] is the local field on spin i
47 std::vector<std::vector<double> > M; // M[i][_j] is M^{(i)}_j
48 std::vector<std::vector<size_t> > kindex; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
49 std::vector<std::vector<std::vector<double> > > cors;
50
51 static const size_t kmax = 31;
52
53 size_t N;
54
55 std::vector<double> Mag;
56
57 public:
58 ENUM2(UpdateType,FULL,LINEAR)
59 ENUM3(InitType,RESPPROP,CLAMPING,EXACT)
60
61 UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
62 InitType Inits() const { return GetPropertyAs<InitType>("inits"); }
63
64 MR( const FactorGraph & fg, const Properties &opts );
65 void init(size_t Nin, double *_w, double *_th);
66 void makekindex();
67 void read_files();
68 void init_cor();
69 double init_cor_resp();
70 void solvemcav();
71 void solveM();
72 double run();
73 Factor belief( const Var &n ) const;
74 Factor belief( const VarSet &/*ns*/ ) const { assert( 0 == 1 ); }
75 std::vector<Factor> beliefs() const;
76 Complex logZ() const { return NAN; }
77 void init() { assert( checkProperties() ); }
78 static const char *Name;
79 std::string identify() const;
80 double _tJ(size_t i, sub_nb A);
81
82 double Omega(size_t i, size_t _j, size_t _l);
83 double T(size_t i, sub_nb A);
84 double T(size_t i, size_t _j);
85 double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
86 double Gamma(size_t i, size_t _l1, size_t _l2);
87
88 double appM(size_t i, sub_nb A);
89 void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
90
91 double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
92 MR* clone() const { assert( 0 == 1 ); }
93
94 bool checkProperties();
95
96 };
97
98
99 // represents a subset of nb[i] as a binary number
100 // the elements of a subset should be thought of as indices in nb[i]
101 class sub_nb {
102 private:
103 size_t s;
104 size_t bits;
105
106 public:
107 // construct full subset containing nr_elmt elements
108 sub_nb(size_t nr_elmt) {
109 #ifdef DAI_DEBUG
110 assert( nr_elmt < sizeof(size_t) / sizeof(char) * 8 );
111 #endif
112 bits = nr_elmt;
113 s = (1U << bits) - 1;
114 }
115
116 // copy constructor
117 sub_nb( const sub_nb & x ) : s(x.s), bits(x.bits) {}
118
119 // assignment operator
120 sub_nb & operator=( const sub_nb &x ) {
121 if( this != &x ) {
122 s = x.s;
123 bits = x.bits;
124 }
125 return *this;
126 }
127
128 // returns number of elements
129 size_t size() {
130 size_t size = 0;
131 for(size_t bit = 0; bit < bits; bit++)
132 if( s & (1U << bit) )
133 size++;
134 return size;
135 }
136
137 // increases s by one (for enumeration in lexicographical order)
138 sub_nb operator++() {
139 s++;
140 if( s >= (1U << bits) )
141 s = 0;
142 return *this;
143 }
144
145 // return i'th element of this subset
146 size_t operator[](size_t i) {
147 size_t bit;
148 for(bit = 0; bit < bits; bit++ )
149 if( s & (1U << bit) ) {
150 if( i == 0 )
151 break;
152 else
153 i--;
154 }
155 #ifdef DAI_DEBUG
156 assert( bit < bits );
157 #endif
158 return bit;
159 }
160
161 // add index _j to this subset
162 sub_nb &operator +=(size_t _j) {
163 s |= (1U << _j);
164 return *this;
165 }
166
167 // return copy with index _j
168 sub_nb operator+(size_t _j) {
169 sub_nb x = *this;
170 x += _j;
171 return x;
172 }
173
174 // delete index _j from this subset
175 sub_nb &operator -=(size_t _j) {
176 s &= ~(1U << _j);
177 return *this;
178 }
179
180 // return copy without index _j
181 sub_nb operator-(size_t _j) {
182 sub_nb x = *this;
183 x -= _j;
184 return x;
185 }
186
187 // empty this subset
188 sub_nb & clear() {
189 s = 0;
190 return *this;
191 }
192
193 // returns true if subset is empty
194 bool empty() { return (s == 0); }
195
196 // return 1 if _j is contained, 0 otherwise ("is element of")
197 size_t operator&(size_t _j) { return s & (1U << _j); }
198
199 friend std::ostream& operator<< (std::ostream& os, const sub_nb x) {
200 if( x.bits == 0 )
201 os << "empty";
202 else {
203 for(size_t bit = x.bits; bit > 0; bit-- )
204 if( x.s & (1U << (bit-1)) )
205 os << "1";
206 else
207 os << "0";
208 }
209 return os;
210 }
211 };
212
213
214 } // end of namespace dai
215
216
217 #endif