1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
4 This file is part of libDAI.
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.
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.
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
22 #ifndef __defined_libdai_mr_h
23 #define __defined_libdai_mr_h
28 #include <dai/factorgraph.h>
29 #include <dai/daialg.h>
31 #include <dai/properties.h>
32 #include <dai/exceptions.h>
41 class MR
: public DAIAlgFG
{
43 bool supported
; // is the underlying factor graph supported?
45 std::vector
<size_t> con
; // con[i] = connectivity of spin i
46 std::vector
<std::vector
<size_t> > nb
; // nb[i] are the neighbours of spin i
47 std::vector
<std::vector
<double> > tJ
; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
48 std::vector
<double> theta
; // theta[i] is the local field on spin i
49 std::vector
<std::vector
<double> > M
; // M[i][_j] is M^{(i)}_j
50 std::vector
<std::vector
<size_t> > kindex
; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
51 std::vector
<std::vector
<std::vector
<double> > > cors
;
53 static const size_t kmax
= 31;
57 std::vector
<double> Mag
;
63 DAI_ENUM(UpdateType
,FULL
,LINEAR
)
64 DAI_ENUM(InitType
,RESPPROP
,CLAMPING
,EXACT
)
72 MR( const FactorGraph
& fg
, const PropertySet
&opts
);
73 void init(size_t Nin
, double *_w
, double *_th
);
77 double init_cor_resp();
81 Factor
belief( const Var
&n
) const;
82 Factor
belief( const VarSet
&/*ns*/ ) const {
83 DAI_THROW(NOT_IMPLEMENTED
);
86 std::vector
<Factor
> beliefs() const;
88 DAI_THROW(NOT_IMPLEMENTED
);
92 /// Clear messages and beliefs corresponding to the nodes in ns
93 virtual void init( const VarSet
&/*ns*/ ) {
94 DAI_THROW(NOT_IMPLEMENTED
);
96 static const char *Name
;
97 std::string
identify() const;
98 double _tJ(size_t i
, sub_nb A
);
100 double Omega(size_t i
, size_t _j
, size_t _l
);
101 double T(size_t i
, sub_nb A
);
102 double T(size_t i
, size_t _j
);
103 double Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
);
104 double Gamma(size_t i
, size_t _l1
, size_t _l2
);
106 double appM(size_t i
, sub_nb A
);
107 void sum_subs(size_t j
, sub_nb A
, double *sum_even
, double *sum_odd
);
109 double sign(double a
) { return (a
>= 0) ? 1.0 : -1.0; }
110 MR
* clone() const { return new MR(*this); }
111 /// Create (virtual constructor)
112 virtual MR
* create() const { return new MR(); }
114 void setProperties( const PropertySet
&opts
);
115 PropertySet
getProperties() const;
116 std::string
printProperties() const;
117 double maxDiff() const { return maxdiff
; }
121 // represents a subset of nb[i] as a binary number
122 // the elements of a subset should be thought of as indices in nb[i]
129 // construct full subset containing nr_elmt elements
130 sub_nb(size_t nr_elmt
) {
132 assert( nr_elmt
< sizeof(size_t) / sizeof(char) * 8 );
135 s
= (1U << bits
) - 1;
139 sub_nb( const sub_nb
& x
) : s(x
.s
), bits(x
.bits
) {}
141 // assignment operator
142 sub_nb
& operator=( const sub_nb
&x
) {
150 // returns number of elements
153 for(size_t bit
= 0; bit
< bits
; bit
++)
154 if( s
& (1U << bit
) )
159 // increases s by one (for enumeration in lexicographical order)
160 sub_nb
operator++() {
162 if( s
>= (1U << bits
) )
167 // return i'th element of this subset
168 size_t operator[](size_t i
) {
170 for(bit
= 0; bit
< bits
; bit
++ )
171 if( s
& (1U << bit
) ) {
178 assert( bit
< bits
);
183 // add index _j to this subset
184 sub_nb
&operator +=(size_t _j
) {
189 // return copy with index _j
190 sub_nb
operator+(size_t _j
) {
196 // delete index _j from this subset
197 sub_nb
&operator -=(size_t _j
) {
202 // return copy without index _j
203 sub_nb
operator-(size_t _j
) {
215 // returns true if subset is empty
216 bool empty() { return (s
== 0); }
218 // return 1 if _j is contained, 0 otherwise ("is element of")
219 size_t operator&(size_t _j
) { return s
& (1U << _j
); }
221 friend std::ostream
& operator<< (std::ostream
& os
, const sub_nb x
) {
225 for(size_t bit
= x
.bits
; bit
> 0; bit
-- )
226 if( x
.s
& (1U << (bit
-1)) )
236 } // end of namespace dai