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 bool LC::checkProperties() {
43 if( !HasProperty("cavity") )
45 if( !HasProperty("updates") )
47 if( !HasProperty("tol") )
49 if (!HasProperty("maxiter") )
51 if (!HasProperty("verbose") )
54 ConvertPropertyTo
<CavityType
>("cavity");
55 ConvertPropertyTo
<UpdateType
>("updates");
56 ConvertPropertyTo
<double>("tol");
57 ConvertPropertyTo
<size_t>("maxiter");
58 ConvertPropertyTo
<size_t>("verbose");
60 if (HasProperty("cavainame") )
61 ConvertPropertyTo
<string
>("cavainame");
62 if (HasProperty("cavaiopts") )
63 ConvertPropertyTo
<Properties
>("cavaiopts");
64 if( HasProperty("reinit") )
65 ConvertPropertyTo
<bool>("reinit");
71 LC::LC(const FactorGraph
& fg
, const Properties
&opts
) : DAIAlgFG(fg
, opts
) {
72 assert( checkProperties() );
75 for( size_t i
=0; i
< nrVars(); i
++ ) {
76 for( _nb_cit I
= nb1(i
).begin(); I
!= nb1(i
).end(); I
++ ) {
81 _iI
.push_back(_iI_entry
);
86 _pancakes
.resize(nrVars());
89 for( size_t i
=0; i
< nrVars(); i
++ )
90 _cavitydists
.push_back(Factor(delta(var(i
))));
93 _phis
.reserve(nr_edges());
94 for( size_t iI
= 0; iI
< nr_edges(); iI
++ ) {
95 size_t i
= edge(iI
).first
;
96 size_t I
= edge(iI
).second
;
97 _phis
.push_back( Factor( factor(I
).vars() / var(i
) ) );
101 for( size_t i
=0; i
< nrVars(); i
++ )
102 _beliefs
.push_back(Factor(var(i
)));
106 string
LC::identify() const {
107 stringstream
result (stringstream::out
);
108 result
<< Name
<< GetProperties();
113 void LC::CalcBelief (size_t i
) {
114 _beliefs
[i
] = _pancakes
[i
].marginal(var(i
));
118 double LC::CalcCavityDist (size_t i
, const string
&name
, const Properties
&opts
) {
123 cout
<< "Initing cavity " << var(i
) << "(" << delta(var(i
)).size() << " vars, " << delta(var(i
)).stateSpace() << " states)" << endl
;
125 if( Cavity() == CavityType::UNIFORM
)
126 Bi
= Factor(delta(var(i
)));
128 InfAlg
*cav
= newInfAlg( name
, *this, opts
);
129 cav
->makeCavity( var(i
) );
131 if( Cavity() == CavityType::FULL
)
132 Bi
= calcMarginal( *cav
, cav
->delta(var(i
)), reInit() );
133 else if( Cavity() == CavityType::PAIR
)
134 Bi
= calcMarginal2ndO( *cav
, cav
->delta(var(i
)), reInit() );
135 else if( Cavity() == CavityType::PAIR2
) {
136 vector
<Factor
> pairbeliefs
= calcPairBeliefsNew( *cav
, cav
->delta(var(i
)), reInit() );
137 for( size_t ij
= 0; ij
< pairbeliefs
.size(); ij
++ )
138 Bi
*= pairbeliefs
[ij
];
139 } else if( Cavity() == CavityType::PAIRINT
) {
140 Bi
= calcMarginal( *cav
, cav
->delta(var(i
)), reInit() );
142 // Set interactions of order > 2 to zero
143 size_t N
= delta(var(i
)).size();
144 Real
*p
= &(*Bi
.p().p().begin());
147 x2x::fill (N
, p
, 2, 0.0);
149 // x2x::logpnorm (N, p);
151 } else if( Cavity() == CavityType::PAIRCUM
) {
152 Bi
= calcMarginal( *cav
, cav
->delta(var(i
)), reInit() );
154 // Set cumulants of order > 2 to zero
155 size_t N
= delta(var(i
)).size();
156 Real
*p
= &(*Bi
.p().p().begin());
159 x2x::fill (N
, p
, 2, 0.0);
163 maxdiff
= cav
->MaxDiff();
166 Bi
.normalize( _normtype
);
167 _cavitydists
[i
] = Bi
;
173 double LC::InitCavityDists (const string
&name
, const Properties
&opts
) {
176 if( Verbose() >= 1 ) {
177 cout
<< "LC::InitCavityDists: ";
178 if( Cavity() == CavityType::UNIFORM
)
179 cout
<< "Using uniform initial cavity distributions" << endl
;
180 else if( Cavity() == CavityType::FULL
)
181 cout
<< "Using full " << name
<< opts
<< "...";
182 else if( Cavity() == CavityType::PAIR
)
183 cout
<< "Using pairwise " << name
<< opts
<< "...";
184 else if( Cavity() == CavityType::PAIR2
)
185 cout
<< "Using pairwise(new) " << name
<< opts
<< "...";
188 double maxdiff
= 0.0;
189 for( size_t i
= 0; i
< nrVars(); i
++ ) {
190 double md
= CalcCavityDist(i
, name
, opts
);
196 if( Verbose() >= 1 ) {
197 cout
<< "used " << toc() - tic
<< " clocks." << endl
;
204 long LC::SetCavityDists (vector
<Factor
> &Q
) {
206 cout
<< "LC::SetCavityDists: Setting initial cavity distributions" << endl
;
207 if( Q
.size() != nrVars() )
209 for( size_t i
= 0; i
< nrVars(); i
++ ) {
210 if( _cavitydists
[i
].vars() != Q
[i
].vars() ) {
213 _cavitydists
[i
] = Q
[i
];
221 for( size_t iI
= 0; iI
< nr_edges(); iI
++ ) {
222 if( Updates() == UpdateType::SEQRND
)
223 _phis
[iI
].randomize();
227 for( size_t i
= 0; i
< nrVars(); i
++ ) {
228 _pancakes
[i
] = _cavitydists
[i
];
230 for( _nb_cit I
= nb1(i
).begin(); I
!= nb1(i
).end(); I
++ ) {
231 _pancakes
[i
] *= factor(*I
);
232 if( Updates() == UpdateType::SEQRND
)
233 _pancakes
[i
] *= _phis
[VV2E(i
,*I
)];
236 _pancakes
[i
].normalize( _normtype
);
243 Factor
LC::NewPancake (size_t iI
, bool & hasNaNs
) {
244 size_t i
= _iI
[iI
].i
;
245 size_t I
= _iI
[iI
].I
;
248 Factor piet
= _pancakes
[i
];
250 // recalculate _pancake[i]
251 VarSet Ivars
= factor(I
).vars();
253 for( VarSet::const_iterator k
= Ivars
.begin(); k
!= Ivars
.end(); k
++ )
255 A_I
*= (_pancakes
[findVar(*k
)] * factor(I
).inverse()).part_sum( Ivars
/ var(i
) );
256 if( Ivars
.size() > 1 )
257 A_I
^= (1.0 / (Ivars
.size() - 1));
258 Factor A_Ii
= (_pancakes
[i
] * factor(I
).inverse() * _phis
[iI
].inverse()).part_sum( Ivars
/ var(i
) );
259 Factor quot
= A_I
.divided_by(A_Ii
);
261 piet
*= quot
.divided_by( _phis
[iI
] ).normalized( _normtype
);
262 _phis
[iI
] = quot
.normalized( _normtype
);
264 piet
.normalize( _normtype
);
266 if( piet
.hasNaNs() ) {
267 cout
<< "LC::NewPancake(" << iI
<< "): has NaNs!" << endl
;
277 cout
<< "Starting " << identify() << "...";
282 Diffs
diffs(nrVars(), 1.0);
284 updateMaxDiff( InitCavityDists(GetPropertyAs
<string
>("cavainame"), GetPropertyAs
<Properties
>("cavaiopts")) );
286 vector
<Factor
> old_beliefs
;
287 for(size_t i
=0; i
< nrVars(); i
++ )
288 old_beliefs
.push_back(belief(i
));
290 bool hasNaNs
= false;
291 for( size_t i
=0; i
< nrVars(); i
++ )
292 if( _pancakes
[i
].hasNaNs() ) {
297 cout
<< "LC::run: initial _pancakes has NaNs!" << endl
;
301 vector
<long> update_seq(nr_iI(),0);
302 for( size_t k
=0; k
< nr_iI(); k
++ )
307 // do several passes over the network until maximum number of iterations has
308 // been reached or until the maximum belief difference is smaller than tolerance
309 for( iter
=0; iter
< MaxIter() && diffs
.max() > Tol(); iter
++ ) {
310 // Sequential updates
311 if( Updates() == UpdateType::SEQRND
)
312 random_shuffle( update_seq
.begin(), update_seq
.end() );
314 for( size_t t
=0; t
< nr_iI(); t
++ ) {
315 long iI
= update_seq
[t
];
317 _pancakes
[i
] = NewPancake(iI
, hasNaNs
);
323 // compare new beliefs with old ones
324 for(size_t i
=0; i
< nrVars(); i
++ ) {
325 diffs
.push( dist( belief(i
), old_beliefs
[i
], Prob::DISTLINF
) );
326 old_beliefs
[i
] = belief(i
);
330 cout
<< "LC::run: maxdiff " << diffs
.max() << " after " << iter
+1 << " passes" << endl
;
333 updateMaxDiff( diffs
.max() );
335 if( Verbose() >= 1 ) {
336 if( diffs
.max() > Tol() ) {
339 cout
<< "LC::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.max() << endl
;
343 cout
<< "converged in " << iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
351 } // end of namespace dai