This repository has been archived on 2023-08-20. You can view files and clone it, but cannot push or open issues or pull requests.
yap-6.3/packages/CLPBN/horus/BeliefProp.cpp

470 lines
12 KiB
C++
Raw Normal View History

2012-05-23 14:56:01 +01:00
#include <cassert>
#include <limits>
#include <algorithm>
#include <iostream>
#include "BeliefProp.h"
2012-05-23 14:56:01 +01:00
#include "FactorGraph.h"
#include "Factor.h"
#include "Indexer.h"
#include "Horus.h"
BeliefProp::BeliefProp (const FactorGraph& fg) : Solver (fg)
2012-05-23 14:56:01 +01:00
{
runned_ = false;
}
BeliefProp::~BeliefProp (void)
2012-05-23 14:56:01 +01:00
{
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < varsI_.size(); i++) {
2012-05-23 14:56:01 +01:00
delete varsI_[i];
}
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < facsI_.size(); i++) {
2012-05-23 14:56:01 +01:00
delete facsI_[i];
}
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-23 14:56:01 +01:00
delete links_[i];
}
}
Params
BeliefProp::solveQuery (VarIds queryVids)
2012-05-23 14:56:01 +01:00
{
assert (queryVids.empty() == false);
2012-05-28 18:26:15 +01:00
return queryVids.size() == 1
? getPosterioriOf (queryVids[0])
: getJointDistributionOf (queryVids);
2012-05-23 14:56:01 +01:00
}
void
BeliefProp::printSolverFlags (void) const
2012-05-23 14:56:01 +01:00
{
stringstream ss;
ss << "belief propagation [" ;
ss << "schedule=" ;
typedef BpOptions::Schedule Sch;
switch (BpOptions::schedule) {
case Sch::SEQ_FIXED: ss << "seq_fixed"; break;
case Sch::SEQ_RANDOM: ss << "seq_random"; break;
case Sch::PARALLEL: ss << "parallel"; break;
case Sch::MAX_RESIDUAL: ss << "max_residual"; break;
}
2012-05-28 18:26:15 +01:00
ss << ",max_iter=" << Util::toString (BpOptions::maxIter);
ss << ",accuracy=" << Util::toString (BpOptions::accuracy);
2012-05-23 14:56:01 +01:00
ss << ",log_domain=" << Util::toString (Globals::logDomain);
ss << "]" ;
cout << ss.str() << endl;
}
Params
BeliefProp::getPosterioriOf (VarId vid)
2012-05-23 14:56:01 +01:00
{
if (runned_ == false) {
runSolver();
}
2012-05-31 21:24:40 +01:00
assert (fg.getVarNode (vid));
VarNode* var = fg.getVarNode (vid);
2012-05-23 14:56:01 +01:00
Params probs;
if (var->hasEvidence()) {
probs.resize (var->range(), LogAware::noEvidence());
probs[var->getEvidence()] = LogAware::withEvidence();
} else {
probs.resize (var->range(), LogAware::multIdenty());
const BpLinks& links = ninf(var)->getLinks();
2012-05-23 14:56:01 +01:00
if (Globals::logDomain) {
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links.size(); i++) {
2012-05-28 19:41:24 +01:00
probs += links[i]->message();
2012-05-23 14:56:01 +01:00
}
LogAware::normalize (probs);
2012-05-24 16:14:13 +01:00
Util::exp (probs);
2012-05-23 14:56:01 +01:00
} else {
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links.size(); i++) {
2012-05-28 19:41:24 +01:00
probs *= links[i]->message();
2012-05-23 14:56:01 +01:00
}
LogAware::normalize (probs);
}
}
return probs;
}
Params
BeliefProp::getJointDistributionOf (const VarIds& jointVarIds)
2012-05-23 14:56:01 +01:00
{
if (runned_ == false) {
runSolver();
}
2012-05-31 21:24:40 +01:00
VarNode* vn = fg.getVarNode (jointVarIds[0]);
2012-05-23 14:56:01 +01:00
const FacNodes& facNodes = vn->neighbors();
2012-05-24 22:55:20 +01:00
size_t idx = facNodes.size();
for (size_t i = 0; i < facNodes.size(); i++) {
2012-05-23 14:56:01 +01:00
if (facNodes[i]->factor().contains (jointVarIds)) {
idx = i;
break;
}
}
2012-05-24 22:55:20 +01:00
if (idx == facNodes.size()) {
2012-05-23 14:56:01 +01:00
return getJointByConditioning (jointVarIds);
}
2012-06-14 11:57:00 +01:00
return getFactorJoint (idx, jointVarIds);
}
Params
BeliefProp::getFactorJoint (
2012-06-14 11:57:00 +01:00
size_t fnIdx,
const VarIds& jointVarIds)
{
if (runned_ == false) {
runSolver();
}
2012-06-14 11:57:00 +01:00
FacNode* fn = fg.facNodes()[fnIdx];
Factor res (fn->factor());
const BpLinks& links = ninf(fn)->getLinks();
for (size_t i = 0; i < links.size(); i++) {
Factor msg ({links[i]->varNode()->varId()},
{links[i]->varNode()->range()},
getVarToFactorMsg (links[i]));
res.multiply (msg);
}
res.sumOutAllExcept (jointVarIds);
res.reorderArguments (jointVarIds);
res.normalize();
Params jointDist = res.params();
if (Globals::logDomain) {
Util::exp (jointDist);
}
return jointDist;
2012-05-23 14:56:01 +01:00
}
void
BeliefProp::runSolver (void)
2012-05-23 14:56:01 +01:00
{
initializeSolver();
nIters_ = 0;
while (!converged() && nIters_ < BpOptions::maxIter) {
nIters_ ++;
if (Globals::verbosity > 1) {
Util::printHeader (string ("Iteration ") + Util::toString (nIters_));
}
switch (BpOptions::schedule) {
case BpOptions::Schedule::SEQ_RANDOM:
2012-05-28 19:41:24 +01:00
std::random_shuffle (links_.begin(), links_.end());
2012-05-23 14:56:01 +01:00
// no break
case BpOptions::Schedule::SEQ_FIXED:
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-23 14:56:01 +01:00
calculateAndUpdateMessage (links_[i]);
}
break;
case BpOptions::Schedule::PARALLEL:
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-23 14:56:01 +01:00
calculateMessage (links_[i]);
}
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-23 14:56:01 +01:00
updateMessage(links_[i]);
}
break;
case BpOptions::Schedule::MAX_RESIDUAL:
maxResidualSchedule();
break;
}
}
if (Globals::verbosity > 0) {
if (nIters_ < BpOptions::maxIter) {
2012-06-14 12:19:51 +01:00
cout << "Belief propagation converged in " ;
2012-05-23 14:56:01 +01:00
cout << nIters_ << " iterations" << endl;
} else {
cout << "The maximum number of iterations was hit, terminating..." ;
cout << endl;
}
cout << endl;
}
runned_ = true;
}
void
BeliefProp::createLinks (void)
2012-05-23 14:56:01 +01:00
{
2012-05-31 21:24:40 +01:00
const FacNodes& facNodes = fg.facNodes();
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < facNodes.size(); i++) {
2012-05-23 14:56:01 +01:00
const VarNodes& neighbors = facNodes[i]->neighbors();
2012-05-24 22:55:20 +01:00
for (size_t j = 0; j < neighbors.size(); j++) {
links_.push_back (new BpLink (facNodes[i], neighbors[j]));
2012-05-23 14:56:01 +01:00
}
}
}
void
BeliefProp::maxResidualSchedule (void)
2012-05-23 14:56:01 +01:00
{
if (nIters_ == 1) {
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-23 14:56:01 +01:00
calculateMessage (links_[i]);
SortedOrder::iterator it = sortedOrder_.insert (links_[i]);
linkMap_.insert (make_pair (links_[i], it));
}
return;
}
2012-05-24 22:55:20 +01:00
for (size_t c = 0; c < links_.size(); c++) {
2012-05-23 14:56:01 +01:00
if (Globals::verbosity > 1) {
cout << "current residuals:" << endl;
for (SortedOrder::iterator it = sortedOrder_.begin();
2012-05-28 14:12:18 +01:00
it != sortedOrder_.end(); ++it) {
2012-05-23 14:56:01 +01:00
cout << " " << setw (30) << left << (*it)->toString();
2012-05-28 19:41:24 +01:00
cout << "residual = " << (*it)->residual() << endl;
2012-05-23 14:56:01 +01:00
}
}
SortedOrder::iterator it = sortedOrder_.begin();
BpLink* link = *it;
2012-05-28 19:41:24 +01:00
if (link->residual() < BpOptions::accuracy) {
2012-05-23 14:56:01 +01:00
return;
}
updateMessage (link);
link->clearResidual();
sortedOrder_.erase (it);
linkMap_.find (link)->second = sortedOrder_.insert (link);
// update the messages that depend on message source --> destin
2012-05-28 19:41:24 +01:00
const FacNodes& factorNeighbors = link->varNode()->neighbors();
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < factorNeighbors.size(); i++) {
2012-05-28 19:41:24 +01:00
if (factorNeighbors[i] != link->facNode()) {
const BpLinks& links = ninf(factorNeighbors[i])->getLinks();
2012-05-24 22:55:20 +01:00
for (size_t j = 0; j < links.size(); j++) {
2012-05-28 19:41:24 +01:00
if (links[j]->varNode() != link->varNode()) {
2012-05-23 14:56:01 +01:00
calculateMessage (links[j]);
BpLinkMap::iterator iter = linkMap_.find (links[j]);
2012-05-23 14:56:01 +01:00
sortedOrder_.erase (iter->second);
iter->second = sortedOrder_.insert (links[j]);
}
}
}
}
if (Globals::verbosity > 1) {
Util::printDashedLine();
}
}
}
void
BeliefProp::calcFactorToVarMsg (BpLink* link)
2012-05-23 14:56:01 +01:00
{
2012-05-28 19:41:24 +01:00
FacNode* src = link->facNode();
const VarNode* dst = link->varNode();
const BpLinks& links = ninf(src)->getLinks();
2012-05-23 14:56:01 +01:00
// calculate the product of messages that were sent
// to factor `src', except from var `dst'
2012-05-28 21:27:52 +01:00
unsigned reps = 1;
unsigned msgSize = Util::sizeExpected (src->factor().ranges());
2012-05-23 14:56:01 +01:00
Params msgProduct (msgSize, LogAware::multIdenty());
if (Globals::logDomain) {
2012-05-28 14:01:05 +01:00
for (size_t i = links.size(); i-- > 0; ) {
2012-05-28 19:41:24 +01:00
if (links[i]->varNode() != dst) {
2012-05-23 14:56:01 +01:00
if (Constants::SHOW_BP_CALCS) {
2012-05-28 19:41:24 +01:00
cout << " message from " << links[i]->varNode()->label();
2012-05-23 14:56:01 +01:00
cout << ": " ;
}
Util::apply_n_times (msgProduct, getVarToFactorMsg (links[i]),
reps, std::plus<double>());
2012-05-23 14:56:01 +01:00
if (Constants::SHOW_BP_CALCS) {
cout << endl;
}
}
2012-05-28 19:47:20 +01:00
reps *= links[i]->varNode()->range();
2012-05-23 14:56:01 +01:00
}
} else {
2012-05-28 14:01:05 +01:00
for (size_t i = links.size(); i-- > 0; ) {
2012-05-28 19:41:24 +01:00
if (links[i]->varNode() != dst) {
2012-05-23 14:56:01 +01:00
if (Constants::SHOW_BP_CALCS) {
2012-05-28 19:41:24 +01:00
cout << " message from " << links[i]->varNode()->label();
2012-05-23 14:56:01 +01:00
cout << ": " ;
}
Util::apply_n_times (msgProduct, getVarToFactorMsg (links[i]),
reps, std::multiplies<double>());
2012-05-23 14:56:01 +01:00
if (Constants::SHOW_BP_CALCS) {
cout << endl;
}
}
2012-05-28 19:47:20 +01:00
reps *= links[i]->varNode()->range();
2012-05-23 14:56:01 +01:00
}
}
Factor result (src->factor().arguments(),
src->factor().ranges(), msgProduct);
result.multiply (src->factor());
if (Constants::SHOW_BP_CALCS) {
cout << " message product: " << msgProduct << endl;
cout << " original factor: " << src->factor().params() << endl;
cout << " factor product: " << result.params() << endl;
}
result.sumOutAllExcept (dst->varId());
if (Constants::SHOW_BP_CALCS) {
cout << " marginalized: " << result.params() << endl;
}
2012-05-28 19:41:24 +01:00
link->nextMessage() = result.params();
LogAware::normalize (link->nextMessage());
2012-05-23 14:56:01 +01:00
if (Constants::SHOW_BP_CALCS) {
2012-05-28 19:41:24 +01:00
cout << " curr msg: " << link->message() << endl;
cout << " next msg: " << link->nextMessage() << endl;
2012-05-23 14:56:01 +01:00
}
}
Params
BeliefProp::getVarToFactorMsg (const BpLink* link) const
2012-05-23 14:56:01 +01:00
{
2012-05-28 19:41:24 +01:00
const VarNode* src = link->varNode();
2012-05-23 14:56:01 +01:00
Params msg;
if (src->hasEvidence()) {
msg.resize (src->range(), LogAware::noEvidence());
msg[src->getEvidence()] = LogAware::withEvidence();
} else {
msg.resize (src->range(), LogAware::one());
}
if (Constants::SHOW_BP_CALCS) {
cout << msg;
}
BpLinks::const_iterator it;
const BpLinks& links = ninf (src)->getLinks();
2012-05-23 14:56:01 +01:00
if (Globals::logDomain) {
2012-05-28 19:10:29 +01:00
for (it = links.begin(); it != links.end(); ++it) {
2012-05-28 19:41:24 +01:00
msg += (*it)->message();
2012-05-23 14:56:01 +01:00
if (Constants::SHOW_BP_CALCS) {
2012-05-28 19:41:24 +01:00
cout << " x " << (*it)->message();
2012-05-23 14:56:01 +01:00
}
}
2012-05-28 19:41:24 +01:00
msg -= link->message();
2012-05-23 14:56:01 +01:00
} else {
2012-05-28 19:10:29 +01:00
for (it = links.begin(); it != links.end(); ++it) {
2012-05-28 19:41:24 +01:00
msg *= (*it)->message();
2012-05-28 19:10:29 +01:00
if (Constants::SHOW_BP_CALCS) {
2012-05-28 19:41:24 +01:00
cout << " x " << (*it)->message();
2012-05-23 14:56:01 +01:00
}
}
2012-05-28 19:41:24 +01:00
msg /= link->message();
2012-05-23 14:56:01 +01:00
}
if (Constants::SHOW_BP_CALCS) {
cout << " = " << msg;
}
return msg;
}
Params
BeliefProp::getJointByConditioning (const VarIds& jointVarIds) const
2012-05-23 14:56:01 +01:00
{
return Solver::getJointByConditioning (GroundSolver::BP, fg, jointVarIds);
2012-05-23 14:56:01 +01:00
}
void
BeliefProp::initializeSolver (void)
2012-05-23 14:56:01 +01:00
{
2012-05-31 21:24:40 +01:00
const VarNodes& varNodes = fg.varNodes();
2012-05-23 14:56:01 +01:00
varsI_.reserve (varNodes.size());
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < varNodes.size(); i++) {
2012-05-23 14:56:01 +01:00
varsI_.push_back (new SPNodeInfo());
}
2012-05-31 21:24:40 +01:00
const FacNodes& facNodes = fg.facNodes();
2012-05-23 14:56:01 +01:00
facsI_.reserve (facNodes.size());
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < facNodes.size(); i++) {
2012-05-23 14:56:01 +01:00
facsI_.push_back (new SPNodeInfo());
}
createLinks();
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-28 19:41:24 +01:00
FacNode* src = links_[i]->facNode();
VarNode* dst = links_[i]->varNode();
ninf (dst)->addBpLink (links_[i]);
ninf (src)->addBpLink (links_[i]);
2012-05-23 14:56:01 +01:00
}
}
bool
BeliefProp::converged (void)
2012-05-23 14:56:01 +01:00
{
if (links_.size() == 0) {
return true;
}
if (nIters_ == 0) {
return false;
}
if (Globals::verbosity > 2) {
cout << endl;
}
if (nIters_ == 1) {
if (Globals::verbosity > 1) {
cout << "no residuals" << endl << endl;
}
return false;
}
bool converged = true;
if (BpOptions::schedule == BpOptions::Schedule::MAX_RESIDUAL) {
2012-05-28 19:41:24 +01:00
double maxResidual = (*(sortedOrder_.begin()))->residual();
2012-05-23 14:56:01 +01:00
if (maxResidual > BpOptions::accuracy) {
converged = false;
} else {
converged = true;
}
} else {
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
2012-05-28 19:41:24 +01:00
double residual = links_[i]->residual();
2012-05-23 14:56:01 +01:00
if (Globals::verbosity > 1) {
cout << links_[i]->toString() + " residual = " << residual << endl;
}
if (residual > BpOptions::accuracy) {
converged = false;
if (Globals::verbosity < 2) {
break;
}
}
}
if (Globals::verbosity > 1) {
cout << endl;
}
}
return converged;
}
void
BeliefProp::printLinkInformation (void) const
2012-05-23 14:56:01 +01:00
{
2012-05-24 22:55:20 +01:00
for (size_t i = 0; i < links_.size(); i++) {
BpLink* l = links_[i];
2012-05-23 14:56:01 +01:00
cout << l->toString() << ":" << endl;
cout << " curr msg = " ;
2012-05-28 19:41:24 +01:00
cout << l->message() << endl;
2012-05-23 14:56:01 +01:00
cout << " next msg = " ;
2012-05-28 19:41:24 +01:00
cout << l->nextMessage() << endl;
cout << " residual = " << l->residual() << endl;
2012-05-23 14:56:01 +01:00
}
}