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
51 if( clamped
->Verbose() >= 2 )
59 logZ0
= clamped
->logZ();
62 // subtract logZ0 to avoid very large numbers
63 Z
= exp(clamped
->logZ() - logZ0
);
64 if( fabs(imag(Z
)) > 1e-5 )
65 cout
<< "Marginal:: WARNING: complex Z (" << Z
<< ")" << endl
;
70 // restore clamped factors
71 clamped
->undoProbs( ns
);
76 return( Pns
.normalized(Prob::NORMPROB
) );
80 vector
<Factor
> calcPairBeliefs( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
81 // convert ns to vector<VarSet>
85 for( VarSet::const_iterator n
= ns
.begin(); n
!= ns
.end(); n
++ )
88 vector
<Factor
> pairbeliefs
;
89 pairbeliefs
.reserve( N
* N
);
90 for( size_t j
= 0; j
< N
; j
++ )
91 for( size_t k
= 0; k
< N
; k
++ )
93 pairbeliefs
.push_back(Factor());
95 pairbeliefs
.push_back(Factor(vns
[j
] | vns
[k
]));
97 InfAlg
*clamped
= obj
.clone();
102 for( size_t j
= 0; j
< N
; j
++ ) {
103 // clamp Var j to its possible values
104 for( size_t j_val
= 0; j_val
< vns
[j
].states(); j_val
++ ) {
105 if( obj
.Verbose() >= 2 )
106 cout
<< j
<< "/" << N
-1 << " (" << j_val
<< "/" << vns
[j
].states() << "): ";
108 // save unclamped factors connected to ns
109 clamped
->saveProbs( ns
);
111 clamped
->clamp( vns
[j
], j_val
);
117 // logZ0 = obj.logZ();
119 if( j
== 0 && j_val
== 0 ) {
120 logZ0
= clamped
->logZ();
122 // subtract logZ0 to avoid very large numbers
123 Complex Z
= exp(clamped
->logZ() - logZ0
);
124 if( fabs(imag(Z
)) > 1e-5 )
125 cout
<< "calcPairBelief:: Warning: complex Z: " << Z
<< endl
;
129 for( size_t k
= 0; k
< N
; k
++ )
131 Factor b_k
= clamped
->belief(vns
[k
]);
132 for( size_t k_val
= 0; k_val
< vns
[k
].states(); k_val
++ )
133 if( vns
[j
].label() < vns
[k
].label() )
134 pairbeliefs
[j
* N
+ k
][j_val
+ (k_val
* vns
[j
].states())] = Z_xj
* b_k
[k_val
];
136 pairbeliefs
[j
* N
+ k
][k_val
+ (j_val
* vns
[k
].states())] = Z_xj
* b_k
[k_val
];
139 // restore clamped factors
140 clamped
->undoProbs( ns
);
146 // Calculate result by taking the geometric average
147 vector
<Factor
> result
;
148 result
.reserve( N
* (N
- 1) / 2 );
149 for( size_t j
= 0; j
< N
; j
++ )
150 for( size_t k
= j
+1; k
< N
; k
++ )
151 result
.push_back( (pairbeliefs
[j
* N
+ k
] * pairbeliefs
[k
* N
+ j
]) ^ 0.5 );
157 Factor
calcMarginal2ndO( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
158 // returns a a probability distribution whose 1st order interactions
159 // are unspecified, whose 2nd order interactions approximate those of
160 // the marginal on ns, and whose higher order interactions are absent.
162 vector
<Factor
> pairbeliefs
= calcPairBeliefs( obj
, ns
, reInit
);
165 for( size_t ij
= 0; ij
< pairbeliefs
.size(); ij
++ )
166 Pns
*= pairbeliefs
[ij
];
168 return( Pns
.normalized(Prob::NORMPROB
) );
172 vector
<Factor
> calcPairBeliefsNew( const InfAlg
& obj
, const VarSet
& ns
, bool reInit
) {
173 vector
<Factor
> result
;
174 result
.reserve( ns
.size() * (ns
.size() - 1) / 2 );
176 InfAlg
*clamped
= obj
.clone();
181 VarSet::const_iterator nj
= ns
.begin();
182 for( long j
= 0; j
< (long)ns
.size() - 1; j
++, nj
++ ) {
184 for( VarSet::const_iterator nk
= nj
; (++nk
) != ns
.end(); k
++ ) {
185 Factor
pairbelief( *nj
| *nk
);
187 // clamp Vars j and k to their possible values
188 for( size_t j_val
= 0; j_val
< nj
->states(); j_val
++ )
189 for( size_t k_val
= 0; k_val
< nk
->states(); k_val
++ ) {
190 // save unclamped factors connected to ns
191 clamped
->saveProbs( ns
);
193 clamped
->clamp( *nj
, j_val
);
194 clamped
->clamp( *nk
, k_val
);
200 if( j_val
== 0 && k_val
== 0 ) {
201 logZ0
= clamped
->logZ();
203 // subtract logZ0 to avoid very large numbers
204 Complex Z
= exp(clamped
->logZ() - logZ0
);
205 if( fabs(imag(Z
)) > 1e-5 )
206 cout
<< "calcPairBelief:: Warning: complex Z: " << Z
<< endl
;
210 // we assume that j.label() < k.label()
211 // i.e. we make an assumption here about the indexing
212 pairbelief
[j_val
+ (k_val
* nj
->states())] = Z_xj
;
214 // restore clamped factors
215 clamped
->undoProbs( ns
);
218 result
.push_back( pairbelief
);
224 assert( result
.size() == (ns
.size() * (ns
.size() - 1) / 2) );
230 } // end of namespace dai