[Patrick Pletscher] Fixed performance issue in FactorGraph::clamp and FactorGraph...
[libdai.git] / src / mr.cpp
index 2c84702..1888b48 100644 (file)
@@ -1,6 +1,10 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
+/*  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
+
+    Copyright (C) 2007  Bastian Wemmenhove
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -27,7 +31,6 @@
 #include <dai/bp.h>
 #include <dai/jtree.h>
 #include <dai/util.h>
-#include <dai/diffs.h>
 
 
 namespace dai {
@@ -39,22 +42,37 @@ using namespace std;
 const char *MR::Name = "MR";
 
 
-bool MR::checkProperties() {
-    if( !HasProperty("updates") )
-        return false;
-    if( !HasProperty("inits") )
-        return false;
-    if( !HasProperty("verbose") )
-        return false;
-    if( !HasProperty("tol") )
-        return false;
+void MR::setProperties( const PropertySet &opts ) {
+    assert( opts.hasKey("tol") );
+    assert( opts.hasKey("verbose") );
+    assert( opts.hasKey("updates") );
+    assert( opts.hasKey("inits") );
     
-    ConvertPropertyTo<UpdateType>("updates");
-    ConvertPropertyTo<InitType>("inits");
-    ConvertPropertyTo<size_t>("verbose");
-    ConvertPropertyTo<double>("tol");
+    props.tol = opts.getStringAs<double>("tol");
+    props.verbose = opts.getStringAs<size_t>("verbose");
+    props.updates = opts.getStringAs<Properties::UpdateType>("updates");
+    props.inits = opts.getStringAs<Properties::InitType>("inits");
+}
+
+
+PropertySet MR::getProperties() const {
+    PropertySet opts;
+    opts.Set( "tol", props.tol );
+    opts.Set( "verbose", props.verbose );
+    opts.Set( "updates", props.updates );
+    opts.Set( "inits", props.inits );
+    return opts;
+}
+
 
-    return true;
+string MR::printProperties() const {
+    stringstream s( stringstream::out );
+    s << "[";
+    s << "tol=" << props.tol << ",";
+    s << "verbose=" << props.verbose << ",";
+    s << "updates=" << props.updates << ",";
+    s << "inits=" << props.inits << "]";
+    return s.str();
 }
 
 
@@ -199,10 +217,10 @@ double MR::init_cor_resp() {
                         }
                     }
                 }
-            } while((md > Tol())&&(runx<runs)); // Precision condition reached -> BP and RP finished
+            } while((md > props.tol)&&(runx<runs)); // Precision condition reached -> BP and RP finished
             if(runx==runs)
-                if( Verbose() >= 2 )
-                    cout << "init_cor_resp: Convergence not reached (md=" << md << ")..." << endl;
+                if( props.verbose >= 2 )
+                    cerr << "init_cor_resp: Convergence not reached (md=" << md << ")..." << endl;
             if(md > maxdev)
                 maxdev = md;
 
@@ -250,31 +268,28 @@ double MR::T(size_t i, sub_nb A) {
     // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
     
     sub_nb _nbi_min_A(con[i]);
-    for( size_t __j = 0; __j < A.size(); __j++ )
-        _nbi_min_A -= A[__j];
+    _nbi_min_A.set();
+    _nbi_min_A &= ~A;
 
     double res = theta[i];
-    for( size_t __j = 0; __j < _nbi_min_A.size(); __j++ ) {
-        size_t _j = _nbi_min_A[__j];
-        res += atanh(tJ[i][_j] * M[i][_j]);
-    }
+    for( size_t _j = 0; _j < _nbi_min_A.size(); _j++ )
+        if( _nbi_min_A.test(_j) )
+            res += atanh(tJ[i][_j] * M[i][_j]);
     return tanh(res);
 }
 
 
 double MR::T(size_t i, size_t _j) {
     sub_nb j(con[i]);
-    j.clear();
-    j += _j;
+    j.set(_j);
     return T(i,j);
 }
 
 
 double MR::Omega(size_t i, size_t _j, size_t _l) {
     sub_nb jl(con[i]);
-    jl.clear();
-    jl += _j;
-    jl += _l;
+    jl.set(_j);
+    jl.set(_l);
     double Tijl = T(i,jl);
     return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl);
 }
@@ -282,11 +297,10 @@ double MR::Omega(size_t i, size_t _j, size_t _l) {
 
 double MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
     sub_nb jll(con[i]);
-    jll.clear();
-    jll += _j;
+    jll.set(_j);
     double Tij = T(i,jll);
-    jll += _l1;
-    jll += _l2;
+    jll.set(_l1);
+    jll.set(_l2);
     double Tijll = T(i,jll);
 
     return (Tijll - Tij) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Tijll + tJ[i][_l2] * M[i][_l2] * Tijll);
@@ -295,10 +309,9 @@ double MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
 
 double MR::Gamma(size_t i, size_t _l1, size_t _l2) {
     sub_nb ll(con[i]);
-    ll.clear();
     double Ti = T(i,ll);
-    ll += _l1;
-    ll += _l2;
+    ll.set(_l1);
+    ll.set(_l2);
     double Till = T(i,ll);
 
     return (Till - Ti) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Till + tJ[i][_l2] * M[i][_l2] * Till);
@@ -311,19 +324,11 @@ double MR::_tJ(size_t i, sub_nb A) {
     //
     // calculate the product of all tJ[i][_j] for _j in A
     
-    size_t Asize = A.size();
-    switch( Asize ) {
-        case 0: return 1.0;
-//      case 1: return tJ[i][A[0]];
-//      case 2: return tJ[i][A[0]] * tJ[i][A[1]];
-//      case 3: return tJ[i][A[0]] * tJ[i][A[1]] * tJ[i][A[2]];
-        default:
-            size_t __j = Asize - 1;
-            size_t _j = A[__j];
-            sub_nb A_j = A - _j;
-            return tJ[i][_j] * _tJ(i, A_j);
-    }
-
+    sub_nb::size_type _j = A.find_first();
+    if( _j == sub_nb::npos )
+        return 1.0;
+    else
+        return tJ[i][_j] * _tJ(i, A.reset(_j));
 }
 
 
@@ -336,24 +341,20 @@ double MR::appM(size_t i, sub_nb A) {
     // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the 
     // entries of the partitions
     
-    size_t Asize = A.size();
-    switch( Asize ) {
-        case 0: return 1.0;
-//      case 1: return M[i][A[0]];
-//      case 2: return M[i][A[0]] * M[i][A[1]] + cors[i][A[0]][A[1]];
-//      case 3: return M[i][A[0]] * M[i][A[1]] * M[i][A[2]] + M[i][A[0]] * cors[i][A[1]][A[2]] + M[i][A[1]] * cors[i][A[0]][A[2]] + M[i][A[2]] * cors[i][A[0]][A[1]];
-        default:
-            size_t __j = Asize - 1;
-            size_t _j = A[__j];
-            sub_nb A_j = A - _j;
-
-            double result = M[i][_j] * appM(i, A_j);
-            for( size_t __k = 0; __k < __j; __k++ ) {
-                size_t _k = A[__k];
-                result += cors[i][_j][_k] * appM(i,A_j - _k);
+    sub_nb::size_type _j = A.find_first();
+    if( _j == sub_nb::npos )
+        return 1.0;
+    else {
+        sub_nb A_j(A); A_j.reset(_j);
+
+        double result = M[i][_j] * appM(i, A_j);
+        for( size_t _k = 0; _k < A_j.size(); _k++ )
+            if( A_j.test(_k) ) {
+                sub_nb A_jk(A_j); A_jk.reset(_k);
+                result += cors[i][_j][_k] * appM(i,A_jk);
             }
 
-            return result;
+        return result;
     }
 }
 
@@ -367,24 +368,25 @@ void MR::sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd) {
     *sum_even = 0.0;
     *sum_odd = 0.0;
 
-    sub_nb S(A.size());
-    S.clear();
-    do { // for all subsets of A
-
-        // construct subset B of A corresponding to S
-        sub_nb B = A;
-        B.clear();
-        size_t Ssize = S.size();
-        for( size_t bit = 0; bit < Ssize; bit++ )
-            B += A[S[bit]];
-        
-        if( S.size() % 2 )
+    sub_nb B(A.size());
+    do {
+        if( B.count() % 2 )
             *sum_odd += _tJ(j,B) * appM(j,B);
         else
             *sum_even += _tJ(j,B) * appM(j,B);
 
-        ++S;
-    } while( !S.empty() );
+        // calc next subset B
+        size_t bit = 0;
+        for( ; bit < A.size(); bit++ )
+            if( A.test(bit) ) {
+                if( B.test(bit) )
+                    B.reset(bit);
+                else {
+                    B.set(bit);
+                    break;
+                }
+            }
+    } while (!B.none());
 }
             
 
@@ -409,14 +411,16 @@ void MR::solvemcav() {
                 assert( nb[j][_i] == i );
 
                 double newM = 0.0;
-                if( Updates() == UpdateType::FULL ) {
+                if( props.updates == Properties::UpdateType::FULL ) {
                     // find indices in nb[j] that do not correspond with i
                     sub_nb _nbj_min_i(con[j]);
-                    _nbj_min_i -= kindex[i][_j];
+                    _nbj_min_i.set();
+                    _nbj_min_i.reset(kindex[i][_j]);
 
                     // find indices in nb[i] that do not correspond with j
                     sub_nb _nbi_min_j(con[i]);
-                    _nbi_min_j -= _j;
+                    _nbi_min_j.set();
+                    _nbi_min_j.reset(_j);
 
                     sum_subs(j, _nbj_min_i, &sum_even, &sum_odd);
                     newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd);
@@ -425,11 +429,13 @@ void MR::solvemcav() {
                     double denom = sum_even + tanh(theta[i]) * sum_odd;
                     double numer = 0.0;
                     for(size_t _k=0; _k<con[i]; _k++) if(_k != _j) {
-                        sum_subs(i, _nbi_min_j - _k, &sum_even, &sum_odd);
+                        sub_nb _nbi_min_jk(_nbi_min_j);
+                        _nbi_min_jk.reset(_k);
+                        sum_subs(i, _nbi_min_jk, &sum_even, &sum_odd);
                         numer += tJ[i][_k] * cors[i][_j][_k] * (tanh(theta[i]) * sum_even + sum_odd);
                     }
                     newM -= numer / denom;
-                } else if( Updates() == UpdateType::LINEAR ) {
+                } else if( props.updates == Properties::UpdateType::LINEAR ) {
                     newM = T(j,_i);
                     for(size_t _l=0; _l<con[i]; _l++) if( _l != _j )
                         newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
@@ -449,22 +455,25 @@ void MR::solvemcav() {
                 M[i][_j] = newM;
             }
         }
-    } while((maxdev>Tol())&&(run<maxruns));
+    } while((maxdev>props.tol)&&(run<maxruns));
 
-    updateMaxDiff( maxdev );
+    _iters = run;
+    if( maxdev > _maxdiff )
+        _maxdiff = maxdev;
 
     if(run==maxruns){
-        if( Verbose() >= 1 )
-            cout << "solve_mcav: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
+        if( props.verbose >= 1 )
+            cerr << "solve_mcav: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
     }
 }
 
  
 void MR::solveM() { 
     for(size_t i=0; i<N; i++) {
-        if( Updates() == UpdateType::FULL ) {
+        if( props.updates == Properties::UpdateType::FULL ) {
             // find indices in nb[i]
             sub_nb _nbi(con[i]);
+            _nbi.set();
 
             // calc numerator1 and denominator1
             double sum_even, sum_odd;
@@ -472,9 +481,8 @@ void MR::solveM() {
 
             Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
 
-        } else if( Updates() == UpdateType::LINEAR ) {
+        } else if( props.updates == Properties::UpdateType::LINEAR ) {
             sub_nb empty(con[i]);
-            empty.clear();
             Mag[i] = T(i,empty);
 
             for(size_t _l1=0; _l1<con[i]; _l1++)
@@ -490,12 +498,12 @@ void MR::solveM() {
 void MR::init_cor() {
     for( size_t i = 0; i < nrVars(); i++ ) {
         vector<Factor> pairq;
-        if( Inits() == InitType::CLAMPING ) {
-            BP bpcav(*this, Properties()("updates",string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("1000UL"))("verbose", string("0UL"))("logdomain", string("0")));
+        if( props.inits == Properties::InitType::CLAMPING ) {
+            BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("10000UL"))("verbose", string("0UL"))("logdomain", string("0")));
             bpcav.makeCavity( i );
             pairq = calcPairBeliefs( bpcav, delta(i), false );
-        } else if( Inits() == InitType::EXACT ) {
-            JTree jtcav(*this, Properties()("updates",string("HUGIN"))("verbose", string("0UL")) );
+        } else if( props.inits == Properties::InitType::EXACT ) {
+            JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", string("0UL")) );
             jtcav.makeCavity( i );
             pairq = calcPairBeliefs( jtcav, delta(i), false );
         }
@@ -503,7 +511,7 @@ void MR::init_cor() {
             VarSet::const_iterator kit = pairq[jk].vars().begin();
             size_t j = findVar( *(kit) );
             size_t k = findVar( *(++kit) );
-            pairq[jk].normalize(Prob::NORMPROB);
+            pairq[jk].normalize();
             double cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
             for( size_t _j = 0; _j < con[i]; _j++ ) if( nb[i][_j] == j )
                 for( size_t _k = 0; _k < con[i]; _k++ ) if( nb[i][_k] == k ) {
@@ -516,16 +524,14 @@ void MR::init_cor() {
 
 
 string MR::identify() const { 
-    stringstream result (stringstream::out);
-    result << Name << GetProperties();
-    return result.str();
+    return string(Name) + printProperties();
 }
 
 
 double MR::run() {
     if( supported ) {
-        if( Verbose() >= 1 )
-            cout << "Starting " << identify() << "...";
+        if( props.verbose >= 1 )
+            cerr << "Starting " << identify() << "...";
 
         double tic = toc();
 //        Diffs diffs(nrVars(), 1.0);
@@ -545,11 +551,13 @@ double MR::run() {
         for(size_t i=0; i<N; i++)
           kindex[i].resize(kmax);
 
-        if( Inits() == InitType::RESPPROP )
-            updateMaxDiff( init_cor_resp() );
-        else if( Inits() == InitType::EXACT )
+        if( props.inits == Properties::InitType::RESPPROP ) {
+            double md = init_cor_resp();
+            if( md > _maxdiff )
+                _maxdiff = md;
+        } else if( props.inits == Properties::InitType::EXACT )
             init_cor(); // FIXME no MaxDiff() calculation
-        else if( Inits() == InitType::CLAMPING )
+        else if( props.inits == Properties::InitType::CLAMPING )
             init_cor(); // FIXME no MaxDiff() calculation
 
         solvemcav();
@@ -557,12 +565,12 @@ double MR::run() {
         Mag.resize(N);
         solveM();
 
-        if( Verbose() >= 1 )
-            cout << "MR needed " << toc() - tic << " clocks." << endl;
+        if( props.verbose >= 1 )
+            cerr << Name << " needed " << toc() - tic << " seconds." << endl;
 
         return 0.0;
     } else
-        return -1.0;
+        return 1.0;
 }
 
 
@@ -601,7 +609,9 @@ vector<Factor> MR::beliefs() const {
 
 
 
-MR::MR( const FactorGraph &fg, const Properties &opts ) : DAIAlgFG(fg, opts), supported(true) {
+MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
+    setProperties( opts );
+
     // check whether all vars in fg are binary
     // check whether connectivity is <= kmax
     for( size_t i = 0; i < fg.nrVars(); i++ )