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
23 #include <dai/daialg.h>
32 /// Calculate the marginal of obj on ns by clamping
33 /// all variables in ns and calculating logZ for each joined state
34 Factor
calcMarginal( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
37 InfAlg
*clamped
= obj
.clone();
42 for( State
s(ns
); s
.valid(); s
++ ) {
43 // save unclamped factors connected to ns
44 clamped
->saveProbs( ns
);
46 // set clamping Factors to delta functions
47 for( VarSet::const_iterator n
= ns
.begin(); n
!= ns
.end(); n
++ )
48 clamped
->clamp( *n
, s(*n
) );
50 // run DAIAlg, calc logZ, store in Pns
57 logZ0
= clamped
->logZ();
60 // subtract logZ0 to avoid very large numbers
61 Z
= exp(clamped
->logZ() - logZ0
);
66 // restore clamped factors
67 clamped
->undoProbs( ns
);
72 return( Pns
.normalized(Prob::NORMPROB
) );
76 vector
<Factor
> calcPairBeliefs( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
77 // convert ns to vector<VarSet>
81 for( VarSet::const_iterator n
= ns
.begin(); n
!= ns
.end(); n
++ )
84 vector
<Factor
> pairbeliefs
;
85 pairbeliefs
.reserve( N
* N
);
86 for( size_t j
= 0; j
< N
; j
++ )
87 for( size_t k
= 0; k
< N
; k
++ )
89 pairbeliefs
.push_back(Factor());
91 pairbeliefs
.push_back(Factor(vns
[j
] | vns
[k
]));
93 InfAlg
*clamped
= obj
.clone();
98 for( size_t j
= 0; j
< N
; j
++ ) {
99 // clamp Var j to its possible values
100 for( size_t j_val
= 0; j_val
< vns
[j
].states(); j_val
++ ) {
101 // save unclamped factors connected to ns
102 clamped
->saveProbs( ns
);
104 clamped
->clamp( vns
[j
], j_val
);
110 // logZ0 = obj.logZ();
112 if( j
== 0 && j_val
== 0 ) {
113 logZ0
= clamped
->logZ();
115 // subtract logZ0 to avoid very large numbers
116 Z_xj
= exp(clamped
->logZ() - logZ0
);
119 for( size_t k
= 0; k
< N
; k
++ )
121 Factor b_k
= clamped
->belief(vns
[k
]);
122 for( size_t k_val
= 0; k_val
< vns
[k
].states(); k_val
++ )
123 if( vns
[j
].label() < vns
[k
].label() )
124 pairbeliefs
[j
* N
+ k
][j_val
+ (k_val
* vns
[j
].states())] = Z_xj
* b_k
[k_val
];
126 pairbeliefs
[j
* N
+ k
][k_val
+ (j_val
* vns
[k
].states())] = Z_xj
* b_k
[k_val
];
129 // restore clamped factors
130 clamped
->undoProbs( ns
);
136 // Calculate result by taking the geometric average
137 vector
<Factor
> result
;
138 result
.reserve( N
* (N
- 1) / 2 );
139 for( size_t j
= 0; j
< N
; j
++ )
140 for( size_t k
= j
+1; k
< N
; k
++ )
141 result
.push_back( (pairbeliefs
[j
* N
+ k
] * pairbeliefs
[k
* N
+ j
]) ^ 0.5 );
147 Factor
calcMarginal2ndO( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
148 // returns a a probability distribution whose 1st order interactions
149 // are unspecified, whose 2nd order interactions approximate those of
150 // the marginal on ns, and whose higher order interactions are absent.
152 vector
<Factor
> pairbeliefs
= calcPairBeliefs( obj
, ns
, reInit
);
155 for( size_t ij
= 0; ij
< pairbeliefs
.size(); ij
++ )
156 Pns
*= pairbeliefs
[ij
];
158 return( Pns
.normalized(Prob::NORMPROB
) );
162 vector
<Factor
> calcPairBeliefsNew( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
163 vector
<Factor
> result
;
164 result
.reserve( ns
.size() * (ns
.size() - 1) / 2 );
166 InfAlg
*clamped
= obj
.clone();
171 VarSet::const_iterator nj
= ns
.begin();
172 for( long j
= 0; j
< (long)ns
.size() - 1; j
++, nj
++ ) {
174 for( VarSet::const_iterator nk
= nj
; (++nk
) != ns
.end(); k
++ ) {
175 Factor
pairbelief( *nj
| *nk
);
177 // clamp Vars j and k to their possible values
178 for( size_t j_val
= 0; j_val
< nj
->states(); j_val
++ )
179 for( size_t k_val
= 0; k_val
< nk
->states(); k_val
++ ) {
180 // save unclamped factors connected to ns
181 clamped
->saveProbs( ns
);
183 clamped
->clamp( *nj
, j_val
);
184 clamped
->clamp( *nk
, k_val
);
190 if( j_val
== 0 && k_val
== 0 ) {
191 logZ0
= clamped
->logZ();
193 // subtract logZ0 to avoid very large numbers
194 Z_xj
= exp(clamped
->logZ() - logZ0
);
197 // we assume that j.label() < k.label()
198 // i.e. we make an assumption here about the indexing
199 pairbelief
[j_val
+ (k_val
* nj
->states())] = Z_xj
;
201 // restore clamped factors
202 clamped
->undoProbs( ns
);
205 result
.push_back( pairbelief
);
211 assert( result
.size() == (ns
.size() * (ns
.size() - 1) / 2) );
217 } // end of namespace dai