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
27 #include <dai/diffs.h>
29 #include <dai/alldai.h>
39 const char *LC::Name
= "LC";
42 void LC::setProperties( const PropertySet
&opts
) {
43 assert( opts
.hasKey("tol") );
44 assert( opts
.hasKey("maxiter") );
45 assert( opts
.hasKey("verbose") );
46 assert( opts
.hasKey("cavity") );
47 assert( opts
.hasKey("updates") );
49 props
.tol
= opts
.getStringAs
<double>("tol");
50 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
51 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
52 props
.cavity
= opts
.getStringAs
<Properties::CavityType
>("cavity");
53 props
.updates
= opts
.getStringAs
<Properties::UpdateType
>("updates");
54 if( opts
.hasKey("cavainame") )
55 props
.cavainame
= opts
.getStringAs
<string
>("cavainame");
56 if( opts
.hasKey("cavaiopts") )
57 props
.cavaiopts
= opts
.getStringAs
<PropertySet
>("cavaiopts");
58 if( opts
.hasKey("reinit") )
59 props
.reinit
= opts
.getStringAs
<bool>("reinit");
63 PropertySet
LC::getProperties() const {
65 opts
.Set( "tol", props
.tol
);
66 opts
.Set( "maxiter", props
.maxiter
);
67 opts
.Set( "verbose", props
.verbose
);
68 opts
.Set( "cavity", props
.cavity
);
69 opts
.Set( "updates", props
.updates
);
70 opts
.Set( "cavainame", props
.cavainame
);
71 opts
.Set( "cavaiopts", props
.cavaiopts
);
72 opts
.Set( "reinit", props
.reinit
);
77 string
LC::printProperties() const {
78 stringstream
s( stringstream::out
);
80 s
<< "tol=" << props
.tol
<< ",";
81 s
<< "maxiter=" << props
.maxiter
<< ",";
82 s
<< "verbose=" << props
.verbose
<< ",";
83 s
<< "cavity=" << props
.cavity
<< ",";
84 s
<< "updates=" << props
.updates
<< ",";
85 s
<< "cavainame=" << props
.cavainame
<< ",";
86 s
<< "cavaiopts=" << props
.cavaiopts
<< ",";
87 s
<< "reinit=" << props
.reinit
<< "]";
92 LC::LC( const FactorGraph
& fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), _pancakes(), _cavitydists(), _phis(), _beliefs(), props(), maxdiff(0.0) {
93 setProperties( opts
);
96 _pancakes
.resize(nrVars());
99 for( size_t i
=0; i
< nrVars(); i
++ )
100 _cavitydists
.push_back(Factor(delta(i
)));
103 _phis
.reserve( nrVars() );
104 for( size_t i
= 0; i
< nrVars(); i
++ ) {
105 _phis
.push_back( vector
<Factor
>() );
106 _phis
[i
].reserve( nbV(i
).size() );
107 foreach( const Neighbor
&I
, nbV(i
) )
108 _phis
[i
].push_back( Factor( factor(I
).vars() / var(i
) ) );
112 for( size_t i
=0; i
< nrVars(); i
++ )
113 _beliefs
.push_back(Factor(var(i
)));
117 string
LC::identify() const {
118 return string(Name
) + printProperties();
122 void LC::CalcBelief (size_t i
) {
123 _beliefs
[i
] = _pancakes
[i
].marginal(var(i
));
127 double LC::CalcCavityDist (size_t i
, const std::string
&name
, const PropertySet
&opts
) {
131 if( props
.verbose
>= 2 )
132 cout
<< "Initing cavity " << var(i
) << "(" << delta(i
).size() << " vars, " << delta(i
).states() << " states)" << endl
;
134 if( props
.cavity
== Properties::CavityType::UNIFORM
)
135 Bi
= Factor(delta(i
));
137 InfAlg
*cav
= newInfAlg( name
, *this, opts
);
138 cav
->makeCavity( i
);
140 if( props
.cavity
== Properties::CavityType::FULL
)
141 Bi
= calcMarginal( *cav
, cav
->fg().delta(i
), props
.reinit
);
142 else if( props
.cavity
== Properties::CavityType::PAIR
)
143 Bi
= calcMarginal2ndO( *cav
, cav
->fg().delta(i
), props
.reinit
);
144 else if( props
.cavity
== Properties::CavityType::PAIR2
) {
145 vector
<Factor
> pairbeliefs
= calcPairBeliefsNew( *cav
, cav
->fg().delta(i
), props
.reinit
);
146 for( size_t ij
= 0; ij
< pairbeliefs
.size(); ij
++ )
147 Bi
*= pairbeliefs
[ij
];
148 } else if( props
.cavity
== Properties::CavityType::PAIRINT
) {
149 Bi
= calcMarginal( *cav
, cav
->fg().delta(i
), props
.reinit
);
151 // Set interactions of order > 2 to zero
152 size_t N
= delta(i
).size();
153 Real
*p
= &(*Bi
.p().p().begin());
156 x2x::fill (N
, p
, 2, 0.0);
158 // x2x::logpnorm (N, p);
160 } else if( props
.cavity
== Properties::CavityType::PAIRCUM
) {
161 Bi
= calcMarginal( *cav
, cav
->fg().delta(i
), props
.reinit
);
163 // Set cumulants of order > 2 to zero
164 size_t N
= delta(i
).size();
165 Real
*p
= &(*Bi
.p().p().begin());
168 x2x::fill (N
, p
, 2, 0.0);
172 maxdiff
= cav
->maxDiff();
175 Bi
.normalize( Prob::NORMPROB
);
176 _cavitydists
[i
] = Bi
;
182 double LC::InitCavityDists( const std::string
&name
, const PropertySet
&opts
) {
185 if( props
.verbose
>= 1 ) {
186 cout
<< "LC::InitCavityDists: ";
187 if( props
.cavity
== Properties::CavityType::UNIFORM
)
188 cout
<< "Using uniform initial cavity distributions" << endl
;
189 else if( props
.cavity
== Properties::CavityType::FULL
)
190 cout
<< "Using full " << name
<< opts
<< "...";
191 else if( props
.cavity
== Properties::CavityType::PAIR
)
192 cout
<< "Using pairwise " << name
<< opts
<< "...";
193 else if( props
.cavity
== Properties::CavityType::PAIR2
)
194 cout
<< "Using pairwise(new) " << name
<< opts
<< "...";
197 double maxdiff
= 0.0;
198 for( size_t i
= 0; i
< nrVars(); i
++ ) {
199 double md
= CalcCavityDist(i
, name
, opts
);
205 if( props
.verbose
>= 1 ) {
206 cout
<< "used " << toc() - tic
<< " clocks." << endl
;
213 long LC::SetCavityDists( std::vector
<Factor
> &Q
) {
214 if( props
.verbose
>= 1 )
215 cout
<< "LC::SetCavityDists: Setting initial cavity distributions" << endl
;
216 if( Q
.size() != nrVars() )
218 for( size_t i
= 0; i
< nrVars(); i
++ ) {
219 if( _cavitydists
[i
].vars() != Q
[i
].vars() ) {
222 _cavitydists
[i
] = Q
[i
];
230 for( size_t i
= 0; i
< nrVars(); ++i
)
231 foreach( const Neighbor
&I
, nbV(i
) )
232 if( props
.updates
== Properties::UpdateType::SEQRND
)
233 _phis
[i
][I
.iter
].randomize();
235 _phis
[i
][I
.iter
].fill(1.0);
236 for( size_t i
= 0; i
< nrVars(); i
++ ) {
237 _pancakes
[i
] = _cavitydists
[i
];
239 foreach( const Neighbor
&I
, nbV(i
) ) {
240 _pancakes
[i
] *= factor(I
);
241 if( props
.updates
== Properties::UpdateType::SEQRND
)
242 _pancakes
[i
] *= _phis
[i
][I
.iter
];
245 _pancakes
[i
].normalize( Prob::NORMPROB
);
252 Factor
LC::NewPancake (size_t i
, size_t _I
, bool & hasNaNs
) {
253 size_t I
= nbV(i
)[_I
];
254 Factor piet
= _pancakes
[i
];
256 // recalculate _pancake[i]
257 VarSet Ivars
= factor(I
).vars();
259 for( VarSet::const_iterator k
= Ivars
.begin(); k
!= Ivars
.end(); k
++ )
261 A_I
*= (_pancakes
[findVar(*k
)] * factor(I
).inverse()).partSum( Ivars
/ var(i
) );
262 if( Ivars
.size() > 1 )
263 A_I
^= (1.0 / (Ivars
.size() - 1));
264 Factor A_Ii
= (_pancakes
[i
] * factor(I
).inverse() * _phis
[i
][_I
].inverse()).partSum( Ivars
/ var(i
) );
265 Factor quot
= A_I
.divided_by(A_Ii
);
267 piet
*= quot
.divided_by( _phis
[i
][_I
] ).normalized( Prob::NORMPROB
);
268 _phis
[i
][_I
] = quot
.normalized( Prob::NORMPROB
);
270 piet
.normalize( Prob::NORMPROB
);
272 if( piet
.hasNaNs() ) {
273 cout
<< "LC::NewPancake(" << i
<< ", " << _I
<< "): has NaNs!" << endl
;
282 if( props
.verbose
>= 1 )
283 cout
<< "Starting " << identify() << "...";
284 if( props
.verbose
>= 2 )
288 Diffs
diffs(nrVars(), 1.0);
290 double md
= InitCavityDists( props
.cavainame
, props
.cavaiopts
);
294 vector
<Factor
> old_beliefs
;
295 for(size_t i
=0; i
< nrVars(); i
++ )
296 old_beliefs
.push_back(belief(i
));
298 bool hasNaNs
= false;
299 for( size_t i
=0; i
< nrVars(); i
++ )
300 if( _pancakes
[i
].hasNaNs() ) {
305 cout
<< "LC::run: initial _pancakes has NaNs!" << endl
;
309 size_t nredges
= nrEdges();
310 vector
<Edge
> update_seq
;
311 update_seq
.reserve( nredges
);
312 for( size_t i
= 0; i
< nrVars(); ++i
)
313 foreach( const Neighbor
&I
, nbV(i
) )
314 update_seq
.push_back( Edge( i
, I
.iter
) );
318 // do several passes over the network until maximum number of iterations has
319 // been reached or until the maximum belief difference is smaller than tolerance
320 for( iter
=0; iter
< props
.maxiter
&& diffs
.maxDiff() > props
.tol
; iter
++ ) {
321 // Sequential updates
322 if( props
.updates
== Properties::UpdateType::SEQRND
)
323 random_shuffle( update_seq
.begin(), update_seq
.end() );
325 for( size_t t
=0; t
< nredges
; t
++ ) {
326 size_t i
= update_seq
[t
].first
;
327 size_t _I
= update_seq
[t
].second
;
328 _pancakes
[i
] = NewPancake( i
, _I
, hasNaNs
);
334 // compare new beliefs with old ones
335 for(size_t i
=0; i
< nrVars(); i
++ ) {
336 diffs
.push( dist( belief(i
), old_beliefs
[i
], Prob::DISTLINF
) );
337 old_beliefs
[i
] = belief(i
);
340 if( props
.verbose
>= 3 )
341 cout
<< "LC::run: maxdiff " << diffs
.maxDiff() << " after " << iter
+1 << " passes" << endl
;
344 if( diffs
.maxDiff() > maxdiff
)
345 maxdiff
= diffs
.maxDiff();
347 if( props
.verbose
>= 1 ) {
348 if( diffs
.maxDiff() > props
.tol
) {
349 if( props
.verbose
== 1 )
351 cout
<< "LC::run: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.maxDiff() << endl
;
353 if( props
.verbose
>= 2 )
355 cout
<< "converged in " << iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
359 return diffs
.maxDiff();
363 } // end of namespace dai