-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
+/* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
+ Radboud University Nijmegen, The Netherlands /
+ Max Planck Institute for Biological Cybernetics, Germany
+
This file is part of libDAI.
libDAI is free software; you can redistribute it and/or modify
#include <map>
#include <set>
#include <dai/lc.h>
-#include <dai/diffs.h>
#include <dai/util.h>
#include <dai/alldai.h>
-#include <dai/x2x.h>
namespace dai {
const char *LC::Name = "LC";
-bool LC::checkProperties() {
- if( !HasProperty("cavity") )
- return false;
- if( !HasProperty("updates") )
- return false;
- if( !HasProperty("tol") )
- return false;
- if (!HasProperty("maxiter") )
- return false;
- if (!HasProperty("verbose") )
- return false;
-
- ConvertPropertyTo<CavityType>("cavity");
- ConvertPropertyTo<UpdateType>("updates");
- ConvertPropertyTo<double>("tol");
- ConvertPropertyTo<size_t>("maxiter");
- ConvertPropertyTo<size_t>("verbose");
-
- if (HasProperty("cavainame") )
- ConvertPropertyTo<string>("cavainame");
- if (HasProperty("cavaiopts") )
- ConvertPropertyTo<Properties>("cavaiopts");
- if( HasProperty("reinit") )
- ConvertPropertyTo<bool>("reinit");
+void LC::setProperties( const PropertySet &opts ) {
+ assert( opts.hasKey("tol") );
+ assert( opts.hasKey("maxiter") );
+ assert( opts.hasKey("verbose") );
+ assert( opts.hasKey("cavity") );
+ assert( opts.hasKey("updates") );
- return true;
+ props.tol = opts.getStringAs<double>("tol");
+ props.maxiter = opts.getStringAs<size_t>("maxiter");
+ props.verbose = opts.getStringAs<size_t>("verbose");
+ props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
+ props.updates = opts.getStringAs<Properties::UpdateType>("updates");
+ if( opts.hasKey("cavainame") )
+ props.cavainame = opts.getStringAs<string>("cavainame");
+ if( opts.hasKey("cavaiopts") )
+ props.cavaiopts = opts.getStringAs<PropertySet>("cavaiopts");
+ if( opts.hasKey("reinit") )
+ props.reinit = opts.getStringAs<bool>("reinit");
+ if( opts.hasKey("damping") )
+ props.damping = opts.getStringAs<double>("damping");
+ else
+ props.damping = 0.0;
+}
+
+
+PropertySet LC::getProperties() const {
+ PropertySet opts;
+ opts.Set( "tol", props.tol );
+ opts.Set( "maxiter", props.maxiter );
+ opts.Set( "verbose", props.verbose );
+ opts.Set( "cavity", props.cavity );
+ opts.Set( "updates", props.updates );
+ opts.Set( "cavainame", props.cavainame );
+ opts.Set( "cavaiopts", props.cavaiopts );
+ opts.Set( "reinit", props.reinit );
+ opts.Set( "damping", props.damping );
+ return opts;
+}
+
+
+string LC::printProperties() const {
+ stringstream s( stringstream::out );
+ s << "[";
+ s << "tol=" << props.tol << ",";
+ s << "maxiter=" << props.maxiter << ",";
+ s << "verbose=" << props.verbose << ",";
+ s << "cavity=" << props.cavity << ",";
+ s << "updates=" << props.updates << ",";
+ s << "cavainame=" << props.cavainame << ",";
+ s << "cavaiopts=" << props.cavaiopts << ",";
+ s << "reinit=" << props.reinit << ",";
+ s << "damping=" << props.damping << "]";
+ return s.str();
}
-LC::LC(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
- assert( checkProperties() );
+LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _pancakes(), _cavitydists(), _phis(), _beliefs(), _maxdiff(0.0), _iters(0), props() {
+ setProperties( opts );
// create pancakes
- _pancakes.resize(nrVars());
+ _pancakes.resize( nrVars() );
// create cavitydists
for( size_t i=0; i < nrVars(); i++ )
- _cavitydists.push_back(Factor(delta(i)));
+ _cavitydists.push_back(Factor( delta(i) ));
// create phis
_phis.reserve( nrVars() );
}
// create beliefs
+ _beliefs.reserve( nrVars() );
for( size_t i=0; i < nrVars(); i++ )
_beliefs.push_back(Factor(var(i)));
}
string LC::identify() const {
- stringstream result (stringstream::out);
- result << Name << GetProperties();
- return result.str();
+ return string(Name) + printProperties();
}
}
-double LC::CalcCavityDist (size_t i, const string &name, const Properties &opts) {
+double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) {
Factor Bi;
double maxdiff = 0;
- if( Verbose() >= 2 )
- cout << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).states() << " states)" << endl;
+ if( props.verbose >= 2 )
+ cout << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).nrStates() << " states)" << endl;
- if( Cavity() == CavityType::UNIFORM )
+ if( props.cavity == Properties::CavityType::UNIFORM )
Bi = Factor(delta(i));
else {
InfAlg *cav = newInfAlg( name, *this, opts );
cav->makeCavity( i );
- if( Cavity() == CavityType::FULL )
- Bi = calcMarginal( *cav, cav->delta(i), reInit() );
- else if( Cavity() == CavityType::PAIR )
- Bi = calcMarginal2ndO( *cav, cav->delta(i), reInit() );
- else if( Cavity() == CavityType::PAIR2 ) {
- vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->delta(i), reInit() );
+ if( props.cavity == Properties::CavityType::FULL )
+ Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
+ else if( props.cavity == Properties::CavityType::PAIR )
+ Bi = calcMarginal2ndO( *cav, cav->fg().delta(i), props.reinit );
+ else if( props.cavity == Properties::CavityType::PAIR2 ) {
+ vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->fg().delta(i), props.reinit );
for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
Bi *= pairbeliefs[ij];
- } else if( Cavity() == CavityType::PAIRINT ) {
- Bi = calcMarginal( *cav, cav->delta(i), reInit() );
-
- // Set interactions of order > 2 to zero
- size_t N = delta(i).size();
- Real *p = &(*Bi.p().p().begin());
- x2x::p2logp (N, p);
- x2x::logp2w (N, p);
- x2x::fill (N, p, 2, 0.0);
- x2x::w2logp (N, p);
-// x2x::logpnorm (N, p);
- x2x::logp2p (N, p);
- } else if( Cavity() == CavityType::PAIRCUM ) {
- Bi = calcMarginal( *cav, cav->delta(i), reInit() );
-
- // Set cumulants of order > 2 to zero
- size_t N = delta(i).size();
- Real *p = &(*Bi.p().p().begin());
- x2x::p2m (N, p);
- x2x::m2c (N, p, N);
- x2x::fill (N, p, 2, 0.0);
- x2x::c2m (N, p, N);
- x2x::m2p (N, p);
}
- maxdiff = cav->MaxDiff();
+ maxdiff = cav->maxDiff();
delete cav;
}
- Bi.normalize( _normtype );
+ Bi.normalize();
_cavitydists[i] = Bi;
return maxdiff;
}
-double LC::InitCavityDists (const string &name, const Properties &opts) {
- clock_t tic = toc();
+double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
+ double tic = toc();
- if( Verbose() >= 1 ) {
- cout << "LC::InitCavityDists: ";
- if( Cavity() == CavityType::UNIFORM )
+ if( props.verbose >= 1 ) {
+ cout << Name << "::InitCavityDists: ";
+ if( props.cavity == Properties::CavityType::UNIFORM )
cout << "Using uniform initial cavity distributions" << endl;
- else if( Cavity() == CavityType::FULL )
+ else if( props.cavity == Properties::CavityType::FULL )
cout << "Using full " << name << opts << "...";
- else if( Cavity() == CavityType::PAIR )
+ else if( props.cavity == Properties::CavityType::PAIR )
cout << "Using pairwise " << name << opts << "...";
- else if( Cavity() == CavityType::PAIR2 )
+ else if( props.cavity == Properties::CavityType::PAIR2 )
cout << "Using pairwise(new) " << name << opts << "...";
}
if( md > maxdiff )
maxdiff = md;
}
- init();
- if( Verbose() >= 1 ) {
- cout << "used " << toc() - tic << " clocks." << endl;
+ if( props.verbose >= 1 ) {
+ cout << Name << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
}
return maxdiff;
}
-long LC::SetCavityDists (vector<Factor> &Q) {
- if( Verbose() >= 1 )
- cout << "LC::SetCavityDists: Setting initial cavity distributions" << endl;
+long LC::SetCavityDists( std::vector<Factor> &Q ) {
+ if( props.verbose >= 1 )
+ cout << Name << "::SetCavityDists: Setting initial cavity distributions" << endl;
if( Q.size() != nrVars() )
return -1;
for( size_t i = 0; i < nrVars(); i++ ) {
} else
_cavitydists[i] = Q[i];
}
- init();
return 0;
}
void LC::init() {
for( size_t i = 0; i < nrVars(); ++i )
foreach( const Neighbor &I, nbV(i) )
- if( Updates() == UpdateType::SEQRND )
+ if( props.updates == Properties::UpdateType::SEQRND )
_phis[i][I.iter].randomize();
else
_phis[i][I.iter].fill(1.0);
- for( size_t i = 0; i < nrVars(); i++ ) {
- _pancakes[i] = _cavitydists[i];
-
- foreach( const Neighbor &I, nbV(i) ) {
- _pancakes[i] *= factor(I);
- if( Updates() == UpdateType::SEQRND )
- _pancakes[i] *= _phis[i][I.iter];
- }
-
- _pancakes[i].normalize( _normtype );
-
- CalcBelief(i);
- }
}
Factor A_I;
for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
if( var(i) != *k )
- A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).part_sum( Ivars / var(i) );
+ A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).partSum( Ivars / var(i) );
if( Ivars.size() > 1 )
A_I ^= (1.0 / (Ivars.size() - 1));
- Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).part_sum( Ivars / var(i) );
+ Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).partSum( Ivars / var(i) );
Factor quot = A_I.divided_by(A_Ii);
+ if( props.damping != 0.0 )
+ quot = (quot^(1.0 - props.damping)) * (_phis[i][_I]^props.damping);
- piet *= quot.divided_by( _phis[i][_I] ).normalized( _normtype );
- _phis[i][_I] = quot.normalized( _normtype );
+ piet *= quot.divided_by( _phis[i][_I] ).normalized();
+ _phis[i][_I] = quot.normalized();
- piet.normalize( _normtype );
+ piet.normalize();
if( piet.hasNaNs() ) {
- cout << "LC::NewPancake(" << i << ", " << _I << "): has NaNs!" << endl;
+ cout << Name << "::NewPancake(" << i << ", " << _I << "): has NaNs!" << endl;
hasNaNs = true;
}
double LC::run() {
- if( Verbose() >= 1 )
+ if( props.verbose >= 1 )
cout << "Starting " << identify() << "...";
- if( Verbose() >= 2 )
+ if( props.verbose >= 2 )
cout << endl;
- clock_t tic = toc();
+ double tic = toc();
Diffs diffs(nrVars(), 1.0);
- updateMaxDiff( InitCavityDists(GetPropertyAs<string>("cavainame"), GetPropertyAs<Properties>("cavaiopts")) );
+ double md = InitCavityDists( props.cavainame, props.cavaiopts );
+ if( md > _maxdiff )
+ _maxdiff = md;
+
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ _pancakes[i] = _cavitydists[i];
+
+ foreach( const Neighbor &I, nbV(i) ) {
+ _pancakes[i] *= factor(I);
+ if( props.updates == Properties::UpdateType::SEQRND )
+ _pancakes[i] *= _phis[i][I.iter];
+ }
+
+ _pancakes[i].normalize();
+
+ CalcBelief(i);
+ }
vector<Factor> old_beliefs;
for(size_t i=0; i < nrVars(); i++ )
break;
}
if( hasNaNs ) {
- cout << "LC::run: initial _pancakes has NaNs!" << endl;
- return NAN;
+ cout << Name << "::run: initial _pancakes has NaNs!" << endl;
+ return 1.0;
}
size_t nredges = nrEdges();
- typedef pair<size_t,size_t> Edge;
vector<Edge> update_seq;
update_seq.reserve( nredges );
for( size_t i = 0; i < nrVars(); ++i )
foreach( const Neighbor &I, nbV(i) )
update_seq.push_back( Edge( i, I.iter ) );
- size_t iter = 0;
-
// do several passes over the network until maximum number of iterations has
// been reached or until the maximum belief difference is smaller than tolerance
- for( iter=0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
+ for( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
// Sequential updates
- if( Updates() == UpdateType::SEQRND )
+ if( props.updates == Properties::UpdateType::SEQRND )
random_shuffle( update_seq.begin(), update_seq.end() );
for( size_t t=0; t < nredges; t++ ) {
size_t _I = update_seq[t].second;
_pancakes[i] = NewPancake( i, _I, hasNaNs);
if( hasNaNs )
- return NAN;
+ return 1.0;
CalcBelief( i );
}
old_beliefs[i] = belief(i);
}
- if( Verbose() >= 3 )
- cout << "LC::run: maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
+ if( props.verbose >= 3 )
+ cout << Name << "::run: maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
}
- updateMaxDiff( diffs.max() );
+ if( diffs.maxDiff() > _maxdiff )
+ _maxdiff = diffs.maxDiff();
- if( Verbose() >= 1 ) {
- if( diffs.max() > Tol() ) {
- if( Verbose() == 1 )
+ if( props.verbose >= 1 ) {
+ if( diffs.maxDiff() > props.tol ) {
+ if( props.verbose == 1 )
cout << endl;
- cout << "LC::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+ cout << Name << "::run: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
} else {
- if( Verbose() >= 2 )
- cout << "LC::run: ";
- cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+ if( props.verbose >= 2 )
+ cout << Name << "::run: ";
+ cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
}
}
- return diffs.max();
+ return diffs.maxDiff();
}
projects / libdai.git / blobdiff