Fixed tabs and trailing whitespaces
[libdai.git] / src / lc.cpp
index e3349ce..6fcf73a 100644 (file)
@@ -1,6 +1,7 @@
-/*  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 {
@@ -39,44 +38,71 @@ using namespace std;
 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");
-    
-    return true;
+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") );
+
+    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;
 }
 
 
-LC::LC(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
-    assert( checkProperties() );
+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 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() );
@@ -88,15 +114,14 @@ LC::LC(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
     }
 
     // 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();
+string LC::identify() const {
+    return string(Name) + printProperties();
 }
 
 
@@ -105,74 +130,51 @@ void LC::CalcBelief (size_t i) {
 }
 
 
-double LC::CalcCavityDist (size_t i, const std::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 )
+        cerr << "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 std::string &name, const Properties &opts) {
+double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
     double tic = toc();
 
-    if( Verbose() >= 1 ) {
-        cout << "LC::InitCavityDists:  ";
-        if( Cavity() == CavityType::UNIFORM )
-            cout << "Using uniform initial cavity distributions" << endl;
-        else if( Cavity() == CavityType::FULL )
-            cout << "Using full " << name << opts << "...";
-        else if( Cavity() == CavityType::PAIR )
-            cout << "Using pairwise " << name << opts << "...";
-        else if( Cavity() == CavityType::PAIR2 )
-            cout << "Using pairwise(new) " << name << opts << "...";
+    if( props.verbose >= 1 ) {
+        cerr << Name << "::InitCavityDists:  ";
+        if( props.cavity == Properties::CavityType::UNIFORM )
+            cerr << "Using uniform initial cavity distributions" << endl;
+        else if( props.cavity == Properties::CavityType::FULL )
+            cerr << "Using full " << name << opts << "...";
+        else if( props.cavity == Properties::CavityType::PAIR )
+            cerr << "Using pairwise " << name << opts << "...";
+        else if( props.cavity == Properties::CavityType::PAIR2 )
+            cerr << "Using pairwise(new) " << name << opts << "...";
     }
 
     double maxdiff = 0.0;
@@ -181,10 +183,9 @@ double LC::InitCavityDists (const std::string &name, const Properties &opts) {
         if( md > maxdiff )
             maxdiff = md;
     }
-    init();
 
-    if( Verbose() >= 1 ) {
-        cout << "used " << toc() - tic << " clocks." << endl;
+    if( props.verbose >= 1 ) {
+        cerr << Name << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
     }
 
     return maxdiff;
@@ -192,8 +193,8 @@ double LC::InitCavityDists (const std::string &name, const Properties &opts) {
 
 
 long LC::SetCavityDists( std::vector<Factor> &Q ) {
-    if( Verbose() >= 1 ) 
-        cout << "LC::SetCavityDists:  Setting initial cavity distributions" << endl;
+    if( props.verbose >= 1 )
+        cerr << Name << "::SetCavityDists:  Setting initial cavity distributions" << endl;
     if( Q.size() != nrVars() )
         return -1;
     for( size_t i = 0; i < nrVars(); i++ ) {
@@ -202,7 +203,6 @@ long LC::SetCavityDists( std::vector<Factor> &Q ) {
         } else
             _cavitydists[i] = Q[i];
     }
-    init();
     return 0;
 }
 
@@ -210,23 +210,10 @@ long LC::SetCavityDists( std::vector<Factor> &Q ) {
 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);
-    }
 }
 
 
@@ -239,19 +226,21 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
     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()).marginal( Ivars / var(i), false );
     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 quot = A_I.divided_by(A_Ii);
+    Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).marginal( Ivars / var(i), false );
+    Factor quot = A_I / 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 / _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;
+        cerr << Name << "::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
         hasNaNs = true;
     }
 
@@ -260,15 +249,31 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
 
 
 double LC::run() {
-    if( Verbose() >= 1 )
-        cout << "Starting " << identify() << "...";
-    if( Verbose() >= 2 )
-        cout << endl;
+    if( props.verbose >= 1 )
+        cerr << "Starting " << identify() << "...";
+    if( props.verbose >= 2 )
+        cerr << endl;
 
     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++ )
@@ -281,8 +286,8 @@ double LC::run() {
             break;
         }
     if( hasNaNs ) {
-        cout << "LC::run:  initial _pancakes has NaNs!" << endl;
-        return -1.0;
+        cerr << Name << "::run:  initial _pancakes has NaNs!" << endl;
+        return 1.0;
     }
 
     size_t nredges = nrEdges();
@@ -292,21 +297,19 @@ double LC::run() {
         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.maxDiff() > 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].first;
             size_t _I = update_seq[t].second;
             _pancakes[i] = NewPancake( i, _I, hasNaNs);
             if( hasNaNs )
-                return -1.0;
+                return 1.0;
             CalcBelief( i );
         }
 
@@ -316,21 +319,22 @@ double LC::run() {
             old_beliefs[i] = belief(i);
         }
 
-        if( Verbose() >= 3 )
-            cout << "LC::run:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
+        if( props.verbose >= 3 )
+            cerr << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
     }
 
-    updateMaxDiff( diffs.maxDiff() );
+    if( diffs.maxDiff() > _maxdiff )
+        _maxdiff = diffs.maxDiff();
 
-    if( Verbose() >= 1 ) {
-        if( diffs.maxDiff() > Tol() ) {
-            if( Verbose() == 1 )
-                cout << endl;
-                cout << "LC::run:  WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+    if( props.verbose >= 1 ) {
+        if( diffs.maxDiff() > props.tol ) {
+            if( props.verbose == 1 )
+                cerr << endl;
+                cerr << 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 )
+                cerr << Name << "::run:  ";
+                cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }