# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+# Directories
+INC = include/dai
+SRC = src
+LIB = lib
+
# We use the BOOST Program Options library
BOOSTFLAGS = -lboost_program_options
CC = g++
# Flags for the C++ compiler
-CCFLAGS = -Wall -fpic -g -DDEBUG -I. -I$(HOME)/include -O3 #-static #-pg #-DVERBOSE
+CCFLAGS = -Wall -W -Wextra -fpic -g -DDEBUG -I./include -Llib -O3 #-static #-pg #-DVERBOSE
# To enable the Matlab interface, define WITH_MATLAB = yes
WITH_MATLAB =
# Replace the following with the extension of compiled MEX files on this platform, e.g. .mexglx for x86
MEXEXT = .mexglx
-HEADERS = bipgraph.h diffs.h index.h var.h factor.h varset.h prob.h daialg.h properties.h alldai.h enum.h x2x.h
+HEADERS = $(INC)/bipgraph.h $(INC)/diffs.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/x2x.h
# target matlabs is disabled by default since it only compiles with a very recent MatLab version
-TARGETS = tests utils libdai.a example testregression
+TARGETS = tests utils $(LIB)/libdai.a example testregression
ifdef WITH_MATLAB
TARGETS := $(TARGETS) matlabs
endif
matlabs : matlab/dai.$(MEXEXT) matlab/dai_readfg.$(MEXEXT) matlab/dai_writefg.$(MEXEXT) matlab/dai_removeshortloops.$(MEXEXT) matlab/dai_potstrength.$(MEXEXT)
-libdai.a : daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
- ar rcs libdai.a daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+$(LIB)/libdai.a : daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+ ar rcs $(LIB)/libdai.a daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
tests : tests/test
echo Testing...this can take a while...
cd tests; ./testregression; cd ..
-doc : *.h *.cpp doxygen.conf
+doc : $(INC)/*.h $(SRC)/*.cpp doxygen.conf
doxygen doxygen.conf
clean :
- rm *.o *.$(MEXEXT) example matlab/*.$(MEXEXT) matlab/*.o tests/test utils/fg2dot utils/createfg utils/remove_short_loops utils/fginfo libdai.a; echo
+ rm *.o *.$(MEXEXT) example matlab/*.$(MEXEXT) matlab/*.o tests/test utils/fg2dot utils/createfg utils/remove_short_loops utils/fginfo $(LIB)/libdai.a; echo
rm -R doc; echo
-daialg.o : daialg.cpp $(HEADERS)
- $(CC) $(CCFLAGS) -c daialg.cpp
+daialg.o : $(SRC)/daialg.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/daialg.cpp
-bp.o : bp.cpp bp.h $(HEADERS)
- $(CC) $(CCFLAGS) -c bp.cpp
+bp.o : $(SRC)/bp.cpp $(INC)/bp.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/bp.cpp
-lc.o : lc.cpp lc.h $(HEADERS)
- $(CC) $(CCFLAGS) -c lc.cpp
+lc.o : $(SRC)/lc.cpp $(INC)/lc.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/lc.cpp
-mf.o : mf.cpp mf.h $(HEADERS)
- $(CC) $(CCFLAGS) -c mf.cpp
+mf.o : $(SRC)/mf.cpp $(INC)/mf.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/mf.cpp
-factorgraph.o : factorgraph.cpp factorgraph.h $(HEADERS)
- $(CC) $(CCFLAGS) -c factorgraph.cpp
+factorgraph.o : $(SRC)/factorgraph.cpp $(INC)/factorgraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/factorgraph.cpp
-util.o : util.cpp util.h $(HEADERS)
- $(CC) $(CCFLAGS) -c util.cpp
+util.o : $(SRC)/util.cpp $(INC)/util.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/util.cpp
-regiongraph.o : regiongraph.cpp regiongraph.h $(HEADERS)
- $(CC) $(CCFLAGS) -c regiongraph.cpp
+regiongraph.o : $(SRC)/regiongraph.cpp $(INC)/regiongraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/regiongraph.cpp
-hak.o : hak.cpp hak.h $(HEADERS) regiongraph.h
- $(CC) $(CCFLAGS) -c hak.cpp
+hak.o : $(SRC)/hak.cpp $(INC)/hak.h $(HEADERS) $(INC)/regiongraph.h
+ $(CC) $(CCFLAGS) -c $(SRC)/hak.cpp
-clustergraph.o : clustergraph.cpp clustergraph.h $(HEADERS)
- $(CC) $(CCFLAGS) -c clustergraph.cpp
+clustergraph.o : $(SRC)/clustergraph.cpp $(INC)/clustergraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/clustergraph.cpp
-jtree.o : jtree.cpp jtree.h $(HEADERS) weightedgraph.h clustergraph.h regiongraph.h
- $(CC) $(CCFLAGS) -c jtree.cpp
+jtree.o : $(SRC)/jtree.cpp $(INC)/jtree.h $(HEADERS) $(INC)/weightedgraph.h $(INC)/clustergraph.h $(INC)/regiongraph.h
+ $(CC) $(CCFLAGS) -c $(SRC)/jtree.cpp
-treeep.o : treeep.cpp treeep.h $(HEADERS) weightedgraph.h clustergraph.h regiongraph.h jtree.h
- $(CC) $(CCFLAGS) -c treeep.cpp
+treeep.o : $(SRC)/treeep.cpp $(INC)/treeep.h $(HEADERS) $(INC)/weightedgraph.h $(INC)/clustergraph.h $(INC)/regiongraph.h $(INC)/jtree.h
+ $(CC) $(CCFLAGS) -c $(SRC)/treeep.cpp
-weightedgraph.o : weightedgraph.cpp weightedgraph.h $(HEADERS)
- $(CC) $(CCFLAGS) -c weightedgraph.cpp
+weightedgraph.o : $(SRC)/weightedgraph.cpp $(INC)/weightedgraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/weightedgraph.cpp
-mr.o : mr.cpp mr.h $(HEADERS)
- $(CC) $(CCFLAGS) -c mr.cpp
+mr.o : $(SRC)/mr.cpp $(INC)/mr.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/mr.cpp
-properties.o : properties.cpp $(HEADERS)
- $(CC) $(CCFLAGS) -c properties.cpp
+properties.o : $(SRC)/properties.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/properties.cpp
-alldai.o : alldai.cpp $(HEADERS)
- $(CC) $(CCFLAGS) -c alldai.cpp
+alldai.o : $(SRC)/alldai.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/alldai.cpp
-x2x.o : x2x.cpp $(HEADERS)
- $(CC) $(CCFLAGS) -c x2x.cpp
+x2x.o : $(SRC)/x2x.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c $(SRC)/x2x.cpp
# EXAMPLE
##########
-example : example.cpp $(HEADERS) libdai.a
- $(CC) $(CCFLAGS) -o example example.cpp -L. libdai.a
+example : $(SRC)/example.cpp $(HEADERS) $(LIB)/libdai.a
+ $(CC) $(CCFLAGS) -o example $(SRC)/example.cpp -ldai
# TESTS
########
-tests/test : tests/test.cpp $(HEADERS) libdai.a
- $(CC) $(CCFLAGS) -o tests/test tests/test.cpp -L. -ldai $(BOOSTFLAGS)
+tests/test : tests/test.cpp $(HEADERS) lib/libdai.a
+ $(CC) $(CCFLAGS) -o tests/test tests/test.cpp -ldai $(BOOSTFLAGS)
# MATLAB INTERFACE
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include "alldai.h"
-
-
-namespace dai {
-
-
-InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const Properties &opts ) {
- if( name == BP::Name )
- return new BP (fg, opts);
- else if( name == MF::Name )
- return new MF (fg, opts);
- else if( name == HAK::Name )
- return new HAK (fg, opts);
- else if( name == LC::Name )
- return new LC (fg, opts);
- else if( name == TreeEP::Name )
- return new TreeEP (fg, opts);
- else if( name == MR::Name )
- return new MR (fg, opts);
- else if( name == JTree::Name )
- return new JTree (fg, opts);
- else
- throw "Unknown inference algorithm";
-}
-
-
-}
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __ALLDAI_H__
-#define __ALLDAI_H__
-
-
-#include "daialg.h"
-#include "bp.h"
-#include "lc.h"
-#include "hak.h"
-#include "mf.h"
-#include "jtree.h"
-#include "treeep.h"
-#include "mr.h"
-
-
-namespace dai {
-
-
-/// newInfAlg constructs a new approximate inference algorithm named name for the
-/// FactorGraph fg with optionts opts and returns a pointer to the new object.
-/// The caller needs to delete it (maybe some sort of smart_ptr might be useful here).
-InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const Properties &opts );
-
-
-/// AINames contains the names of all approximate inference algorithms
-static const char* DAINames[] = {BP::Name, MF::Name, HAK::Name, LC::Name, TreeEP::Name, MR::Name, JTree::Name};
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __BIPGRAPH_H__
-#define __BIPGRAPH_H__
-
-
-#include <vector>
-#include <algorithm>
-#include <boost/numeric/ublas/vector_sparse.hpp>
-
-
-namespace dai {
-
-
-/// A BipartiteGraph represents a graph with two types of nodes
-template <class T1, class T2> class BipartiteGraph {
- public:
- typedef std::pair<size_t,size_t> _edge_t;
- typedef std::vector<size_t> _nb_t;
- typedef _nb_t::const_iterator _nb_cit;
-
-
- private:
- /// Vertices of first type
- std::vector<T1> _V1;
-
- /// Vertices of second type
- std::vector<T2> _V2;
-
- /// Edges, which are pairs (v1,v2) with v1 in _V1 and v2 in _V2
- std::vector<_edge_t> _E12;
-
- /// Conversion matrix that computes which index of _E12 corresponds to a vertex index pair (v1,v2)
- std::vector<boost::numeric::ublas::mapped_vector<size_t> > _E12ind;
-
- /// Neighbour indices of vertices of first type
- std::vector<_nb_t> _nb1;
-
- /// Neighbour indices of vertices of second type
- std::vector<_nb_t> _nb2;
-
-
- public:
- /// Default constructor
- BipartiteGraph<T1,T2> () {};
-
- /// Copy constructor
- BipartiteGraph<T1,T2> ( const BipartiteGraph<T1,T2> & x ) : _V1(x._V1), _V2(x._V2), _E12(x._E12), _E12ind(x._E12ind), _nb1(x._nb1), _nb2(x._nb2) {};
-
- /// Assignment operator
- BipartiteGraph<T1,T2> & operator=(const BipartiteGraph<T1,T2> & x) {
- if( this != &x ) {
- _V1 = x._V1;
- _V2 = x._V2;
- _E12 = x._E12;
- _E12ind = x._E12ind;
- _nb1 = x._nb1;
- _nb2 = x._nb2;
- }
- return *this;
- }
-
- /// Provides read access to node of first type
- const T1 & V1( size_t i1 ) const { return _V1[i1]; }
- /// Provides full access to node of first type
- T1 & V1( size_t i1 ) { return _V1[i1]; }
- /// Provides read access to all nodes of first type
- const std::vector<T1> & V1s() const { return _V1; }
- /// Provides full access to all nodes of first type
- std::vector<T1> & V1s() { return _V1; }
-
- /// Provides read access to node of second type
- const T2 & V2( size_t i2 ) const { return _V2[i2]; }
- /// Provides full access to node of second type
- T2 & V2( size_t i2 ) { return _V2[i2]; }
- /// Provides read access to all nodes of second type
- const std::vector<T2> & V2s() const { return _V2; }
- /// Provides full access to all nodes of second type
- std::vector<T2> & V2s() { return _V2; }
-
- /// Provides read access to edge
- const _edge_t & edge(size_t ind) const { return _E12[ind]; }
- /// Provides full access to edge
- _edge_t & edge(size_t ind) { return _E12[ind]; }
- /// Provides read access to all edges
- const std::vector<_edge_t> & edges() const { return _E12; }
- /// Provides full access to all edges
- std::vector<_edge_t> & edges() { return _E12; }
- /// Returns number of edges
- size_t nr_edges() const { return _E12.size(); }
-
- /// Provides read access to neighbours of node of first type
- const _nb_t & nb1( size_t i1 ) const { return _nb1[i1]; }
- /// Provides full access to neighbours of node of first type
- _nb_t & nb1( size_t i1 ) { return _nb1[i1]; }
-
- /// Provides read access to neighbours of node of second type
- const _nb_t & nb2( size_t i2 ) const { return _nb2[i2]; }
- /// Provides full access to neighbours of node of second type
- _nb_t & nb2( size_t i2 ) { return _nb2[i2]; }
-
- /// Converts the pair of indices (i1,i2) to the corresponding edge index
- size_t VV2E( const size_t i1, const size_t i2 ) const { return _E12ind[i1](i2); }
-
-
- /// Regenerates internal structures (_E12ind, _nb1 and _nb2) based upon _V1, _V2 and _E12
- void Regenerate() {
- // Calculate _nb1 and _nb2
-
- // Start with empty vectors
- _nb1.clear();
- _nb1.resize(_V1.size());
- // Start with empty vectors
- _nb2.clear();
- _nb2.resize(_V2.size());
- // Each edge yields a neighbour pair
- for( std::vector<_edge_t>::const_iterator e = _E12.begin(); e != _E12.end(); e++ ) {
- _nb1[e->first].push_back(e->second);
- _nb2[e->second].push_back(e->first);
- }
- // Remove duplicates from _nb1
- for( size_t i1 = 0; i1 < _V1.size(); i1++ ) {
- _nb_t::iterator new_end = unique(_nb1[i1].begin(), _nb1[i1].end());
- _nb1[i1].erase( new_end, _nb1[i1].end() );
- }
- // Remove duplicates from _nb2
- for( size_t i2 = 0; i2 < _V2.size(); i2++ ) {
- _nb_t::iterator new_end = unique(_nb2[i2].begin(), _nb2[i2].end());
- _nb2[i2].erase( new_end, _nb2[i2].end() );
- }
-
- // Calculate _E12ind
- _E12ind = std::vector<boost::numeric::ublas::mapped_vector<size_t> >(_V1.size(), boost::numeric::ublas::mapped_vector<size_t>(_V2.size()) );
- for( size_t k = 0; k < _E12.size(); k++ )
- _E12ind[_E12[k].first](_E12[k].second) = k;
-
- }
-};
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include <iostream>
-#include <sstream>
-#include <map>
-#include <set>
-#include <algorithm>
-#include "bp.h"
-#include "diffs.h"
-#include "util.h"
-#include "properties.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-const char *BP::Name = "BP";
-
-
-bool BP::checkProperties() {
- if( !HasProperty("updates") )
- return false;
- if( !HasProperty("tol") )
- return false;
- if (!HasProperty("maxiter") )
- return false;
- if (!HasProperty("verbose") )
- return false;
-
- ConvertPropertyTo<double>("tol");
- ConvertPropertyTo<size_t>("maxiter");
- ConvertPropertyTo<size_t>("verbose");
- ConvertPropertyTo<UpdateType>("updates");
-
- return true;
-}
-
-
-void BP::Regenerate() {
- DAIAlgFG::Regenerate();
-
- // clear messages
- _messages.clear();
- _messages.reserve(nr_edges());
-
- // clear indices
- _indices.clear();
- _indices.reserve(nr_edges());
-
- // create messages and indices
- for( vector<_edge_t>::const_iterator iI=edges().begin(); iI!=edges().end(); iI++ ) {
- _messages.push_back( Prob( var(iI->first).states() ) );
-
- vector<size_t> ind( factor(iI->second).stateSpace(), 0 );
- Index i (var(iI->first), factor(iI->second).vars() );
- for( size_t j = 0; i >= 0; ++i,++j )
- ind[j] = i;
- _indices.push_back( ind );
- }
-
- // create new_messages
- _newmessages = _messages;
-}
-
-
-void BP::init() {
- assert( checkProperties() );
- for( vector<Prob>::iterator mij = _messages.begin(); mij != _messages.end(); mij++ )
- mij->fill(1.0 / mij->size());
- _newmessages = _messages;
-}
-
-
-void BP::calcNewMessage (size_t iI) {
- // calculate updated message I->i
- size_t i = edge(iI).first;
- size_t I = edge(iI).second;
-
-/* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
-
- Factor prod( factor( I ) );
- for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ )
- if( *j != i ) { // for all j in I \ i
- for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); J++ )
- if( *J != I ) { // for all J in nb(j) \ I
- prod *= Factor( *j, message(*j,*J) );
- Factor marg = prod.marginal(var(i));
-*/
-
- Prob prod( factor(I).p() );
-
- // Calculate product of incoming messages and factor I
- for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ )
- if( *j != i ) { // for all j in I \ i
- // ind is the precalculated Index(j,I) i.e. to x_I == k corresponds x_j == ind[k]
- _ind_t* ind = &(index(*j,I));
-
- // prod_j will be the product of messages coming into j
- Prob prod_j( var(*j).states() );
- for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); J++ )
- if( *J != I ) // for all J in nb(j) \ I
- prod_j *= message(*j,*J);
-
- // multiply prod with prod_j
- for( size_t r = 0; r < prod.size(); r++ )
- prod[r] *= prod_j[(*ind)[r]];
- }
-
- // Marginalize onto i
- Prob marg( var(i).states(), 0.0 );
- // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
- _ind_t* ind = &(index(i,I));
- for( size_t r = 0; r < prod.size(); r++ )
- marg[(*ind)[r]] += prod[r];
- marg.normalize( _normtype );
-
- // Store result
- _newmessages[iI] = marg;
-}
-
-
-// BP::run does not check for NANs for performance reasons
-// Somehow NaNs do not often occur in BP...
-double BP::run() {
- if( Verbose() >= 1 )
- cout << "Starting " << identify() << "...";
- if( Verbose() >= 3)
- cout << endl;
-
- clock_t tic = toc();
- Diffs diffs(nrVars(), 1.0);
-
- vector<size_t> edge_seq;
- vector<double> residuals;
-
- vector<Factor> old_beliefs;
- old_beliefs.reserve( nrVars() );
- for( size_t i = 0; i < nrVars(); i++ )
- old_beliefs.push_back(belief1(i));
-
- size_t iter = 0;
-
- if( Updates() == UpdateType::SEQMAX ) {
- // do the first pass
- for(size_t iI = 0; iI < nr_edges(); iI++ )
- calcNewMessage(iI);
-
- // calculate initial residuals
- residuals.reserve(nr_edges());
- for( size_t iI = 0; iI < nr_edges(); iI++ )
- residuals.push_back( dist( _newmessages[iI], _messages[iI], Prob::DISTLINF ) );
- } else {
- edge_seq.reserve( nr_edges() );
- for( size_t i = 0; i < nr_edges(); i++ )
- edge_seq.push_back( i );
- }
-
- // 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++ ) {
- if( Updates() == UpdateType::SEQMAX ) {
- // Residuals-BP by Koller et al.
- for( size_t t = 0; t < nr_edges(); t++ ) {
- // update the message with the largest residual
- size_t iI = max_element(residuals.begin(), residuals.end()) - residuals.begin();
- _messages[iI] = _newmessages[iI];
- residuals[iI] = 0;
-
- // I->i has been updated, which means that residuals for all
- // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
- size_t i = edge(iI).first;
- size_t I = edge(iI).second;
- for( _nb_cit J = nb1(i).begin(); J != nb1(i).end(); J++ )
- if( *J != I )
- for( _nb_cit j = nb2(*J).begin(); j != nb2(*J).end(); j++ )
- if( *j != i ) {
- size_t jJ = VV2E(*j,*J);
- calcNewMessage(jJ);
- residuals[jJ] = dist( _newmessages[jJ], _messages[jJ], Prob::DISTLINF );
- }
- }
- } else if( Updates() == UpdateType::PARALL ) {
- // Parallel updates
- for( size_t t = 0; t < nr_edges(); t++ )
- calcNewMessage(t);
-
- for( size_t t = 0; t < nr_edges(); t++ )
- _messages[t] = _newmessages[t];
- } else {
- // Sequential updates
- if( Updates() == UpdateType::SEQRND )
- random_shuffle( edge_seq.begin(), edge_seq.end() );
-
- for( size_t t = 0; t < nr_edges(); t++ ) {
- size_t k = edge_seq[t];
- calcNewMessage(k);
- _messages[k] = _newmessages[k];
- }
- }
-
- // calculate new beliefs and compare with old ones
- for( size_t i = 0; i < nrVars(); i++ ) {
- Factor nb( belief1(i) );
- diffs.push( dist( nb, old_beliefs[i], Prob::DISTLINF ) );
- old_beliefs[i] = nb;
- }
-
- if( Verbose() >= 3 )
- cout << "BP::run: maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
- }
-
- updateMaxDiff( diffs.max() );
-
- if( Verbose() >= 1 ) {
- if( diffs.max() > Tol() ) {
- if( Verbose() == 1 )
- cout << endl;
- cout << "BP::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
- } else {
- if( Verbose() >= 3 )
- cout << "BP::run: ";
- cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
- }
- }
-
- return diffs.max();
-}
-
-
-Factor BP::belief1( size_t i ) const {
- Prob prod( var(i).states() );
- for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ )
- prod *= newMessage(i,*I);
-
- prod.normalize( Prob::NORMPROB );
- return( Factor( var(i), prod ) );
-}
-
-
-Factor BP::belief (const Var &n) const {
- return( belief1( findVar( n ) ) );
-}
-
-
-vector<Factor> BP::beliefs() const {
- vector<Factor> result;
- for( size_t i = 0; i < nrVars(); i++ )
- result.push_back( belief1(i) );
- for( size_t I = 0; I < nrFactors(); I++ )
- result.push_back( belief2(I) );
- return result;
-}
-
-
-Factor BP::belief( const VarSet &ns ) const {
- if( ns.size() == 1 )
- return belief( *(ns.begin()) );
- else {
- size_t I;
- for( I = 0; I < nrFactors(); I++ )
- if( factor(I).vars() >> ns )
- break;
- assert( I != nrFactors() );
- return belief2(I).marginal(ns);
- }
-}
-
-
-Factor BP::belief2 (size_t I) const {
- Prob prod( factor(I).p() );
-
- for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ ) {
- // ind is the precalculated Index(j,I) i.e. to x_I == k corresponds x_j == ind[k]
- const _ind_t *ind = &(index(*j, I));
-
- // prod_j will be the product of messages coming into j
- Prob prod_j( var(*j).states() );
- for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); J++ )
- if( *J != I ) // for all J in nb(j) \ I
- prod_j *= newMessage(*j,*J);
-
- // multiply prod with prod_j
- for( size_t r = 0; r < prod.size(); r++ )
- prod[r] *= prod_j[(*ind)[r]];
- }
-
- Factor result( factor(I).vars(), prod );
- result.normalize( Prob::NORMPROB );
-
- return( result );
-
-/* UNOPTIMIZED VERSION
-
- Factor prod( factor(I) );
- for( _nb_cit i = nb2(I).begin(); i != nb2(I).end(); i++ ) {
- for( _nb_cit J = nb1(*i).begin(); J != nb1(*i).end(); J++ )
- if( *J != I )
- prod *= Factor( var(*i), newMessage(*i,*J)) );
- }
- return prod.normalize( Prob::NORMPROB );*/
-}
-
-
-Complex BP::logZ() const {
- Complex sum = 0.0;
- for(size_t i = 0; i < nrVars(); i++ )
- sum += Complex(1.0 - nb1(i).size()) * belief1(i).entropy();
- for( size_t I = 0; I < nrFactors(); I++ )
- sum -= KL_dist( belief2(I), factor(I) );
- return sum;
-}
-
-
-string BP::identify() const {
- stringstream result (stringstream::out);
- result << Name << GetProperties();
- return result.str();
-}
-
-
-void BP::init( const VarSet &ns ) {
- for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) {
- size_t ni = findVar( *n );
- for( _nb_cit I = nb1(ni).begin(); I != nb1(ni).end(); I++ )
- message(ni,*I).fill( 1.0 );
- }
-}
-
-
-}
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __BP_H__
-#define __BP_H__
-
-
-#include "daialg.h"
-#include "factorgraph.h"
-#include "enum.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-class BP : public DAIAlgFG {
- protected:
- typedef vector<size_t> _ind_t;
-
- vector<_ind_t> _indices;
- vector<Prob> _messages, _newmessages;
-
- public:
- ENUM4(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
- UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
-
- // default constructor
- BP() : DAIAlgFG() {};
- // copy constructor
- BP(const BP & x) : DAIAlgFG(x), _indices(x._indices), _messages(x._messages), _newmessages(x._newmessages) {};
- BP* clone() const { return new BP(*this); }
- // construct BP object from FactorGraph
- BP(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
- assert( checkProperties() );
- Regenerate();
- }
- // assignment operator
- BP & operator=(const BP & x) {
- if(this!=&x) {
- DAIAlgFG::operator=(x);
- _messages = x._messages;
- _newmessages = x._newmessages;
- _indices = x._indices;
- }
- return *this;
- }
-
- static const char *Name;
-
- Prob & message(size_t i1, size_t i2) { return( _messages[VV2E(i1,i2)] ); }
- const Prob & message(size_t i1, size_t i2) const { return( _messages[VV2E(i1,i2)] ); }
- Prob & newMessage(size_t i1, size_t i2) { return( _newmessages[VV2E(i1,i2)] ); }
- const Prob & newMessage(size_t i1, size_t i2) const { return( _newmessages[VV2E(i1,i2)] ); }
- _ind_t & index(size_t i1, size_t i2) { return( _indices[VV2E(i1,i2)] ); }
- const _ind_t & index(size_t i1, size_t i2) const { return( _indices[VV2E(i1,i2)] ); }
-
- string identify() const;
- void Regenerate();
- void init();
- void calcNewMessage(size_t iI);
- double run();
- Factor belief1 (size_t i) const;
- Factor belief2 (size_t I) const;
- Factor belief (const Var &n) const;
- Factor belief (const VarSet &n) const;
- vector<Factor> beliefs() const;
- Complex logZ() const;
-
- void init( const VarSet &ns );
- void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
- bool checkProperties();
-};
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include <set>
-#include <vector>
-#include <iostream>
-#include "varset.h"
-#include "clustergraph.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-ClusterGraph ClusterGraph::VarElim( const vector<Var> & ElimSeq ) const {
- const long verbose = 0;
-
- // Make a copy
- ClusterGraph _Cl(*this);
-
- ClusterGraph result;
- _Cl.eraseNonMaximal();
-
- // Do variable elimination
- for( vector<Var>::const_iterator n = ElimSeq.begin(); n != ElimSeq.end(); n++ ) {
- assert( _Cl.vars() && *n );
-
- if( verbose >= 1 )
- cout << "Cost of eliminating " << *n << ": " << _Cl.eliminationCost( *n ) << " new edges" << endl;
-
- result.insert( _Cl.Delta(*n) );
-
- if( verbose >= 1 )
- cout << "_Cl = " << _Cl << endl;
-
- if( verbose >= 1 )
- cout << "After inserting " << _Cl.delta(*n) << ", _Cl = ";
- _Cl.insert( _Cl.delta(*n) );
- if( verbose >= 1 )
- cout << _Cl << endl;
-
- if( verbose >= 1 )
- cout << "After erasing clusters that contain " << *n << ", _Cl = ";
- _Cl.eraseSubsuming( *n );
- if( verbose >= 1 )
- cout << _Cl << endl;
-
- if( verbose >= 1 )
- cout << "After erasing nonmaximal clusters, _Cl = ";
- _Cl.eraseNonMaximal();
- if( verbose >= 1 )
- cout << _Cl << endl;
- }
-
- return result;
-}
-
-
-ClusterGraph ClusterGraph::VarElim_MinFill() const {
- const long verbose = 0;
-
- // Make a copy
- ClusterGraph _Cl(*this);
- VarSet _vars( vars() );
-
- ClusterGraph result;
- _Cl.eraseNonMaximal();
-
- // Do variable elimination
- while( !_vars.empty() ) {
- if( verbose >= 1 )
- cout << "Var Eliminiation cost" << endl;
- VarSet::iterator lowest = _vars.end();
- size_t lowest_cost = -1UL;
- for( VarSet::const_iterator n = _vars.begin(); n != _vars.end(); n++ ) {
- size_t cost = _Cl.eliminationCost( *n );
- if( verbose >= 1 )
- cout << *n << " " << cost << endl;
- if( lowest == _vars.end() || lowest_cost > cost ) {
- lowest = n;
- lowest_cost = cost;
- }
- }
- Var n = *lowest;
-
- if( verbose >= 1 )
- cout << "Lowest: " << n << " (" << lowest_cost << ")" << endl;
-
- result.insert( _Cl.Delta(n) );
-
- _Cl.insert( _Cl.delta(n) );
- _Cl.eraseSubsuming( n );
- _Cl.eraseNonMaximal();
- _vars /= n;
-
- }
-
- return result;
-}
-
-
-}
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __CLUSTERGRAPH_H__
-#define __CLUSTERGRAPH_H__
-
-
-#include <set>
-#include <vector>
-#include "varset.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-/// A ClusterGraph is a hypergraph with VarSets as nodes.
-/// It is implemented as a set<VarSet> in which the adjacency
-/// relationship is computed realtime. It may be better to
-/// implement it as a bipartitegraph, though. The additional
-/// functionality compared to a simple set<VarSet> is
-/// finding maximal clusters, finding cliques, etc...
-class ClusterGraph : public set<VarSet> {
- public:
- /// Default constructor
- ClusterGraph() : set<VarSet>() {}
-
- /// Construct from vector<VarSet>
- ClusterGraph( const vector<VarSet> & cls ) {
- insert( cls.begin(), cls.end() );
- }
-
- /// Copy constructor
- ClusterGraph(const ClusterGraph &x) : set<VarSet>(x) {}
-
- /// Assignment operator
- ClusterGraph& operator=( const ClusterGraph &x ) {
- if( this != &x ) {
- set<VarSet>::operator=( x );
- }
- return *this;
- }
-
- /// Returns true if ns is a maximal member of *this under inclusion (VarSet::operator<<)
- bool isMaximal( const VarSet &ns ) const {
- if( count( ns ) ) {
- // ns is a member
- bool maximal = true;
- for( const_iterator x = begin(); x != end() && maximal; x++ )
- if( (ns << *x) && (ns != *x) )
- maximal = false;
- return maximal;
- } else
- return false;
- }
-
- /// Erase all VarSets that are not maximal
- ClusterGraph& eraseNonMaximal() {
- for( iterator x = begin(); x != end(); )
- if( !isMaximal(*x) )
- erase(x++);
- else
- x++;
- return *this;
- }
-
- /// Return union of all members
- VarSet vars() const {
- VarSet result;
- for( iterator x = begin(); x != end(); x++ )
- result |= *x;
- return result;
- }
-
- /// Returns true if n1 and n2 are adjacent, i.e.\ by
- /// definition, are both contained in some cluster in *this
- bool adj( const Var& n1, const Var& n2 ) {
- bool result = false;
- for( iterator x = begin(); (x != end()) && (!result); x++ )
- if( (*x && n1) && (*x && n2) )
- result = true;
- return result;
- }
-
- /// Returns union of clusters that contain n, minus n
- VarSet Delta( const Var& n ) const {
- VarSet result;
- for( iterator x = begin(); x != end(); x++ )
- if( (*x && n) )
- result |= *x;
- return result;
- }
-
- /// Returns union of clusters that contain n, minus n
- VarSet delta( const Var& n ) const {
- return Delta( n ) / n;
- }
-
- /// Erases all members that contain n
- ClusterGraph& eraseSubsuming( const Var& n ) {
- for( iterator x = begin(); x != end(); )
- if( (*x && n) )
- erase(x++);
- else
- x++;
- return *this;
- }
-
- /// Send to output stream
- friend std::ostream & operator << ( std::ostream & os, const ClusterGraph & cl ) {
- os << "{";
- ClusterGraph::const_iterator x = cl.begin();
- if( x != cl.end() )
- os << *(x++);
- for( ; x != cl.end(); x++ )
- os << ", " << *x;
- os << "}";
- return os;
- }
-
- /// Convert to vector<VarSet>
- vector<VarSet> toVector() const {
- vector<VarSet> result;
- result.reserve( size() );
- for( const_iterator x = begin(); x != end(); x++ )
- result.push_back( *x );
- return result;
- }
-
- /// Calculate cost of eliminating variable n
- /// using as a measure "number of added edges in the adjacency graph"
- /// where the adjacency graph has the variables as its nodes and
- /// connects nodes i1 and i2 iff i1 and i2 occur in some common factor I
- size_t eliminationCost( const Var& n ) {
- VarSet d_n = delta( n );
- size_t cost = 0;
-
- // for each unordered pair {i1,i2} adjacent to n
- for( VarSet::const_iterator i1 = d_n.begin(); i1 != d_n.end(); i1++ ) {
- VarSet d_i1 = delta( *i1 );
- for( VarSet::const_iterator i2 = i1; (++i2) != d_n.end(); ) {
- // if i1 and i2 are not adjacent, eliminating n would make them adjacent
- if( !adj(*i1, *i2) )
- cost++;
- }
- }
- return cost;
- }
-
- /// Perform Variable Elimination without Probs, i.e. only keeping track of
- /// the interactions that are created along the way.
- /// Input: a set of outer clusters and an elimination sequence
- /// Output: a set of elimination "cliques"
- ClusterGraph VarElim( const vector<Var> &ElimSeq ) const;
-
- /// As Taylan does it
- ClusterGraph VarElim_MinFill() const;
-};
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include "daialg.h"
-
-
-namespace dai {
-
-
-/// Calculate the marginal of obj on ns by clamping
-/// all variables in ns and calculating logZ for each joined state
-Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) {
- Factor Pns (ns);
-
- multind mi( ns );
-
- InfAlg *clamped = obj.clone();
- if( !reInit )
- clamped->init();
-
- Complex logZ0;
- for( size_t j = 0; j < mi.max(); j++ ) {
- // save unclamped factors connected to ns
- clamped->saveProbs( ns );
-
- // set clamping Factors to delta functions
- vector<size_t> vi = mi.vi( j );
- size_t k = 0;
- for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++, k++ )
- clamped->clamp( *n, vi[k] );
-
- // run DAIAlg, calc logZ, store in Pns
- if( clamped->Verbose() >= 2 )
- cout << j << ": ";
- if( reInit )
- clamped->init();
- clamped->run();
-
- Complex Z;
- if( j == 0 ) {
- logZ0 = clamped->logZ();
- Z = 1.0;
- } else {
- // subtract logZ0 to avoid very large numbers
- Z = exp(clamped->logZ() - logZ0);
- if( fabs(imag(Z)) > 1e-5 )
- cout << "Marginal:: WARNING: complex Z (" << Z << ")" << endl;
- }
-
- Pns[j] = real(Z);
-
- // restore clamped factors
- clamped->undoProbs( ns );
- }
-
- delete clamped;
-
- return( Pns.normalized(Prob::NORMPROB) );
-}
-
-
-vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit ) {
- // convert ns to vector<VarSet>
- size_t N = ns.size();
- vector<Var> vns;
- vns.reserve( N );
- for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
- vns.push_back( *n );
-
- vector<Factor> pairbeliefs;
- pairbeliefs.reserve( N * N );
- for( size_t j = 0; j < N; j++ )
- for( size_t k = 0; k < N; k++ )
- if( j == k )
- pairbeliefs.push_back(Factor());
- else
- pairbeliefs.push_back(Factor(vns[j] | vns[k]));
-
- InfAlg *clamped = obj.clone();
- if( !reInit )
- clamped->init();
-
- Complex logZ0;
- for( size_t j = 0; j < N; j++ ) {
- // clamp Var j to its possible values
- for( size_t j_val = 0; j_val < vns[j].states(); j_val++ ) {
- if( obj.Verbose() >= 2 )
- cout << j << "/" << N-1 << " (" << j_val << "/" << vns[j].states() << "): ";
-
- // save unclamped factors connected to ns
- clamped->saveProbs( ns );
-
- clamped->clamp( vns[j], j_val );
- if( reInit )
- clamped->init();
- clamped->run();
-
- //if( j == 0 )
- // logZ0 = obj.logZ();
- double Z_xj = 1.0;
- if( j == 0 && j_val == 0 ) {
- logZ0 = clamped->logZ();
- } else {
- // subtract logZ0 to avoid very large numbers
- Complex Z = exp(clamped->logZ() - logZ0);
- if( fabs(imag(Z)) > 1e-5 )
- cout << "calcPairBelief:: Warning: complex Z: " << Z << endl;
- Z_xj = real(Z);
- }
-
- for( size_t k = 0; k < N; k++ )
- if( k != j ) {
- Factor b_k = clamped->belief(vns[k]);
- for( size_t k_val = 0; k_val < vns[k].states(); k_val++ )
- if( vns[j].label() < vns[k].label() )
- pairbeliefs[j * N + k][j_val + (k_val * vns[j].states())] = Z_xj * b_k[k_val];
- else
- pairbeliefs[j * N + k][k_val + (j_val * vns[k].states())] = Z_xj * b_k[k_val];
- }
-
- // restore clamped factors
- clamped->undoProbs( ns );
- }
- }
-
- delete clamped;
-
- // Calculate result by taking the geometric average
- vector<Factor> result;
- result.reserve( N * (N - 1) / 2 );
- for( size_t j = 0; j < N; j++ )
- for( size_t k = j+1; k < N; k++ )
- result.push_back( (pairbeliefs[j * N + k] * pairbeliefs[k * N + j]) ^ 0.5 );
-
- return result;
-}
-
-
-Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit ) {
- // returns a a probability distribution whose 1st order interactions
- // are unspecified, whose 2nd order interactions approximate those of
- // the marginal on ns, and whose higher order interactions are absent.
-
- vector<Factor> pairbeliefs = calcPairBeliefs( obj, ns, reInit );
-
- Factor Pns (ns);
- for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
- Pns *= pairbeliefs[ij];
-
- return( Pns.normalized(Prob::NORMPROB) );
-}
-
-
-vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit ) {
- vector<Factor> result;
- result.reserve( ns.size() * (ns.size() - 1) / 2 );
-
- InfAlg *clamped = obj.clone();
- if( !reInit )
- clamped->init();
-
- Complex logZ0;
- VarSet::const_iterator nj = ns.begin();
- for( long j = 0; j < (long)ns.size() - 1; j++, nj++ ) {
- size_t k = 0;
- for( VarSet::const_iterator nk = nj; (++nk) != ns.end(); k++ ) {
- Factor pairbelief( *nj | *nk );
-
- // clamp Vars j and k to their possible values
- for( size_t j_val = 0; j_val < nj->states(); j_val++ )
- for( size_t k_val = 0; k_val < nk->states(); k_val++ ) {
- // save unclamped factors connected to ns
- clamped->saveProbs( ns );
-
- clamped->clamp( *nj, j_val );
- clamped->clamp( *nk, k_val );
- if( reInit )
- clamped->init();
- clamped->run();
-
- double Z_xj = 1.0;
- if( j_val == 0 && k_val == 0 ) {
- logZ0 = clamped->logZ();
- } else {
- // subtract logZ0 to avoid very large numbers
- Complex Z = exp(clamped->logZ() - logZ0);
- if( fabs(imag(Z)) > 1e-5 )
- cout << "calcPairBelief:: Warning: complex Z: " << Z << endl;
- Z_xj = real(Z);
- }
-
- // we assume that j.label() < k.label()
- // i.e. we make an assumption here about the indexing
- pairbelief[j_val + (k_val * nj->states())] = Z_xj;
-
- // restore clamped factors
- clamped->undoProbs( ns );
- }
-
- result.push_back( pairbelief );
- }
- }
-
- delete clamped;
-
- assert( result.size() == (ns.size() * (ns.size() - 1) / 2) );
-
- return result;
-}
-
-
-}
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __DAIALG_H__
-#define __DAIALG_H__
-
-
-#include <string>
-#include <iostream>
-#include <vector>
-#include "factorgraph.h"
-#include "regiongraph.h"
-#include "properties.h"
-
-
-namespace dai {
-
-
-/// The InfAlg class is the common denominator of the various approximate inference algorithms.
-/// A InfAlg object represents a discrete factorized probability distribution over multiple variables
-/// together with an inference algorithm.
-class InfAlg {
- private:
- /// Properties of the algorithm (replaces _tol, _maxiter, _verbose)
- Properties _properties;
-
- /// Maximum difference encountered so far
- double _maxdiff;
-
-
- public:
- /// Default constructor
- InfAlg() : _properties(), _maxdiff(0.0) {}
-
- /// Constructor with options
- InfAlg( const Properties &opts ) : _properties(opts), _maxdiff(0.0) {}
-
- /// Copy constructor
- InfAlg( const InfAlg & x ) : _properties(x._properties), _maxdiff(x._maxdiff) {}
-
- /// Clone (virtual copy constructor)
- virtual InfAlg* clone() const = 0;
-
- /// Assignment operator
- InfAlg & operator=( const InfAlg & x ) {
- if( this != &x ) {
- _properties = x._properties;
- _maxdiff = x._maxdiff;
- }
- return *this;
- }
-
- /// Virtual desctructor
- // (this is needed because this class contains virtual functions)
- virtual ~InfAlg() {}
-
- /// Returns true if a property with the given key is present
- bool HasProperty(const PropertyKey &key) const { return _properties.hasKey(key); }
-
- /// Gets a property
- const PropertyValue & GetProperty(const PropertyKey &key) const { return _properties.Get(key); }
-
- /// Gets a property, casted as ValueType
- template<typename ValueType>
- ValueType GetPropertyAs(const PropertyKey &key) const { return _properties.GetAs<ValueType>(key); }
-
- /// Sets a property
- void SetProperty(const PropertyKey &key, const PropertyValue &val) { _properties[key] = val; }
-
- /// Converts a property from string to ValueType, if necessary
- template<typename ValueType>
- void ConvertPropertyTo(const PropertyKey &key) { _properties.ConvertTo<ValueType>(key); }
-
- /// Gets all properties
- const Properties & GetProperties() const { return _properties; }
-
- /// Sets properties
- void SetProperties(const Properties &p) { _properties = p; }
-
- /// Sets tolerance
- void Tol( double tol ) { SetProperty("tol", tol); }
- /// Gets tolerance
- double Tol() const { return GetPropertyAs<double>("tol"); }
-
- /// Sets maximum number of iterations
- void MaxIter( size_t maxiter ) { SetProperty("maxiter", maxiter); }
- /// Gets maximum number of iterations
- size_t MaxIter() const { return GetPropertyAs<size_t>("maxiter"); }
-
- /// Sets verbosity
- void Verbose( size_t verbose ) { SetProperty("verbose", verbose); }
- /// Gets verbosity
- size_t Verbose() const { return GetPropertyAs<size_t>("verbose"); }
-
- /// Sets maximum difference encountered so far
- void MaxDiff( double maxdiff ) { _maxdiff = maxdiff; }
- /// Gets maximum difference encountered so far
- double MaxDiff() const { return _maxdiff; }
- /// Updates maximum difference encountered so far
- void updateMaxDiff( double maxdiff ) { if( maxdiff > _maxdiff ) _maxdiff = maxdiff; }
- /// Sets maximum difference encountered so far to zero
- void clearMaxDiff() { _maxdiff = 0.0; }
-
- /// Identifies itself for logging purposes
- virtual std::string identify() const = 0;
-
- /// Get single node belief
- virtual Factor belief( const Var &n ) const = 0;
-
- /// Get general belief
- virtual Factor belief( const VarSet &n ) const = 0;
-
- /// Get all beliefs
- virtual vector<Factor> beliefs() const = 0;
-
- /// Get log partition sum
- virtual Complex logZ() const = 0;
-
- /// Clear messages and beliefs
- virtual void init() = 0;
-
- /// The actual approximate inference algorithm
- virtual double run() = 0;
-
- /// Save factor I
- virtual void saveProb( size_t I ) = 0;
- /// Save Factors involving ns
- virtual void saveProbs( const VarSet &ns ) = 0;
-
- /// Restore factor I
- virtual void undoProb( size_t I ) = 0;
- /// Restore Factors involving ns
- virtual void undoProbs( const VarSet &ns ) = 0;
-
- /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
- virtual void clamp( const Var & n, size_t i ) = 0;
-
- /// Return all variables that interact with n
- virtual VarSet delta( const Var & n ) const = 0;
-
- /// Set all factors interacting with n to 1
- virtual void makeCavity( const Var & n ) = 0;
-
- /// Set factor I to 1
- virtual void makeFactorCavity( size_t I ) = 0;
-
- /// Get index of variable n
- virtual size_t findVar( const Var & n ) const = 0;
-
- /// Get index of first factor involving ns
- virtual size_t findFactor( const VarSet &ns ) const = 0;
-
- /// Get number of variables
- virtual size_t nrVars() const = 0;
-
- /// Get number of factors
- virtual size_t nrFactors() const = 0;
-
- /// Get const reference to variable i
- virtual const Var & var(size_t i) const = 0;
-
- /// Get reference to variable i
- virtual Var & var(size_t i) = 0;
-
- /// Get const reference to factor I
- virtual const Factor & factor( size_t I ) const = 0;
-
- /// Get reference to factor I
- virtual Factor & factor( size_t I ) = 0;
-
- /// Factor I has been updated
- virtual void updatedFactor( size_t I ) = 0;
-
- /// Checks whether all necessary properties have been set
- /// and casts string-valued properties to other values if necessary
- virtual bool checkProperties() = 0;
-};
-
-
-template <class T>
-class DAIAlg : public InfAlg, public T {
- public:
- /// Default constructor
- DAIAlg() : InfAlg(), T() {}
-
- /// Construct DAIAlg with empty T but using the specified properties
- DAIAlg( const Properties &opts ) : InfAlg( opts ), T() {}
-
- /// Construct DAIAlg using the specified properties
- DAIAlg( const T & t, const Properties &opts ) : InfAlg( opts ), T(t) {}
-
- /// Copy constructor
- DAIAlg( const DAIAlg & x ) : InfAlg(x), T(x) {}
-
- /// Save factor I (using T::saveProb)
- void saveProb( size_t I ) { T::saveProb( I ); }
- /// Save Factors involving ns (using T::saveProbs)
- void saveProbs( const VarSet &ns ) { T::saveProbs( ns ); }
-
- /// Restore factor I (using T::undoProb)
- void undoProb( size_t I ) { T::undoProb( I ); }
- /// Restore Factors involving ns (using T::undoProbs)
- void undoProbs( const VarSet &ns ) { T::undoProbs( ns ); }
-
- /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$) (using T::clamp)
- void clamp( const Var & n, size_t i ) { T::clamp( n, i ); }
-
- /// Return all variables that interact with n (using T::delta)
- VarSet delta( const Var & n ) const { return T::delta(n); }
-
- /// Set all factors interacting with n to 1 (using T::makeCavity)
- void makeCavity( const Var & n ) { T::makeCavity(n); }
-
- /// Set factor I to 1 (using T::makeFactorCavity)
- void makeFactorCavity( size_t I ) { T::makeFactorCavity(I); }
-
- /// Get index of variable n (using T::findVar)
- size_t findVar( const Var & n ) const { return T::findVar(n); }
-
- /// Get index of first factor involving ns (using T::findFactor)
- size_t findFactor( const VarSet &ns ) const { return T::findFactor(ns); }
-
- /// Get number of variables (using T::nrFactors)
- size_t nrVars() const { return T::nrVars(); }
-
- /// Get number of factors (using T::nrFactors)
- size_t nrFactors() const { return T::nrFactors(); }
-
- /// Get const reference to variable i (using T::var)
- const Var & var(size_t i) const { return T::var(i); }
-
- /// Get reference to variable i (using T::var)
- Var & var(size_t i) { return T::var(i); }
-
- /// Get const reference to factor I (using T::factor)
- const Factor & factor( size_t I ) const { return T::factor(I); }
-
- /// Get reference to factor I (using T::factor)
- Factor & factor( size_t I ) { return T::factor(I); }
-
- /// Factor I has been updated (using T::updatedFactor)
- void updatedFactor( size_t I ) { T::updatedFactor(I); }
-};
-
-
-typedef DAIAlg<FactorGraph> DAIAlgFG;
-typedef DAIAlg<RegionGraph> DAIAlgRG;
-
-
-/// Calculate the marginal of obj on ns by clamping
-/// all variables in ns and calculating logZ for each joined state
-Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit );
-
-
-/// Calculate beliefs of all pairs in ns (by clamping
-/// nodes in ns and calculating logZ and the beliefs for each state)
-vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit );
-
-
-/// Calculate beliefs of all pairs in ns (by clamping
-/// pairs in ns and calculating logZ for each joined state)
-vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit );
-
-
-/// Calculate 2nd order interactions of the marginal of obj on ns
-Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit );
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __DIFFS_H__
-#define __DIFFS_H__
-
-
-#include <vector>
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-class Diffs : public vector<double> {
- private:
- size_t _maxsize;
- double _def;
- vector<double>::iterator _pos;
- vector<double>::iterator _maxpos;
- public:
- Diffs(long maxsize, double def) : vector<double>(), _maxsize(maxsize), _def(def) {
- this->reserve(_maxsize);
- _pos = begin();
- _maxpos = begin();
- };
- double max() {
- if( size() < _maxsize )
- return _def;
- else
- return( *_maxpos );
- }
- void push(double x) {
- if( size() < _maxsize ) {
- push_back(x);
- _pos = end();
- _maxpos = max_element(begin(),end());
- }
- else {
- if( _pos == end() )
- _pos = begin();
- if( _maxpos == _pos ) {
- *_pos++ = x;
- _maxpos = max_element(begin(),end());
- } else {
- if( x > *_maxpos )
- _maxpos = _pos;
- *_pos++ = x;
- }
- }
- }
-};
-
-
-}
-
-
-#endif
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
-INPUT =
+INPUT = src include/dai
# If the value of the INPUT tag contains directories, you can use the
# FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __ENUM_H
-#define __ENUM_H
-
-
-#include <cstring>
-#include <iostream>
-
-
-namespace dai {
-
-
-// C++ enums are too limited for my purposes. This defines wrapper classes
-// that provide much more functionality than a simple enum. The only
-// disadvantage is that one wrapper class needs to be written for each
-// number of values an enum can take... a better solution is needed.
-
-
-#define ENUM2(x,a,b) class x {\
- public:\
- enum value {a, b};\
-\
- x() : v(a) {}\
-\
- x(value w) : v(w) {}\
-\
- x(char const *w) {\
- static char const* labels[] = {#a, #b};\
- size_t i = 0;\
- for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
- if( strcmp( w, labels[i] ) == 0 ) {\
- v = (value)i;\
- break;\
- }\
- if( i == sizeof(labels) / sizeof(char const *) )\
- throw "Unknown " #x " value";\
- }\
-\
- operator value () const { return v; }\
-\
- operator size_t () const { return (size_t)v; }\
-\
- operator char const* () const {\
- static char const* labels[] = {#a, #b};\
- return labels[v];\
- }\
-\
- friend std::istream& operator >> (std::istream& is, x& y) {\
- std::string s;\
- is >> s;\
- y = x(s.c_str());\
- return is;\
- }\
-\
- friend std::ostream& operator << (std::ostream& os, const x& y) {\
- os << (const char *)y;\
- return os;\
- }\
-\
- private:\
- value v;\
-};
-
-
-#define ENUM3(x,a,b,c) class x {\
- public:\
- enum value {a, b, c};\
-\
- x() : v(a) {}\
-\
- x(value w) : v(w) {}\
-\
- x(char const *w) {\
- static char const* labels[] = {#a, #b, #c};\
- size_t i = 0;\
- for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
- if( strcmp( w, labels[i] ) == 0 ) {\
- v = (value)i;\
- break;\
- }\
- if( i == sizeof(labels) / sizeof(char const *) )\
- throw "Unknown " #x " value";\
- }\
-\
- operator value () const { return v; }\
-\
- operator size_t () const { return (size_t)v; }\
-\
- operator char const* () const {\
- static char const* labels[] = {#a, #b, #c};\
- return labels[v];\
- }\
-\
- friend std::istream& operator >> (std::istream& is, x& y) {\
- std::string s;\
- is >> s;\
- y = x(s.c_str());\
- return is;\
- }\
-\
- friend std::ostream& operator << (std::ostream& os, const x& y) {\
- os << (const char *)y;\
- return os;\
- }\
-\
- private:\
- value v;\
-};
-
-
-#define ENUM4(x,a,b,c,d) class x {\
- public:\
- enum value {a, b, c, d};\
-\
- x() : v(a) {}\
-\
- x(value w) : v(w) {}\
-\
- x(char const *w) {\
- static char const* labels[] = {#a, #b, #c, #d};\
- size_t i = 0;\
- for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
- if( strcmp( w, labels[i] ) == 0 ) {\
- v = (value)i;\
- break;\
- }\
- if( i == sizeof(labels) / sizeof(char const *) )\
- throw "Unknown " #x " value";\
- }\
-\
- operator value () const { return v; }\
-\
- operator size_t () const { return (size_t)v; }\
-\
- operator char const* () const {\
- static char const* labels[] = {#a, #b, #c, #d};\
- return labels[v];\
- }\
-\
- friend std::istream& operator >> (std::istream& is, x& y) {\
- std::string s;\
- is >> s;\
- y = x(s.c_str());\
- return is;\
- }\
-\
- friend std::ostream& operator << (std::ostream& os, const x& y) {\
- os << (const char *)y;\
- return os;\
- }\
-\
- private:\
- value v;\
-};
-
-
-#define ENUM5(x,a,b,c,d,e) class x {\
- public:\
- enum value {a, b, c, d, e};\
-\
- x() : v(a) {}\
-\
- x(value w) : v(w) {}\
-\
- x(char const *w) {\
- static char const* labels[] = {#a, #b, #c, #d, #e};\
- size_t i = 0;\
- for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
- if( strcmp( w, labels[i] ) == 0 ) {\
- v = (value)i;\
- break;\
- }\
- if( i == sizeof(labels) / sizeof(char const *) )\
- throw "Unknown " #x " value";\
- }\
-\
- operator value () const { return v; }\
-\
- operator size_t () const { return (size_t)v; }\
-\
- operator char const* () const {\
- static char const* labels[] = {#a, #b, #c, #d, #e};\
- return labels[v];\
- }\
-\
- friend std::istream& operator >> (std::istream& is, x& y) {\
- std::string s;\
- is >> s;\
- y = x(s.c_str());\
- return is;\
- }\
-\
- friend std::ostream& operator << (std::ostream& os, const x& y) {\
- os << (const char *)y;\
- return os;\
- }\
-\
- private:\
- value v;\
-};
-
-
-#define ENUM6(x,a,b,c,d,e,f) class x {\
- public:\
- enum value {a, b, c, d, e, f};\
-\
- x() : v(a) {}\
-\
- x(value w) : v(w) {}\
-\
- x(char const *w) {\
- static char const* labels[] = {#a, #b, #c, #d, #e, #f};\
- size_t i = 0;\
- for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
- if( strcmp( w, labels[i] ) == 0 ) {\
- v = (value)i;\
- break;\
- }\
- if( i == sizeof(labels) / sizeof(char const *) )\
- throw "Unknown " #x " value";\
- }\
-\
- operator value () const { return v; }\
-\
- operator size_t () const { return (size_t)v; }\
-\
- operator char const* () const {\
- static char const* labels[] = {#a, #b, #c, #d, #e, #f};\
- return labels[v];\
- }\
-\
- friend std::istream& operator >> (std::istream& is, x& y) {\
- std::string s;\
- is >> s;\
- y = x(s.c_str());\
- return is;\
- }\
-\
- friend std::ostream& operator << (std::ostream& os, const x& y) {\
- os << (const char *)y;\
- return os;\
- }\
-\
- private:\
- value v;\
-};
-
-
-}
-
-
-#endif
#include <iostream>
-#include "alldai.h"
+#include <dai/alldai.h>
using namespace dai;
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __FACTOR_H__
-#define __FACTOR_H__
-
-
-#include <iostream>
-#include <cmath>
-#include "prob.h"
-#include "varset.h"
-#include "index.h"
-
-
-namespace dai {
-
-
-template<typename T> class TFactor;
-typedef TFactor<Real> Factor;
-typedef TFactor<Complex> CFactor;
-
-
-// predefine friends
-template<typename T> Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt );
-template<typename T> Complex KL_dist( const TFactor<T> & p, const TFactor<T> & q );
-template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P);
-
-
-// T should be castable from and to double and to complex
-template <typename T> class TFactor {
- protected:
- VarSet _vs;
- TProb<T> _p;
-
- public:
- // Default constructor
- TFactor () : _vs(), _p(1,1.0) {}
-
- // Construct Factor from VarSet
- TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.stateSpace()) {}
-
- // Construct Factor from VarSet and initial value
- TFactor( const VarSet& ns, Real p ) : _vs(ns), _p(_vs.stateSpace(),p) {}
-
- // Construct Factor from VarSet and initial array
- TFactor( const VarSet& ns, const Real* p ) : _vs(ns), _p(_vs.stateSpace(),p) {}
-
- // Construct Factor from VarSet and TProb<T>
- TFactor( const VarSet& ns, const TProb<T> p ) : _vs(ns), _p(p) {
-#ifdef DEBUG
- assert( _vs.stateSpace() == _p.size() );
-#endif
- }
-
- // Construct Factor from Var
- TFactor( const Var& n ) : _vs(n), _p(n.states()) {}
-
- // Copy constructor
- TFactor( const TFactor<T> &x ) : _vs(x._vs), _p(x._p) {}
-
- // Assignment operator
- TFactor<T> & operator= (const TFactor<T> &x) {
- if( this != &x ) {
- _vs = x._vs;
- _p = x._p;
- }
- return *this;
- }
-
- const TProb<T> & p() const { return _p; }
- TProb<T> & p() { return _p; }
- const VarSet & vars() const { return _vs; }
- size_t stateSpace() const {
-#ifdef DEBUG
- assert( _vs.stateSpace() == _p.size() );
-#endif
- return _p.size();
- }
-
- T operator[] (size_t i) const { return _p[i]; }
- T& operator[] (size_t i) { return _p[i]; }
- TFactor<T> & fill (T p)
- { _p.fill( p ); return(*this); }
- TFactor<T> & randomize ()
- { _p.randomize(); return(*this); }
- TFactor<T> operator* (T x) const {
- Factor result = *this;
- result.p() *= x;
- return result;
- }
- TFactor<T>& operator*= (T x) {
- _p *= x;
- return *this;
- }
- TFactor<T> operator/ (T x) const {
- Factor result = *this;
- result.p() /= x;
- return result;
- }
- TFactor<T>& operator/= (T x) {
- _p /= x;
- return *this;
- }
- TFactor<T> operator* (const TFactor<T>& Q) const;
- TFactor<T>& operator*= (const TFactor<T>& Q) { return( *this = (*this * Q) ); }
- TFactor<T> operator+ (const TFactor<T>& Q) const {
-#ifdef DEBUG
- assert( Q._vs == _vs );
-#endif
- TFactor<T> sum(*this);
- sum._p += Q._p;
- return sum;
- }
- TFactor<T> operator- (const TFactor<T>& Q) const {
-#ifdef DEBUG
- assert( Q._vs == _vs );
-#endif
- TFactor<T> sum(*this);
- sum._p -= Q._p;
- return sum;
- }
- TFactor<T>& operator+= (const TFactor<T>& Q) {
-#ifdef DEBUG
- assert( Q._vs == _vs );
-#endif
- _p += Q._p;
- return *this;
- }
- TFactor<T>& operator-= (const TFactor<T>& Q) {
-#ifdef DEBUG
- assert( Q._vs == _vs );
-#endif
- _p -= Q._p;
- return *this;
- }
-
- TFactor<T> operator^ (Real a) const { TFactor<T> x; x._vs = _vs; x._p = _p^a; return x; }
- TFactor<T>& operator^= (Real a) { _p ^= a; return *this; }
-
- TFactor<T>& makeZero( Real epsilon ) {
- _p.makeZero( epsilon );
- return *this;
- }
-
- TFactor<T> inverse() const {
- TFactor<T> inv;
- inv._vs = _vs;
- inv._p = _p.inverse(true); // FIXME
- return inv;
- }
-
- TFactor<T> divided_by( const TFactor<T>& denom ) const {
-#ifdef DEBUG
- assert( denom._vs == _vs );
-#endif
- TFactor<T> quot(*this);
- quot._p /= denom._p;
- return quot;
- }
-
- TFactor<T>& divide( const TFactor<T>& denom ) {
-#ifdef DEBUG
- assert( denom._vs == _vs );
-#endif
- _p /= denom._p;
- return *this;
- }
-
- TFactor<T> exp() const {
- TFactor<T> e;
- e._vs = _vs;
- e._p = _p.exp();
- return e;
- }
-
- TFactor<T> log() const {
- TFactor<T> l;
- l._vs = _vs;
- l._p = _p.log();
- return l;
- }
-
- TFactor<T> log0() const {
- TFactor<T> l0;
- l0._vs = _vs;
- l0._p = _p.log0();
- return l0;
- }
-
- CFactor clog0() const {
- CFactor l0;
- l0._vs = _vs;
- l0._p = _p.clog0();
- return l0;
- }
-
- T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); }
- TFactor<T> normalized( typename Prob::NormType norm ) const {
- TFactor<T> result;
- result._vs = _vs;
- result._p = _p.normalized( norm );
- return result;
- }
-
- // returns slice of this factor where the subset ns is in state ns_state
- Factor slice( const VarSet & ns, size_t ns_state ) const {
- assert( ns << _vs );
- VarSet nsrem = _vs / ns;
- Factor result( nsrem, 0.0 );
-
- // OPTIMIZE ME
- Index i_ns (ns, _vs);
- Index i_nsrem (nsrem, _vs);
- for( size_t i = 0; i < stateSpace(); i++, ++i_ns, ++i_nsrem )
- if( (size_t)i_ns == ns_state )
- result._p[i_nsrem] = _p[i];
-
- return result;
- }
-
- // returns unnormalized marginal
- TFactor<T> part_sum(const VarSet & ns) const;
- // returns normalized marginal
- TFactor<T> marginal(const VarSet & ns) const { return part_sum(ns).normalized( Prob::NORMPROB ); }
-
- bool hasNaNs() const { return _p.hasNaNs(); }
- bool hasNegatives() const { return _p.hasNegatives(); }
- T totalSum() const { return _p.totalSum(); }
- T maxAbs() const { return _p.maxAbs(); }
- T max() const { return _p.max(); }
- Complex entropy() const { return _p.entropy(); }
- T strength( const Var &i, const Var &j ) const;
-
- friend Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt ) {
- if( x._vs.empty() || y._vs.empty() )
- return -1;
- else {
-#ifdef DEBUG
- assert( x._vs == y._vs );
-#endif
- return dist( x._p, y._p, dt );
- }
- }
- friend Complex KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
- template<class U> friend std::ostream& operator<< (std::ostream& os, const TFactor<U>& P);
-};
-
-
-template<typename T> TFactor<T> TFactor<T>::part_sum(const VarSet & ns) const {
-#ifdef DEBUG
- assert( ns << _vs );
-#endif
-
- TFactor<T> res( ns, 0.0 );
-
- Index i_res( ns, _vs );
- for( size_t i = 0; i < _p.size(); i++, ++i_res )
- res._p[i_res] += _p[i];
-
- return res;
-}
-
-
-template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P) {
- os << "(" << P.vars() << " <";
- for( size_t i = 0; i < P._p.size(); i++ )
- os << P._p[i] << " ";
- os << ">)";
- return os;
-}
-
-
-template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& Q) const {
- TFactor<T> prod( _vs | Q._vs, 0.0 );
-
- Index i1(_vs, prod._vs);
- Index i2(Q._vs, prod._vs);
-
- for( size_t i = 0; i < prod._p.size(); i++, ++i1, ++i2 )
- prod._p[i] += _p[i1] * Q._p[i2];
-
- return prod;
-}
-
-
-template<typename T> Complex KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
- if( P._vs.empty() || Q._vs.empty() )
- return -1;
- else {
-#ifdef DEBUG
- assert( P._vs == Q._vs );
-#endif
- return KL_dist( P._p, Q._p );
- }
-}
-
-
-// calculate N(psi, i, j)
-template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
-#ifdef DEBUG
- assert( _vs && i );
- assert( _vs && j );
- assert( i != j );
-#endif
- VarSet ij = i | j;
-
- T max = 0.0;
- for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
- for( size_t alpha2 = 0; alpha2 < i.states(); alpha2++ )
- if( alpha2 != alpha1 )
- for( size_t beta1 = 0; beta1 < j.states(); beta1++ )
- for( size_t beta2 = 0; beta2 < j.states(); beta2++ )
- if( beta2 != beta1 ) {
- size_t as = 1, bs = 1;
- if( i < j )
- bs = i.states();
- else
- as = j.states();
- T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).max();
- T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).max();
- T f = f1 * f2;
- if( f > max )
- max = f;
- }
-
- return std::tanh( 0.25 * std::log( max ) );
-}
-
-
-template<typename T> TFactor<T> RemoveFirstOrderInteractions( const TFactor<T> & psi ) {
- TFactor<T> result = psi;
-
- VarSet vars = psi.vars();
- for( size_t iter = 0; iter < 100; iter++ ) {
- for( VarSet::const_iterator n = vars.begin(); n != vars.end(); n++ )
- result = result * result.part_sum(*n).inverse();
- result.normalize( Prob::NORMPROB );
- }
-
- return result;
-}
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include <iostream>
-#include <iterator>
-#include <map>
-#include <set>
-#include <fstream>
-#include <string>
-#include <algorithm>
-#include <functional>
-#include <tr1/unordered_map>
-#include "factorgraph.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-FactorGraph::FactorGraph( const vector<Factor> &P ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {
- // add factors, obtain variables
- set<Var> _vars;
- V2s().reserve( P.size() );
- size_t nrEdges = 0;
- for( vector<Factor>::const_iterator p2 = P.begin(); p2 != P.end(); p2++ ) {
- V2s().push_back( *p2 );
- copy( p2->vars().begin(), p2->vars().end(), inserter( _vars, _vars.begin() ) );
- nrEdges += p2->vars().size();
- }
-
- // add _vars
- V1s().reserve( _vars.size() );
- for( set<Var>::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
- V1s().push_back( *p1 );
-
- // create graph structure
- createGraph( nrEdges );
-}
-
-
-/// Part of constructors (creates edges, neighbours and adjacency matrix)
-void FactorGraph::createGraph( size_t nrEdges ) {
- // create a mapping for indices
- std::tr1::unordered_map<size_t, size_t> hashmap;
-
- for( size_t i = 0; i < vars().size(); i++ )
- hashmap[vars()[i].label()] = i;
-
- // create edges
- edges().reserve( nrEdges );
- for( size_t i2 = 0; i2 < nrFactors(); i2++ ) {
- const VarSet& ns = factor(i2).vars();
- for( VarSet::const_iterator q = ns.begin(); q != ns.end(); q++ )
- edges().push_back(_edge_t(hashmap[q->label()], i2));
- }
-
- // calc neighbours and adjacency matrix
- Regenerate();
-}
-
-
-/*FactorGraph& FactorGraph::addFactor( const Factor &I ) {
- // add Factor
- _V2.push_back( I );
-
- // add new vars in Factor
- for( VarSet::const_iterator i = I.vars().begin(); i != I.vars().end(); i++ ) {
- size_t i_ind = find(vars().begin(), vars().end(), *i) - vars().begin();
- if( i_ind == vars().size() )
- _V1.push_back( *i );
- _E12.push_back( _edge_t( i_ind, nrFactors() - 1 ) );
- }
-
- Regenerate();
- return(*this);
-}*/
-
-
-ostream& operator << (ostream& os, const FactorGraph& fg) {
- os << fg.nrFactors() << endl;
-
- for( size_t I = 0; I < fg.nrFactors(); I++ ) {
- os << endl;
- os << fg.factor(I).vars().size() << endl;
- for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
- os << i->label() << " ";
- os << endl;
- for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
- os << i->states() << " ";
- os << endl;
- size_t nr_nonzeros = 0;
- for( size_t k = 0; k < fg.factor(I).stateSpace(); k++ )
- if( fg.factor(I)[k] != 0.0 )
- nr_nonzeros++;
- os << nr_nonzeros << endl;
- for( size_t k = 0; k < fg.factor(I).stateSpace(); k++ )
- if( fg.factor(I)[k] != 0.0 ) {
- char buf[20];
- sprintf(buf,"%18.14g", fg.factor(I)[k]);
- os << k << " " << buf << endl;
- }
- }
-
- return(os);
-}
-
-
-istream& operator >> (istream& is, FactorGraph& fg) {
- long verbose = 0;
-
- try {
- vector<Factor> factors;
- size_t nr_f;
- string line;
-
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> nr_f;
- if( is.fail() )
- throw "ReadFromFile: unable to read number of Factors";
- if( verbose >= 2 )
- cout << "Reading " << nr_f << " factors..." << endl;
-
- getline (is,line);
- if( is.fail() )
- throw "ReadFromFile: empty line expected";
-
- for( size_t I = 0; I < nr_f; I++ ) {
- if( verbose >= 3 )
- cout << "Reading factor " << I << "..." << endl;
- size_t nr_members;
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> nr_members;
- if( verbose >= 3 )
- cout << " nr_members: " << nr_members << endl;
-
- vector<long> labels;
- for( size_t mi = 0; mi < nr_members; mi++ ) {
- long mi_label;
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> mi_label;
- labels.push_back(mi_label);
- }
- if( verbose >= 3 ) {
- cout << " labels: ";
- copy (labels.begin(), labels.end(), ostream_iterator<int>(cout, " "));
- cout << endl;
- }
-
- vector<size_t> dims;
- for( size_t mi = 0; mi < nr_members; mi++ ) {
- size_t mi_dim;
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> mi_dim;
- dims.push_back(mi_dim);
- }
- if( verbose >= 3 ) {
- cout << " dimensions: ";
- copy (dims.begin(), dims.end(), ostream_iterator<int>(cout, " "));
- cout << endl;
- }
-
- // add the Factor
- VarSet I_vars;
- for( size_t mi = 0; mi < nr_members; mi++ )
- I_vars.insert( Var(labels[mi], dims[mi]) );
- factors.push_back(Factor(I_vars,0.0));
-
- // calculate permutation sigma (internally, members are sorted)
- vector<long> sigma(nr_members,0);
- VarSet::iterator j = I_vars.begin();
- for( size_t mi = 0; mi < nr_members; mi++,j++ ) {
- long search_for = j->label();
- vector<long>::iterator j_loc = find(labels.begin(),labels.end(),search_for);
- sigma[mi] = j_loc - labels.begin();
- }
- if( verbose >= 3 ) {
- cout << " sigma: ";
- copy( sigma.begin(), sigma.end(), ostream_iterator<int>(cout," "));
- cout << endl;
- }
-
- // calculate multindices
- vector<size_t> sdims(nr_members,0);
- for( size_t k = 0; k < nr_members; k++ ) {
- sdims[k] = dims[sigma[k]];
- }
- multind mi(dims);
- multind smi(sdims);
- if( verbose >= 3 ) {
- cout << " mi.max(): " << mi.max() << endl;
- cout << " ";
- for( size_t k=0; k < nr_members; k++ )
- cout << labels[k] << " ";
- cout << " ";
- for( size_t k=0; k < nr_members; k++ )
- cout << labels[sigma[k]] << " ";
- cout << endl;
- }
-
- // read values
- size_t nr_nonzeros;
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> nr_nonzeros;
- if( verbose >= 3 )
- cout << " nonzeroes: " << nr_nonzeros << endl;
- for( size_t k = 0; k < nr_nonzeros; k++ ) {
- size_t li;
- double val;
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> li;
- while( (is.peek()) == '#' )
- getline(is,line);
- is >> val;
-
- vector<size_t> vi = mi.vi(li);
- vector<size_t> svi(vi.size(),0);
- for( size_t k = 0; k < vi.size(); k++ )
- svi[k] = vi[sigma[k]];
- size_t sli = smi.li(svi);
- if( verbose >= 3 ) {
- cout << " " << li << ": ";
- copy( vi.begin(), vi.end(), ostream_iterator<size_t>(cout," "));
- cout << "-> ";
- copy( svi.begin(), svi.end(), ostream_iterator<size_t>(cout," "));
- cout << ": " << sli << endl;
- }
- factors.back()[sli] = val;
- }
- }
-
- if( verbose >= 3 ) {
- cout << "factors:" << endl;
- copy(factors.begin(), factors.end(), ostream_iterator<Factor>(cout,"\n"));
- }
-
- fg = FactorGraph(factors);
- } catch (char *e) {
- cout << e << endl;
- }
-
- return is;
-}
-
-
-VarSet FactorGraph::delta(const Var & n) const {
- // calculate Markov Blanket
- size_t i = findVar( n );
-
- VarSet del;
- for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ )
- for( _nb_cit j = nb2(*I).begin(); j != nb2(*I).end(); j++ )
- if( *j != i )
- del |= var(*j);
-
- return del;
-}
-
-
-VarSet FactorGraph::Delta(const Var & n) const {
- return( delta(n) | n );
-}
-
-
-void FactorGraph::makeFactorCavity(size_t I) {
- // fill Factor I with ones
- factor(I).fill(1.0);
-}
-
-
-void FactorGraph::makeCavity(const Var & n) {
- // fills all Factors that include Var n with ones
- size_t i = findVar( n );
-
- for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ )
- factor(*I).fill(1.0);
-}
-
-
-bool FactorGraph::hasNegatives() const {
- bool result = false;
- for( size_t I = 0; I < nrFactors() && !result; I++ )
- if( factor(I).hasNegatives() )
- result = true;
- return result;
-}
-
-
-/*FactorGraph & FactorGraph::DeleteFactor(size_t I) {
- // Go through all edges
- for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
- if( edge->second >= I ) {
- if( edge->second == I )
- edge->second = -1UL;
- else
- (edge->second)--;
- }
- // Remove all edges containing I
- for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
- if( edge->second == -1UL )
- edge = _E12.erase( edge );
-// vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select2nd<_edge_t>() ) );
-// _E12.erase( new_end, _E12.end() );
-
- // Erase the factor
- _V2.erase( _V2.begin() + I );
-
- Regenerate();
-
- return *this;
-}
-
-
-FactorGraph & FactorGraph::DeleteVar(size_t i) {
- // Go through all edges
- for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
- if( edge->first >= i ) {
- if( edge->first == i )
- edge->first = -1UL;
- else
- (edge->first)--;
- }
- // Remove all edges containing i
- for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
- if( edge->first == -1UL )
- edge = _E12.erase( edge );
-
-// vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select1st<_edge_t>() ) );
-// _E12.erase( new_end, _E12.end() );
-
- // Erase the variable
- _V1.erase( _V1.begin() + i );
-
- Regenerate();
-
- return *this;
-}*/
-
-
-long FactorGraph::ReadFromFile(const char *filename) {
- ifstream infile;
- infile.open (filename);
- if (infile.is_open()) {
- infile >> *this;
- infile.close();
- return 0;
- } else {
- cout << "ERROR OPENING FILE" << endl;
- return 1;
- }
-}
-
-
-long FactorGraph::WriteToFile(const char *filename) const {
- ofstream outfile;
- outfile.open (filename);
- if (outfile.is_open()) {
- try {
- outfile << *this;
- } catch (char *e) {
- cout << e << endl;
- return 1;
- }
- outfile.close();
- return 0;
- } else {
- cout << "ERROR OPENING FILE" << endl;
- return 1;
- }
-}
-
-
-long FactorGraph::WriteToDotFile(const char *filename) const {
- ofstream outfile;
- outfile.open (filename);
- if (outfile.is_open()) {
- try {
- outfile << "graph G {" << endl;
- outfile << "graph[size=\"9,9\"];" << endl;
- outfile << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
- for( size_t i = 0; i < nrVars(); i++ )
- outfile << "\tx" << var(i).label() << ";" << endl;
- outfile << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
- for( size_t I = 0; I < nrFactors(); I++ )
- outfile << "\tp" << I << ";" << endl;
- for( size_t iI = 0; iI < nr_edges(); iI++ )
- outfile << "\tx" << var(edge(iI).first).label() << " -- p" << edge(iI).second << ";" << endl;
- outfile << "}" << endl;
- } catch (char *e) {
- cout << e << endl;
- return 1;
- }
- outfile.close();
- return 0;
- } else {
- cout << "ERROR OPENING FILE" << endl;
- return 1;
- }
-}
-
-
-bool hasShortLoops( const vector<Factor> &P ) {
- bool found = false;
- vector<Factor>::const_iterator I, J;
- for( I = P.begin(); I != P.end(); I++ ) {
- J = I;
- J++;
- for( ; J != P.end(); J++ )
- if( (I->vars() & J->vars()).size() >= 2 ) {
- found = true;
- break;
- }
- if( found )
- break;
- }
- return found;
-}
-
-
-void RemoveShortLoops(vector<Factor> &P) {
- bool found = true;
- while( found ) {
- found = false;
- vector<Factor>::iterator I, J;
- for( I = P.begin(); I != P.end(); I++ ) {
- J = I;
- J++;
- for( ; J != P.end(); J++ )
- if( (I->vars() & J->vars()).size() >= 2 ) {
- found = true;
- break;
- }
- if( found )
- break;
- }
- if( found ) {
- cout << "Merging factors " << I->vars() << " and " << J->vars() << endl;
- *I *= *J;
- P.erase(J);
- }
- }
-}
-
-
-Factor FactorGraph::ExactMarginal(const VarSet & x) const {
- Factor P;
- for( size_t I = 0; I < nrFactors(); I++ )
- P *= factor(I);
- return P.marginal(x);
-}
-
-
-Real FactorGraph::ExactlogZ() const {
- Factor P;
- for( size_t I = 0; I < nrFactors(); I++ )
- P *= factor(I);
- return std::log(P.totalSum());
-}
-
-
-vector<VarSet> FactorGraph::Cliques() const {
- vector<VarSet> result;
-
- for( size_t I = 0; I < nrFactors(); I++ ) {
- bool maximal = true;
- for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
- if( (factor(J).vars() >> factor(I).vars()) && !(factor(J).vars() == factor(I).vars()) )
- maximal = false;
-
- if( maximal )
- result.push_back( factor(I).vars() );
- }
-
- return result;
-}
-
-
-void FactorGraph::clamp( const Var & n, size_t i ) {
- assert( i <= n.states() );
-
-/* if( do_surgery ) {
- size_t ni = find( vars().begin(), vars().end(), n) - vars().begin();
-
- if( ni != nrVars() ) {
- for( _nb_cit I = nb1(ni).begin(); I != nb1(ni).end(); I++ ) {
- if( factor(*I).size() == 1 )
- // Remove this single-variable factor
- // I = (_V2.erase(I))--;
- _E12.erase( _E12.begin() + VV2E(ni, *I) );
- else {
- // Replace it by the slice
- Index ind_I_min_n( factor(*I), factor(*I) / n );
- Index ind_n( factor(*I), n );
- Factor slice_I( factor(*I) / n );
- for( size_t ind_I = 0; ind_I < factor(*I).stateSpace(); ++ind_I, ++ind_I_min_n, ++ind_n )
- if( ind_n == i )
- slice_I[ind_I_min_n] = factor(*I)[ind_I];
- factor(*I) = slice_I;
-
- // Remove the edge between n and I
- _E12.erase( _E12.begin() + VV2E(ni, *I) );
- }
- }
-
- Regenerate();
-
- // remove all unconnected factors
- for( size_t I = 0; I < nrFactors(); I++ )
- if( nb2(I).size() == 0 )
- DeleteFactor(I--);
-
- DeleteVar( ni );
-
- // FIXME
- }
- } */
-
- // The cheap solution (at least in terms of coding time) is to multiply every factor
- // that contains the variable with a delta function
-
- Factor delta_n_i(n,0.0);
- delta_n_i[i] = 1.0;
-
- // For all factors that contain n
- for( size_t I = 0; I < nrFactors(); I++ )
- if( factor(I).vars() && n )
- // Multiply it with a delta function
- factor(I) *= delta_n_i;
-
- return;
-}
-
-
-void FactorGraph::saveProb( size_t I ) {
- map<size_t,Prob>::iterator it = _undoProbs.find( I );
- if( it != _undoProbs.end() )
- cout << "FactorGraph::saveProb: WARNING: _undoProbs[I] already defined!" << endl;
- _undoProbs[I] = factor(I).p();
-}
-
-
-void FactorGraph::undoProb( size_t I ) {
- map<size_t,Prob>::iterator it = _undoProbs.find( I );
- if( it != _undoProbs.end() ) {
- factor(I).p() = (*it).second;
- _undoProbs.erase(it);
- }
-}
-
-
-void FactorGraph::saveProbs( const VarSet &ns ) {
- if( !_undoProbs.empty() )
- cout << "FactorGraph::saveProbs: WARNING: _undoProbs not empy!" << endl;
- for( size_t I = 0; I < nrFactors(); I++ )
- if( factor(I).vars() && ns )
- _undoProbs[I] = factor(I).p();
-}
-
-
-void FactorGraph::undoProbs( const VarSet &ns ) {
- for( map<size_t,Prob>::iterator uI = _undoProbs.begin(); uI != _undoProbs.end(); ) {
- if( factor((*uI).first).vars() && ns ) {
-// cout << "undoing " << factor((*uI).first).vars() << endl;
-// cout << "from " << factor((*uI).first).p() << " to " << (*uI).second << endl;
- factor((*uI).first).p() = (*uI).second;
- _undoProbs.erase(uI++);
- } else
- uI++;
- }
-}
-
-
-bool FactorGraph::isConnected() const {
- if( nrVars() == 0 )
- return false;
- else {
- Var n = var( 0 );
-
- VarSet component = n;
-
- VarSet remaining;
- for( size_t i = 1; i < nrVars(); i++ )
- remaining |= var(i);
-
- bool found_new_vars = true;
- while( found_new_vars ) {
- VarSet new_vars;
- for( VarSet::const_iterator m = remaining.begin(); m != remaining.end(); m++ )
- if( delta(*m) && component )
- new_vars |= *m;
-
- if( new_vars.empty() )
- found_new_vars = false;
- else
- found_new_vars = true;
-
- component |= new_vars;
- remaining /= new_vars;
- };
- return remaining.empty();
- }
-}
-
-
-}
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __FACTORGRAPH_H__
-#define __FACTORGRAPH_H__
-
-
-#include <iostream>
-#include <map>
-#include "bipgraph.h"
-#include "factor.h"
-
-#include <tr1/unordered_map>
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-bool hasShortLoops(const vector<Factor> &P);
-void RemoveShortLoops(vector<Factor> &P);
-
-class FactorGraph : public BipartiteGraph<Var,Factor> {
- protected:
- map<size_t,Prob> _undoProbs;
- Prob::NormType _normtype;
-
- public:
- /// Default constructor
- FactorGraph() : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {};
- /// Copy constructor
- FactorGraph(const FactorGraph & x) : BipartiteGraph<Var,Factor>(x), _undoProbs(), _normtype(x._normtype) {};
- /// Construct FactorGraph from vector of Factors
- FactorGraph(const vector<Factor> &P);
- // Construct a FactorGraph from given factor and variable iterators
- template<typename FactorInputIterator, typename VarInputIterator>
- FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint = 0, size_t nr_var_hint = 0 );
-
- /// Assignment operator
- FactorGraph & operator=(const FactorGraph & x) {
- if(this!=&x) {
- BipartiteGraph<Var,Factor>::operator=(x);
- _undoProbs = x._undoProbs;
- _normtype = x._normtype;
- }
- return *this;
- }
- virtual ~FactorGraph() {}
-
- // aliases
- Var & var(size_t i) { return V1(i); }
- const Var & var(size_t i) const { return V1(i); }
- const vector<Var> & vars() const { return V1s(); }
- vector<Var> & vars() { return V1s(); }
- size_t nrVars() const { return V1s().size(); }
- Factor & factor(size_t I) { return V2(I); }
- const Factor & factor(size_t I) const { return V2(I); }
- const vector<Factor> & factors() const { return V2s(); }
- vector<Factor> & factors() { return V2s(); }
- size_t nrFactors() const { return V2s().size(); }
-
- /// Provides read access to neighbours of variable
- const _nb_t & nbV( size_t i1 ) const { return nb1(i1); }
- /// Provides full access to neighbours of variable
- _nb_t & nbV( size_t i1 ) { return nb1(i1); }
- /// Provides read access to neighbours of factor
- const _nb_t & nbF( size_t i2 ) const { return nb2(i2); }
- /// Provides full access to neighbours of factor
- _nb_t & nbF( size_t i2 ) { return nb2(i2); }
-
- size_t findVar(const Var & n) const {
- size_t i = find( vars().begin(), vars().end(), n ) - vars().begin();
- assert( i != nrVars() );
- return i;
- }
- size_t findFactor(const VarSet &ns) const {
- size_t I;
- for( I = 0; I < nrFactors(); I++ )
- if( factor(I).vars() == ns )
- break;
- assert( I != nrFactors() );
- return I;
- }
-
- friend ostream& operator << (ostream& os, const FactorGraph& fg);
- friend istream& operator >> (istream& is, FactorGraph& fg);
-
- VarSet delta(const Var & n) const;
- VarSet Delta(const Var & n) const;
- virtual void makeFactorCavity(size_t I);
- virtual void makeCavity(const Var & n);
-
- long ReadFromFile(const char *filename);
- long WriteToFile(const char *filename) const;
- long WriteToDotFile(const char *filename) const;
-
- Factor ExactMarginal(const VarSet & x) const;
- Real ExactlogZ() const;
-
- virtual void clamp( const Var & n, size_t i );
-
- bool hasNegatives() const;
- Prob::NormType NormType() const { return _normtype; }
-
- vector<VarSet> Cliques() const;
-
- virtual void undoProbs( const VarSet &ns );
- void saveProbs( const VarSet &ns );
- virtual void undoProb( size_t I );
- void saveProb( size_t I );
-
- bool isConnected() const;
-
- virtual void updatedFactor( size_t I ) {};
-
- private:
- /// Part of constructors (creates edges, neighbours and adjacency matrix)
- void createGraph( size_t nrEdges );
-};
-
-
-// assumes that the set of variables in [var_begin,var_end) is the union of the variables in the factors in [fact_begin, fact_end)
-template<typename FactorInputIterator, typename VarInputIterator>
-FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {
- // add factors
- size_t nrEdges = 0;
- V2s().reserve( nr_fact_hint );
- for( FactorInputIterator p2 = fact_begin; p2 != fact_end; ++p2 ) {
- V2s().push_back( *p2 );
- nrEdges += p2->vars().size();
- }
-
- // add variables
- V1s().reserve( nr_var_hint );
- for( VarInputIterator p1 = var_begin; p1 != var_end; ++p1 )
- V1s().push_back( *p1 );
-
- // create graph structure
- createGraph( nrEdges );
-}
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include <map>
-#include "hak.h"
-#include "util.h"
-#include "diffs.h"
-
-
-namespace dai {
-
-
-const char *HAK::Name = "HAK";
-
-
-bool HAK::checkProperties() {
- if( !HasProperty("tol") )
- return false;
- if (!HasProperty("maxiter") )
- return false;
- if (!HasProperty("verbose") )
- return false;
- if( !HasProperty("doubleloop") )
- return false;
- if( !HasProperty("clusters") )
- return false;
-
- ConvertPropertyTo<double>("tol");
- ConvertPropertyTo<size_t>("maxiter");
- ConvertPropertyTo<size_t>("verbose");
- ConvertPropertyTo<bool>("doubleloop");
- ConvertPropertyTo<ClustersType>("clusters");
-
- if( HasProperty("loopdepth") )
- ConvertPropertyTo<size_t>("loopdepth");
- else if( Clusters() == ClustersType::LOOP )
- return false;
-
- return true;
-}
-
-
-void HAK::constructMessages() {
- // Create outer beliefs
- _Qa.clear();
- _Qa.reserve(nr_ORs());
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
- _Qa.push_back( Factor( OR(alpha).vars() ) );
-
- // Create inner beliefs
- _Qb.clear();
- _Qb.reserve(nr_IRs());
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- _Qb.push_back( Factor( IR(beta) ) );
-
- // Create messages
- _muab.clear();
- _muab.reserve(nr_Redges());
- _muba.clear();
- _muba.reserve(nr_Redges());
- for( vector<R_edge_t>::const_iterator ab = Redges().begin(); ab != Redges().end(); ab++ ) {
- _muab.push_back( Factor( IR(ab->second) ) );
- _muba.push_back( Factor( IR(ab->second) ) );
- }
-}
-
-
-HAK::HAK(const RegionGraph & rg, const Properties &opts) : DAIAlgRG(rg, opts) {
- assert( checkProperties() );
-
- constructMessages();
-}
-
-
-void HAK::findLoopClusters( const FactorGraph & fg, set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
- for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
- VarSet ind = fg.delta( *in );
- if( (newcl.size()) >= 2 && (ind >> root) ) {
- allcl.insert( newcl | *in );
- }
- else if( length > 1 )
- findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
- }
-}
-
-
-HAK::HAK(const FactorGraph & fg, const Properties &opts) : DAIAlgRG(opts) {
- assert( checkProperties() );
-
- vector<VarSet> cl;
- if( Clusters() == ClustersType::MIN ) {
- cl = fg.Cliques();
- } else if( Clusters() == ClustersType::DELTA ) {
- for( size_t i = 0; i < fg.nrVars(); i++ )
- cl.push_back(fg.Delta(fg.var(i)));
- } else if( Clusters() == ClustersType::LOOP ) {
- cl = fg.Cliques();
- set<VarSet> scl;
- for( vector<Var>::const_iterator i0 = fg.vars().begin(); i0 != fg.vars().end(); i0++ ) {
- VarSet i0d = fg.delta(*i0);
- if( LoopDepth() > 1 )
- findLoopClusters( fg, scl, *i0, *i0, LoopDepth() - 1, fg.delta(*i0) );
- }
- for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
- cl.push_back(*c);
- if( Verbose() >= 3 ) {
- cout << "HAK uses the following clusters: " << endl;
- for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
- cout << *cli << endl;
- }
- } else
- throw "Invalid Clusters type";
-
- RegionGraph rg(fg,cl);
- RegionGraph::operator=(rg);
- constructMessages();
-
- if( Verbose() >= 3 )
- cout << "HAK regiongraph: " << *this << endl;
-}
-
-
-string HAK::identify() const {
- stringstream result (stringstream::out);
- result << Name << GetProperties();
- return result.str();
-}
-
-
-void HAK::init( const VarSet &ns ) {
- for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
- if( alpha->vars() && ns )
- alpha->fill( 1.0 / alpha->stateSpace() );
-
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- if( IR(beta) && ns ) {
- _Qb[beta].fill( 1.0 / IR(beta).stateSpace() );
- for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ ) {
- muab(*alpha,beta).fill( 1.0 / IR(beta).stateSpace() );
- muba(beta,*alpha).fill( 1.0 / IR(beta).stateSpace() );
- }
- }
-}
-
-
-void HAK::init() {
- assert( checkProperties() );
-
- for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
- alpha->fill( 1.0 / alpha->stateSpace() );
-
- for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
- beta->fill( 1.0 / beta->stateSpace() );
-
- for( size_t ab = 0; ab < nr_Redges(); ab++ ) {
- _muab[ab].fill( 1.0 / _muab[ab].stateSpace() );
- _muba[ab].fill( 1.0 / _muba[ab].stateSpace() );
- }
-}
-
-
-double HAK::doGBP() {
- if( Verbose() >= 1 )
- cout << "Starting " << identify() << "...";
- if( Verbose() >= 3)
- cout << endl;
-
- clock_t tic = toc();
-
- // Check whether counting numbers won't lead to problems
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- assert( nbIR(beta).size() + IR(beta).c() != 0.0 );
-
- // Keep old beliefs to check convergence
- vector<Factor> old_beliefs;
- old_beliefs.reserve( nrVars() );
- for( size_t i = 0; i < nrVars(); i++ )
- old_beliefs.push_back( belief( var(i) ) );
-
- // Differences in single node beliefs
- Diffs diffs(nrVars(), 1.0);
-
- 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( size_t beta = 0; beta < nr_IRs(); beta++ ) {
- for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
- muab(*alpha,beta) = _Qa[*alpha].marginal(IR(beta)).divided_by( muba(beta,*alpha) );
-
- Factor Qb_new;
- for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
- Qb_new *= muab(*alpha,beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
- Qb_new.normalize( _normtype );
- if( Qb_new.hasNaNs() ) {
- cout << "HAK::doGBP: Qb_new has NaNs!" << endl;
- return NAN;
- }
-// _Qb[beta] = Qb_new.makeZero(1e-100); // damping?
- _Qb[beta] = Qb_new;
-
- for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ ) {
- muba(beta,*alpha) = _Qb[beta].divided_by( muab(*alpha,beta) );
-
- Factor Qa_new = OR(*alpha);
- for( R_nb_cit gamma = nbOR(*alpha).begin(); gamma != nbOR(*alpha).end(); gamma++ )
- Qa_new *= muba(*gamma,*alpha);
- Qa_new ^= (1.0 / OR(*alpha).c());
- Qa_new.normalize( _normtype );
- if( Qa_new.hasNaNs() ) {
- cout << "HAK::doGBP: Qa_new has NaNs!" << endl;
- return NAN;
- }
-// _Qa[*alpha] = Qa_new.makeZero(1e-100); // damping?
- _Qa[*alpha] = Qa_new;
- }
- }
-
- // Calculate new single variable beliefs and compare with old ones
- for( size_t i = 0; i < nrVars(); i++ ) {
- Factor new_belief = belief( var( i ) );
- diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
- old_beliefs[i] = new_belief;
- }
-
- if( Verbose() >= 3 )
- cout << "HAK::doGBP: maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
- }
-
- updateMaxDiff( diffs.max() );
-
- if( Verbose() >= 1 ) {
- if( diffs.max() > Tol() ) {
- if( Verbose() == 1 )
- cout << endl;
- cout << "HAK::doGBP: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
- } else {
- if( Verbose() >= 2 )
- cout << "HAK::doGBP: ";
- cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
- }
- }
-
- return diffs.max();
-}
-
-
-double HAK::doDoubleLoop() {
- if( Verbose() >= 1 )
- cout << "Starting " << identify() << "...";
- if( Verbose() >= 3)
- cout << endl;
-
- clock_t tic = toc();
-
- // Save original outer regions
- vector<FRegion> org_ORs = ORs();
-
- // Save original inner counting numbers and set negative counting numbers to zero
- vector<double> org_IR_cs( nr_IRs(), 0.0 );
- for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
- org_IR_cs[beta] = IR(beta).c();
- if( IR(beta).c() < 0.0 )
- IR(beta).c() = 0.0;
- }
-
- // Keep old beliefs to check convergence
- vector<Factor> old_beliefs;
- old_beliefs.reserve( nrVars() );
- for( size_t i = 0; i < nrVars(); i++ )
- old_beliefs.push_back( belief( var(i) ) );
-
- // Differences in single node beliefs
- Diffs diffs(nrVars(), 1.0);
-
- size_t outer_maxiter = MaxIter();
- double outer_tol = Tol();
- size_t outer_verbose = Verbose();
- double org_maxdiff = MaxDiff();
-
- // Set parameters for inner loop
- MaxIter( 5 );
- Verbose( outer_verbose ? outer_verbose - 1 : 0 );
-
- size_t outer_iter = 0;
- for( outer_iter = 0; outer_iter < outer_maxiter && diffs.max() > outer_tol; outer_iter++ ) {
- // Calculate new outer regions
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
- OR(alpha) = org_ORs[alpha];
- for( R_nb_cit beta = nbOR(alpha).begin(); beta != nbOR(alpha).end(); beta++ )
- OR(alpha) *= _Qb[*beta] ^ ((IR(*beta).c() - org_IR_cs[*beta]) / nbIR(*beta).size());
- }
-
- // Inner loop
- if( isnan( doGBP() ) )
- return NAN;
-
- // Calculate new single variable beliefs and compare with old ones
- for( size_t i = 0; i < nrVars(); i++ ) {
- Factor new_belief = belief( var( i ) );
- diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
- old_beliefs[i] = new_belief;
- }
-
- if( Verbose() >= 3 )
- cout << "HAK::doDoubleLoop: maxdiff " << diffs.max() << " after " << outer_iter+1 << " passes" << endl;
- }
-
- // restore _maxiter, _verbose and _maxdiff
- MaxIter( outer_maxiter );
- Verbose( outer_verbose );
- MaxDiff( org_maxdiff );
-
- updateMaxDiff( diffs.max() );
-
- // Restore original outer regions
- ORs() = org_ORs;
-
- // Restore original inner counting numbers
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- IR(beta).c() = org_IR_cs[beta];
-
- if( Verbose() >= 1 ) {
- if( diffs.max() > Tol() ) {
- if( Verbose() == 1 )
- cout << endl;
- cout << "HAK::doDoubleLoop: WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
- } else {
- if( Verbose() >= 3 )
- cout << "HAK::doDoubleLoop: ";
- cout << "converged in " << outer_iter << " passes (" << toc() - tic << " clocks)." << endl;
- }
- }
-
- return diffs.max();
-}
-
-
-double HAK::run() {
- if( DoubleLoop() )
- return doDoubleLoop();
- else
- return doGBP();
-}
-
-
-Factor HAK::belief( const VarSet &ns ) const {
- vector<Factor>::const_iterator beta;
- for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
- if( beta->vars() >> ns )
- break;
- if( beta != _Qb.end() )
- return( beta->marginal(ns) );
- else {
- vector<Factor>::const_iterator alpha;
- for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
- if( alpha->vars() >> ns )
- break;
- assert( alpha != _Qa.end() );
- return( alpha->marginal(ns) );
- }
-}
-
-
-Factor HAK::belief( const Var &n ) const {
- return belief( (VarSet)n );
-}
-
-
-vector<Factor> HAK::beliefs() const {
- vector<Factor> result;
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- result.push_back( Qb(beta) );
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
- result.push_back( Qa(alpha) );
- return result;
-}
-
-
-Complex HAK::logZ() const {
- Complex sum = 0.0;
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- sum += Complex(IR(beta).c()) * Qb(beta).entropy();
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
- sum += Complex(OR(alpha).c()) * Qa(alpha).entropy();
- sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
- }
- return sum;
-}
-
-
-}
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __HAK_H__
-#define __HAK_H__
-
-
-#include "daialg.h"
-#include "regiongraph.h"
-#include "enum.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-/// HAK provides an implementation of the single and double-loop algorithms by Heskes, Albers and Kappen
-class HAK : public DAIAlgRG {
- protected:
- vector<Factor> _Qa;
- vector<Factor> _Qb;
- vector<Factor> _muab;
- vector<Factor> _muba;
-
- public:
- /// Default constructor
- HAK() : DAIAlgRG() {};
-
- /// Copy constructor
- HAK(const HAK & x) : DAIAlgRG(x), _Qa(x._Qa), _Qb(x._Qb), _muab(x._muab), _muba(x._muba) {};
-
- /// Clone function
- HAK* clone() const { return new HAK(*this); }
-
- /// Construct from RegionGraph
- HAK(const RegionGraph & rg, const Properties &opts);
-
- /// Construct from RactorGraph using "clusters" option
- HAK(const FactorGraph & fg, const Properties &opts);
-
- /// Assignment operator
- HAK & operator=(const HAK & x) {
- if( this != &x ) {
- DAIAlgRG::operator=(x);
- _Qa = x._Qa;
- _Qb = x._Qb;
- _muab = x._muab;
- _muba = x._muba;
- }
- return *this;
- }
-
- static const char *Name;
-
- ENUM3(ClustersType,MIN,DELTA,LOOP)
- ClustersType Clusters() const { return GetPropertyAs<ClustersType>("clusters"); }
- bool DoubleLoop() { return GetPropertyAs<bool>("doubleloop"); }
- size_t LoopDepth() { return GetPropertyAs<size_t>("loopdepth"); }
-
- Factor & muab( size_t alpha, size_t beta ) { return _muab[ORIR2E(alpha,beta)]; }
- Factor & muba( size_t beta, size_t alpha ) { return _muba[ORIR2E(alpha,beta)]; }
- const Factor& Qa( size_t alpha ) const { return _Qa[alpha]; };
- const Factor& Qb( size_t beta ) const { return _Qb[beta]; };
-
-// void Regenerate();
- double doGBP();
- double doDoubleLoop();
- double run();
- void init();
- string identify() const;
- Factor belief( const Var &n ) const;
- Factor belief( const VarSet &ns ) const;
- vector<Factor> beliefs() const;
- Complex logZ () const;
-
- void init( const VarSet &ns );
- void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
- bool checkProperties();
-
- private:
- void constructMessages();
- void findLoopClusters( const FactorGraph &fg, set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars );
-};
-
-
-}
-
-
-#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_alldai_h
+#define __defined_libdai_alldai_h
+
+
+#include <string>
+#include <dai/daialg.h>
+#include <dai/bp.h>
+#include <dai/lc.h>
+#include <dai/hak.h>
+#include <dai/mf.h>
+#include <dai/jtree.h>
+#include <dai/treeep.h>
+#include <dai/mr.h>
+
+
+namespace dai {
+
+
+/// newInfAlg constructs a new approximate inference algorithm named name for the
+/// FactorGraph fg with optionts opts and returns a pointer to the new object.
+/// The caller needs to delete it (maybe some sort of smart_ptr might be useful here).
+InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const Properties &opts );
+
+
+/// AINames contains the names of all approximate inference algorithms
+static const char* DAINames[] = {BP::Name, MF::Name, HAK::Name, LC::Name, TreeEP::Name, MR::Name, JTree::Name};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_bipgraph_h
+#define __defined_libdai_bipgraph_h
+
+
+#include <vector>
+#include <algorithm>
+#include <boost/numeric/ublas/vector_sparse.hpp>
+
+
+namespace dai {
+
+
+/// A BipartiteGraph represents a graph with two types of nodes
+template <class T1, class T2> class BipartiteGraph {
+ public:
+ typedef std::pair<size_t,size_t> _edge_t;
+ typedef std::vector<size_t> _nb_t;
+ typedef _nb_t::const_iterator _nb_cit;
+
+
+ private:
+ /// Vertices of first type
+ std::vector<T1> _V1;
+
+ /// Vertices of second type
+ std::vector<T2> _V2;
+
+ /// Edges, which are pairs (v1,v2) with v1 in _V1 and v2 in _V2
+ std::vector<_edge_t> _E12;
+
+ /// Conversion matrix that computes which index of _E12 corresponds to a vertex index pair (v1,v2)
+ std::vector<boost::numeric::ublas::mapped_vector<size_t> > _E12ind;
+
+ /// Neighbour indices of vertices of first type
+ std::vector<_nb_t> _nb1;
+
+ /// Neighbour indices of vertices of second type
+ std::vector<_nb_t> _nb2;
+
+
+ public:
+ /// Default constructor
+ BipartiteGraph<T1,T2> () {};
+
+ /// Copy constructor
+ BipartiteGraph<T1,T2> ( const BipartiteGraph<T1,T2> & x ) : _V1(x._V1), _V2(x._V2), _E12(x._E12), _E12ind(x._E12ind), _nb1(x._nb1), _nb2(x._nb2) {};
+
+ /// Assignment operator
+ BipartiteGraph<T1,T2> & operator=(const BipartiteGraph<T1,T2> & x) {
+ if( this != &x ) {
+ _V1 = x._V1;
+ _V2 = x._V2;
+ _E12 = x._E12;
+ _E12ind = x._E12ind;
+ _nb1 = x._nb1;
+ _nb2 = x._nb2;
+ }
+ return *this;
+ }
+
+ /// Provides read access to node of first type
+ const T1 & V1( size_t i1 ) const { return _V1[i1]; }
+ /// Provides full access to node of first type
+ T1 & V1( size_t i1 ) { return _V1[i1]; }
+ /// Provides read access to all nodes of first type
+ const std::vector<T1> & V1s() const { return _V1; }
+ /// Provides full access to all nodes of first type
+ std::vector<T1> & V1s() { return _V1; }
+
+ /// Provides read access to node of second type
+ const T2 & V2( size_t i2 ) const { return _V2[i2]; }
+ /// Provides full access to node of second type
+ T2 & V2( size_t i2 ) { return _V2[i2]; }
+ /// Provides read access to all nodes of second type
+ const std::vector<T2> & V2s() const { return _V2; }
+ /// Provides full access to all nodes of second type
+ std::vector<T2> & V2s() { return _V2; }
+
+ /// Provides read access to edge
+ const _edge_t & edge(size_t ind) const { return _E12[ind]; }
+ /// Provides full access to edge
+ _edge_t & edge(size_t ind) { return _E12[ind]; }
+ /// Provides read access to all edges
+ const std::vector<_edge_t> & edges() const { return _E12; }
+ /// Provides full access to all edges
+ std::vector<_edge_t> & edges() { return _E12; }
+ /// Returns number of edges
+ size_t nr_edges() const { return _E12.size(); }
+
+ /// Provides read access to neighbours of node of first type
+ const _nb_t & nb1( size_t i1 ) const { return _nb1[i1]; }
+ /// Provides full access to neighbours of node of first type
+ _nb_t & nb1( size_t i1 ) { return _nb1[i1]; }
+
+ /// Provides read access to neighbours of node of second type
+ const _nb_t & nb2( size_t i2 ) const { return _nb2[i2]; }
+ /// Provides full access to neighbours of node of second type
+ _nb_t & nb2( size_t i2 ) { return _nb2[i2]; }
+
+ /// Converts the pair of indices (i1,i2) to the corresponding edge index
+ size_t VV2E( const size_t i1, const size_t i2 ) const { return _E12ind[i1](i2); }
+
+
+ /// Regenerates internal structures (_E12ind, _nb1 and _nb2) based upon _V1, _V2 and _E12
+ void Regenerate() {
+ // Calculate _nb1 and _nb2
+
+ // Start with empty vectors
+ _nb1.clear();
+ _nb1.resize(_V1.size());
+ // Start with empty vectors
+ _nb2.clear();
+ _nb2.resize(_V2.size());
+ // Each edge yields a neighbour pair
+ for( std::vector<_edge_t>::const_iterator e = _E12.begin(); e != _E12.end(); e++ ) {
+ _nb1[e->first].push_back(e->second);
+ _nb2[e->second].push_back(e->first);
+ }
+ // Remove duplicates from _nb1
+ for( size_t i1 = 0; i1 < _V1.size(); i1++ ) {
+ _nb_t::iterator new_end = unique(_nb1[i1].begin(), _nb1[i1].end());
+ _nb1[i1].erase( new_end, _nb1[i1].end() );
+ }
+ // Remove duplicates from _nb2
+ for( size_t i2 = 0; i2 < _V2.size(); i2++ ) {
+ _nb_t::iterator new_end = unique(_nb2[i2].begin(), _nb2[i2].end());
+ _nb2[i2].erase( new_end, _nb2[i2].end() );
+ }
+
+ // Calculate _E12ind
+ _E12ind = std::vector<boost::numeric::ublas::mapped_vector<size_t> >(_V1.size(), boost::numeric::ublas::mapped_vector<size_t>(_V2.size()) );
+ for( size_t k = 0; k < _E12.size(); k++ )
+ _E12ind[_E12[k].first](_E12[k].second) = k;
+
+ }
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_bp_h
+#define __defined_libdai_bp_h
+
+
+#include <string>
+#include <dai/daialg.h>
+#include <dai/factorgraph.h>
+#include <dai/enum.h>
+
+
+namespace dai {
+
+
+class BP : public DAIAlgFG {
+ protected:
+ typedef std::vector<size_t> _ind_t;
+
+ std::vector<_ind_t> _indices;
+ std::vector<Prob> _messages, _newmessages;
+
+ public:
+ ENUM4(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+
+ // default constructor
+ BP() : DAIAlgFG() {};
+ // copy constructor
+ BP(const BP & x) : DAIAlgFG(x), _indices(x._indices), _messages(x._messages), _newmessages(x._newmessages) {};
+ BP* clone() const { return new BP(*this); }
+ // construct BP object from FactorGraph
+ BP(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
+ assert( checkProperties() );
+ Regenerate();
+ }
+ // assignment operator
+ BP & operator=(const BP & x) {
+ if(this!=&x) {
+ DAIAlgFG::operator=(x);
+ _messages = x._messages;
+ _newmessages = x._newmessages;
+ _indices = x._indices;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+
+ Prob & message(size_t i1, size_t i2) { return( _messages[VV2E(i1,i2)] ); }
+ const Prob & message(size_t i1, size_t i2) const { return( _messages[VV2E(i1,i2)] ); }
+ Prob & newMessage(size_t i1, size_t i2) { return( _newmessages[VV2E(i1,i2)] ); }
+ const Prob & newMessage(size_t i1, size_t i2) const { return( _newmessages[VV2E(i1,i2)] ); }
+ _ind_t & index(size_t i1, size_t i2) { return( _indices[VV2E(i1,i2)] ); }
+ const _ind_t & index(size_t i1, size_t i2) const { return( _indices[VV2E(i1,i2)] ); }
+
+ std::string identify() const;
+ void Regenerate();
+ void init();
+ void calcNewMessage(size_t iI);
+ double run();
+ Factor belief1 (size_t i) const;
+ Factor belief2 (size_t I) const;
+ Factor belief (const Var &n) const;
+ Factor belief (const VarSet &n) const;
+ std::vector<Factor> beliefs() const;
+ Complex logZ() const;
+
+ void init( const VarSet &ns );
+ void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
+ bool checkProperties();
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_clustergraph_h
+#define __defined_libdai_clustergraph_h
+
+
+#include <set>
+#include <vector>
+#include <dai/varset.h>
+
+
+namespace dai {
+
+
+/// A ClusterGraph is a hypergraph with VarSets as nodes.
+/// It is implemented as a set<VarSet> in which the adjacency
+/// relationship is computed realtime. It may be better to
+/// implement it as a bipartitegraph, though. The additional
+/// functionality compared to a simple set<VarSet> is
+/// finding maximal clusters, finding cliques, etc...
+class ClusterGraph : public std::set<VarSet> {
+ public:
+ /// Default constructor
+ ClusterGraph() : std::set<VarSet>() {}
+
+ /// Construct from vector<VarSet>
+ ClusterGraph( const std::vector<VarSet> & cls ) {
+ insert( cls.begin(), cls.end() );
+ }
+
+ /// Copy constructor
+ ClusterGraph(const ClusterGraph &x) : std::set<VarSet>(x) {}
+
+ /// Assignment operator
+ ClusterGraph& operator=( const ClusterGraph &x ) {
+ if( this != &x ) {
+ std::set<VarSet>::operator=( x );
+ }
+ return *this;
+ }
+
+ /// Returns true if ns is a maximal member of *this under inclusion (VarSet::operator<<)
+ bool isMaximal( const VarSet &ns ) const {
+ if( count( ns ) ) {
+ // ns is a member
+ bool maximal = true;
+ for( const_iterator x = begin(); x != end() && maximal; x++ )
+ if( (ns << *x) && (ns != *x) )
+ maximal = false;
+ return maximal;
+ } else
+ return false;
+ }
+
+ /// Erase all VarSets that are not maximal
+ ClusterGraph& eraseNonMaximal() {
+ for( iterator x = begin(); x != end(); )
+ if( !isMaximal(*x) )
+ erase(x++);
+ else
+ x++;
+ return *this;
+ }
+
+ /// Return union of all members
+ VarSet vars() const {
+ VarSet result;
+ for( iterator x = begin(); x != end(); x++ )
+ result |= *x;
+ return result;
+ }
+
+ /// Returns true if n1 and n2 are adjacent, i.e.\ by
+ /// definition, are both contained in some cluster in *this
+ bool adj( const Var& n1, const Var& n2 ) {
+ bool result = false;
+ for( iterator x = begin(); (x != end()) && (!result); x++ )
+ if( (*x && n1) && (*x && n2) )
+ result = true;
+ return result;
+ }
+
+ /// Returns union of clusters that contain n, minus n
+ VarSet Delta( const Var& n ) const {
+ VarSet result;
+ for( iterator x = begin(); x != end(); x++ )
+ if( (*x && n) )
+ result |= *x;
+ return result;
+ }
+
+ /// Returns union of clusters that contain n, minus n
+ VarSet delta( const Var& n ) const {
+ return Delta( n ) / n;
+ }
+
+ /// Erases all members that contain n
+ ClusterGraph& eraseSubsuming( const Var& n ) {
+ for( iterator x = begin(); x != end(); )
+ if( (*x && n) )
+ erase(x++);
+ else
+ x++;
+ return *this;
+ }
+
+ /// Send to output stream
+ friend std::ostream & operator << ( std::ostream & os, const ClusterGraph & cl ) {
+ os << "{";
+ ClusterGraph::const_iterator x = cl.begin();
+ if( x != cl.end() )
+ os << *(x++);
+ for( ; x != cl.end(); x++ )
+ os << ", " << *x;
+ os << "}";
+ return os;
+ }
+
+ /// Convert to vector<VarSet>
+ std::vector<VarSet> toVector() const {
+ std::vector<VarSet> result;
+ result.reserve( size() );
+ for( const_iterator x = begin(); x != end(); x++ )
+ result.push_back( *x );
+ return result;
+ }
+
+ /// Calculate cost of eliminating variable n
+ /// using as a measure "number of added edges in the adjacency graph"
+ /// where the adjacency graph has the variables as its nodes and
+ /// connects nodes i1 and i2 iff i1 and i2 occur in some common factor I
+ size_t eliminationCost( const Var& n ) {
+ VarSet d_n = delta( n );
+ size_t cost = 0;
+
+ // for each unordered pair {i1,i2} adjacent to n
+ for( VarSet::const_iterator i1 = d_n.begin(); i1 != d_n.end(); i1++ ) {
+ VarSet d_i1 = delta( *i1 );
+ for( VarSet::const_iterator i2 = i1; (++i2) != d_n.end(); ) {
+ // if i1 and i2 are not adjacent, eliminating n would make them adjacent
+ if( !adj(*i1, *i2) )
+ cost++;
+ }
+ }
+ return cost;
+ }
+
+ /// Perform Variable Elimination without Probs, i.e. only keeping track of
+ /// the interactions that are created along the way.
+ /// Input: a set of outer clusters and an elimination sequence
+ /// Output: a set of elimination "cliques"
+ ClusterGraph VarElim( const std::vector<Var> &ElimSeq ) const;
+
+ /// As Taylan does it
+ ClusterGraph VarElim_MinFill() const;
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_daialg_h
+#define __defined_libdai_daialg_h
+
+
+#include <string>
+#include <iostream>
+#include <vector>
+#include <dai/factorgraph.h>
+#include <dai/regiongraph.h>
+#include <dai/properties.h>
+
+
+namespace dai {
+
+
+/// The InfAlg class is the common denominator of the various approximate inference algorithms.
+/// A InfAlg object represents a discrete factorized probability distribution over multiple variables
+/// together with an inference algorithm.
+class InfAlg {
+ private:
+ /// Properties of the algorithm (replaces _tol, _maxiter, _verbose)
+ Properties _properties;
+
+ /// Maximum difference encountered so far
+ double _maxdiff;
+
+
+ public:
+ /// Default constructor
+ InfAlg() : _properties(), _maxdiff(0.0) {}
+
+ /// Constructor with options
+ InfAlg( const Properties &opts ) : _properties(opts), _maxdiff(0.0) {}
+
+ /// Copy constructor
+ InfAlg( const InfAlg & x ) : _properties(x._properties), _maxdiff(x._maxdiff) {}
+
+ /// Clone (virtual copy constructor)
+ virtual InfAlg* clone() const = 0;
+
+ /// Assignment operator
+ InfAlg & operator=( const InfAlg & x ) {
+ if( this != &x ) {
+ _properties = x._properties;
+ _maxdiff = x._maxdiff;
+ }
+ return *this;
+ }
+
+ /// Virtual desctructor
+ // (this is needed because this class contains virtual functions)
+ virtual ~InfAlg() {}
+
+ /// Returns true if a property with the given key is present
+ bool HasProperty(const PropertyKey &key) const { return _properties.hasKey(key); }
+
+ /// Gets a property
+ const PropertyValue & GetProperty(const PropertyKey &key) const { return _properties.Get(key); }
+
+ /// Gets a property, casted as ValueType
+ template<typename ValueType>
+ ValueType GetPropertyAs(const PropertyKey &key) const { return _properties.GetAs<ValueType>(key); }
+
+ /// Sets a property
+ void SetProperty(const PropertyKey &key, const PropertyValue &val) { _properties[key] = val; }
+
+ /// Converts a property from string to ValueType, if necessary
+ template<typename ValueType>
+ void ConvertPropertyTo(const PropertyKey &key) { _properties.ConvertTo<ValueType>(key); }
+
+ /// Gets all properties
+ const Properties & GetProperties() const { return _properties; }
+
+ /// Sets properties
+ void SetProperties(const Properties &p) { _properties = p; }
+
+ /// Sets tolerance
+ void Tol( double tol ) { SetProperty("tol", tol); }
+ /// Gets tolerance
+ double Tol() const { return GetPropertyAs<double>("tol"); }
+
+ /// Sets maximum number of iterations
+ void MaxIter( size_t maxiter ) { SetProperty("maxiter", maxiter); }
+ /// Gets maximum number of iterations
+ size_t MaxIter() const { return GetPropertyAs<size_t>("maxiter"); }
+
+ /// Sets verbosity
+ void Verbose( size_t verbose ) { SetProperty("verbose", verbose); }
+ /// Gets verbosity
+ size_t Verbose() const { return GetPropertyAs<size_t>("verbose"); }
+
+ /// Sets maximum difference encountered so far
+ void MaxDiff( double maxdiff ) { _maxdiff = maxdiff; }
+ /// Gets maximum difference encountered so far
+ double MaxDiff() const { return _maxdiff; }
+ /// Updates maximum difference encountered so far
+ void updateMaxDiff( double maxdiff ) { if( maxdiff > _maxdiff ) _maxdiff = maxdiff; }
+ /// Sets maximum difference encountered so far to zero
+ void clearMaxDiff() { _maxdiff = 0.0; }
+
+ /// Identifies itself for logging purposes
+ virtual std::string identify() const = 0;
+
+ /// Get single node belief
+ virtual Factor belief( const Var &n ) const = 0;
+
+ /// Get general belief
+ virtual Factor belief( const VarSet &n ) const = 0;
+
+ /// Get all beliefs
+ virtual std::vector<Factor> beliefs() const = 0;
+
+ /// Get log partition sum
+ virtual Complex logZ() const = 0;
+
+ /// Clear messages and beliefs
+ virtual void init() = 0;
+
+ /// The actual approximate inference algorithm
+ virtual double run() = 0;
+
+ /// Save factor I
+ virtual void saveProb( size_t I ) = 0;
+ /// Save Factors involving ns
+ virtual void saveProbs( const VarSet &ns ) = 0;
+
+ /// Restore factor I
+ virtual void undoProb( size_t I ) = 0;
+ /// Restore Factors involving ns
+ virtual void undoProbs( const VarSet &ns ) = 0;
+
+ /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
+ virtual void clamp( const Var & n, size_t i ) = 0;
+
+ /// Return all variables that interact with n
+ virtual VarSet delta( const Var & n ) const = 0;
+
+ /// Set all factors interacting with n to 1
+ virtual void makeCavity( const Var & n ) = 0;
+
+ /// Set factor I to 1
+ virtual void makeFactorCavity( size_t I ) = 0;
+
+ /// Get index of variable n
+ virtual size_t findVar( const Var & n ) const = 0;
+
+ /// Get index of first factor involving ns
+ virtual size_t findFactor( const VarSet &ns ) const = 0;
+
+ /// Get number of variables
+ virtual size_t nrVars() const = 0;
+
+ /// Get number of factors
+ virtual size_t nrFactors() const = 0;
+
+ /// Get const reference to variable i
+ virtual const Var & var(size_t i) const = 0;
+
+ /// Get reference to variable i
+ virtual Var & var(size_t i) = 0;
+
+ /// Get const reference to factor I
+ virtual const Factor & factor( size_t I ) const = 0;
+
+ /// Get reference to factor I
+ virtual Factor & factor( size_t I ) = 0;
+
+ /// Factor I has been updated
+ virtual void updatedFactor( size_t I ) = 0;
+
+ /// Checks whether all necessary properties have been set
+ /// and casts string-valued properties to other values if necessary
+ virtual bool checkProperties() = 0;
+};
+
+
+template <class T>
+class DAIAlg : public InfAlg, public T {
+ public:
+ /// Default constructor
+ DAIAlg() : InfAlg(), T() {}
+
+ /// Construct DAIAlg with empty T but using the specified properties
+ DAIAlg( const Properties &opts ) : InfAlg( opts ), T() {}
+
+ /// Construct DAIAlg using the specified properties
+ DAIAlg( const T & t, const Properties &opts ) : InfAlg( opts ), T(t) {}
+
+ /// Copy constructor
+ DAIAlg( const DAIAlg & x ) : InfAlg(x), T(x) {}
+
+ /// Save factor I (using T::saveProb)
+ void saveProb( size_t I ) { T::saveProb( I ); }
+ /// Save Factors involving ns (using T::saveProbs)
+ void saveProbs( const VarSet &ns ) { T::saveProbs( ns ); }
+
+ /// Restore factor I (using T::undoProb)
+ void undoProb( size_t I ) { T::undoProb( I ); }
+ /// Restore Factors involving ns (using T::undoProbs)
+ void undoProbs( const VarSet &ns ) { T::undoProbs( ns ); }
+
+ /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$) (using T::clamp)
+ void clamp( const Var & n, size_t i ) { T::clamp( n, i ); }
+
+ /// Return all variables that interact with n (using T::delta)
+ VarSet delta( const Var & n ) const { return T::delta(n); }
+
+ /// Set all factors interacting with n to 1 (using T::makeCavity)
+ void makeCavity( const Var & n ) { T::makeCavity(n); }
+
+ /// Set factor I to 1 (using T::makeFactorCavity)
+ void makeFactorCavity( size_t I ) { T::makeFactorCavity(I); }
+
+ /// Get index of variable n (using T::findVar)
+ size_t findVar( const Var & n ) const { return T::findVar(n); }
+
+ /// Get index of first factor involving ns (using T::findFactor)
+ size_t findFactor( const VarSet &ns ) const { return T::findFactor(ns); }
+
+ /// Get number of variables (using T::nrFactors)
+ size_t nrVars() const { return T::nrVars(); }
+
+ /// Get number of factors (using T::nrFactors)
+ size_t nrFactors() const { return T::nrFactors(); }
+
+ /// Get const reference to variable i (using T::var)
+ const Var & var( size_t i ) const { return T::var(i); }
+
+ /// Get reference to variable i (using T::var)
+ Var & var(size_t i) { return T::var(i); }
+
+ /// Get const reference to factor I (using T::factor)
+ const Factor & factor( size_t I ) const { return T::factor(I); }
+
+ /// Get reference to factor I (using T::factor)
+ Factor & factor( size_t I ) { return T::factor(I); }
+
+ /// Factor I has been updated (using T::updatedFactor)
+ void updatedFactor( size_t I ) { T::updatedFactor(I); }
+};
+
+
+typedef DAIAlg<FactorGraph> DAIAlgFG;
+typedef DAIAlg<RegionGraph> DAIAlgRG;
+
+
+/// Calculate the marginal of obj on ns by clamping
+/// all variables in ns and calculating logZ for each joined state
+Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit );
+
+
+/// Calculate beliefs of all pairs in ns (by clamping
+/// nodes in ns and calculating logZ and the beliefs for each state)
+std::vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit );
+
+
+/// Calculate beliefs of all pairs in ns (by clamping
+/// pairs in ns and calculating logZ for each joined state)
+std::vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit );
+
+
+/// Calculate 2nd order interactions of the marginal of obj on ns
+Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit );
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_diffs_h
+#define __defined_libdai_diffs_h
+
+
+#include <vector>
+
+
+namespace dai {
+
+
+class Diffs : public std::vector<double> {
+ private:
+ size_t _maxsize;
+ double _def;
+ std::vector<double>::iterator _pos;
+ std::vector<double>::iterator _maxpos;
+ public:
+ Diffs(long maxsize, double def) : std::vector<double>(), _maxsize(maxsize), _def(def) {
+ this->reserve(_maxsize);
+ _pos = begin();
+ _maxpos = begin();
+ };
+ double max() {
+ if( size() < _maxsize )
+ return _def;
+ else
+ return( *_maxpos );
+ }
+ void push(double x) {
+ if( size() < _maxsize ) {
+ push_back(x);
+ _pos = end();
+ _maxpos = max_element(begin(),end());
+ }
+ else {
+ if( _pos == end() )
+ _pos = begin();
+ if( _maxpos == _pos ) {
+ *_pos++ = x;
+ _maxpos = max_element(begin(),end());
+ } else {
+ if( x > *_maxpos )
+ _maxpos = _pos;
+ *_pos++ = x;
+ }
+ }
+ }
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_enum_h
+#define __defined_libdai_enum_h
+
+
+#include <cstring>
+#include <iostream>
+
+
+namespace dai {
+
+
+// C++ enums are too limited for my purposes. This defines wrapper classes
+// that provide much more functionality than a simple enum. The only
+// disadvantage is that one wrapper class needs to be written for each
+// number of values an enum can take... a better solution is needed.
+
+
+#define ENUM2(x,a,b) class x {\
+ public:\
+ enum value {a, b};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM3(x,a,b,c) class x {\
+ public:\
+ enum value {a, b, c};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM4(x,a,b,c,d) class x {\
+ public:\
+ enum value {a, b, c, d};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c, #d};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c, #d};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM5(x,a,b,c,d,e) class x {\
+ public:\
+ enum value {a, b, c, d, e};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c, #d, #e};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c, #d, #e};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM6(x,a,b,c,d,e,f) class x {\
+ public:\
+ enum value {a, b, c, d, e, f};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c, #d, #e, #f};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c, #d, #e, #f};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_factor_h
+#define __defined_libdai_factor_h
+
+
+#include <iostream>
+#include <cmath>
+#include <dai/prob.h>
+#include <dai/varset.h>
+#include <dai/index.h>
+
+
+namespace dai {
+
+
+template<typename T> class TFactor;
+typedef TFactor<Real> Factor;
+typedef TFactor<Complex> CFactor;
+
+
+// predefine friends
+template<typename T> Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt );
+template<typename T> Complex KL_dist( const TFactor<T> & p, const TFactor<T> & q );
+template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P);
+
+
+// T should be castable from and to double and to complex
+template <typename T> class TFactor {
+ protected:
+ VarSet _vs;
+ TProb<T> _p;
+
+ public:
+ // Default constructor
+ TFactor () : _vs(), _p(1,1.0) {}
+
+ // Construct Factor from VarSet
+ TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.stateSpace()) {}
+
+ // Construct Factor from VarSet and initial value
+ TFactor( const VarSet& ns, Real p ) : _vs(ns), _p(_vs.stateSpace(),p) {}
+
+ // Construct Factor from VarSet and initial array
+ TFactor( const VarSet& ns, const Real* p ) : _vs(ns), _p(_vs.stateSpace(),p) {}
+
+ // Construct Factor from VarSet and TProb<T>
+ TFactor( const VarSet& ns, const TProb<T> p ) : _vs(ns), _p(p) {
+#ifdef DEBUG
+ assert( _vs.stateSpace() == _p.size() );
+#endif
+ }
+
+ // Construct Factor from Var
+ TFactor( const Var& n ) : _vs(n), _p(n.states()) {}
+
+ // Copy constructor
+ TFactor( const TFactor<T> &x ) : _vs(x._vs), _p(x._p) {}
+
+ // Assignment operator
+ TFactor<T> & operator= (const TFactor<T> &x) {
+ if( this != &x ) {
+ _vs = x._vs;
+ _p = x._p;
+ }
+ return *this;
+ }
+
+ const TProb<T> & p() const { return _p; }
+ TProb<T> & p() { return _p; }
+ const VarSet & vars() const { return _vs; }
+ size_t stateSpace() const {
+#ifdef DEBUG
+ assert( _vs.stateSpace() == _p.size() );
+#endif
+ return _p.size();
+ }
+
+ T operator[] (size_t i) const { return _p[i]; }
+ T& operator[] (size_t i) { return _p[i]; }
+ TFactor<T> & fill (T p)
+ { _p.fill( p ); return(*this); }
+ TFactor<T> & randomize ()
+ { _p.randomize(); return(*this); }
+ TFactor<T> operator* (T x) const {
+ Factor result = *this;
+ result.p() *= x;
+ return result;
+ }
+ TFactor<T>& operator*= (T x) {
+ _p *= x;
+ return *this;
+ }
+ TFactor<T> operator/ (T x) const {
+ Factor result = *this;
+ result.p() /= x;
+ return result;
+ }
+ TFactor<T>& operator/= (T x) {
+ _p /= x;
+ return *this;
+ }
+ TFactor<T> operator* (const TFactor<T>& Q) const;
+ TFactor<T>& operator*= (const TFactor<T>& Q) { return( *this = (*this * Q) ); }
+ TFactor<T> operator+ (const TFactor<T>& Q) const {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ TFactor<T> sum(*this);
+ sum._p += Q._p;
+ return sum;
+ }
+ TFactor<T> operator- (const TFactor<T>& Q) const {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ TFactor<T> sum(*this);
+ sum._p -= Q._p;
+ return sum;
+ }
+ TFactor<T>& operator+= (const TFactor<T>& Q) {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ _p += Q._p;
+ return *this;
+ }
+ TFactor<T>& operator-= (const TFactor<T>& Q) {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ _p -= Q._p;
+ return *this;
+ }
+
+ TFactor<T> operator^ (Real a) const { TFactor<T> x; x._vs = _vs; x._p = _p^a; return x; }
+ TFactor<T>& operator^= (Real a) { _p ^= a; return *this; }
+
+ TFactor<T>& makeZero( Real epsilon ) {
+ _p.makeZero( epsilon );
+ return *this;
+ }
+
+ TFactor<T> inverse() const {
+ TFactor<T> inv;
+ inv._vs = _vs;
+ inv._p = _p.inverse(true); // FIXME
+ return inv;
+ }
+
+ TFactor<T> divided_by( const TFactor<T>& denom ) const {
+#ifdef DEBUG
+ assert( denom._vs == _vs );
+#endif
+ TFactor<T> quot(*this);
+ quot._p /= denom._p;
+ return quot;
+ }
+
+ TFactor<T>& divide( const TFactor<T>& denom ) {
+#ifdef DEBUG
+ assert( denom._vs == _vs );
+#endif
+ _p /= denom._p;
+ return *this;
+ }
+
+ TFactor<T> exp() const {
+ TFactor<T> e;
+ e._vs = _vs;
+ e._p = _p.exp();
+ return e;
+ }
+
+ TFactor<T> log() const {
+ TFactor<T> l;
+ l._vs = _vs;
+ l._p = _p.log();
+ return l;
+ }
+
+ TFactor<T> log0() const {
+ TFactor<T> l0;
+ l0._vs = _vs;
+ l0._p = _p.log0();
+ return l0;
+ }
+
+ CFactor clog0() const {
+ CFactor l0;
+ l0._vs = _vs;
+ l0._p = _p.clog0();
+ return l0;
+ }
+
+ T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); }
+ TFactor<T> normalized( typename Prob::NormType norm ) const {
+ TFactor<T> result;
+ result._vs = _vs;
+ result._p = _p.normalized( norm );
+ return result;
+ }
+
+ // returns slice of this factor where the subset ns is in state ns_state
+ Factor slice( const VarSet & ns, size_t ns_state ) const {
+ assert( ns << _vs );
+ VarSet nsrem = _vs / ns;
+ Factor result( nsrem, 0.0 );
+
+ // OPTIMIZE ME
+ Index i_ns (ns, _vs);
+ Index i_nsrem (nsrem, _vs);
+ for( size_t i = 0; i < stateSpace(); i++, ++i_ns, ++i_nsrem )
+ if( (size_t)i_ns == ns_state )
+ result._p[i_nsrem] = _p[i];
+
+ return result;
+ }
+
+ // returns unnormalized marginal
+ TFactor<T> part_sum(const VarSet & ns) const;
+ // returns normalized marginal
+ TFactor<T> marginal(const VarSet & ns) const { return part_sum(ns).normalized( Prob::NORMPROB ); }
+
+ bool hasNaNs() const { return _p.hasNaNs(); }
+ bool hasNegatives() const { return _p.hasNegatives(); }
+ T totalSum() const { return _p.totalSum(); }
+ T maxAbs() const { return _p.maxAbs(); }
+ T max() const { return _p.max(); }
+ Complex entropy() const { return _p.entropy(); }
+ T strength( const Var &i, const Var &j ) const;
+
+ friend Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt ) {
+ if( x._vs.empty() || y._vs.empty() )
+ return -1;
+ else {
+#ifdef DEBUG
+ assert( x._vs == y._vs );
+#endif
+ return dist( x._p, y._p, dt );
+ }
+ }
+ friend Complex KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
+ template<class U> friend std::ostream& operator<< (std::ostream& os, const TFactor<U>& P);
+};
+
+
+template<typename T> TFactor<T> TFactor<T>::part_sum(const VarSet & ns) const {
+#ifdef DEBUG
+ assert( ns << _vs );
+#endif
+
+ TFactor<T> res( ns, 0.0 );
+
+ Index i_res( ns, _vs );
+ for( size_t i = 0; i < _p.size(); i++, ++i_res )
+ res._p[i_res] += _p[i];
+
+ return res;
+}
+
+
+template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P) {
+ os << "(" << P.vars() << " <";
+ for( size_t i = 0; i < P._p.size(); i++ )
+ os << P._p[i] << " ";
+ os << ">)";
+ return os;
+}
+
+
+template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& Q) const {
+ TFactor<T> prod( _vs | Q._vs, 0.0 );
+
+ Index i1(_vs, prod._vs);
+ Index i2(Q._vs, prod._vs);
+
+ for( size_t i = 0; i < prod._p.size(); i++, ++i1, ++i2 )
+ prod._p[i] += _p[i1] * Q._p[i2];
+
+ return prod;
+}
+
+
+template<typename T> Complex KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
+ if( P._vs.empty() || Q._vs.empty() )
+ return -1;
+ else {
+#ifdef DEBUG
+ assert( P._vs == Q._vs );
+#endif
+ return KL_dist( P._p, Q._p );
+ }
+}
+
+
+// calculate N(psi, i, j)
+template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
+#ifdef DEBUG
+ assert( _vs && i );
+ assert( _vs && j );
+ assert( i != j );
+#endif
+ VarSet ij = i | j;
+
+ T max = 0.0;
+ for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
+ for( size_t alpha2 = 0; alpha2 < i.states(); alpha2++ )
+ if( alpha2 != alpha1 )
+ for( size_t beta1 = 0; beta1 < j.states(); beta1++ )
+ for( size_t beta2 = 0; beta2 < j.states(); beta2++ )
+ if( beta2 != beta1 ) {
+ size_t as = 1, bs = 1;
+ if( i < j )
+ bs = i.states();
+ else
+ as = j.states();
+ T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).max();
+ T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).max();
+ T f = f1 * f2;
+ if( f > max )
+ max = f;
+ }
+
+ return std::tanh( 0.25 * std::log( max ) );
+}
+
+
+template<typename T> TFactor<T> RemoveFirstOrderInteractions( const TFactor<T> & psi ) {
+ TFactor<T> result = psi;
+
+ VarSet vars = psi.vars();
+ for( size_t iter = 0; iter < 100; iter++ ) {
+ for( VarSet::const_iterator n = vars.begin(); n != vars.end(); n++ )
+ result = result * result.part_sum(*n).inverse();
+ result.normalize( Prob::NORMPROB );
+ }
+
+ return result;
+}
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_factorgraph_h
+#define __defined_libdai_factorgraph_h
+
+
+#include <iostream>
+#include <map>
+#include <tr1/unordered_map>
+#include <dai/bipgraph.h>
+#include <dai/factor.h>
+
+
+namespace dai {
+
+
+bool hasShortLoops( const std::vector<Factor> &P );
+void RemoveShortLoops( std::vector<Factor> &P );
+
+
+class FactorGraph : public BipartiteGraph<Var,Factor> {
+ protected:
+ std::map<size_t,Prob> _undoProbs;
+ Prob::NormType _normtype;
+
+ public:
+ /// Default constructor
+ FactorGraph() : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {};
+ /// Copy constructor
+ FactorGraph(const FactorGraph & x) : BipartiteGraph<Var,Factor>(x), _undoProbs(), _normtype(x._normtype) {};
+ /// Construct FactorGraph from vector of Factors
+ FactorGraph(const std::vector<Factor> &P);
+ // Construct a FactorGraph from given factor and variable iterators
+ template<typename FactorInputIterator, typename VarInputIterator>
+ FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint = 0, size_t nr_var_hint = 0 );
+
+ /// Assignment operator
+ FactorGraph & operator=(const FactorGraph & x) {
+ if(this!=&x) {
+ BipartiteGraph<Var,Factor>::operator=(x);
+ _undoProbs = x._undoProbs;
+ _normtype = x._normtype;
+ }
+ return *this;
+ }
+ virtual ~FactorGraph() {}
+
+ // aliases
+ Var & var(size_t i) { return V1(i); }
+ const Var & var(size_t i) const { return V1(i); }
+ const std::vector<Var> & vars() const { return V1s(); }
+ std::vector<Var> & vars() { return V1s(); }
+ size_t nrVars() const { return V1s().size(); }
+ Factor & factor(size_t I) { return V2(I); }
+ const Factor & factor(size_t I) const { return V2(I); }
+ const std::vector<Factor> & factors() const { return V2s(); }
+ std::vector<Factor> & factors() { return V2s(); }
+ size_t nrFactors() const { return V2s().size(); }
+
+ /// Provides read access to neighbours of variable
+ const _nb_t & nbV( size_t i1 ) const { return nb1(i1); }
+ /// Provides full access to neighbours of variable
+ _nb_t & nbV( size_t i1 ) { return nb1(i1); }
+ /// Provides read access to neighbours of factor
+ const _nb_t & nbF( size_t i2 ) const { return nb2(i2); }
+ /// Provides full access to neighbours of factor
+ _nb_t & nbF( size_t i2 ) { return nb2(i2); }
+
+ size_t findVar(const Var & n) const {
+ size_t i = find( vars().begin(), vars().end(), n ) - vars().begin();
+ assert( i != nrVars() );
+ return i;
+ }
+ size_t findFactor(const VarSet &ns) const {
+ size_t I;
+ for( I = 0; I < nrFactors(); I++ )
+ if( factor(I).vars() == ns )
+ break;
+ assert( I != nrFactors() );
+ return I;
+ }
+
+ friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
+ friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
+
+ VarSet delta(const Var & n) const;
+ VarSet Delta(const Var & n) const;
+ virtual void makeFactorCavity(size_t I);
+ virtual void makeCavity(const Var & n);
+
+ long ReadFromFile(const char *filename);
+ long WriteToFile(const char *filename) const;
+ long WriteToDotFile(const char *filename) const;
+
+ Factor ExactMarginal(const VarSet & x) const;
+ Real ExactlogZ() const;
+
+ virtual void clamp( const Var & n, size_t i );
+
+ bool hasNegatives() const;
+ Prob::NormType NormType() const { return _normtype; }
+
+ std::vector<VarSet> Cliques() const;
+
+ virtual void undoProbs( const VarSet &ns );
+ void saveProbs( const VarSet &ns );
+ virtual void undoProb( size_t I );
+ void saveProb( size_t I );
+
+ bool isConnected() const;
+
+ virtual void updatedFactor( size_t /*I*/ ) {};
+
+ private:
+ /// Part of constructors (creates edges, neighbours and adjacency matrix)
+ void createGraph( size_t nrEdges );
+};
+
+
+// assumes that the set of variables in [var_begin,var_end) is the union of the variables in the factors in [fact_begin, fact_end)
+template<typename FactorInputIterator, typename VarInputIterator>
+FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {
+ // add factors
+ size_t nrEdges = 0;
+ V2s().reserve( nr_fact_hint );
+ for( FactorInputIterator p2 = fact_begin; p2 != fact_end; ++p2 ) {
+ V2s().push_back( *p2 );
+ nrEdges += p2->vars().size();
+ }
+
+ // add variables
+ V1s().reserve( nr_var_hint );
+ for( VarInputIterator p1 = var_begin; p1 != var_end; ++p1 )
+ V1s().push_back( *p1 );
+
+ // create graph structure
+ createGraph( nrEdges );
+}
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_hak_h
+#define __defined_libdai_hak_h
+
+
+#include <string>
+#include <dai/daialg.h>
+#include <dai/regiongraph.h>
+#include <dai/enum.h>
+
+
+namespace dai {
+
+
+/// HAK provides an implementation of the single and double-loop algorithms by Heskes, Albers and Kappen
+class HAK : public DAIAlgRG {
+ protected:
+ std::vector<Factor> _Qa;
+ std::vector<Factor> _Qb;
+ std::vector<Factor> _muab;
+ std::vector<Factor> _muba;
+
+ public:
+ /// Default constructor
+ HAK() : DAIAlgRG() {};
+
+ /// Copy constructor
+ HAK(const HAK & x) : DAIAlgRG(x), _Qa(x._Qa), _Qb(x._Qb), _muab(x._muab), _muba(x._muba) {};
+
+ /// Clone function
+ HAK* clone() const { return new HAK(*this); }
+
+ /// Construct from RegionGraph
+ HAK(const RegionGraph & rg, const Properties &opts);
+
+ /// Construct from RactorGraph using "clusters" option
+ HAK(const FactorGraph & fg, const Properties &opts);
+
+ /// Assignment operator
+ HAK & operator=(const HAK & x) {
+ if( this != &x ) {
+ DAIAlgRG::operator=(x);
+ _Qa = x._Qa;
+ _Qb = x._Qb;
+ _muab = x._muab;
+ _muba = x._muba;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+
+ ENUM3(ClustersType,MIN,DELTA,LOOP)
+ ClustersType Clusters() const { return GetPropertyAs<ClustersType>("clusters"); }
+ bool DoubleLoop() { return GetPropertyAs<bool>("doubleloop"); }
+ size_t LoopDepth() { return GetPropertyAs<size_t>("loopdepth"); }
+
+ Factor & muab( size_t alpha, size_t beta ) { return _muab[ORIR2E(alpha,beta)]; }
+ Factor & muba( size_t beta, size_t alpha ) { return _muba[ORIR2E(alpha,beta)]; }
+ const Factor& Qa( size_t alpha ) const { return _Qa[alpha]; };
+ const Factor& Qb( size_t beta ) const { return _Qb[beta]; };
+
+// void Regenerate();
+ double doGBP();
+ double doDoubleLoop();
+ double run();
+ void init();
+ std::string identify() const;
+ Factor belief( const Var &n ) const;
+ Factor belief( const VarSet &ns ) const;
+ std::vector<Factor> beliefs() const;
+ Complex logZ () const;
+
+ void init( const VarSet &ns );
+ void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
+ bool checkProperties();
+
+ private:
+ void constructMessages();
+ void findLoopClusters( const FactorGraph &fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars );
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_index_h
+#define __defined_libdai_index_h
+
+
+#include <vector>
+#include <dai/varset.h>
+
+
+namespace dai {
+
+
+/* Example:
+ *
+ * Index i ({s_j_1,s_j_2,...,s_j_m}, {s_1,...,s_N}); // j_k in {1,...,N}
+ * for( ; i>=0; ++i ) {
+ * // loops over all states of (s_1,...,s_N)
+ * // i is linear index of corresponding state of (s_j_1, ..., s_j_m)
+ * }
+ */
+
+
+class Index
+{
+private:
+ long _index;
+ std::vector<int> _count,_max,_sum;
+public:
+ Index () { _index=-1; };
+ Index (const VarSet& P, const VarSet& ns)
+ {
+ long sum=1;
+ VarSet::const_iterator j=ns.begin();
+ for(VarSet::const_iterator i=P.begin();i!=P.end();++i)
+ {
+ for(;j!=ns.end()&&j->label()<=i->label();++j)
+ {
+ _count.push_back(0);
+ _max.push_back(j->states());
+ _sum.push_back((i->label()==j->label())?sum:0);
+ };
+ sum*=i->states();
+ };
+ for(;j!=ns.end();++j)
+ {
+ _count.push_back(0);
+ _max.push_back(j->states());
+ _sum.push_back(0);
+ };
+ _index=0;
+ };
+ Index (const Index & ind) : _index(ind._index), _count(ind._count), _max(ind._max), _sum(ind._sum) {};
+ Index & operator=(const Index & ind) {
+ if(this!=&ind) {
+ _index = ind._index;
+ _count = ind._count;
+ _max = ind._max;
+ _sum = ind._sum;
+ }
+ return *this;
+ }
+ Index& clear ()
+ {
+ for(unsigned i=0;i!=_count.size();++i) _count[i]=0;
+ _index=0;
+ return(*this);
+ };
+ operator long () const { return(_index); };
+ Index& operator ++ ()
+ {
+ if(_index>=0)
+ {
+ unsigned i;
+ for(i=0;(i<_count.size())
+ &&(_index+=_sum[i],++_count[i]==_max[i]);++i)
+ {
+ _index-=_sum[i]*_max[i];
+ _count[i]=0;
+ };
+ if(i==_count.size()) _index=-1;
+ };
+ return(*this);
+ };
+};
+
+
+class multind {
+ private:
+ std::vector<size_t> _dims; // dimensions
+ std::vector<size_t> _pdims; // products of dimensions
+
+ public:
+ multind(const std::vector<size_t> di) {
+ _dims = di;
+ size_t prod = 1;
+ for( std::vector<size_t>::const_iterator i=di.begin(); i!=di.end(); i++ ) {
+ _pdims.push_back(prod);
+ prod = prod * (*i);
+ }
+ _pdims.push_back(prod);
+ }
+ multind(const VarSet& ns) {
+ _dims.reserve( ns.size() );
+ _pdims.reserve( ns.size() + 1 );
+ size_t prod = 1;
+ for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) {
+ _pdims.push_back( prod );
+ prod *= n->states();
+ _dims.push_back( n->states() );
+ }
+ _pdims.push_back( prod );
+ }
+ std::vector<size_t> vi(size_t li) const { // linear index to vector index
+ std::vector<size_t> v(_dims.size(),0);
+ assert(li < _pdims.back());
+ for( long j = v.size()-1; j >= 0; j-- ) {
+ size_t q = li / _pdims[j];
+ v[j] = q;
+ li = li - q * _pdims[j];
+ }
+ return v;
+ }
+ size_t li(const std::vector<size_t> vi) const { // linear index
+ size_t s = 0;
+ assert(vi.size() == _dims.size());
+ for( size_t j = 0; j < vi.size(); j++ )
+ s += vi[j] * _pdims[j];
+ return s;
+ }
+ size_t max() const { return( _pdims.back() ); };
+
+ // FIXME add an iterator, which increases a vector index just using addition
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_jtree_h
+#define __defined_libdai_jtree_h
+
+
+#include <vector>
+#include <string>
+#include <dai/daialg.h>
+#include <dai/varset.h>
+#include <dai/regiongraph.h>
+#include <dai/factorgraph.h>
+#include <dai/clustergraph.h>
+#include <dai/weightedgraph.h>
+#include <dai/enum.h>
+
+
+namespace dai {
+
+
+class JTree : public DAIAlgRG {
+ protected:
+ DEdgeVec _RTree; // rooted tree
+ std::vector<Factor> _Qa;
+ std::vector<Factor> _Qb;
+ std::vector<Factor> _mes;
+ double _logZ;
+
+
+ public:
+ ENUM2(UpdateType,HUGIN,SHSH)
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+
+ JTree() : DAIAlgRG(), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {};
+ JTree( const JTree& x ) : DAIAlgRG(x), _RTree(x._RTree), _Qa(x._Qa), _Qb(x._Qb), _mes(x._mes), _logZ(x._logZ) {};
+ JTree* clone() const { return new JTree(*this); }
+ JTree & operator=( const JTree& x ) {
+ if( this != &x ) {
+ DAIAlgRG::operator=(x);
+ _RTree = x._RTree;
+ _Qa = x._Qa;
+ _Qb = x._Qb;
+ _mes = x._mes;
+ _logZ = x._logZ;
+ }
+ return *this;
+ }
+ JTree( const FactorGraph &fg, const Properties &opts, bool automatic=true );
+ void GenerateJT( const std::vector<VarSet> &Cliques );
+
+ Factor & message(size_t i1, size_t i2) { return( _mes[ORIR2E(i1,i2)] ); }
+ const Factor & message(size_t i1, size_t i2) const { return( _mes[ORIR2E(i1,i2)] ); }
+
+ static const char *Name;
+ std::string identify() const;
+// void Regenerate();
+ void init() {
+ assert( checkProperties() );
+ }
+ void runHUGIN();
+ void runShaferShenoy();
+ double run();
+ Factor belief( const Var &n ) const;
+ Factor belief( const VarSet &ns ) const;
+ std::vector<Factor> beliefs() const;
+ Complex logZ() const;
+
+ void init( const VarSet &/*ns*/ ) {}
+ void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
+
+ size_t findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t PreviousRoot=(size_t)-1 ) const;
+ Factor calcMarginal( const VarSet& ns );
+ bool checkProperties();
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_lc_h
+#define __defined_libdai_lc_h
+
+
+#include <string>
+#include <dai/daialg.h>
+#include <dai/enum.h>
+#include <dai/factorgraph.h>
+
+
+namespace dai {
+
+
+class LC : public DAIAlgFG {
+ protected:
+ typedef struct { size_t i; size_t I; } _iI_type;
+
+ std::vector<Factor> _pancakes; // used by all LC types (psi_I is stored in the pancake)
+ std::vector<Factor> _cavitydists; // used by all LC types to store the approximate cavity distribution
+ /// _phis[VV2E(i,I)] corresponds to \f$ \phi^{\setminus i}_I(x_{I \setminus i}) \f$
+ std::vector<Factor> _phis;
+
+ /// Single variable beliefs
+ std::vector<Factor> _beliefs;
+
+ /// For each pair (i,j) with j in delta(i), store i and the common factor I
+ std::vector<_iI_type> _iI;
+
+ public:
+ ENUM6(CavityType,FULL,PAIR,PAIR2,PAIRINT,PAIRCUM,UNIFORM)
+ ENUM3(UpdateType,SEQFIX,SEQRND,NONE)
+
+ CavityType Cavity() const { return GetPropertyAs<CavityType>("cavity"); }
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+ bool reInit() const { return GetPropertyAs<bool>("reinit"); }
+
+ /// Default constructor
+ LC() : DAIAlgFG() {};
+ /// Copy constructor
+ LC(const LC & x) : DAIAlgFG(x), _pancakes(x._pancakes), _cavitydists(x._cavitydists), _phis(x._phis), _beliefs(x._beliefs), _iI(x._iI) {};
+ /// Clone function
+ LC* clone() const { return new LC(*this); }
+ /// Construct LC object from a FactorGraph and parameters
+ LC(const FactorGraph & fg, const Properties &opts);
+ /// Assignment operator
+ LC& operator=(const LC & x) {
+ if( this != &x ) {
+ DAIAlgFG::operator=(x);
+ _pancakes = x._pancakes;
+ _cavitydists = x._cavitydists;
+ _phis = x._phis;
+ _beliefs = x._beliefs;
+ _iI = x._iI;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+ double CalcCavityDist( size_t i, const std::string &name, const Properties &opts );
+ double InitCavityDists( const std::string &name, const Properties &opts );
+ long SetCavityDists( std::vector<Factor> &Q );
+
+ void init();
+ Factor NewPancake (size_t iI, bool & hasNaNs);
+ double run();
+
+ std::string identify() const;
+ Factor belief (const Var &n) const { return( _beliefs[findVar(n)] ); }
+ Factor belief (const VarSet &/*ns*/) const { assert( 0 == 1 ); }
+ std::vector<Factor> beliefs() const { return _beliefs; }
+ Complex logZ() const { return NAN; }
+ void CalcBelief (size_t i);
+ const Factor & belief (size_t i) const { return _beliefs[i]; };
+ const Factor & pancake (size_t i) const { return _pancakes[i]; };
+ const Factor & cavitydist (size_t i) const { return _cavitydists[i]; };
+ size_t nr_iI() const { return _iI.size(); };
+
+ void clamp( const Var &/*n*/, size_t /*i*/ ) { assert( 0 == 1 ); }
+ void undoProbs( const VarSet &/*ns*/ ) { assert( 0 == 1 ); }
+ void saveProbs( const VarSet &/*ns*/ ) { assert( 0 == 1 ); }
+ void makeFactorCavity(size_t /*I*/) { assert( 0 == 1 ); }
+ virtual void makeCavity(const Var & /*n*/) { assert( 0 == 1 ); }
+ bool checkProperties();
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_mf_h
+#define __defined_libdai_mf_h
+
+
+#include <string>
+#include <dai/daialg.h>
+#include <dai/factorgraph.h>
+
+
+namespace dai {
+
+
+class MF : public DAIAlgFG {
+ protected:
+ std::vector<Factor> _beliefs;
+
+ public:
+ // default constructor
+ MF() : DAIAlgFG(), _beliefs() {};
+ // copy constructor
+ MF(const MF & x) : DAIAlgFG(x), _beliefs(x._beliefs) {};
+ MF* clone() const { return new MF(*this); }
+ // construct MF object from FactorGraph
+ MF(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
+ assert( checkProperties() );
+ Regenerate();
+ }
+ // assignment operator
+ MF & operator=(const MF & x) {
+ if(this!=&x) {
+ DAIAlgFG::operator=(x);
+ _beliefs = x._beliefs;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+ std::string identify() const;
+ void Regenerate();
+ void init();
+ double run();
+ Factor belief1 (size_t i) const;
+ Factor belief (const Var &n) const;
+ Factor belief (const VarSet &ns) const;
+ std::vector<Factor> beliefs() const;
+ Complex logZ() const;
+
+ void init( const VarSet &ns );
+ void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
+ bool checkProperties();
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_mr_h
+#define __defined_libdai_mr_h
+
+
+#include <vector>
+#include <string>
+#include <dai/factorgraph.h>
+#include <dai/daialg.h>
+#include <dai/enum.h>
+
+
+namespace dai {
+
+
+class sub_nb;
+
+
+class MR : public DAIAlgFG {
+ private:
+ bool supported; // is the underlying factor graph supported?
+
+ std::vector<size_t> con; // con[i] = connectivity of spin i
+ std::vector<std::vector<size_t> > nb; // nb[i] are the neighbours of spin i
+ std::vector<std::vector<double> > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
+ std::vector<double> theta; // theta[i] is the local field on spin i
+ std::vector<std::vector<double> > M; // M[i][_j] is M^{(i)}_j
+ std::vector<std::vector<size_t> > kindex; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
+ std::vector<std::vector<std::vector<double> > > cors;
+
+ static const size_t kmax = 31;
+
+ size_t N;
+
+ std::vector<double> Mag;
+
+ public:
+ ENUM2(UpdateType,FULL,LINEAR)
+ ENUM3(InitType,RESPPROP,CLAMPING,EXACT)
+
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+ InitType Inits() const { return GetPropertyAs<InitType>("inits"); }
+
+ MR( const FactorGraph & fg, const Properties &opts );
+ void init(size_t _N, double *_w, double *_th);
+ void makekindex();
+ void read_files();
+ void init_cor();
+ double init_cor_resp();
+ void solvemcav();
+ void solveM();
+ double run();
+ Factor belief( const Var &n ) const;
+ Factor belief( const VarSet &/*ns*/ ) const { assert( 0 == 1 ); }
+ std::vector<Factor> beliefs() const;
+ Complex logZ() const { return NAN; }
+ void init() { assert( checkProperties() ); }
+ static const char *Name;
+ std::string identify() const;
+ double _tJ(size_t i, sub_nb A);
+
+ double Omega(size_t i, size_t _j, size_t _l);
+ double T(size_t i, sub_nb A);
+ double T(size_t i, size_t _j);
+ double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
+ double Gamma(size_t i, size_t _l1, size_t _l2);
+
+ double appM(size_t i, sub_nb A);
+ void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
+
+ double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
+ MR* clone() const { assert( 0 == 1 ); }
+
+ bool checkProperties();
+
+};
+
+
+// represents a subset of nb[i] as a binary number
+// the elements of a subset should be thought of as indices in nb[i]
+class sub_nb {
+ private:
+ size_t s;
+ size_t bits;
+
+ public:
+ // construct full subset containing nr_elmt elements
+ sub_nb(size_t nr_elmt) {
+#ifdef DEBUG
+ assert( nr_elmt < sizeof(size_t) / sizeof(char) * 8 );
+#endif
+ bits = nr_elmt;
+ s = (1U << bits) - 1;
+ }
+
+ // copy constructor
+ sub_nb( const sub_nb & x ) : s(x.s), bits(x.bits) {}
+
+ // assignment operator
+ sub_nb & operator=( const sub_nb &x ) {
+ if( this != &x ) {
+ s = x.s;
+ bits = x.bits;
+ }
+ return *this;
+ }
+
+ // returns number of elements
+ size_t size() {
+ size_t size = 0;
+ for(size_t bit = 0; bit < bits; bit++)
+ if( s & (1U << bit) )
+ size++;
+ return size;
+ }
+
+ // increases s by one (for enumeration in lexicographical order)
+ sub_nb operator++() {
+ s++;
+ if( s >= (1U << bits) )
+ s = 0;
+ return *this;
+ }
+
+ // return i'th element of this subset
+ size_t operator[](size_t i) {
+ size_t bit;
+ for(bit = 0; bit < bits; bit++ )
+ if( s & (1U << bit) ) {
+ if( i == 0 )
+ break;
+ else
+ i--;
+ }
+#ifdef DEBUG
+ assert( bit < bits );
+#endif
+ return bit;
+ }
+
+ // add index _j to this subset
+ sub_nb &operator +=(size_t _j) {
+ s |= (1U << _j);
+ return *this;
+ }
+
+ // return copy with index _j
+ sub_nb operator+(size_t _j) {
+ sub_nb x = *this;
+ x += _j;
+ return x;
+ }
+
+ // delete index _j from this subset
+ sub_nb &operator -=(size_t _j) {
+ s &= ~(1U << _j);
+ return *this;
+ }
+
+ // return copy without index _j
+ sub_nb operator-(size_t _j) {
+ sub_nb x = *this;
+ x -= _j;
+ return x;
+ }
+
+ // empty this subset
+ sub_nb & clear() {
+ s = 0;
+ return *this;
+ }
+
+ // returns true if subset is empty
+ bool empty() { return (s == 0); }
+
+ // return 1 if _j is contained, 0 otherwise ("is element of")
+ size_t operator&(size_t _j) { return s & (1U << _j); }
+
+ friend std::ostream& operator<< (std::ostream& os, const sub_nb x) {
+ if( x.bits == 0 )
+ os << "empty";
+ else {
+ for(size_t bit = x.bits; bit > 0; bit-- )
+ if( x.s & (1U << (bit-1)) )
+ os << "1";
+ else
+ os << "0";
+ }
+ return os;
+ }
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_prob_h
+#define __defined_libdai_prob_h
+
+
+#include <complex>
+#include <cmath>
+#include <vector>
+#include <iostream>
+#include <cassert>
+#include <dai/util.h>
+
+
+namespace dai {
+
+
+typedef double Real;
+typedef std::complex<double> Complex;
+
+template<typename T> class TProb;
+typedef TProb<Real> Prob;
+typedef TProb<Complex> CProb;
+
+
+/// TProb<T> implements a probability vector of type T.
+/// T should be castable from and to double and to complex.
+template <typename T> class TProb {
+ protected:
+ /// The entries
+ std::vector<T> _p;
+
+ private:
+ /// Calculate x times log(x), or 0 if x == 0
+ Complex xlogx( Real x ) const { return( x == 0.0 ? 0.0 : Complex(x) * std::log(Complex(x))); }
+
+ public:
+ /// NORMPROB means that the sum of all entries should be 1
+ /// NORMLINF means that the maximum absolute value of all entries should be 1
+ typedef enum { NORMPROB, NORMLINF } NormType;
+ /// DISTL1 is the L-1 distance (sum of absolute values of pointwise difference)
+ /// DISTLINF is the L-inf distance (maximum absolute value of pointwise difference)
+ /// DISTTV is the Total Variation distance
+ typedef enum { DISTL1, DISTLINF, DISTTV } DistType;
+
+ /// Default constructor
+ TProb() : _p() {}
+
+ /// Construct uniform distribution of given length
+ TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
+
+ /// Construct with given length and initial value
+ TProb( size_t n, Real p ) : _p(std::vector<T>(n,(T)p)) {}
+
+ /// Construct with given length and initial array
+ TProb( size_t n, const Real* p ) {
+ // Reserve-push_back is faster than resize-copy
+ _p.reserve( n );
+ for( size_t i = 0; i < n; i++ )
+ _p.push_back( p[i] );
+ }
+
+ /// Copy constructor
+ TProb( const TProb<T> & x ) : _p(x._p) {}
+
+ /// Assignment operator
+ TProb<T> & operator=( const TProb<T> &x ) {
+ if( this != &x ) {
+ _p = x._p;
+ }
+ return *this;
+ }
+
+ /// Provide read access to _p
+ const std::vector<T> & p() const { return _p; }
+
+ /// Provide full access to _p
+ std::vector<T> & p() { return _p; }
+
+ /// Provide read access to ith element of _p
+ T operator[]( size_t i ) const { return _p[i]; }
+
+ /// Provide full access to ith element of _p
+ T& operator[]( size_t i ) { return _p[i]; }
+
+ /// Set all elements to x
+ TProb<T> & fill(T x) {
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] = x;
+ return *this;
+ }
+
+ /// Set all elements to iid random numbers from uniform(0,1) distribution
+ TProb<T> & randomize() {
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] = rnd_uniform();
+ return *this;
+ }
+
+ /// Return size
+ size_t size() const {
+ return _p.size();
+ }
+
+ /// Make entries zero if (Real) absolute value smaller than epsilon
+ TProb<T>& makeZero (Real epsilon) {
+ for( size_t i = 0; i < size(); i++ )
+ if( fabs((Real)_p[i]) < epsilon )
+ _p[i] = 0;
+ return *this;
+ }
+
+ /// Multiplication with T x
+ TProb<T>& operator*= (T x) {
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] *= x;
+ return *this;
+ }
+
+ /// Return product of *this with T x
+ TProb<T> operator* (T x) const {
+ TProb<T> prod( *this );
+ prod *= x;
+ return prod;
+ }
+
+ /// Division by T x
+ TProb<T>& operator/= (T x) {
+#ifdef DEBUG
+ assert( x != 0.0 );
+#endif
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] /= x;
+ return *this;
+ }
+
+ /// Return quotient of *this and T x
+ TProb<T> operator/ (T x) const {
+ TProb<T> prod( *this );
+ prod /= x;
+ return prod;
+ }
+
+ /// Pointwise multiplication with q
+ TProb<T>& operator*= (const TProb<T> & q) {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] *= q[i];
+ return *this;
+ }
+
+ /// Return product of *this with q
+ TProb<T> operator* (const TProb<T> & q) const {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ TProb<T> prod( *this );
+ prod *= q;
+ return prod;
+ }
+
+ /// Pointwise addition with q
+ TProb<T>& operator+= (const TProb<T> & q) {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] += q[i];
+ return *this;
+ }
+
+ /// Pointwise subtraction of q
+ TProb<T>& operator-= (const TProb<T> & q) {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] -= q[i];
+ return *this;
+ }
+
+ /// Return sum of *this and q
+ TProb<T> operator+ (const TProb<T> & q) const {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ TProb<T> sum( *this );
+ sum += q;
+ return sum;
+ }
+
+ /// Return *this minus q
+ TProb<T> operator- (const TProb<T> & q) const {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ TProb<T> sum( *this );
+ sum -= q;
+ return sum;
+ }
+
+ /// Pointwise division by q
+ TProb<T>& operator/= (const TProb<T> & q) {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ for( size_t i = 0; i < size(); i++ ) {
+#ifdef DEBUG
+// assert( q[i] != 0.0 );
+#endif
+ if( q[i] == 0.0 ) // FIXME
+ _p[i] = 0.0;
+ else
+ _p[i] /= q[i];
+ }
+ return *this;
+ }
+
+ /// Pointwise division by q, division by zero yields infinity
+ TProb<T>& divide (const TProb<T> & q) {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] /= q[i];
+ return *this;
+ }
+
+ /// Return quotient of *this with q
+ TProb<T> operator/ (const TProb<T> & q) const {
+#ifdef DEBUG
+ assert( size() == q.size() );
+#endif
+ TProb<T> quot( *this );
+ quot /= q;
+ return quot;
+ }
+
+ /// Return pointwise inverse
+ TProb<T> inverse(bool zero = false) const {
+ TProb<T> inv;
+ inv._p.reserve( size() );
+ if( zero )
+ for( size_t i = 0; i < size(); i++ )
+ inv._p.push_back( _p[i] == 0.0 ? 0.0 : 1.0 / _p[i] );
+ else
+ for( size_t i = 0; i < size(); i++ ) {
+#ifdef DEBUG
+ assert( _p[i] != 0.0 );
+#endif
+ inv._p.push_back( 1.0 / _p[i] );
+ }
+ return inv;
+ }
+
+ /// Return *this to the power of a (pointwise)
+ TProb<T>& operator^= (Real a) {
+ if( a != 1.0 ) {
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] = std::pow( _p[i], a );
+ }
+ return *this;
+ }
+
+ /// Pointwise power of a
+ TProb<T> operator^ (Real a) const {
+ TProb<T> power;
+ if( a != 1.0 ) {
+ power._p.reserve( size() );
+ for( size_t i = 0; i < size(); i++ )
+ power._p.push_back( std::pow( _p[i], a ) );
+ } else
+ power = *this;
+ return power;
+ }
+
+ /// Pointwise exp
+ TProb<T> exp() const {
+ TProb<T> e;
+ e._p.reserve( size() );
+ for( size_t i = 0; i < size(); i++ )
+ e._p.push_back( std::exp( _p[i] ) );
+ return e;
+ }
+
+ /// Pointwise log
+ TProb<T> log() const {
+ TProb<T> l;
+ l._p.reserve( size() );
+ for( size_t i = 0; i < size(); i++ )
+ l._p.push_back( std::log( _p[i] ) );
+ return l;
+ }
+
+ /// Pointwise log (or 0 if == 0)
+ TProb<T> log0() const {
+ TProb<T> l0;
+ l0._p.reserve( size() );
+ for( size_t i = 0; i < size(); i++ )
+ l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( _p[i] ) );
+ return l0;
+ }
+
+ /// Pointwise (complex) log (or 0 if == 0)
+/* CProb clog0() const {
+ CProb l0;
+ l0._p.reserve( size() );
+ for( size_t i = 0; i < size(); i++ )
+ l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( Complex( _p[i] ) ) );
+ return l0;
+ }*/
+
+ /// Return distance of p and q
+ friend Real dist( const TProb<T> & p, const TProb<T> & q, DistType dt ) {
+#ifdef DEBUG
+ assert( p.size() == q.size() );
+#endif
+ Real result = 0.0;
+ switch( dt ) {
+ case DISTL1:
+ for( size_t i = 0; i < p.size(); i++ )
+ result += fabs((Real)p[i] - (Real)q[i]);
+ break;
+
+ case DISTLINF:
+ for( size_t i = 0; i < p.size(); i++ ) {
+ Real z = fabs((Real)p[i] - (Real)q[i]);
+ if( z > result )
+ result = z;
+ }
+ break;
+
+ case DISTTV:
+ for( size_t i = 0; i < p.size(); i++ )
+ result += fabs((Real)p[i] - (Real)q[i]);
+ result *= 0.5;
+ break;
+ }
+ return result;
+ }
+
+ /// Return (complex) Kullback-Leibler distance with q
+ friend Complex KL_dist( const TProb<T> & p, const TProb<T> & q ) {
+#ifdef DEBUG
+ assert( p.size() == q.size() );
+#endif
+ Complex result = 0.0;
+ for( size_t i = 0; i < p.size(); i++ ) {
+ if( (Real) p[i] != 0.0 ) {
+ Complex p_i = p[i];
+ Complex q_i = q[i];
+ result += p_i * (std::log(p_i) - std::log(q_i));
+ }
+ }
+ return result;
+ }
+
+ /// Return sum of all entries
+ T totalSum() const {
+ T Z = 0.0;
+ for( size_t i = 0; i < size(); i++ )
+ Z += _p[i];
+ return Z;
+ }
+
+ /// Converts entries to Real and returns maximum absolute value
+ T maxAbs() const {
+ T Z = 0.0;
+ for( size_t i = 0; i < size(); i++ ) {
+ Real mag = fabs( (Real) _p[i] );
+ if( mag > Z )
+ Z = mag;
+ }
+ return Z;
+ }
+
+ /// Returns maximum value
+ T max() const {
+ T Z = 0.0;
+ for( size_t i = 0; i < size(); i++ ) {
+ if( _p[i] > Z )
+ Z = _p[i];
+ }
+ return Z;
+ }
+
+ /// Normalize, using the specified norm
+ T normalize( NormType norm ) {
+ T Z = 0.0;
+ if( norm == NORMPROB )
+ Z = totalSum();
+ else if( norm == NORMLINF )
+ Z = maxAbs();
+#ifdef DEBUG
+ assert( Z != 0.0 );
+#endif
+ T Zi = 1.0 / Z;
+ for( size_t i = 0; i < size(); i++ )
+ _p[i] *= Zi;
+ return Z;
+ }
+
+ /// Return normalized copy of *this, using the specified norm
+ TProb<T> normalized( NormType norm ) const {
+ T Z = 0.0;
+ if( norm == NORMPROB )
+ Z = totalSum();
+ else if( norm == NORMLINF )
+ Z = maxAbs();
+#ifdef DEBUG
+ assert( Z != 0.0 );
+#endif
+ Z = 1.0 / Z;
+
+ TProb<T> result;
+ result._p.reserve( size() );
+ for( size_t i = 0; i < size(); i++ )
+ result._p.push_back( _p[i] * Z );
+ return result;
+ }
+
+ /// Returns true if one or more entries are NaN
+ bool hasNaNs() const {
+ bool NaNs = false;
+ for( size_t i = 0; i < size() && !NaNs; i++ )
+ if( isnan( _p[i] ) )
+ NaNs = true;
+ return NaNs;
+ }
+
+ /// Returns true if one or more entries are negative
+ bool hasNegatives() const {
+ bool Negatives = false;
+ for( size_t i = 0; i < size() && !Negatives; i++ )
+ if( _p[i] < 0.0 )
+ Negatives = true;
+ return Negatives;
+ }
+
+ /// Returns (complex) entropy
+ Complex entropy() const {
+ Complex S = 0.0;
+ for( size_t i = 0; i < size(); i++ )
+ S -= xlogx(_p[i]);
+ return S;
+ }
+
+ friend std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
+ for( size_t i = 0; i < P.size(); i++ )
+ os << P._p[i] << " ";
+ os << std::endl;
+ return os;
+ }
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_properties_h
+#define __defined_libdai_properties_h
+
+
+#include <iostream>
+#include <sstream>
+#include <boost/any.hpp>
+#include <map>
+#include <cassert>
+#include <typeinfo>
+
+
+namespace dai {
+
+
+typedef std::string PropertyKey;
+typedef boost::any PropertyValue;
+typedef std::pair<PropertyKey, PropertyValue> Property;
+
+
+/// Sends a Properties object to an output stream
+std::ostream& operator<< (std::ostream & os, const Property & p);
+
+
+/// The Properties class represents a set of properties
+class Properties : public std::map<PropertyKey, PropertyValue> {
+ public:
+ /// Gets a property
+ const PropertyValue & Get(const PropertyKey &key) const {
+ Properties::const_iterator x = find(key);
+#ifdef DEBUG
+ if( x == this->end() )
+ std::cerr << "Get cannot find property " << key << std::endl;
+#endif
+ assert( x != this->end() );
+ return x->second;
+ }
+
+ /// Sets a property
+ Properties & Set(const PropertyKey &key, const PropertyValue &val) { this->operator[](key) = val; return *this; }
+
+ /// Gets a property, casted as ValueType
+ template<typename ValueType>
+ ValueType GetAs(const PropertyKey &key) const {
+ try {
+ return boost::any_cast<ValueType>(Get(key));
+ } catch( const boost::bad_any_cast & ) {
+ std::cerr << "Cannot cast property " << key << " to ";
+ std::cerr << typeid(ValueType).name() << std::endl;
+ return boost::any_cast<ValueType>(Get(key));
+ }
+ }
+
+ /// Converts a property from string to ValueType, if necessary
+ template<typename ValueType>
+ void ConvertTo(const PropertyKey &key) {
+ PropertyValue val = Get(key);
+ if( val.type() != typeid(ValueType) ) {
+ assert( val.type() == typeid(std::string) );
+
+ std::stringstream ss;
+ ss << GetAs<std::string>(key);
+ ValueType result;
+ ss >> result;
+
+ Set(key, result);
+ }
+ }
+
+ /// Shorthand for (temporarily) adding properties, e.g. Properties p()("method","BP")("verbose",1)("tol",1e-9)
+ Properties operator()(const PropertyKey &key, const PropertyValue &val) const { Properties copy = *this; return copy.Set(key,val); }
+
+ /// Check if a property with given key exists
+ bool hasKey(const PropertyKey &key) const { Properties::const_iterator x = find(key); return (x != this->end()); }
+
+ /// Sends a Properties object to an output stream
+ friend std::ostream& operator<< (std::ostream & os, const Properties & ps);
+
+ /// Reads a Properties object from an input stream
+ friend std::istream& operator >> (std::istream& is, Properties & ps);
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_regiongraph_h
+#define __defined_libdai_regiongraph_h
+
+
+#include <iostream>
+#include <dai/bipgraph.h>
+#include <dai/factorgraph.h>
+#include <dai/weightedgraph.h>
+
+
+namespace dai {
+
+
+/// A Region is a set of variables with a counting number
+class Region : public VarSet {
+ protected:
+ /// Counting number
+ double _c;
+
+ public:
+ /// Default constructor
+ Region() : VarSet(), _c(1.0) {}
+
+ /// Construct Region from a VarSet and a counting number
+ Region(const VarSet & x, double c) : VarSet(x), _c(c) {}
+
+ /// Copy constructor
+ Region(const Region & x) : VarSet(x), _c(x._c) {}
+
+ /// Assignment operator
+ Region & operator=(const Region & x) {
+ if( this != &x ) {
+ VarSet::operator=(x);
+ _c = x._c;
+ }
+ return *this;
+ }
+
+ /// Provide read access to counting number
+ const double & c() const { return _c; }
+ /// Provide full access to counting number
+ double & c() { return _c; }
+};
+
+
+/// A FRegion is a factor with a counting number
+class FRegion : public Factor {
+ protected:
+ /// Counting number
+ double _c;
+
+ public:
+ /// Default constructor
+ FRegion() : Factor(), _c(1.0) {}
+
+ /// Constructs FRegion from a Factor and a counting number
+ FRegion( const Factor & x, double c ) : Factor(x), _c(c) {}
+
+ /// Copy constructor
+ FRegion( const FRegion & x ) : Factor(x), _c(x._c) {}
+
+ /// Assignment operator
+ FRegion & operator=(const FRegion & x) {
+ if( this != &x ) {
+ Factor::operator=(x);
+ _c = x._c;
+ }
+ return *this;
+ }
+
+ /// Provide read access to counting number
+ const double & c() const { return _c; }
+ /// Provide full access to counting number
+ double & c() { return _c; }
+};
+
+
+typedef BipartiteGraph<FRegion,Region> BipRegGraph;
+
+
+/// A RegionGraph is a bipartite graph consisting of outer regions (type FRegion) and inner regions (type Region)
+class RegionGraph : public FactorGraph, BipRegGraph {
+ public:
+ typedef BipRegGraph::_nb_t R_nb_t;
+ typedef R_nb_t::const_iterator R_nb_cit;
+ typedef BipRegGraph::_edge_t R_edge_t;
+
+
+ protected:
+ /// Give back the OR index that corresponds to a factor index
+ typedef std::map<size_t,size_t>::const_iterator fac2OR_cit;
+ std::map<size_t,size_t> _fac2OR;
+
+
+ public:
+ /// Default constructor
+ RegionGraph() : FactorGraph(), BipRegGraph(), _fac2OR() {}
+
+ /// Constructs a RegionGraph from a FactorGraph
+ RegionGraph(const FactorGraph & fg) : FactorGraph(fg), BipRegGraph(), _fac2OR() {}
+
+ /// Constructs a RegionGraph from a FactorGraph, a vector of outer regions, a vector of inner regions and a vector of edges
+ RegionGraph(const FactorGraph & fg, const std::vector<Region> & ors, const std::vector<Region> & irs, const std::vector<R_edge_t> & edges);
+
+ /// Constructs a RegionGraph from a FactorGraph and a vector of outer VarSets (CVM style)
+ RegionGraph(const FactorGraph & fg, const std::vector<VarSet> & cl);
+
+ /// Copy constructor
+ RegionGraph(const RegionGraph & x) : FactorGraph(x), BipRegGraph(x), _fac2OR(x._fac2OR) {}
+
+ /// Assignment operator
+ RegionGraph & operator=(const RegionGraph & x) {
+ if( this != &x ) {
+ FactorGraph::operator=(x);
+ BipRegGraph::operator=(x);
+ _fac2OR = x._fac2OR;
+ }
+ return *this;
+ }
+
+ /// Provides read access to outer region
+ const FRegion & OR(long alpha) const { return BipRegGraph::V1(alpha); }
+ /// Provides access to outer region
+ FRegion & OR(long alpha) { return BipRegGraph::V1(alpha); }
+ /// Provides read access to all outer regions
+ const std::vector<FRegion> & ORs() const { return BipRegGraph::V1s(); }
+ /// Provides access to all outer regions
+ std::vector<FRegion> &ORs() { return BipRegGraph::V1s(); }
+ /// Returns number of outer regions
+ size_t nr_ORs() const { return BipRegGraph::V1s().size(); }
+
+ /// Provides read access to inner region
+ const Region & IR(long beta) const { return BipRegGraph::V2(beta); }
+ /// Provides access to inner region
+ Region & IR(long beta) { return BipRegGraph::V2(beta); }
+ /// Provides read access to all inner regions
+ const std::vector<Region> & IRs() const { return BipRegGraph::V2s(); }
+ /// Provides access to all inner regions
+ std::vector<Region> & IRs() { return BipRegGraph::V2s(); }
+ /// Returns number of inner regions
+ size_t nr_IRs() const { return BipRegGraph::V2s().size(); }
+
+ /// Provides read access to edge
+ const R_edge_t & Redge(size_t ind) const { return BipRegGraph::edge(ind); }
+ /// Provides full access to edge
+ R_edge_t & Redge(size_t ind) { return BipRegGraph::edge(ind); }
+ /// Provides read access to all edges
+ const std::vector<R_edge_t> & Redges() const { return BipRegGraph::edges(); }
+ /// Provides full access to all edges
+ std::vector<R_edge_t> & Redges() { return BipRegGraph::edges(); }
+ /// Returns number of edges
+ size_t nr_Redges() const { return BipRegGraph::edges().size(); }
+
+ /// Provides read access to neighbours of outer region
+ const R_nb_t & nbOR( size_t i1 ) const { return BipRegGraph::nb1(i1); }
+ /// Provides full access to neighbours of outer region
+ R_nb_t & nbOR( size_t i1 ) { return BipRegGraph::nb1(i1); }
+
+ /// Provides read access to neighbours of inner region
+ const R_nb_t & nbIR( size_t i2 ) const { return BipRegGraph::nb2(i2); }
+ /// Provides full access to neighbours of inner region
+ R_nb_t & nbIR( size_t i2 ) { return BipRegGraph::nb2(i2); }
+
+ /// Converts the pair of outer/inner region indices (i1,i2) to the corresponding edge index
+ size_t ORIR2E( const size_t i1, const size_t i2 ) const { return BipRegGraph::VV2E(i1, i2); }
+
+ void Regenerate() { BipRegGraph::Regenerate(); }
+
+
+ /// Calculates counting numbers of inner regions based upon counting numbers of outer regions
+ void Calc_Counting_Numbers();
+ /// Check whether the counting numbers are valid
+ bool Check_Counting_Numbers();
+
+ /// Recompute all outer regions
+ void RecomputeORs();
+
+ /// Recompute all outer regions involving the variables in ns
+ void RecomputeORs( const VarSet & ns );
+
+ /// Recompute all outer regions involving factor I
+ void RecomputeOR( size_t I );
+
+ /// We have to overload FactorGraph::clamp because the corresponding outer regions have to be recomputed
+ void clamp( const Var &n, size_t i ) { FactorGraph::clamp( n, i ); RecomputeORs( n ); }
+
+ /// We have to overload FactorGraph::makeCavity because the corresponding outer regions have to be recomputed
+ void makeCavity( const Var &n ) { FactorGraph::makeCavity( n ); RecomputeORs( n ); }
+
+ /// We have to overload FactorGraph::makeFactorCavity because the corresponding outer regions have to be recomputed
+ void makeFactorCavity( size_t I ) { FactorGraph::makeFactorCavity( I ); RecomputeOR( I ); }
+
+ /// We have to overload FactorGraph::undoProbs because the corresponding outer regions have to be recomputed
+ void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs( ns ); RecomputeORs( ns ); }
+
+ /// We have to overload FactorGraph::undoProb because the corresponding outer regions have to be recomputed
+ void undoProb( size_t I ) { FactorGraph::undoProb( I ); RecomputeOR( I ); }
+
+ /// If updateFactor is called, we know that factor I has been changed and we should recompute the outer regions involving I
+ void updatedFactor( size_t I ) { RecomputeOR( I ); }
+
+ /// Send RegionGraph to output stream
+ friend std::ostream & operator << ( std::ostream & os, const RegionGraph & rg );
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_treeep_h
+#define __defined_libdai_treeep_h
+
+
+#include <vector>
+#include <string>
+#include <dai/daialg.h>
+#include <dai/varset.h>
+#include <dai/regiongraph.h>
+#include <dai/factorgraph.h>
+#include <dai/clustergraph.h>
+#include <dai/weightedgraph.h>
+#include <dai/jtree.h>
+#include <dai/enum.h>
+
+
+namespace dai {
+
+
+class TreeEPSubTree {
+ protected:
+ std::vector<Factor> _Qa;
+ std::vector<Factor> _Qb;
+ DEdgeVec _RTree;
+ std::vector<size_t> _a; // _Qa[alpha] <-> superTree._Qa[_a[alpha]]
+ std::vector<size_t> _b; // _Qb[beta] <-> superTree._Qb[_b[beta]]
+ // _Qb[beta] <-> _RTree[beta]
+ const Factor * _I;
+ VarSet _ns;
+ VarSet _nsrem;
+ double _logZ;
+
+
+ public:
+ TreeEPSubTree() : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(NULL), _ns(), _nsrem(), _logZ(0.0) {}
+ TreeEPSubTree( const TreeEPSubTree &x) : _Qa(x._Qa), _Qb(x._Qb), _RTree(x._RTree), _a(x._a), _b(x._b), _I(x._I), _ns(x._ns), _nsrem(x._nsrem), _logZ(x._logZ) {}
+ TreeEPSubTree & operator=( const TreeEPSubTree& x ) {
+ if( this != &x ) {
+ _Qa = x._Qa;
+ _Qb = x._Qb;
+ _RTree = x._RTree;
+ _a = x._a;
+ _b = x._b;
+ _I = x._I;
+ _ns = x._ns;
+ _nsrem = x._nsrem;
+ _logZ = x._logZ;
+ }
+ return *this;
+ }
+
+ TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I );
+ void init();
+ void InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb );
+ void HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb );
+ double logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const;
+ const Factor *& I() { return _I; }
+};
+
+
+class TreeEP : public JTree {
+ protected:
+ std::map<size_t, TreeEPSubTree> _Q;
+
+ public:
+ ENUM2(TypeType,ORG,ALT)
+ TypeType Type() const { return GetPropertyAs<TypeType>("type"); }
+
+ TreeEP() : JTree(), _Q() {};
+ TreeEP( const TreeEP& x ) : JTree(x), _Q(x._Q) {
+ for( size_t I = 0; I < nrFactors(); I++ )
+ if( offtree( I ) )
+ _Q[I].I() = &factor(I);
+ }
+ TreeEP* clone() const { return new TreeEP(*this); }
+ TreeEP & operator=( const TreeEP& x ) {
+ if( this != &x ) {
+ JTree::operator=(x);
+ _Q = x._Q;
+ for( size_t I = 0; I < nrFactors(); I++ )
+ if( offtree( I ) )
+ _Q[I].I() = &factor(I);
+ }
+ return *this;
+ }
+ TreeEP( const FactorGraph &fg, const Properties &opts );
+ void ConstructRG( const DEdgeVec &tree );
+
+ static const char *Name;
+ std::string identify() const;
+ void init();
+ double run();
+ Complex logZ() const;
+
+ bool offtree(size_t I) const { return !_fac2OR.count(I); }
+
+ void init( const VarSet &/*ns*/ ) { init(); }
+ void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
+ bool checkProperties();
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_util_h
+#define __defined_libdai_util_h
+
+
+#include <sys/times.h>
+#include <cstdio>
+
+
+namespace dai {
+
+
+clock_t toc();
+void rnd_seed( int seed );
+double rnd_uniform();
+double rnd_stdnormal();
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_var_h
+#define __defined_libdai_var_h
+
+
+#include <iostream>
+
+
+namespace dai {
+
+
+/// Represents a discrete variable
+class Var {
+ private:
+ /// Internal label of the variable
+ long _label;
+
+ /// Number of possible values
+ size_t _states;
+
+ public:
+ /// Default constructor
+ Var() : _label(-1), _states(0) {};
+ /// Constructor
+ Var(long label, size_t states) : _label(label), _states(states) {};
+
+ /// Read access to label
+ long label() const { return _label; };
+ /// Access to label
+ long & label() { return _label; };
+
+ /// Read access to states
+ size_t states () const { return _states; };
+ /// Access to states
+ size_t& states () { return _states; };
+
+ /// Smaller-than operator (compares labels)
+ bool operator < (const Var& n) const { return( _label < n._label ); };
+ /// Larger-than operator (compares labels)
+ bool operator > (const Var& n) const { return( _label > n._label ); };
+ /// Smaller-than-or-equal-to operator (compares labels)
+ bool operator <= (const Var& n) const { return( _label <= n._label ); };
+ /// Larger-than-or-equal-to operator (compares labels)
+ bool operator >= (const Var& n) const { return( _label >= n._label ); };
+ /// Not-equal-to operator (compares labels)
+ bool operator != (const Var& n) const { return( _label != n._label ); };
+ /// Equal-to operator (compares labels)
+ bool operator == (const Var& n) const { return( _label == n._label ); };
+
+ /// Stream output operator
+ friend std::ostream& operator << (std::ostream& os, const Var& n) {
+ return( os << "[" << n.label() << "]" );
+ };
+};
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_varset_h
+#define __defined_libdai_varset_h
+
+
+#include <set>
+#include <algorithm>
+#include <iostream>
+#include <cassert>
+#include <dai/var.h>
+
+
+namespace dai {
+
+
+/// VarSet represents a set of variables and is a descendant of set<Var>.
+/// In addition, it provides an easy interface for set-theoretic operations
+/// by operator overloading.
+class VarSet : private std::set<Var> {
+ protected:
+ /// Product of number of states of all contained variables
+ size_t _statespace;
+
+ /// Check whether ns is a subset
+ bool includes( const VarSet& ns ) const {
+ return std::includes( begin(), end(), ns.begin(), ns.end() );
+ }
+
+ /// Calculate statespace
+ size_t calcStateSpace() {
+ _statespace = 1;
+ for( const_iterator i = begin(); i != end(); ++i )
+ _statespace *= i->states();
+ return _statespace;
+ }
+
+
+ public:
+ /// Default constructor
+ VarSet() : _statespace(0) {};
+
+ /// Construct a VarSet with one variable
+ VarSet( const Var &n ) : _statespace( n.states() ) {
+ insert( n );
+ }
+
+ /// Construct a VarSet with two variables
+ VarSet( const Var &n1, const Var &n2 ) {
+ insert( n1 );
+ insert( n2 );
+ calcStateSpace();
+ }
+
+ /// Construct from a set<Var>
+ VarSet( const std::set<Var> &ns ) {
+ std::set<Var>::operator=( ns );
+ calcStateSpace();
+ }
+
+ /// Copy constructor
+ VarSet( const VarSet &x ) : std::set<Var>( x ), _statespace( x._statespace ) {}
+
+ /// Assignment operator
+ VarSet & operator=( const VarSet &x ) {
+ if( this != &x ) {
+ std::set<Var>::operator=( x );
+ _statespace = x._statespace;
+ }
+ return *this;
+ }
+
+
+ /// Return statespace, i.e. the product of the number of states of each variable
+ size_t stateSpace() const {
+#ifdef DEBUG
+ size_t x = 1;
+ for( const_iterator i = begin(); i != end(); ++i )
+ x *= i->states();
+ assert( x == _statespace );
+#endif
+ return _statespace;
+ }
+
+
+ /// Erase one variable
+ VarSet& operator/= (const Var& n) {
+ erase( n );
+ calcStateSpace();
+ return *this;
+ }
+
+ /// Add one variable
+ VarSet& operator|= (const Var& n) {
+ insert( n );
+ calcStateSpace();
+ return *this;
+ }
+
+ /// Setminus operator (result contains all variables except those in ns)
+ VarSet operator/ (const VarSet& ns) const {
+ VarSet res;
+ std::set_difference( begin(), end(), ns.begin(), ns.end(), inserter( res, res.begin() ) );
+ res.calcStateSpace();
+ return res;
+ }
+
+ /// Set-union operator (result contains all variables plus those in ns)
+ VarSet operator| (const VarSet& ns) const {
+ VarSet res;
+ std::set_union( begin(), end(), ns.begin(), ns.end(), inserter( res, res.begin() ) );
+ res.calcStateSpace();
+ return res;
+ }
+
+ /// Set-intersection operator (result contains all variables that are also contained in ns)
+ VarSet operator& (const VarSet& ns) const {
+ VarSet res;
+ std::set_intersection( begin(), end(), ns.begin(), ns.end(), inserter( res, res.begin() ) );
+ res.calcStateSpace();
+ return res;
+ }
+
+ /// Erases from *this all variables in ns
+ VarSet& operator/= (const VarSet& ns) {
+ return (*this = (*this / ns));
+ }
+
+ /// Adds to *this all variables in ns
+ VarSet& operator|= (const VarSet& ns) {
+ return (*this = (*this | ns));
+ }
+
+ /// Erases from *this all variables not in ns
+ VarSet& operator&= (const VarSet& ns) {
+ return (*this = (*this & ns));
+ }
+
+
+ /// Returns false if both *this and ns are empty
+ bool operator|| (const VarSet& ns) const {
+ return !( this->empty() && ns.empty() );
+ }
+
+ /// Returns true if *this and ns contain common variables
+ bool operator&& (const VarSet& ns) const {
+ return !( (*this & ns).empty() );
+ }
+
+ /// Returns true if *this is a subset of ns
+ bool operator<< (const VarSet& ns) const {
+ return ns.includes( *this );
+ }
+
+ /// Returns true if ns is a subset of *this
+ bool operator>> (const VarSet& ns) const {
+ return includes( ns );
+ }
+
+ /// Returns true if *this contains the variable n
+ bool operator&& (const Var& n) const {
+ return( find( n ) == end() ? false : true );
+ }
+
+
+ /// Sends a VarSet to an output stream
+ friend std::ostream& operator<< (std::ostream & os, const VarSet& ns) {
+ for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++)
+ os << *n;
+ return( os );
+ }
+
+
+/* The following makes part of the public interface of set<Var> available.
+ * It is important to note that insert functions have to be overloaded,
+ * because they have to recalculate the statespace. A different approach
+ * would be to publicly inherit from set<Var> and only overload the insert
+ * methods.
+ */
+
+ // Additional interface from set<Var> that has to be provided
+ using std::set<Var>::const_iterator;
+ using std::set<Var>::iterator;
+ using std::set<Var>::const_reference;
+ using std::set<Var>::begin;
+ using std::set<Var>::end;
+ using std::set<Var>::size;
+ using std::set<Var>::empty;
+
+ /// Copy of set<Var>::insert which additionally calculates the new statespace
+ std::pair<iterator, bool> insert( const Var& x ) {
+ std::pair<iterator, bool> result = std::set<Var>::insert( x );
+ calcStateSpace();
+ return result;
+ }
+
+ /// Copy of set<Var>::insert which additionally calculates the new statespace
+ iterator insert( iterator pos, const value_type& x ) {
+ iterator result = std::set<Var>::insert( pos, x );
+ calcStateSpace();
+ return result;
+ }
+
+ /// Test for equality (ignore _statespace member)
+ friend bool operator==( const VarSet &a, const VarSet &b ) {
+ return operator==( (std::set<Var>)a, (std::set<Var>)b );
+ }
+
+ /// Test for inequality (ignore _statespace member)
+ friend bool operator!=( const VarSet &a, const VarSet &b ) {
+ return operator!=( (std::set<Var>)a, (std::set<Var>)b );
+ }
+
+ friend bool operator<( const VarSet &a, const VarSet &b ) {
+ return operator<( (std::set<Var>)a, (std::set<Var>)b );
+ }
+};
+
+
+/// For two Vars n1 and n2, the expression n1 | n2 gives the Varset containing n1 and n2
+inline VarSet operator| (const Var& n1, const Var& n2) {
+ return( VarSet(n1, n2) );
+}
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_weightedgraph_h
+#define __defined_libdai_weightedgraph_h
+
+
+#include <vector>
+#include <map>
+#include <iostream>
+#include <set>
+
+
+namespace dai {
+
+
+/// Directed edge
+class DEdge {
+ public:
+ size_t n1, n2;
+
+ DEdge() {}
+ DEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
+ bool operator==( const DEdge &x ) const {
+ return ((n1 == x.n1) && (n2 == x.n2));
+ }
+ bool operator!=( const DEdge &x ) const {
+ return !(*this == x);
+ }
+ bool operator<( const DEdge &x ) const {
+ return( (n1 < x.n1) || ((n1 == x.n1) && (n2 < x.n2)) );
+ }
+ friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
+ os << "(" << e.n1 << "," << e.n2 << ")";
+ return os;
+ }
+};
+
+
+/// Undirected edge
+class UEdge {
+ public:
+ size_t n1, n2;
+
+ UEdge() {}
+ UEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
+ UEdge( const DEdge & e ) : n1(e.n1), n2(e.n2) {}
+ bool operator==( const UEdge &x ) {
+ return ((n1 == x.n1) && (n2 == x.n2)) || ((n1 == x.n2) && (n2 == x.n1));
+ }
+ bool operator<( const UEdge &x ) const {
+ size_t s = n1, l = n2;
+ if( s > l )
+ std::swap( s, l );
+ size_t xs = x.n1, xl = x.n2;
+ if( xs > xl )
+ std::swap( xs, xl );
+ return( (s < xs) || ((s == xs) && (l < xl)) );
+ }
+ friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
+ if( e.n1 < e.n2 )
+ os << "{" << e.n1 << "," << e.n2 << "}";
+ else
+ os << "{" << e.n2 << "," << e.n1 << "}";
+ return os;
+ }
+};
+
+
+typedef std::vector<UEdge> UEdgeVec;
+typedef std::vector<DEdge> DEdgeVec;
+std::ostream & operator << (std::ostream & os, const DEdgeVec & rt);
+template<class T> class WeightedGraph : public std::map<UEdge, T> {};
+typedef std::set<UEdge> Graph;
+
+
+/// Use Prim's algorithm to construct a maximal spanning tree from the weighted graph Graph
+template<typename T> DEdgeVec MaxSpanningTreePrim( const WeightedGraph<T> & Graph ) {
+ const long verbose = 0;
+
+ DEdgeVec result;
+ if( Graph.size() == 0 )
+ return result;
+ else {
+ // Make a copy
+ WeightedGraph<T> Gr = Graph;
+
+ // Nodes in the tree
+ std::set<size_t> treeV;
+
+ // Start with one node
+ treeV.insert( Gr.begin()->first.n1 );
+
+ // Perform Prim's algorithm
+ while( Gr.size() ) {
+ typename WeightedGraph<T>::iterator largest = Gr.end();
+
+ for( typename WeightedGraph<T>::iterator e = Gr.begin(); e != Gr.end(); ) {
+ if( verbose >= 1 )
+ std::cout << "considering edge " << e->first << "...";
+ bool e1_in_treeV = treeV.count( e->first.n1 );
+ bool e2_in_treeV = treeV.count( e->first.n2 );
+ if( e1_in_treeV && e2_in_treeV ) {
+ if( verbose >= 1 )
+ std::cout << "red";
+ Gr.erase( e++ ); // Nice trick!
+ } else if( e1_in_treeV || e2_in_treeV ) {
+ if( verbose >= 1 )
+ std::cout << e->second;
+ if( (largest == Gr.end()) || (e->second > largest->second) ) {
+ largest = e; // largest edge connected to the tree (until now)
+ if( verbose >= 1 )
+ std::cout << " and largest!";
+ }
+ e++;
+ } else {
+ if( verbose >= 1 )
+ std::cout << "out of reach";
+ e++;
+ }
+ if( verbose >= 1 )
+ std::cout << std::endl;
+ }
+
+ if( largest != Gr.end() ) {
+ if( verbose >= 1 )
+ std::cout << "largest = " << largest->first << std::endl;
+ // Add directed edge, pointing away from the root
+ if( treeV.count( largest->first.n1 ) ) {
+ result.push_back( DEdge( largest->first.n1, largest->first.n2 ) );
+ treeV.insert( largest->first.n2 );
+ } else {
+ result.push_back( DEdge( largest->first.n2, largest->first.n1 ) );
+ treeV.insert( largest->first.n1 );
+ }
+ Gr.erase( largest );
+ }
+ }
+
+ return result;
+ }
+}
+
+
+/// Calculate rooted tree from a tree T and a root
+DEdgeVec GrowRootedTree( const Graph & T, size_t Root );
+
+
+UEdgeVec RandomDRegularGraph( size_t N, size_t d );
+
+
+} // end of namespace dai
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __defined_libdai_x2x_h
+#define __defined_libdai_x2x_h
+
+
+#include <cmath>
+#include <cstring>
+
+
+namespace x2x {
+
+ // Probability tables store -1 first, then +1
+
+ /// Convert moments to cumulants upto order k
+ void m2c (int N, double *x, int k);
+
+ /// Convert cumulants to moments upto order k
+ void c2m (int N, double *x, int k);
+
+ /// Convert (generalized) weights to log probability or energy
+ void w2logp (int N, double *x);
+
+ /// Convert log probability or energy to (generalized) weights
+ void logp2w (int N, double *x);
+
+ /// Convert probability to moments
+ void p2m (int N, double *x);
+
+ /// Convert moments to probability
+ void m2p (int N, double *x);
+
+ /// Convert log probability to probability
+ void logp2p (int N, double *x);
+
+ /// Convert probability to log probability
+ void p2logp (int N, double *x);
+
+ /// Normalize a log probability table
+ void logpnorm (int N, double *x);
+
+ /// Normalize a probability table, use logpnorm whenever possible
+ void pnorm (int N, double *x);
+
+ /// Fills table with v for all entries with more than k indices
+ /// used for example when cumulants or moments are converted upto some order
+ void fill (int N, double *x, int k, double v);
+
+} // end of namespace x2x
+
+
+#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#ifndef __INDEX_H__
-#define __INDEX_H__
-
-
-#include <vector>
-#include "varset.h"
-
-
-namespace dai {
-
-
-/* Example:
- *
- * Index i ({s_j_1,s_j_2,...,s_j_m}, {s_1,...,s_N}); // j_k in {1,...,N}
- * for( ; i>=0; ++i ) {
- * // loops over all states of (s_1,...,s_N)
- * // i is linear index of corresponding state of (s_j_1, ..., s_j_m)
- * }
- */
-
-
-class Index
-{
-private:
- long _index;
- std::vector<int> _count,_max,_sum;
-public:
- Index () { _index=-1; };
- Index (const VarSet& P, const VarSet& ns)
- {
- long sum=1;
- VarSet::const_iterator j=ns.begin();
- for(VarSet::const_iterator i=P.begin();i!=P.end();++i)
- {
- for(;j!=ns.end()&&j->label()<=i->label();++j)
- {
- _count.push_back(0);
- _max.push_back(j->states());
- _sum.push_back((i->label()==j->label())?sum:0);
- };
- sum*=i->states();
- };
- for(;j!=ns.end();++j)
- {
- _count.push_back(0);
- _max.push_back(j->states());
- _sum.push_back(0);
- };
- _index=0;
- };
- Index (const Index & ind) : _index(ind._index), _count(ind._count), _max(ind._max), _sum(ind._sum) {};
- Index & operator=(const Index & ind) {
- if(this!=&ind) {
- _index = ind._index;
- _count = ind._count;
- _max = ind._max;
- _sum = ind._sum;
- }
- return *this;
- }
- Index& clear ()
- {
- for(unsigned i=0;i!=_count.size();++i) _count[i]=0;
- _index=0;
- return(*this);
- };
- operator long () const { return(_index); };
- Index& operator ++ ()
- {
- if(_index>=0)
- {
- unsigned i;
- for(i=0;(i<_count.size())
- &&(_index+=_sum[i],++_count[i]==_max[i]);++i)
- {
- _index-=_sum[i]*_max[i];
- _count[i]=0;
- };
- if(i==_count.size()) _index=-1;
- };
- return(*this);
- };
-};
-
-
-class multind {
- private:
- std::vector<size_t> _dims; // dimensions
- std::vector<size_t> _pdims; // products of dimensions
-
- public:
- multind(const std::vector<size_t> di) {
- _dims = di;
- size_t prod = 1;
- for( std::vector<size_t>::const_iterator i=di.begin(); i!=di.end(); i++ ) {
- _pdims.push_back(prod);
- prod = prod * (*i);
- }
- _pdims.push_back(prod);
- }
- multind(const VarSet& ns) {
- _dims.reserve( ns.size() );
- _pdims.reserve( ns.size() + 1 );
- size_t prod = 1;
- for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) {
- _pdims.push_back( prod );
- prod *= n->states();
- _dims.push_back( n->states() );
- }
- _pdims.push_back( prod );
- }
- std::vector<size_t> vi(size_t li) const { // linear index to vector index
- std::vector<size_t> v(_dims.size(),0);
- assert(li < _pdims.back());
- for( long j = v.size()-1; j >= 0; j-- ) {
- size_t q = li / _pdims[j];
- v[j] = q;
- li = li - q * _pdims[j];
- }
- return v;
- }
- size_t li(const std::vector<size_t> vi) const { // linear index
- size_t s = 0;
- assert(vi.size() == _dims.size());
- for( size_t j = 0; j < vi.size(); j++ )
- s += vi[j] * _pdims[j];
- return s;
- }
- size_t max() const { return( _pdims.back() ); };
-
- // FIXME add an iterator, which increases a vector index just using addition
-};
-
-
-}
-
-
-#endif
+++ /dev/null
-/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
- Radboud University Nijmegen, The Netherlands
-
- This file is part of libDAI.
-
- libDAI is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- libDAI is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with libDAI; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-*/
-
-
-#include <iostream>
-#include "jtree.h"
-
-
-namespace dai {
-
-
-using namespace std;
-
-
-const char *JTree::Name = "JTREE";
-
-
-bool JTree::checkProperties() {
- if (!HasProperty("verbose") )
- return false;
- if( !HasProperty("updates") )
- return false;
-
- ConvertPropertyTo<size_t>("verbose");
- ConvertPropertyTo<UpdateType>("updates");
-
- return true;
-}
-
-
-JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic) : DAIAlgRG(fg, opts), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {
- assert( checkProperties() );
-
- if( automatic ) {
- ClusterGraph _cg;
-
- // Copy factors
- for( size_t I = 0; I < nrFactors(); I++ )
- _cg.insert( factor(I).vars() );
- if( Verbose() >= 3 )
- cout << "Initial clusters: " << _cg << endl;
-
- // Retain only maximal clusters
- _cg.eraseNonMaximal();
- if( Verbose() >= 3 )
- cout << "Maximal clusters: " << _cg << endl;
-
- vector<VarSet> ElimVec = _cg.VarElim_MinFill().eraseNonMaximal().toVector();
- if( Verbose() >= 3 ) {
- cout << "VarElim_MinFill result: {" << endl;
- for( size_t i = 0; i < ElimVec.size(); i++ ) {
- if( i != 0 )
- cout << ", ";
- cout << ElimVec[i];
- }
- cout << "}" << endl;
- }
-
- GenerateJT( ElimVec );
- }
-}
-
-
-void JTree::GenerateJT( const vector<VarSet> &Cliques ) {
- // Construct a weighted graph (each edge is weighted with the cardinality
- // of the intersection of the nodes, where the nodes are the elements of
- // Cliques).
- WeightedGraph<int> JuncGraph;
- for( size_t i = 0; i < Cliques.size(); i++ )
- for( size_t j = i+1; j < Cliques.size(); j++ ) {
- size_t w = (Cliques[i] & Cliques[j]).size();
- JuncGraph[UEdge(i,j)] = w;
- }
-
- // Construct maximal spanning tree using Prim's algorithm
- _RTree = MaxSpanningTreePrim( JuncGraph );
-
- // Construct corresponding region graph
-
- // Create outer regions
- ORs().reserve( Cliques.size() );
- for( size_t i = 0; i < Cliques.size(); i++ )
- ORs().push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
-
- // For each factor, find an outer region that subsumes that factor.
- // Then, multiply the outer region with that factor.
- for( size_t I = 0; I < nrFactors(); I++ ) {
- size_t alpha;
- for( alpha = 0; alpha < nr_ORs(); alpha++ )
- if( OR(alpha).vars() >> factor(I).vars() ) {
-// OR(alpha) *= factor(I);
- _fac2OR[I] = alpha;
- break;
- }
- assert( alpha != nr_ORs() );
- }
- RecomputeORs();
-
- // Create inner regions and edges
- IRs().reserve( _RTree.size() );
- Redges().reserve( 2 * _RTree.size() );
- for( size_t i = 0; i < _RTree.size(); i++ ) {
- Redges().push_back( R_edge_t( _RTree[i].n1, IRs().size() ) );
- Redges().push_back( R_edge_t( _RTree[i].n2, IRs().size() ) );
- // inner clusters have counting number -1
- IRs().push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
- }
-
- // Regenerate BipartiteGraph internals
- Regenerate();
-
- // Create messages and beliefs
- _Qa.clear();
- _Qa.reserve( nr_ORs() );
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
- _Qa.push_back( OR(alpha) );
-
- _Qb.clear();
- _Qb.reserve( nr_IRs() );
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- _Qb.push_back( Factor( IR(beta), 1.0 ) );
-
- _mes.clear();
- _mes.reserve( nr_Redges() );
- for( size_t e = 0; e < nr_Redges(); e++ )
- _mes.push_back( Factor( IR(Redge(e).second), 1.0 ) );
-
- // Check counting numbers
- Check_Counting_Numbers();
-
- if( Verbose() >= 3 ) {
- cout << "Resulting regiongraph: " << *this << endl;
- }
-}
-
-
-string JTree::identify() const {
- stringstream result (stringstream::out);
- result << Name << GetProperties();
- return result.str();
-}
-
-
-Factor JTree::belief( const VarSet &ns ) const {
- vector<Factor>::const_iterator beta;
- for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
- if( beta->vars() >> ns )
- break;
- if( beta != _Qb.end() )
- return( beta->marginal(ns) );
- else {
- vector<Factor>::const_iterator alpha;
- for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
- if( alpha->vars() >> ns )
- break;
- assert( alpha != _Qa.end() );
- return( alpha->marginal(ns) );
- }
-}
-
-
-vector<Factor> JTree::beliefs() const {
- vector<Factor> result;
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- result.push_back( _Qb[beta] );
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
- result.push_back( _Qa[alpha] );
- return result;
-}
-
-
-Factor JTree::belief( const Var &n ) const {
- return belief( (VarSet)n );
-}
-
-
-// Needs no init
-void JTree::runHUGIN() {
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
- _Qa[alpha] = OR(alpha);
-
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- _Qb[beta].fill( 1.0 );
-
- // CollectEvidence
- _logZ = 0.0;
- for( size_t i = _RTree.size(); (i--) != 0; ) {
-// Make outer region _RTree[i].n1 consistent with outer region _RTree[i].n2
-// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
- Factor new_Qb = _Qa[_RTree[i].n2].part_sum( IR( i ) );
- _logZ += log(new_Qb.normalize( Prob::NORMPROB ));
- _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] );
- _Qb[i] = new_Qb;
- }
- if( _RTree.empty() )
- _logZ += log(_Qa[0].normalize( Prob::NORMPROB ) );
- else
- _logZ += log(_Qa[_RTree[0].n1].normalize( Prob::NORMPROB ));
-
- // DistributeEvidence
- for( size_t i = 0; i < _RTree.size(); i++ ) {
-// Make outer region _RTree[i].n2 consistent with outer region _RTree[i].n1
-// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
- Factor new_Qb = _Qa[_RTree[i].n1].marginal( IR( i ) );
- _Qa[_RTree[i].n2] *= new_Qb.divided_by( _Qb[i] );
- _Qb[i] = new_Qb;
- }
-
- // Normalize
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
- _Qa[alpha].normalize( Prob::NORMPROB );
-}
-
-
-// Really needs no init! Initial messages can be anything.
-void JTree::runShaferShenoy() {
- // First pass
- _logZ = 0.0;
- for( size_t e = _RTree.size(); (e--) != 0; ) {
- // send a message from _RTree[e].n2 to _RTree[e].n1
- // or, actually, from the seperator IR(e) to _RTree[e].n1
-
- size_t i = _RTree[e].n2;
- size_t j = _RTree[e].n1;
-
- Factor piet = OR(i);
- for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
- if( *k != e )
- piet *= message( i, *k );
- message( j, e ) = piet.part_sum( IR(e) );
- _logZ += log( message(j,e).normalize( Prob::NORMPROB ) );
- }
-
- // Second pass
- for( size_t e = 0; e < _RTree.size(); e++ ) {
- size_t i = _RTree[e].n1;
- size_t j = _RTree[e].n2;
-
- Factor piet = OR(i);
- for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
- if( *k != e )
- piet *= message( i, *k );
- message( j, e ) = piet.marginal( IR(e) );
- }
-
- // Calculate beliefs
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
- Factor piet = OR(alpha);
- for( R_nb_cit k = nbOR(alpha).begin(); k != nbOR(alpha).end(); k++ )
- piet *= message( alpha, *k );
- if( _RTree.empty() ) {
- _logZ += log( piet.normalize( Prob::NORMPROB ) );
- _Qa[alpha] = piet;
- } else if( alpha == _RTree[0].n1 ) {
- _logZ += log( piet.normalize( Prob::NORMPROB ) );
- _Qa[alpha] = piet;
- } else
- _Qa[alpha] = piet.normalized( Prob::NORMPROB );
- }
-
- // Only for logZ (and for belief)...
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- _Qb[beta] = _Qa[nbIR(beta)[0]].marginal( IR(beta) );
-}
-
-
-double JTree::run() {
- if( Updates() == UpdateType::HUGIN )
- runHUGIN();
- else if( Updates() == UpdateType::SHSH )
- runShaferShenoy();
- return 0.0;
-}
-
-
-Complex JTree::logZ() const {
- Complex sum = 0.0;
- for( size_t beta = 0; beta < nr_IRs(); beta++ )
- sum += Complex(IR(beta).c()) * _Qb[beta].entropy();
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
- sum += Complex(OR(alpha).c()) * _Qa[alpha].entropy();
- sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
- }
- return sum;
-}
-
-
-
-size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t PreviousRoot ) const {
- // find new root clique (the one with maximal statespace overlap with ns)
- size_t maxval = 0, maxalpha = 0;
- for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
- size_t val = (ns & OR(alpha).vars()).stateSpace();
- if( val > maxval ) {
- maxval = val;
- maxalpha = alpha;
- }
- }
-
-// for( size_t e = 0; e < _RTree.size(); e++ )
-// cout << OR(_RTree[e].n1).vars() << "->" << OR(_RTree[e].n2).vars() << ", ";
-// cout << endl;
- // grow new tree
- Graph oldTree;
- for( DEdgeVec::const_iterator e = _RTree.begin(); e != _RTree.end(); e++ )
- oldTree.insert( UEdge(e->n1, e->n2) );
- DEdgeVec newTree = GrowRootedTree( oldTree, maxalpha );
-// cout << ns << ": ";
-// for( size_t e = 0; e < newTree.size(); e++ )
-// cout << OR(newTree[e].n1).vars() << "->" << OR(newTree[e].n2).vars() << ", ";
-// cout << endl;
-
- // identify subtree that contains variables of ns which are not in the new root
- VarSet nsrem = ns / OR(maxalpha).vars();
-// cout << "nsrem:" << nsrem << endl;
- set<DEdge> subTree;
- // for each variable in ns that is not in the root clique
- for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ ) {
- // find first occurence of *n in the tree, which is closest to the root
- size_t e = 0;
- for( ; e != newTree.size(); e++ ) {
- if( OR(newTree[e].n2).vars() && *n )
- break;
- }
- assert( e != newTree.size() );
-
- // track-back path to root and add edges to subTree
- subTree.insert( newTree[e] );
- size_t pos = newTree[e].n1;
- for( ; e > 0; e-- )
- if( newTree[e-1].n2 == pos ) {
- subTree.insert( newTree[e-1] );
- pos = newTree[e-1].n1;
- }
- }
- if( PreviousRoot != (size_t)-1 && PreviousRoot != maxalpha) {
- // find first occurence of PreviousRoot in the tree, which is closest to the new root
- size_t e = 0;
- for( ; e != newTree.size(); e++ ) {
- if( newTree[e].n2 == PreviousRoot )
- break;
- }
- assert( e != newTree.size() );
-
- // track-back path to root and add edges to subTree
- subTree.insert( newTree[e] );
- size_t pos = newTree[e].n1;
- for( ; e > 0; e-- )
- if( newTree[e-1].n2 == pos ) {
- subTree.insert( newTree[e-1] );
- pos = newTree[e-1].n1;
- }
- }
-// cout << "subTree: " << endl;
-// for( set<DEdge>::const_iterator sTi = subTree.begin(); sTi != subTree.end(); sTi++ )
-// cout << OR(sTi->n1).vars() << "->" << OR(sTi->n2).vars() << ", ";
-// cout << endl;
-
- // Resulting Tree is a reordered copy of newTree
- // First add edges in subTree to Tree
- Tree.clear();
- for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
- if( subTree.count( *e ) ) {
- Tree.push_back( *e );
-// cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ", ";
- }
-// cout << endl;
- // Then add edges pointing away from nsrem
- // FIXME
-/* for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
- for( set<DEdge>::const_iterator sTi = subTree.begin(); sTi != subTree.end(); sTi++ )
- if( *e != *sTi ) {
- if( e->n1 == sTi->n1 || e->n1 == sTi->n2 ||
- e->n2 == sTi->n1 || e->n2 == sTi->n2 ) {
- Tree.push_back( *e );
-// cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ", ";
- }
- }*/
- // FIXME
-/* for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
- if( find( Tree.begin(), Tree.end(), *e) == Tree.end() ) {
- bool found = false;
- for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ )
- if( (OR(e->n1).vars() && *n) ) {
- found = true;
- break;
- }
- if( found ) {
- Tree.push_back( *e );
- cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ", ";
- }
- }
- cout << endl;*/
- size_t subTreeSize = Tree.size();
- // Then add remaining edges
- for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
- if( find( Tree.begin(), Tree.end(), *e ) == Tree.end() )
- Tree.push_back( *e );
-
- return subTreeSize;
-}
-
-
-// Cutset conditioning
-// assumes that run() has been called already
-Factor JTree::calcMarginal( const VarSet& ns ) {
- vector<Factor>::const_iterator beta;
- for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
- if( beta->vars() >> ns )
- break;
- if( beta != _Qb.end() )
- return( beta->marginal(ns) );
- else {
- vector<Factor>::const_iterator alpha;
- for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
- if( alpha->vars() >> ns )
- break;
- if( alpha != _Qa.end() )
- return( alpha->marginal(ns) );
- else {
- // Find subtree to do efficient inference
- DEdgeVec T;
- size_t Tsize = findEfficientTree( ns, T );
-
- // Find remaining variables (which are not in the new root)
- VarSet nsrem = ns / OR(T.front().n1).vars();
- Factor Pns (ns, 0.0);
-
- multind mi( nsrem );
-
- // Save _Qa and _Qb on the subtree
- map<size_t,Factor> _Qa_old;
- map<size_t,Factor> _Qb_old;
- vector<size_t> b(Tsize, 0);
- for( size_t i = Tsize; (i--) != 0; ) {
- size_t alpha1 = T[i].n1;
- size_t alpha2 = T[i].n2;
- size_t beta;
- for( beta = 0; beta < nr_IRs(); beta++ )
- if( UEdge( _RTree[beta].n1, _RTree[beta].n2 ) == UEdge( alpha1, alpha2 ) )
- break;
- assert( beta != nr_IRs() );
- b[i] = beta;
-
- if( !_Qa_old.count( alpha1 ) )
- _Qa_old[alpha1] = _Qa[alpha1];
- if( !_Qa_old.count( alpha2 ) )
- _Qa_old[alpha2] = _Qa[alpha2];
- if( !_Qb_old.count( beta ) )
- _Qb_old[beta] = _Qb[beta];
- }
-
- // For all states of nsrem
- for( size_t j = 0; j < mi.max(); j++ ) {
- vector<size_t> vi = mi.vi( j );
-
- // CollectEvidence
- double logZ = 0.0;
- for( size_t i = Tsize; (i--) != 0; ) {
- // Make outer region T[i].n1 consistent with outer region T[i].n2
- // IR(i) = seperator OR(T[i].n1) && OR(T[i].n2)
-
- size_t k = 0;
- for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++, k++ )
- if( _Qa[T[i].n2].vars() >> *n ) {
- Factor piet( *n, 0.0 );
- piet[vi[k]] = 1.0;
- _Qa[T[i].n2] *= piet;
- }
-
- Factor new_Qb = _Qa[T[i].n2].part_sum( IR( b[i] ) );
- logZ += log(new_Qb.normalize( Prob::NORMPROB ));
- _Qa[T[i].n1] *= new_Qb.divided_by( _Qb[b[i]] );
- _Qb[b[i]] = new_Qb;
- }
- logZ += log(_Qa[T[0].n1].normalize( Prob::NORMPROB ));
-
- Factor piet( nsrem, 0.0 );
- piet[j] = exp(logZ);
- Pns += piet * _Qa[T[0].n1].part_sum( ns / nsrem ); // OPTIMIZE ME
-
- // Restore clamped beliefs
- for( map<size_t,Factor>::const_iterator alpha = _Qa_old.begin(); alpha != _Qa_old.end(); alpha++ )
- _Qa[alpha->first] = alpha->second;
- for( map<size_t,Factor>::const_iterator beta = _Qb_old.begin(); beta != _Qb_old.end(); beta++ )
- _Qb[beta->first] = beta->second;
- }
-
- return( Pns.normalized(Prob::NORMPROB) );
- }
- }
-}
-
-
-}
projects / libdai.git / commitdiff