-/* 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 <dai/bp.h>
#include <dai/jtree.h>
#include <dai/util.h>
-#include <dai/diffs.h>
namespace dai {
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();
}
}
}
}
- } 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 )
+ if( props.verbose >= 2 )
cout << "init_cor_resp: Convergence not reached (md=" << md << ")..." << endl;
if(md > maxdev)
maxdev = md;
// 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);
}
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);
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);
//
// 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));
}
// 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;
}
}
*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());
}
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);
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];
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 )
+ if( props.verbose >= 1 )
cout << "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;
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++)
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 );
}
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 ) {
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 )
+ if( props.verbose >= 1 )
cout << "Starting " << identify() << "...";
- clock_t tic = toc();
+ double tic = toc();
// Diffs diffs(nrVars(), 1.0);
M.resize(N);
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();
Mag.resize(N);
solveM();
- if( Verbose() >= 1 )
- cout << "MR needed " << toc() - tic << " clocks." << endl;
+ if( props.verbose >= 1 )
+ cout << Name << " needed " << toc() - tic << " seconds." << endl;
return 0.0;
} else
- return NAN;
+ return 1.0;
}
-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++ )
projects / libdai.git / blobdiff