diogo/yap-6.3
Archived
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changes by Tiago

Browse Source
This commit is contained in:
Vítor Santos Costa 2011-12-12 15:30:24 +00:00
parent 33bf0bc0f5
commit 7df0fd90f5
30 changed files with 498 additions and 8818 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

149
packages/CLPBN/clpbn/bp/BPNodeInfo.cpp
View File

@ -1,149 +0,0 @@
#include <cassert>
#include <cmath>
#include <iostream>
#include "BPNodeInfo.h"
#include "BPSolver.h"
BPNodeInfo::BPNodeInfo (BayesNode* node)
{
node_ = node;
ds_ = node->getDomainSize();
piValsCalc_ = false;
ldValsCalc_ = false;
nPiMsgsRcv_ = 0;
nLdMsgsRcv_ = 0;
piVals_.resize (ds_, 1);
ldVals_.resize (ds_, 1);
const BnNodeSet& childs = node->getChilds();
for (unsigned i = 0; i < childs.size(); i++) {
cmsgs_.insert (make_pair (childs[i], false));
}
const BnNodeSet& parents = node->getParents();
for (unsigned i = 0; i < parents.size(); i++) {
pmsgs_.insert (make_pair (parents[i], false));
}
}
ParamSet
BPNodeInfo::getBeliefs (void) const
{
double sum = 0.0;
ParamSet beliefs (ds_);
for (unsigned xi = 0; xi < ds_; xi++) {
double prod = piVals_[xi] * ldVals_[xi];
beliefs[xi] = prod;
sum += prod;
}
assert (sum);
//normalize the beliefs
for (unsigned xi = 0; xi < ds_; xi++) {
beliefs[xi] /= sum;
}
return beliefs;
}
bool
BPNodeInfo::readyToSendPiMsgTo (const BayesNode* child) const
{
for (unsigned i = 0; i < inChildLinks_.size(); i++) {
if (inChildLinks_[i]->getSource() != child
&& !inChildLinks_[i]->messageWasSended()) {
return false;
}
}
return true;
}
bool
BPNodeInfo::readyToSendLambdaMsgTo (const BayesNode* parent) const
{
for (unsigned i = 0; i < inParentLinks_.size(); i++) {
if (inParentLinks_[i]->getSource() != parent
&& !inParentLinks_[i]->messageWasSended()) {
return false;
}
}
return true;
}
double
BPNodeInfo::getPiValue (unsigned idx) const
{
assert (idx >=0 && idx < ds_);
return piVals_[idx];
}
void
BPNodeInfo::setPiValue (unsigned idx, Param value)
{
assert (idx >=0 && idx < ds_);
piVals_[idx] = value;
}
double
BPNodeInfo::getLambdaValue (unsigned idx) const
{
assert (idx >=0 && idx < ds_);
return ldVals_[idx];
}
void
BPNodeInfo::setLambdaValue (unsigned idx, Param value)
{
assert (idx >=0 && idx < ds_);
ldVals_[idx] = value;
}
double
BPNodeInfo::getBeliefChange (void)
{
double change = 0.0;
if (oldBeliefs_.size() == 0) {
oldBeliefs_ = getBeliefs();
change = 9999999999.0;
} else {
ParamSet currentBeliefs = getBeliefs();
for (unsigned xi = 0; xi < ds_; xi++) {
change += abs (currentBeliefs[xi] - oldBeliefs_[xi]);
}
oldBeliefs_ = currentBeliefs;
}
return change;
}
bool
BPNodeInfo::receivedBottomInfluence (void) const
{
// if all lambda values are equal, then neither
// this node neither its descendents have evidence,
// we can use this to don't send lambda messages his parents
bool childInfluenced = false;
for (unsigned xi = 1; xi < ds_; xi++) {
if (ldVals_[xi] != ldVals_[0]) {
childInfluenced = true;
break;
}
}
return childInfluenced;
}

82
packages/CLPBN/clpbn/bp/BPNodeInfo.h
View File

@ -1,82 +0,0 @@
#ifndef BP_BP_NODE_H
#define BP_BP_NODE_H
#include <vector>
#include <map>
#include "BPSolver.h"
#include "BayesNode.h"
#include "Shared.h"
//class Edge;
using namespace std;
class BPNodeInfo
{
public:
BPNodeInfo (int);
BPNodeInfo (BayesNode*);
ParamSet getBeliefs (void) const;
double getPiValue (unsigned) const;
void setPiValue (unsigned, Param);
double getLambdaValue (unsigned) const;
void setLambdaValue (unsigned, Param);
double getBeliefChange (void);
bool receivedBottomInfluence (void) const;
ParamSet& getPiValues (void) { return piVals_; }
ParamSet& getLambdaValues (void) { return ldVals_; }
bool arePiValuesCalculated (void) { return piValsCalc_; }
bool areLambdaValuesCalculated (void) { return ldValsCalc_; }
void markPiValuesAsCalculated (void) { piValsCalc_ = true; }
void markLambdaValuesAsCalculated (void) { ldValsCalc_ = true; }
void incNumPiMsgsRcv (void) { nPiMsgsRcv_ ++; }
void incNumLambdaMsgsRcv (void) { nLdMsgsRcv_ ++; }
bool receivedAllPiMessages (void)
{
return node_->getParents().size() == nPiMsgsRcv_;
}
bool receivedAllLambdaMessages (void)
{
return node_->getChilds().size() == nLdMsgsRcv_;
}
bool readyToSendPiMsgTo (const BayesNode*) const ;
bool readyToSendLambdaMsgTo (const BayesNode*) const;
CEdgeSet getIncomingParentLinks (void) { return inParentLinks_; }
CEdgeSet getIncomingChildLinks (void) { return inChildLinks_; }
CEdgeSet getOutcomingParentLinks (void) { return outParentLinks_; }
CEdgeSet getOutcomingChildLinks (void) { return outChildLinks_; }
void addIncomingParentLink (Edge* l) { inParentLinks_.push_back (l); }
void addIncomingChildLink (Edge* l) { inChildLinks_.push_back (l); }
void addOutcomingParentLink (Edge* l) { outParentLinks_.push_back (l); }
void addOutcomingChildLink (Edge* l) { outChildLinks_.push_back (l); }
private:
DISALLOW_COPY_AND_ASSIGN (BPNodeInfo);
ParamSet piVals_; // pi values
ParamSet ldVals_; // lambda values
ParamSet oldBeliefs_;
unsigned nPiMsgsRcv_;
unsigned nLdMsgsRcv_;
bool piValsCalc_;
bool ldValsCalc_;
EdgeSet inParentLinks_;
EdgeSet inChildLinks_;
EdgeSet outParentLinks_;
EdgeSet outChildLinks_;
unsigned ds_;
const BayesNode* node_;
map<const BayesNode*, bool> pmsgs_;
map<const BayesNode*, bool> cmsgs_;
};
#endif //BP_BP_NODE_H

905
packages/CLPBN/clpbn/bp/BPSolver.cpp
View File

@ -1,905 +0,0 @@
#include <cstdlib>
#include <limits>
#include <time.h>
#include <iostream>
#include <sstream>
#include <iomanip>
#include "BPSolver.h"
BPSolver::BPSolver (const BayesNet& bn) : Solver (&bn)
{
bn_ = &bn;
useAlwaysLoopySolver_ = false;
//jointCalcType_ = CHAIN_RULE;
jointCalcType_ = JUNCTION_NODE;
}
BPSolver::~BPSolver (void)
{
for (unsigned i = 0; i < nodesI_.size(); i++) {
delete nodesI_[i];
}
for (unsigned i = 0; i < links_.size(); i++) {
delete links_[i];
}
}
void
BPSolver::runSolver (void)
{
clock_t start_ = clock();
unsigned size = bn_->getNumberOfNodes();
unsigned nIters = 0;
initializeSolver();
if (bn_->isSingleConnected() && !useAlwaysLoopySolver_) {
runPolyTreeSolver();
Statistics::numSolvedPolyTrees ++;
} else {
runLoopySolver();
Statistics::numSolvedLoopyNets ++;
if (nIter_ >= SolverOptions::maxIter) {
Statistics::numUnconvergedRuns ++;
} else {
nIters = nIter_;
}
if (DL >= 2) {
cout << endl;
if (nIter_ < SolverOptions::maxIter) {
cout << "Belief propagation converged in " ;
cout << nIter_ << " iterations" << endl;
} else {
cout << "The maximum number of iterations was hit, terminating..." ;
cout << endl;
}
}
}
double time = (double (clock() - start_)) / CLOCKS_PER_SEC;
Statistics::updateStats (size, nIters, time);
if (EXPORT_TO_DOT && size > EXPORT_MIN_SIZE) {
stringstream ss;
ss << size << "." ;
ss << Statistics::getCounting (size) << ".dot" ;
bn_->exportToDotFormat (ss.str().c_str());
}
}
ParamSet
BPSolver::getPosterioriOf (Vid vid) const
{
BayesNode* node = bn_->getBayesNode (vid);
assert (node);
return nodesI_[node->getIndex()]->getBeliefs();
}
ParamSet
BPSolver::getJointDistributionOf (const VidSet& jointVids)
{
if (DL >= 2) {
cout << "calculating joint distribution on: " ;
for (unsigned i = 0; i < jointVids.size(); i++) {
Variable* var = bn_->getBayesNode (jointVids[i]);
cout << var->getLabel() << " " ;
}
cout << endl;
}
if (jointCalcType_ == JUNCTION_NODE) {
return getJointByJunctionNode (jointVids);
} else {
return getJointByChainRule (jointVids);
}
}
void
BPSolver::initializeSolver (void)
{
if (DL >= 2) {
if (!useAlwaysLoopySolver_) {
cout << "-> solver type = polytree solver" << endl;
cout << "-> schedule = n/a";
} else {
cout << "-> solver = loopy solver" << endl;
cout << "-> schedule = ";
switch (SolverOptions::schedule) {
case SolverOptions::S_SEQ_FIXED: cout << "sequential fixed" ; break;
case SolverOptions::S_SEQ_RANDOM: cout << "sequential random" ; break;
case SolverOptions::S_PARALLEL: cout << "parallel" ; break;
case SolverOptions::S_MAX_RESIDUAL: cout << "max residual" ; break;
}
}
cout << endl;
cout << "-> joint method = " ;
if (jointCalcType_ == JUNCTION_NODE) {
cout << "junction node" << endl;
} else {
cout << "chain rule " << endl;
}
cout << "-> max iters = " << SolverOptions::maxIter << endl;
cout << "-> accuracy = " << SolverOptions::accuracy << endl;
cout << endl;
}
CBnNodeSet nodes = bn_->getBayesNodes();
for (unsigned i = 0; i < nodesI_.size(); i++) {
delete nodesI_[i];
}
nodesI_.clear();
nodesI_.reserve (nodes.size());
links_.clear();
sortedOrder_.clear();
edgeMap_.clear();
for (unsigned i = 0; i < nodes.size(); i++) {
nodesI_.push_back (new BPNodeInfo (nodes[i]));
}
BnNodeSet roots = bn_->getRootNodes();
for (unsigned i = 0; i < roots.size(); i++) {
const ParamSet& params = roots[i]->getParameters();
ParamSet& piVals = M(roots[i])->getPiValues();
for (unsigned ri = 0; ri < roots[i]->getDomainSize(); ri++) {
piVals[ri] = params[ri];
}
}
for (unsigned i = 0; i < nodes.size(); i++) {
CBnNodeSet parents = nodes[i]->getParents();
for (unsigned j = 0; j < parents.size(); j++) {
Edge* newLink = new Edge (parents[j], nodes[i], PI_MSG);
links_.push_back (newLink);
M(nodes[i])->addIncomingParentLink (newLink);
M(parents[j])->addOutcomingChildLink (newLink);
}
CBnNodeSet childs = nodes[i]->getChilds();
for (unsigned j = 0; j < childs.size(); j++) {
Edge* newLink = new Edge (childs[j], nodes[i], LAMBDA_MSG);
links_.push_back (newLink);
M(nodes[i])->addIncomingChildLink (newLink);
M(childs[j])->addOutcomingParentLink (newLink);
}
}
for (unsigned i = 0; i < nodes.size(); i++) {
if (nodes[i]->hasEvidence()) {
ParamSet& piVals = M(nodes[i])->getPiValues();
ParamSet& ldVals = M(nodes[i])->getLambdaValues();
for (unsigned xi = 0; xi < nodes[i]->getDomainSize(); xi++) {
piVals[xi] = 0.0;
ldVals[xi] = 0.0;
}
piVals[nodes[i]->getEvidence()] = 1.0;
ldVals[nodes[i]->getEvidence()] = 1.0;
}
}
}
void
BPSolver::runPolyTreeSolver (void)
{
CBnNodeSet nodes = bn_->getBayesNodes();
for (unsigned i = 0; i < nodes.size(); i++) {
if (nodes[i]->isRoot()) {
M(nodes[i])->markPiValuesAsCalculated();
}
if (nodes[i]->isLeaf()) {
M(nodes[i])->markLambdaValuesAsCalculated();
}
}
bool finish = false;
while (!finish) {
finish = true;
for (unsigned i = 0; i < nodes.size(); i++) {
if (M(nodes[i])->arePiValuesCalculated() == false
&& M(nodes[i])->receivedAllPiMessages()) {
if (!nodes[i]->hasEvidence()) {
updatePiValues (nodes[i]);
}
M(nodes[i])->markPiValuesAsCalculated();
finish = false;
}
if (M(nodes[i])->areLambdaValuesCalculated() == false
&& M(nodes[i])->receivedAllLambdaMessages()) {
if (!nodes[i]->hasEvidence()) {
updateLambdaValues (nodes[i]);
}
M(nodes[i])->markLambdaValuesAsCalculated();
finish = false;
}
if (M(nodes[i])->arePiValuesCalculated()) {
CEdgeSet outChildLinks = M(nodes[i])->getOutcomingChildLinks();
for (unsigned j = 0; j < outChildLinks.size(); j++) {
BayesNode* child = outChildLinks[j]->getDestination();
if (!outChildLinks[j]->messageWasSended()) {
if (M(nodes[i])->readyToSendPiMsgTo (child)) {
outChildLinks[j]->setNextMessage (getMessage (outChildLinks[j]));
outChildLinks[j]->updateMessage();
M(child)->incNumPiMsgsRcv();
}
finish = false;
}
}
}
if (M(nodes[i])->areLambdaValuesCalculated()) {
CEdgeSet outParentLinks = M(nodes[i])->getOutcomingParentLinks();
for (unsigned j = 0; j < outParentLinks.size(); j++) {
BayesNode* parent = outParentLinks[j]->getDestination();
if (!outParentLinks[j]->messageWasSended()) {
if (M(nodes[i])->readyToSendLambdaMsgTo (parent)) {
outParentLinks[j]->setNextMessage (getMessage (outParentLinks[j]));
outParentLinks[j]->updateMessage();
M(parent)->incNumLambdaMsgsRcv();
}
finish = false;
}
}
}
}
}
}
void
BPSolver::runLoopySolver()
{
nIter_ = 0;
while (!converged() && nIter_ < SolverOptions::maxIter) {
nIter_++;
if (DL >= 2) {
cout << endl;
cout << "****************************************" ;
cout << "****************************************" ;
cout << endl;
cout << " Iteration " << nIter_ << endl;
cout << "****************************************" ;
cout << "****************************************" ;
cout << endl;
}
switch (SolverOptions::schedule) {
case SolverOptions::S_SEQ_RANDOM:
random_shuffle (links_.begin(), links_.end());
// no break
case SolverOptions::S_SEQ_FIXED:
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->setNextMessage (getMessage (links_[i]));
links_[i]->updateMessage();
updateValues (links_[i]);
}
break;
case SolverOptions::S_PARALLEL:
for (unsigned i = 0; i < links_.size(); i++) {
//calculateNextMessage (links_[i]);
}
for (unsigned i = 0; i < links_.size(); i++) {
//updateMessage (links_[i]);
//updateValues (links_[i]);
}
break;
case SolverOptions::S_MAX_RESIDUAL:
maxResidualSchedule();
break;
}
}
}
bool
BPSolver::converged (void) const
{
// this can happen if the graph is fully disconnected
if (links_.size() == 0) {
return true;
}
if (nIter_ == 0 || nIter_ == 1) {
return false;
}
bool converged = true;
if (SolverOptions::schedule == SolverOptions::S_MAX_RESIDUAL) {
Param maxResidual = (*(sortedOrder_.begin()))->getResidual();
if (maxResidual < SolverOptions::accuracy) {
converged = true;
} else {
converged = false;
}
} else {
CBnNodeSet nodes = bn_->getBayesNodes();
for (unsigned i = 0; i < nodes.size(); i++) {
if (!nodes[i]->hasEvidence()) {
double change = M(nodes[i])->getBeliefChange();
if (DL >= 2) {
cout << nodes[i]->getLabel() + " belief change = " ;
cout << change << endl;
}
if (change > SolverOptions::accuracy) {
converged = false;
if (DL == 0) break;
}
}
}
}
return converged;
}
void
BPSolver::maxResidualSchedule (void)
{
if (nIter_ == 1) {
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->setNextMessage (getMessage (links_[i]));
links_[i]->updateResidual();
SortedOrder::iterator it = sortedOrder_.insert (links_[i]);
edgeMap_.insert (make_pair (links_[i], it));
if (DL >= 2) {
cout << "calculating " << links_[i]->toString() << endl;
}
}
return;
}
for (unsigned c = 0; c < sortedOrder_.size(); c++) {
if (DL >= 2) {
cout << endl << "current residuals:" << endl;
for (SortedOrder::iterator it = sortedOrder_.begin();
it != sortedOrder_.end(); it ++) {
cout << " " << setw (30) << left << (*it)->toString();
cout << "residual = " << (*it)->getResidual() << endl;
}
}
SortedOrder::iterator it = sortedOrder_.begin();
Edge* edge = *it;
if (DL >= 2) {
cout << "updating " << edge->toString() << endl;
}
if (edge->getResidual() < SolverOptions::accuracy) {
return;
}
edge->updateMessage();
updateValues (edge);
edge->clearResidual();
sortedOrder_.erase (it);
edgeMap_.find (edge)->second = sortedOrder_.insert (edge);
// update the messages that depend on message source --> destin
CEdgeSet outChildLinks =
M(edge->getDestination())->getOutcomingChildLinks();
for (unsigned i = 0; i < outChildLinks.size(); i++) {
if (outChildLinks[i]->getDestination() != edge->getSource()) {
if (DL >= 2) {
cout << " calculating " << outChildLinks[i]->toString() << endl;
}
outChildLinks[i]->setNextMessage (getMessage (outChildLinks[i]));
outChildLinks[i]->updateResidual();
EdgeMap::iterator iter = edgeMap_.find (outChildLinks[i]);
sortedOrder_.erase (iter->second);
iter->second = sortedOrder_.insert (outChildLinks[i]);
}
}
CEdgeSet outParentLinks =
M(edge->getDestination())->getOutcomingParentLinks();
for (unsigned i = 0; i < outParentLinks.size(); i++) {
if (outParentLinks[i]->getDestination() != edge->getSource()) {
//&& !outParentLinks[i]->getDestination()->hasEvidence()) FIXME{
if (DL >= 2) {
cout << " calculating " << outParentLinks[i]->toString() << endl;
}
outParentLinks[i]->setNextMessage (getMessage (outParentLinks[i]));
outParentLinks[i]->updateResidual();
EdgeMap::iterator iter = edgeMap_.find (outParentLinks[i]);
sortedOrder_.erase (iter->second);
iter->second = sortedOrder_.insert (outParentLinks[i]);
}
}
}
}
void
BPSolver::updatePiValues (BayesNode* x)
{
// π(Xi)
if (DL >= 3) {
cout << "updating " << PI << " values for " << x->getLabel() << endl;
}
CEdgeSet parentLinks = M(x)->getIncomingParentLinks();
assert (x->getParents() == parentLinks.size());
const vector<CptEntry>& entries = x->getCptEntries();
stringstream* calcs1 = 0;
stringstream* calcs2 = 0;
ParamSet messageProducts (entries.size());
for (unsigned k = 0; k < entries.size(); k++) {
if (DL >= 5) {
calcs1 = new stringstream;
calcs2 = new stringstream;
}
double messageProduct = 1.0;
const DConf& conf = entries[k].getDomainConfiguration();
for (unsigned i = 0; i < parentLinks.size(); i++) {
assert (parentLinks[i]->getSource() == parents[i]);
assert (parentLinks[i]->getDestination() == x);
messageProduct *= parentLinks[i]->getMessage()[conf[i]];
if (DL >= 5) {
if (i != 0) *calcs1 << "." ;
if (i != 0) *calcs2 << "*" ;
*calcs1 << PI << "(" << parentLinks[i]->getSource()->getLabel();
*calcs1 << " --> " << x->getLabel() << ")" ;
*calcs1 << "[" ;
*calcs1 << parentLinks[i]->getSource()->getDomain()[conf[i]];
*calcs1 << "]";
*calcs2 << parentLinks[i]->getMessage()[conf[i]];
}
}
messageProducts[k] = messageProduct;
if (DL >= 5) {
cout << " mp" << k;
cout << " = " << (*calcs1).str();
if (parentLinks.size() == 1) {
cout << " = " << messageProduct << endl;
} else {
cout << " = " << (*calcs2).str();
cout << " = " << messageProduct << endl;
}
delete calcs1;
delete calcs2;
}
}
for (unsigned xi = 0; xi < x->getDomainSize(); xi++) {
double sum = 0.0;
if (DL >= 5) {
calcs1 = new stringstream;
calcs2 = new stringstream;
}
for (unsigned k = 0; k < entries.size(); k++) {
sum += x->getProbability (xi, entries[k]) * messageProducts[k];
if (DL >= 5) {
if (k != 0) *calcs1 << " + " ;
if (k != 0) *calcs2 << " + " ;
*calcs1 << x->cptEntryToString (xi, entries[k]);
*calcs1 << ".mp" << k;
*calcs2 << x->getProbability (xi, entries[k]);
*calcs2 << "*" << messageProducts[k];
}
}
M(x)->setPiValue (xi, sum);
if (DL >= 5) {
cout << " " << PI << "(" << x->getLabel() << ")" ;
cout << "[" << x->getDomain()[xi] << "]" ;
cout << " = " << (*calcs1).str();
cout << " = " << (*calcs2).str();
cout << " = " << sum << endl;
delete calcs1;
delete calcs2;
}
}
}
void
BPSolver::updateLambdaValues (BayesNode* x)
{
// λ(Xi)
if (DL >= 3) {
cout << "updating " << LD << " values for " << x->getLabel() << endl;
}
CEdgeSet childLinks = M(x)->getIncomingChildLinks();
stringstream* calcs1 = 0;
stringstream* calcs2 = 0;
for (unsigned xi = 0; xi < x->getDomainSize(); xi++) {
double product = 1.0;
if (DL >= 5) {
calcs1 = new stringstream;
calcs2 = new stringstream;
}
for (unsigned i = 0; i < childLinks.size(); i++) {
assert (childLinks[i]->getDestination() == x);
product *= childLinks[i]->getMessage()[xi];
if (DL >= 5) {
if (i != 0) *calcs1 << "." ;
if (i != 0) *calcs2 << "*" ;
*calcs1 << LD << "(" << childLinks[i]->getSource()->getLabel();
*calcs1 << "-->" << x->getLabel() << ")" ;
*calcs1 << "[" << x->getDomain()[xi] << "]" ;
*calcs2 << childLinks[i]->getMessage()[xi];
}
}
M(x)->setLambdaValue (xi, product);
if (DL >= 5) {
cout << " " << LD << "(" << x->getLabel() << ")" ;
cout << "[" << x->getDomain()[xi] << "]" ;
cout << " = " << (*calcs1).str();
if (childLinks.size() == 1) {
cout << " = " << product << endl;
} else {
cout << " = " << (*calcs2).str();
cout << " = " << product << endl;
}
delete calcs1;
delete calcs2;
}
}
}
ParamSet
BPSolver::calculateNextPiMessage (Edge* edge)
{
// πX(Zi)
BayesNode* z = edge->getSource();
BayesNode* x = edge->getDestination();
ParamSet zxPiNextMessage (z->getDomainSize());
CEdgeSet zChildLinks = M(z)->getIncomingChildLinks();
stringstream* calcs1 = 0;
stringstream* calcs2 = 0;
for (unsigned zi = 0; zi < z->getDomainSize(); zi++) {
double product = M(z)->getPiValue (zi);
if (DL >= 5) {
calcs1 = new stringstream;
calcs2 = new stringstream;
*calcs1 << PI << "(" << z->getLabel() << ")";
*calcs1 << "[" << z->getDomain()[zi] << "]" ;
*calcs2 << product;
}
for (unsigned i = 0; i < zChildLinks.size(); i++) {
assert (zChildLinks[i]->getDestination() == z);
if (zChildLinks[i]->getSource() != x) {
product *= zChildLinks[i]->getMessage()[zi];
if (DL >= 5) {
*calcs1 << "." << LD << "(" << zChildLinks[i]->getSource()->getLabel();
*calcs1 << "-->" << z->getLabel() << ")";
*calcs1 << "[" << z->getDomain()[zi] + "]" ;
*calcs2 << " * " << zChildLinks[i]->getMessage()[zi];
}
}
}
zxPiNextMessage[zi] = product;
if (DL >= 5) {
cout << " " << PI << "(" << z->getLabel();
cout << "-->" << x->getLabel() << ")" ;
cout << "[" << z->getDomain()[zi] << "]" ;
cout << " = " << (*calcs1).str();
if (zChildLinks.size() == 1) {
cout << " = " << product << endl;
} else {
cout << " = " << (*calcs2).str();
cout << " = " << product << endl;
}
delete calcs1;
delete calcs2;
}
}
return zxPiNextMessage;
}
ParamSet
BPSolver::calculateNextLambdaMessage (Edge* edge)
{
// λY(Xi)
BayesNode* y = edge->getSource();
BayesNode* x = edge->getDestination();
if (!M(y)->receivedBottomInfluence()) {
//cout << "returning 1" << endl;
//return edge->getMessage();
}
if (x->hasEvidence()) {
//cout << "returning 2" << endl;
//return edge->getMessage();
}
ParamSet yxLambdaNextMessage (x->getDomainSize());
CEdgeSet yParentLinks = M(y)->getIncomingParentLinks();
const vector<CptEntry>& allEntries = y->getCptEntries();
int parentIndex = y->getIndexOfParent (x);
stringstream* calcs1 = 0;
stringstream* calcs2 = 0;
vector<CptEntry> entries;
DConstraint constr = make_pair (parentIndex, 0);
for (unsigned i = 0; i < allEntries.size(); i++) {
if (allEntries[i].matchConstraints(constr)) {
entries.push_back (allEntries[i]);
}
}
ParamSet messageProducts (entries.size());
for (unsigned k = 0; k < entries.size(); k++) {
if (DL >= 5) {
calcs1 = new stringstream;
calcs2 = new stringstream;
}
double messageProduct = 1.0;
const DConf& conf = entries[k].getDomainConfiguration();
for (unsigned i = 0; i < yParentLinks.size(); i++) {
assert (yParentLinks[i]->getDestination() == y);
if (yParentLinks[i]->getSource() != x) {
if (DL >= 5) {
if (messageProduct != 1.0) *calcs1 << "*" ;
if (messageProduct != 1.0) *calcs2 << "*" ;
*calcs1 << PI << "(" << yParentLinks[i]->getSource()->getLabel();
*calcs1 << "-->" << y->getLabel() << ")" ;
*calcs1 << "[" ;
*calcs1 << yParentLinks[i]->getSource()->getDomain()[conf[i]];
*calcs1 << "]" ;
*calcs2 << yParentLinks[i]->getMessage()[conf[i]];
}
messageProduct *= yParentLinks[i]->getMessage()[conf[i]];
}
}
messageProducts[k] = messageProduct;
if (DL >= 5) {
cout << " mp" << k;
cout << " = " << (*calcs1).str();
if (yParentLinks.size() == 1) {
cout << 1 << endl;
} else if (yParentLinks.size() == 2) {
cout << " = " << messageProduct << endl;
} else {
cout << " = " << (*calcs2).str();
cout << " = " << messageProduct << endl;
}
delete calcs1;
delete calcs2;
}
}
for (unsigned xi = 0; xi < x->getDomainSize(); xi++) {
if (DL >= 5) {
calcs1 = new stringstream;
calcs2 = new stringstream;
}
vector<CptEntry> entries;
DConstraint constr = make_pair (parentIndex, xi);
for (unsigned i = 0; i < allEntries.size(); i++) {
if (allEntries[i].matchConstraints(constr)) {
entries.push_back (allEntries[i]);
}
}
double outerSum = 0.0;
for (unsigned yi = 0; yi < y->getDomainSize(); yi++) {
if (DL >= 5) {
(yi != 0) ? *calcs1 << " + {" : *calcs1 << "{" ;
(yi != 0) ? *calcs2 << " + {" : *calcs2 << "{" ;
}
double innerSum = 0.0;
for (unsigned k = 0; k < entries.size(); k++) {
if (DL >= 5) {
if (k != 0) *calcs1 << " + " ;
if (k != 0) *calcs2 << " + " ;
*calcs1 << y->cptEntryToString (yi, entries[k]);
*calcs1 << ".mp" << k;
*calcs2 << y->getProbability (yi, entries[k]);
*calcs2 << "*" << messageProducts[k];
}
innerSum += y->getProbability (yi, entries[k]) * messageProducts[k];
}
outerSum += innerSum * M(y)->getLambdaValue (yi);
if (DL >= 5) {
*calcs1 << "}." << LD << "(" << y->getLabel() << ")" ;
*calcs1 << "[" << y->getDomain()[yi] << "]";
*calcs2 << "}*" << M(y)->getLambdaValue (yi);
}
}
yxLambdaNextMessage[xi] = outerSum;
if (DL >= 5) {
cout << " " << LD << "(" << y->getLabel();
cout << "-->" << x->getLabel() << ")" ;
cout << "[" << x->getDomain()[xi] << "]" ;
cout << " = " << (*calcs1).str();
cout << " = " << (*calcs2).str();
cout << " = " << outerSum << endl;
delete calcs1;
delete calcs2;
}
}
return yxLambdaNextMessage;
}
ParamSet
BPSolver::getJointByJunctionNode (const VidSet& jointVids) const
{
BnNodeSet jointVars;
for (unsigned i = 0; i < jointVids.size(); i++) {
jointVars.push_back (bn_->getBayesNode (jointVids[i]));
}
BayesNet* mrn = bn_->getMinimalRequesiteNetwork (jointVids);
BnNodeSet parents;
unsigned dsize = 1;
for (unsigned i = 0; i < jointVars.size(); i++) {
parents.push_back (mrn->getBayesNode (jointVids[i]));
dsize *= jointVars[i]->getDomainSize();
}
unsigned nParams = dsize * dsize;
ParamSet params (nParams);
for (unsigned i = 0; i < nParams; i++) {
unsigned row = i / dsize;
unsigned col = i % dsize;
if (row == col) {
params[i] = 1;
} else {
params[i] = 0;
}
}
unsigned maxVid = std::numeric_limits<unsigned>::max();
Distribution* dist = new Distribution (params);
mrn->addNode (maxVid, dsize, NO_EVIDENCE, parents, dist);
mrn->setIndexes();
BPSolver solver (*mrn);
solver.runSolver();
const ParamSet& results = solver.getPosterioriOf (maxVid);
delete mrn;
delete dist;
return results;
}
ParamSet
BPSolver::getJointByChainRule (const VidSet& jointVids) const
{
BnNodeSet jointVars;
for (unsigned i = 0; i < jointVids.size(); i++) {
jointVars.push_back (bn_->getBayesNode (jointVids[i]));
}
BayesNet* mrn = bn_->getMinimalRequesiteNetwork (jointVids[0]);
BPSolver solver (*mrn);
solver.runSolver();
ParamSet prevBeliefs = solver.getPosterioriOf (jointVids[0]);
delete mrn;
VarSet observedVars = {jointVars[0]};
for (unsigned i = 1; i < jointVids.size(); i++) {
mrn = bn_->getMinimalRequesiteNetwork (jointVids[i]);
ParamSet newBeliefs;
vector<DConf> confs =
Util::getDomainConfigurations (observedVars);
for (unsigned j = 0; j < confs.size(); j++) {
for (unsigned k = 0; k < observedVars.size(); k++) {
if (!observedVars[k]->hasEvidence()) {
BayesNode* node = mrn->getBayesNode (observedVars[k]->getVarId());
if (node) {
node->setEvidence (confs[j][k]);
}
}
}
BPSolver solver (*mrn);
solver.runSolver();
ParamSet beliefs = solver.getPosterioriOf (jointVids[i]);
for (unsigned k = 0; k < beliefs.size(); k++) {
newBeliefs.push_back (beliefs[k]);
}
}
int count = -1;
for (unsigned j = 0; j < newBeliefs.size(); j++) {
if (j % jointVars[i]->getDomainSize() == 0) {
count ++;
}
newBeliefs[j] *= prevBeliefs[count];
}
prevBeliefs = newBeliefs;
observedVars.push_back (jointVars[i]);
delete mrn;
}
return prevBeliefs;
}
void
BPSolver::printMessageStatusOf (const BayesNode* var) const
{
cout << left;
cout << setw (10) << "domain" ;
cout << setw (20) << PI << "(" + var->getLabel() + ")" ;
cout << setw (20) << LD << "(" + var->getLabel() + ")" ;
cout << setw (16) << "belief" ;
cout << endl;
cout << "--------------------------------" ;
cout << "--------------------------------" ;
cout << endl;
BPNodeInfo* x = M(var);
ParamSet& piVals = x->getPiValues();
ParamSet& ldVals = x->getLambdaValues();
ParamSet beliefs = x->getBeliefs();
const Domain& domain = var->getDomain();
CBnNodeSet& childs = var->getChilds();
for (unsigned xi = 0; xi < var->getDomainSize(); xi++) {
cout << setw (10) << domain[xi];
cout << setw (19) << piVals[xi];
cout << setw (19) << ldVals[xi];
cout.precision (PRECISION);
cout << setw (16) << beliefs[xi];
cout << endl;
}
cout << endl;
if (childs.size() > 0) {
string s = "(" + var->getLabel() + ")" ;
for (unsigned j = 0; j < childs.size(); j++) {
cout << setw (10) << "domain" ;
cout << setw (28) << PI + childs[j]->getLabel() + s;
cout << setw (28) << LD + childs[j]->getLabel() + s;
cout << endl;
cout << "--------------------------------" ;
cout << "--------------------------------" ;
cout << endl;
/* FIXME
const ParamSet& piMessage = x->getPiMessage (childs[j]);
const ParamSet& lambdaMessage = x->getLambdaMessage (childs[j]);
for (unsigned xi = 0; xi < var->getDomainSize(); xi++) {
cout << setw (10) << domain[xi];
cout.precision (PRECISION);
cout << setw (27) << piMessage[xi];
cout.precision (PRECISION);
cout << setw (27) << lambdaMessage[xi];
cout << endl;
}
cout << endl;
*/
}
}
}
void
BPSolver::printAllMessageStatus (void) const
{
CBnNodeSet nodes = bn_->getBayesNodes();
for (unsigned i = 0; i < nodes.size(); i++) {
printMessageStatusOf (nodes[i]);
}
}

192
packages/CLPBN/clpbn/bp/BPSolver.h
View File

@ -1,192 +0,0 @@
#ifndef BP_BP_SOLVER_H
#define BP_BP_SOLVER_H
#include <vector>
#include <set>
#include "Solver.h"
#include "BayesNet.h"
#include "BPNodeInfo.h"
#include "Shared.h"
using namespace std;
class BPNodeInfo;
static const string PI = "pi" ;
static const string LD = "ld" ;
enum MessageType {PI_MSG, LAMBDA_MSG};
enum JointCalcType {CHAIN_RULE, JUNCTION_NODE};
class Edge
{
public:
Edge (BayesNode* s, BayesNode* d, MessageType t)
{
source_ = s;
destin_ = d;
type_ = t;
if (type_ == PI_MSG) {
currMsg_.resize (s->getDomainSize(), 1);
nextMsg_.resize (s->getDomainSize(), 1);
} else {
currMsg_.resize (d->getDomainSize(), 1);
nextMsg_.resize (d->getDomainSize(), 1);
}
msgSended_ = false;
residual_ = 0.0;
}
//void setMessage (ParamSet msg)
//{
// Util::normalize (msg);
// residual_ = Util::getMaxNorm (currMsg_, msg);
// currMsg_ = msg;
//}
void setNextMessage (CParamSet msg)
{
nextMsg_ = msg;
Util::normalize (nextMsg_);
residual_ = Util::getMaxNorm (currMsg_, nextMsg_);
}
void updateMessage (void)
{
currMsg_ = nextMsg_;
if (DL >= 3) {
cout << "updating " << toString() << endl;
}
msgSended_ = true;
}
void updateResidual (void)
{
residual_ = Util::getMaxNorm (currMsg_, nextMsg_);
}
string toString (void) const
{
stringstream ss;
if (type_ == PI_MSG) {
ss << PI;
} else if (type_ == LAMBDA_MSG) {
ss << LD;
} else {
abort();
}
ss << "(" << source_->getLabel();
ss << " --> " << destin_->getLabel() << ")" ;
return ss.str();
}
BayesNode* getSource (void) const { return source_; }
BayesNode* getDestination (void) const { return destin_; }
MessageType getMessageType (void) const { return type_; }
CParamSet getMessage (void) const { return currMsg_; }
bool messageWasSended (void) const { return msgSended_; }
double getResidual (void) const { return residual_; }
void clearResidual (void) { residual_ = 0.0; }
private:
BayesNode* source_;
BayesNode* destin_;
MessageType type_;
ParamSet currMsg_;
ParamSet nextMsg_;
bool msgSended_;
double residual_;
};
class BPSolver : public Solver
{
public:
BPSolver (const BayesNet&);
~BPSolver (void);
void runSolver (void);
ParamSet getPosterioriOf (Vid) const;
ParamSet getJointDistributionOf (const VidSet&);
private:
DISALLOW_COPY_AND_ASSIGN (BPSolver);
void initializeSolver (void);
void runPolyTreeSolver (void);
void runLoopySolver (void);
void maxResidualSchedule (void);
bool converged (void) const;
void updatePiValues (BayesNode*);
void updateLambdaValues (BayesNode*);
ParamSet calculateNextLambdaMessage (Edge* edge);
ParamSet calculateNextPiMessage (Edge* edge);
ParamSet getJointByJunctionNode (const VidSet&) const;
ParamSet getJointByChainRule (const VidSet&) const;
void printMessageStatusOf (const BayesNode*) const;
void printAllMessageStatus (void) const;
ParamSet getMessage (Edge* edge)
{
if (DL >= 3) {
cout << " calculating " << edge->toString() << endl;
}
if (edge->getMessageType() == PI_MSG) {
return calculateNextPiMessage (edge);
} else if (edge->getMessageType() == LAMBDA_MSG) {
return calculateNextLambdaMessage (edge);
} else {
abort();
}
return ParamSet();
}
void updateValues (Edge* edge)
{
if (!edge->getDestination()->hasEvidence()) {
if (edge->getMessageType() == PI_MSG) {
updatePiValues (edge->getDestination());
} else if (edge->getMessageType() == LAMBDA_MSG) {
updateLambdaValues (edge->getDestination());
} else {
abort();
}
}
}
BPNodeInfo* M (const BayesNode* node) const
{
assert (node);
assert (node == bn_->getBayesNode (node->getVarId()));
assert (node->getIndex() < nodesI_.size());
return nodesI_[node->getIndex()];
}
const BayesNet* bn_;
vector<BPNodeInfo*> nodesI_;
unsigned nIter_;
vector<Edge*> links_;
bool useAlwaysLoopySolver_;
JointCalcType jointCalcType_;
struct compare
{
inline bool operator() (const Edge* e1, const Edge* e2)
{
return e1->getResidual() > e2->getResidual();
}
};
typedef multiset<Edge*, compare> SortedOrder;
SortedOrder sortedOrder_;
typedef map<Edge*, SortedOrder::iterator> EdgeMap;
EdgeMap edgeMap_;
};
#endif //BP_BP_SOLVER_H

811
packages/CLPBN/clpbn/bp/BpNetwork.cpp
View File

@ -1,811 +0,0 @@
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <sstream>
#include "BpNetwork.h"
#include "BpNode.h"
#include "CptEntry.h"
BpNetwork::BpNetwork (void)
{
schedule_ = SEQUENTIAL_SCHEDULE;
maxIter_ = 150;
stableThreashold_ = 0.00000000000000000001;
}
BpNetwork::~BpNetwork (void)
{
for (unsigned int i = 0; i < nodes_.size(); i++) {
delete nodes_[i];
}
nodes_.clear();
}
void
BpNetwork::setSolverParameters (Schedule schedule,
int maxIter,
double stableThreashold)
{
if (maxIter <= 0) {
cerr << "error: maxIter must be greater or equal to 1" << endl;
abort();
}
if (stableThreashold <= 0.0 || stableThreashold >= 1.0) {
cerr << "error: stableThreashold must be greater than 0.0 " ;
cerr << "and lesser than 1.0" << endl;
abort();
}
schedule_ = schedule;
maxIter_ = maxIter;
stableThreashold_ = stableThreashold;
}
void
BpNetwork::runSolver (BayesianNode* queryVar)
{
vector<BayesianNode*> queryVars;
queryVars.push_back (queryVar);
runSolver (queryVars);
}
void
BpNetwork::runSolver (vector<BayesianNode*> queryVars)
{
if (queryVars.size() > 1) {
addJunctionNode (queryVars);
}
else {
string varName = queryVars[0]->getVariableName();
queryNode_ = static_cast<BpNode*> (getNode (varName));
}
if (!isPolyTree()) {
if (DL_ >= 1) {
cout << "The graph is not single connected. " ;
cout << "Iterative belief propagation will be used." ;
cout << endl << endl;
}
schedule_ = PARALLEL_SCHEDULE;
}
if (schedule_ == SEQUENTIAL_SCHEDULE) {
initializeSolver (queryVars);
runNeapolitanSolver();
for (unsigned int i = 0; i < nodes_.size(); i++) {
if (nodes_[i]->hasEvidence()) {
BpNode* v = static_cast<BpNode*> (nodes_[i]);
addEvidence (v);
vector<BpNode*> parents = cast (v->getParents());
for (unsigned int i = 0; i < parents.size(); i++) {
if (!parents[i]->hasEvidence()) {
sendLambdaMessage (v, parents[i]);
}
}
vector<BpNode*> childs = cast (v->getChilds());
for (unsigned int i = 0; i < childs.size(); i++) {
sendPiMessage (v, childs[i]);
}
}
}
} else if (schedule_ == PARALLEL_SCHEDULE) {
BpNode::enableParallelSchedule();
initializeSolver (queryVars);
for (unsigned int i = 0; i < nodes_.size(); i++) {
if (nodes_[i]->hasEvidence()) {
addEvidence (static_cast<BpNode*> (nodes_[i]));
}
}
runIterativeBpSolver();
}
}
void
BpNetwork::printCurrentStatus (void)
{
for (unsigned int i = 0; i < nodes_.size(); i++) {
printCurrentStatusOf (static_cast<BpNode*> (nodes_[i]));
}
}
void
BpNetwork::printCurrentStatusOf (BpNode* x)
{
vector<BpNode*> childs = cast (x->getChilds());
vector<string> domain = x->getDomain();
cout << left;
cout << setw (10) << "domain" ;
cout << setw (20) << "π(" + x->getVariableName() + ")" ;
cout << setw (20) << "λ(" + x->getVariableName() + ")" ;
cout << setw (16) << "belief" ;
cout << endl;
cout << "--------------------------------" ;
cout << "--------------------------------" ;
cout << endl;
double* piValues = x->getPiValues();
double* lambdaValues = x->getLambdaValues();
double* beliefs = x->getBeliefs();
for (int xi = 0; xi < x->getDomainSize(); xi++) {
cout << setw (10) << domain[xi];
cout << setw (19) << piValues[xi];
cout << setw (19) << lambdaValues[xi];
cout.precision (PRECISION_);
cout << setw (16) << beliefs[xi];
cout << endl;
}
cout << endl;
if (childs.size() > 0) {
string s = "(" + x->getVariableName() + ")" ;
for (unsigned int j = 0; j < childs.size(); j++) {
cout << setw (10) << "domain" ;
cout << setw (28) << "π" + childs[j]->getVariableName() + s;
cout << setw (28) << "λ" + childs[j]->getVariableName() + s;
cout << endl;
cout << "--------------------------------" ;
cout << "--------------------------------" ;
cout << endl;
for (int xi = 0; xi < x->getDomainSize(); xi++) {
cout << setw (10) << domain[xi];
cout.precision (PRECISION_);
cout << setw (27) << x->getPiMessage(childs[j], xi);
cout.precision (PRECISION_);
cout << setw (27) << x->getLambdaMessage(childs[j], xi);
cout << endl;
}
cout << endl;
}
}
}
void
BpNetwork::printBeliefs (void)
{
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*> (nodes_[i]);
vector<string> domain = x->getDomain();
cout << setw (20) << left << x->getVariableName() ;
cout << setw (26) << "belief" ;
cout << endl;
cout << "--------------------------------------" ;
cout << endl;
double* beliefs = x->getBeliefs();
for (int xi = 0; xi < x->getDomainSize(); xi++) {
cout << setw (20) << domain[xi];
cout.precision (PRECISION_);
cout << setw (26) << beliefs[xi];
cout << endl;
}
cout << endl;
}
}
vector<double>
BpNetwork::getBeliefs (void)
{
return getBeliefs (queryNode_);
}
vector<double>
BpNetwork::getBeliefs (BpNode* x)
{
double* beliefs = x->getBeliefs();
vector<double> beliefsVec;
for (int xi = 0; xi < x->getDomainSize(); xi++) {
beliefsVec.push_back (beliefs[xi]);
}
return beliefsVec;
}
void
BpNetwork::initializeSolver (vector<BayesianNode*> queryVars)
{
if (DL_ >= 1) {
cout << "Initializing solver" << endl;
if (schedule_ == SEQUENTIAL_SCHEDULE) {
cout << "-> schedule = sequential" << endl;
} else {
cout << "-> schedule = parallel" << endl;
}
cout << "-> max iters = " << maxIter_ << endl;
cout << "-> stable threashold = " << stableThreashold_ << endl;
cout << "-> query vars = " ;
for (unsigned int i = 0; i < queryVars.size(); i++) {
cout << queryVars[i]->getVariableName() << " " ;
}
cout << endl;
}
nIter_ = 0;
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* node = static_cast<BpNode*> (nodes_[i]);
node->allocateMemory();
}
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*> (nodes_[i]);
double* piValues = x->getPiValues();
double* lambdaValues = x->getLambdaValues();
for (int xi = 0; xi < x->getDomainSize(); xi++) {
piValues[xi] = 1.0;
lambdaValues[xi] = 1.0;
}
vector<BpNode*> xChilds = cast (x->getChilds());
for (unsigned int j = 0; j < xChilds.size(); j++) {
double* piMessages = x->getPiMessages (xChilds[j]);
for (int xi = 0; xi < x->getDomainSize(); xi++) {
piMessages[xi] = 1.0;
//x->setPiMessage (xChilds[j], xi, 1.0);
}
}
vector<BpNode*> xParents = cast (x->getParents());
for (unsigned int j = 0; j < xParents.size(); j++) {
double* lambdaMessages = xParents[j]->getLambdaMessages (x);
for (int xi = 0; xi < xParents[j]->getDomainSize(); xi++) {
lambdaMessages[xi] = 1.0;
//xParents[j]->setLambdaMessage (x, xi, 1.0);
}
}
}
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*> (nodes_[i]);
x->normalizeMessages();
}
printCurrentStatus();
vector<BpNode*> roots = cast (getRootNodes());
for (unsigned int i = 0; i < roots.size(); i++) {
double* params = roots[i]->getParameters();
double* piValues = roots[i]->getPiValues();
for (int ri = 0; ri < roots[i]->getDomainSize(); ri++) {
piValues[ri] = params[ri];
}
}
}
void
BpNetwork::addJunctionNode (vector<BayesianNode*> queryVars)
{
const string VAR_NAME = "_Jn";
int nStates = 1;
vector<BayesianNode*> parents;
vector<string> domain;
for (unsigned int i = 0; i < queryVars.size(); i++) {
parents.push_back (queryVars[i]);
nStates *= queryVars[i]->getDomainSize();
}
for (int i = 0; i < nStates; i++) {
stringstream ss;
ss << "_jn" << i;
domain.push_back (ss.str()); // FIXME make domain optional
}
int nParams = nStates * nStates;
double* params = new double [nParams];
for (int i = 0; i < nParams; i++) {
int row = i / nStates;
int col = i % nStates;
if (row == col) {
params[i] = 1;
} else {
params[i] = 0;
}
}
addNode (VAR_NAME, parents, params, nParams, domain);
queryNode_ = static_cast<BpNode*> (getNode (VAR_NAME));
printNetwork();
}
void
BpNetwork::addEvidence (BpNode* v)
{
if (DL_ >= 1) {
cout << "Adding evidence: node " ;
cout << "`" << v->getVariableName() << "' was instantiated as " ;
cout << "`" << v->getDomain()[v->getEvidence()] << "'" ;
cout << endl;
}
double* piValues = v->getPiValues();
double* lambdaValues = v->getLambdaValues();
for (int vi = 0; vi < v->getDomainSize(); vi++) {
if (vi == v->getEvidence()) {
piValues[vi] = 1.0;
lambdaValues[vi] = 1.0;
} else {
piValues[vi] = 0.0;
lambdaValues[vi] = 0.0;
}
}
}
void
BpNetwork::runNeapolitanSolver (void)
{
vector<BpNode*> roots = cast (getRootNodes());
for (unsigned int i = 0; i < roots.size(); i++) {
vector<BpNode*> childs = cast (roots[i]->getChilds());
for (unsigned int j = 0; j < childs.size(); j++) {
sendPiMessage (roots[i], static_cast<BpNode*> (childs[j]));
}
}
}
void
BpNetwork::sendPiMessage (BpNode* z, BpNode* x)
{
nIter_ ++;
if (!(maxIter_ == -1 || nIter_ < maxIter_)) {
cout << "the maximum number of iterations was achieved, terminating..." ;
cout << endl;
return;
}
if (DL_ >= 1) {
cout << "π message " << z->getVariableName();
cout << " --> " << x->getVariableName() << endl;
}
updatePiMessages(z, x);
if (!x->hasEvidence()) {
updatePiValues (x);
vector<BpNode*> xChilds = cast (x->getChilds());
for (unsigned int i = 0; i < xChilds.size(); i++) {
sendPiMessage (x, xChilds[i]);
}
}
bool isAllOnes = true;
double* lambdaValues = x->getLambdaValues();
for (int xi = 0; xi < x->getDomainSize(); xi++) {
if (lambdaValues[xi] != 1.0) {
isAllOnes = false;
break;
}
}
if (!isAllOnes) {
vector<BpNode*> xParents = cast (x->getParents());
for (unsigned int i = 0; i < xParents.size(); i++) {
if (xParents[i] != z && !xParents[i]->hasEvidence()) {
sendLambdaMessage (x, xParents[i]);
}
}
}
}
void
BpNetwork::sendLambdaMessage (BpNode* y, BpNode* x)
{
nIter_ ++;
if (!(maxIter_ == -1 || nIter_ < maxIter_)) {
cout << "the maximum number of iterations was achieved, terminating..." ;
cout << endl;
return;
}
if (DL_ >= 1) {
cout << "λ message " << y->getVariableName();
cout << " --> " << x->getVariableName() << endl;
}
updateLambdaMessages (x, y);
updateLambdaValues (x);
vector<BpNode*> xParents = cast (x->getParents());
for (unsigned int i = 0; i < xParents.size(); i++) {
if (!xParents[i]->hasEvidence()) {
sendLambdaMessage (x, xParents[i]);
}
}
vector<BpNode*> xChilds = cast (x->getChilds());
for (unsigned int i = 0; i < xChilds.size(); i++) {
if (xChilds[i] != y) {
sendPiMessage (x, xChilds[i]);
}
}
}
void
BpNetwork::updatePiValues (BpNode* x)
{
// π(Xi)
vector<BpNode*> parents = cast (x->getParents());
for (int xi = 0; xi < x->getDomainSize(); xi++) {
stringstream calcs1;
stringstream calcs2;
if (DL_ >= 2) {
calcs1 << "π("<< x->getDomain()[xi] << ")" << endl << "= " ;
}
double sum = 0.0;
vector<pair<int, int> > constraints;
vector<CptEntry> entries = x->getCptEntriesOfRow (xi);
for (unsigned int k = 0; k < entries.size(); k++) {
double prod = x->getProbability (entries[k]);
if (DL_ >= 2) {
if (k != 0) calcs1 << endl << "+ " ;
calcs1 << x->entryToString (entries[k]);
if (DL_ >= 3) {
(k == 0) ? calcs2 << "(" << prod : calcs2 << endl << "+ (" << prod;
}
}
vector<int> insts = entries[k].getDomainInstantiations();
for (unsigned int i = 0; i < parents.size(); i++) {
double value = parents[i]->getPiMessage (x, insts[i + 1]);
prod *= value;
if (DL_ >= 2) {
calcs1 << ".π" << x->getVariableName();
calcs1 << "(" << parents[i]->getDomain()[insts[i + 1]] << ")";
if (DL_ >= 3) calcs2 << "x" << value;
}
}
sum += prod;
if (DL_ >= 3) calcs2 << ")";
}
x->setPiValue (xi, sum);
if (DL_ >= 2) {
cout << calcs1.str();
if (DL_ >= 3) cout << endl << "= " << calcs2.str();
cout << " = " << sum << endl;
}
}
}
void
BpNetwork::updatePiMessages (BpNode* z, BpNode* x)
{
// πX(Zi)
vector<BpNode*> zChilds = cast (z->getChilds());
for (int zi = 0; zi < z->getDomainSize(); zi++) {
stringstream calcs1;
stringstream calcs2;
if (DL_ >= 2) {
calcs1 << "π" << x->getVariableName();
calcs1 << "(" << z->getDomain()[zi] << ") = " ;
}
double prod = z->getPiValue (zi);
if (DL_ >= 2) {
calcs1 << "π(" << z->getDomain()[zi] << ")" ;
if (DL_ >= 3) calcs2 << prod;
}
for (unsigned int i = 0; i < zChilds.size(); i++) {
if (zChilds[i] != x) {
double value = z->getLambdaMessage (zChilds[i], zi);
prod *= value;
if (DL_ >= 2) {
calcs1 << ".λ" << zChilds[i]->getVariableName();
calcs1 << "(" << z->getDomain()[zi] + ")" ;
if (DL_ >= 3) calcs2 << " x " << value;
}
}
}
z->setPiMessage (x, zi, prod);
if (DL_ >= 2) {
cout << calcs1.str();
if (DL_ >= 3) cout << " = " << calcs2.str();
cout << " = " << prod << endl;
}
}
}
void
BpNetwork::updateLambdaValues (BpNode* x)
{
// λ(Xi)
vector<BpNode*> childs = cast (x->getChilds());
for (int xi = 0; xi < x->getDomainSize(); xi++) {
stringstream calcs1;
stringstream calcs2;
if (DL_ >= 2) {
calcs1 << "λ" << "(" << x->getDomain()[xi] << ") = " ;
}
double prod = 1.0;
for (unsigned int i = 0; i < childs.size(); i++) {
double val = x->getLambdaMessage (childs[i], xi);
prod *= val;
if (DL_ >= 2) {
if (i != 0) calcs1 << "." ;
calcs1 << "λ" << childs[i]->getVariableName();
calcs1 << "(" << x->getDomain()[xi] + ")" ;
if (DL_ >= 3) (i == 0) ? calcs2 << val : calcs2 << " x " << val;
}
}
x->setLambdaValue (xi, prod);
if (DL_ >= 2) {
cout << calcs1.str();
if (childs.size() == 0) {
cout << 1 << endl;
} else {
if (DL_ >= 3) cout << " = " << calcs2.str();
cout << " = " << prod << endl;
}
}
}
}
void
BpNetwork::updateLambdaMessages (BpNode* x, BpNode* y)
{
// λY(Xi)
int parentIndex = y->getIndexOfParent (x) + 1;
vector<BpNode*> yParents = cast (y->getParents());
for (int xi = 0; xi < x->getDomainSize(); xi++) {
stringstream calcs1;
stringstream calcs2;
if (DL_ >= 2) {
calcs1 << "λ" << y->getVariableName() ;
calcs1 << "(" << x->getDomain()[xi] << ")" << endl << "= " ;
}
double outer_sum = 0.0;
for (int yi = 0; yi < y->getDomainSize(); yi++) {
if (DL_ >= 2) {
(yi == 0) ? calcs1 << "[" : calcs1 << endl << "+ [" ;
if (DL_ >= 3) {
(yi == 0) ? calcs2 << "[" : calcs2 << endl << "+ [" ;
}
}
double inner_sum = 0.0;
vector<pair<int, int> > constraints;
constraints.push_back (make_pair (0, yi));
constraints.push_back (make_pair (parentIndex, xi));
vector<CptEntry> entries = y->getCptEntries (constraints);
for (unsigned int k = 0; k < entries.size(); k++) {
double prod = y->getProbability (entries[k]);
if (DL_ >= 2) {
if (k != 0) calcs1 << " + " ;
calcs1 << y->entryToString (entries[k]);
if (DL_ >= 3) {
(k == 0) ? calcs2 << "(" << prod : calcs2 << " + (" << prod;
}
}
vector<int> insts = entries[k].getDomainInstantiations();
for (unsigned int i = 0; i < yParents.size(); i++) {
if (yParents[i] != x) {
double val = yParents[i]->getPiMessage (y, insts[i + 1]);
prod *= val;
if (DL_ >= 2) {
calcs1 << ".π" << y->getVariableName();
calcs1 << "(" << yParents[i]->getDomain()[insts[i + 1]] << ")" ;
if (DL_ >= 3) calcs2 << "x" << val;
}
}
}
inner_sum += prod;
if (DL_ >= 3) {
calcs2 << ")" ;
}
}
outer_sum += inner_sum * y->getLambdaValue (yi);
if (DL_ >= 2) {
calcs1 << "].λ(" << y->getDomain()[yi] << ")";
if (DL_ >= 3) calcs2 << "]x" << y->getLambdaValue (yi);
}
}
x->setLambdaMessage (y, xi, outer_sum);
if (DL_ >= 2) {
cout << calcs1.str();
if (DL_ >= 3) cout << endl << "= " << calcs2.str();
cout << " = " << outer_sum << endl;
}
}
}
void
BpNetwork::runIterativeBpSolver()
{
int nIter = 0;
maxIter_ = 100;
bool converged = false;
while (nIter < maxIter_ && !converged) {
if (DL_ >= 1) {
cout << endl << endl;
cout << "****************************************" ;
cout << "****************************************" ;
cout << endl;
cout << " Iteration " << nIter + 1 << endl;
cout << "****************************************" ;
cout << "****************************************" ;
}
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*>(nodes_[i]);
vector<BpNode*> xParents = cast (x->getParents());
for (unsigned int j = 0; j < xParents.size(); j++) {
//if (!xParents[j]->hasEvidence()) {
if (DL_ >= 1) {
cout << endl << "λ message " << x->getVariableName();
cout << " --> " << xParents[j]->getVariableName() << endl;
}
updateLambdaMessages (xParents[j], x);
//}
}
}
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*>(nodes_[i]);
vector<BpNode*> xChilds = cast (x->getChilds());
for (unsigned int j = 0; j < xChilds.size(); j++) {
if (DL_ >= 1) {
cout << endl << "π message " << x->getVariableName();
cout << " --> " << xChilds[j]->getVariableName() << endl;
}
updatePiMessages (x, xChilds[j]);
}
}
/*
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*>(nodes_[i]);
vector<BpNode*> xChilds = cast (x->getChilds());
for (unsigned int j = 0; j < xChilds.size(); j++) {
if (DL_ >= 1) {
cout << "π message " << x->getVariableName();
cout << " --> " << xChilds[j]->getVariableName() << endl;
}
updatePiMessages (x, xChilds[j]);
}
vector<BpNode*> xParents = cast (x->getParents());
for (unsigned int j = 0; j < xParents.size(); j++) {
//if (!xParents[j]->hasEvidence()) {
if (DL_ >= 1) {
cout << "λ message " << x->getVariableName();
cout << " --> " << xParents[j]->getVariableName() << endl;
}
updateLambdaMessages (xParents[j], x);
//}
}
}
*/
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*> (nodes_[i]);
//cout << endl << "SWAPING MESSAGES FOR " << x->getVariableName() << ":" ;
//cout << endl << endl;
//printCurrentStatusOf (x);
x->swapMessages();
x->normalizeMessages();
//cout << endl << "messages swaped " << endl;
//printCurrentStatusOf (x);
}
converged = true;
for (unsigned int i = 0; i < nodes_.size(); i++) {
BpNode* x = static_cast<BpNode*>(nodes_[i]);
if (DL_ >= 1) {
cout << endl << "var " << x->getVariableName() << ":" << endl;
}
//if (!x->hasEvidence()) {
updatePiValues (x);
updateLambdaValues (x);
double change = x->getBeliefChange();
if (DL_ >= 1) {
cout << "belief change = " << change << endl;
}
if (change > stableThreashold_) {
converged = false;
}
//}
}
if (converged) {
// converged = false;
}
if (DL_ >= 2) {
cout << endl;
printCurrentStatus();
}
nIter++;
}
if (DL_ >= 1) {
cout << endl;
if (converged) {
cout << "Iterative belief propagation converged in " ;
cout << nIter << " iterations" << endl;
} else {
cout << "Iterative belief propagation converged didn't converge" ;
cout << endl;
}
if (DL_ == 1) {
cout << endl;
printBeliefs();
}
cout << endl;
}
}
void
BpNetwork::addNode (string varName,
vector<BayesianNode*> parents,
int evidence,
int distId)
{
for (unsigned int i = 0; i < dists_.size(); i++) {
if (dists_[i]->id == distId) {
BpNode* node = new BpNode (varName, parents, dists_[i], evidence);
nodes_.push_back (node);
break;
}
}
}
void
BpNetwork::addNode (string varName,
vector<BayesianNode*> parents,
double* params,
int nParams,
vector<string> domain)
{
Distribution* dist = new Distribution (params, nParams, domain);
BpNode* node = new BpNode (varName, parents, dist);
dists_.push_back (dist);
nodes_.push_back (node);
}
vector<BpNode*>
BpNetwork::cast (vector<BayesianNode*> nodes)
{
vector<BpNode*> castedNodes (nodes.size());
for (unsigned int i = 0; i < nodes.size(); i++) {
castedNodes[i] = static_cast<BpNode*> (nodes[i]);
}
return castedNodes;
}

66
packages/CLPBN/clpbn/bp/BpNetwork.h
View File

@ -1,66 +0,0 @@
#ifndef BP_BP_NETWORK_H
#define BP_BP_NETWORK_H
#include <vector>
#include <string>
#include "BayesianNetwork.h"
using namespace std;
class BpNode;
enum Schedule
{
SEQUENTIAL_SCHEDULE,
PARALLEL_SCHEDULE
};
class BpNetwork : public BayesianNetwork
{
public:
// constructs
BpNetwork (void);
// destruct
~BpNetwork (void);
// methods
void setSolverParameters (Schedule, int, double);
void runSolver (BayesianNode* queryVar);
void runSolver (vector<BayesianNode*>);
void printCurrentStatus (void);
void printCurrentStatusOf (BpNode*);
void printBeliefs (void);
vector<double> getBeliefs (void);
vector<double> getBeliefs (BpNode*);
private:
BpNetwork (const BpNetwork&); // disallow copy
void operator= (const BpNetwork&); // disallow assign
// methods
void initializeSolver (vector<BayesianNode*>);
void addJunctionNode (vector<BayesianNode*>);
void addEvidence (BpNode*);
void runNeapolitanSolver (void);
void sendLambdaMessage (BpNode*, BpNode*);
void sendPiMessage (BpNode*, BpNode*);
void updatePiValues (BpNode*);
void updatePiMessages (BpNode*, BpNode*);
void updateLambdaValues (BpNode*);
void updateLambdaMessages (BpNode*, BpNode*);
void runIterativeBpSolver (void);
void addNode (string, vector<BayesianNode*>, int, int);
void addNode (string, vector<BayesianNode*>,
double*, int, vector<string>);
vector<BpNode*> cast (vector<BayesianNode*>);
// members
Schedule schedule_;
int nIter_;
int maxIter_;
double stableThreashold_;
BpNode* queryNode_;
static const int DL_ = 3;
static const int PRECISION_ = 10;
};
#endif // BP_BP_NETWORK_H

250
packages/CLPBN/clpbn/bp/BpNode.cpp
View File

@ -1,250 +0,0 @@
#include <iostream>
#include <cassert>
#include <cmath>
#include "BpNode.h"
bool BpNode::calculateMessageResidual_ = true;
BpNode::BpNode (BayesNode* node)
{
ds_ = node->getDomainSize();
const NodeSet& childs = node->getChilds();
piVals_.resize (ds_, 1);
ldVals_.resize (ds_, 1);
if (calculateMessageResidual_) {
piResiduals_.resize (childs.size(), 0.0);
ldResiduals_.resize (childs.size(), 0.0);
}
childs_ = &childs;
for (unsigned i = 0; i < childs.size(); i++) {
//indexMap_.insert (make_pair (childs[i]->getVarId(), i));
currPiMsgs_.push_back (ParamSet (ds_, 1));
currLdMsgs_.push_back (ParamSet (ds_, 1));
nextPiMsgs_.push_back (ParamSet (ds_, 1));
nextLdMsgs_.push_back (ParamSet (ds_, 1));
}
}
ParamSet
BpNode::getBeliefs (void) const
{
double sum = 0.0;
ParamSet beliefs (ds_);
for (int xi = 0; xi < ds_; xi++) {
double prod = piVals_[xi] * ldVals_[xi];
beliefs[xi] = prod;
sum += prod;
}
assert (sum);
//normalize the beliefs
for (int xi = 0; xi < ds_; xi++) {
beliefs[xi] /= sum;
}
return beliefs;
}
double
BpNode::getPiValue (int idx) const
{
assert (idx >=0 && idx < ds_);
return piVals_[idx];
}
void
BpNode::setPiValue (int idx, double value)
{
assert (idx >=0 && idx < ds_);
piVals_[idx] = value;
}
double
BpNode::getLambdaValue (int idx) const
{
assert (idx >=0 && idx < ds_);
return ldVals_[idx];
}
void
BpNode::setLambdaValue (int idx, double value)
{
assert (idx >=0 && idx < ds_);
ldVals_[idx] = value;
}
ParamSet&
BpNode::getPiValues (void)
{
return piVals_;
}
ParamSet&
BpNode::getLambdaValues (void)
{
return ldVals_;
}
double
BpNode::getPiMessageValue (const BayesNode* destination, int idx) const
{
assert (idx >=0 && idx < ds_);
return currPiMsgs_[getIndex(destination)][idx];
}
double
BpNode::getLambdaMessageValue (const BayesNode* source, int idx) const
{
assert (idx >=0 && idx < ds_);
return currLdMsgs_[getIndex(source)][idx];
}
const ParamSet&
BpNode::getPiMessage (const BayesNode* destination) const
{
return currPiMsgs_[getIndex(destination)];
}
const ParamSet&
BpNode::getLambdaMessage (const BayesNode* source) const
{
return currLdMsgs_[getIndex(source)];
}
ParamSet&
BpNode::piNextMessageReference (const BayesNode* destination)
{
return nextPiMsgs_[getIndex(destination)];
}
ParamSet&
BpNode::lambdaNextMessageReference (const BayesNode* source)
{
return nextLdMsgs_[getIndex(source)];
}
void
BpNode::updatePiMessage (const BayesNode* destination)
{
int idx = getIndex (destination);
currPiMsgs_[idx] = nextPiMsgs_[idx];
Util::normalize (currPiMsgs_[idx]);
}
void
BpNode::updateLambdaMessage (const BayesNode* source)
{
int idx = getIndex (source);
currLdMsgs_[idx] = nextLdMsgs_[idx];
Util::normalize (currLdMsgs_[idx]);
}
double
BpNode::getBeliefChange (void)
{
double change = 0.0;
if (oldBeliefs_.size() == 0) {
oldBeliefs_ = getBeliefs();
change = 9999999999.0;
} else {
ParamSet currentBeliefs = getBeliefs();
for (int xi = 0; xi < ds_; xi++) {
change += abs (currentBeliefs[xi] - oldBeliefs_[xi]);
}
oldBeliefs_ = currentBeliefs;
}
return change;
}
void
BpNode::updatePiResidual (const BayesNode* destination)
{
int idx = getIndex (destination);
Util::normalize (nextPiMsgs_[idx]);
//piResiduals_[idx] = Util::getL1dist (
// currPiMsgs_[idx], nextPiMsgs_[idx]);
piResiduals_[idx] = Util::getMaxNorm (
currPiMsgs_[idx], nextPiMsgs_[idx]);
}
void
BpNode::updateLambdaResidual (const BayesNode* source)
{
int idx = getIndex (source);
Util::normalize (nextLdMsgs_[idx]);
//ldResiduals_[idx] = Util::getL1dist (
// currLdMsgs_[idx], nextLdMsgs_[idx]);
ldResiduals_[idx] = Util::getMaxNorm (
currLdMsgs_[idx], nextLdMsgs_[idx]);
}
void
BpNode::clearPiResidual (const BayesNode* destination)
{
piResiduals_[getIndex(destination)] = 0;
}
void
BpNode::clearLambdaResidual (const BayesNode* source)
{
ldResiduals_[getIndex(source)] = 0;
}
bool
BpNode::hasReceivedChildInfluence (void) const
{
// if all lambda values are equal, then neither
// this node neither its descendents have evidence,
// we can use this to don't send lambda messages his parents
bool childInfluenced = false;
for (int xi = 1; xi < ds_; xi++) {
if (ldVals_[xi] != ldVals_[0]) {
childInfluenced = true;
break;
}
}
return childInfluenced;
}

99
packages/CLPBN/clpbn/bp/BpNode.h
View File

@ -1,99 +0,0 @@
#ifndef BP_BPNODE_H
#define BP_BPNODE_H
#include <vector>
#include <map>
#include <string>
#include <unordered_map>
#include "BayesNode.h"
#include "Shared.h"
using namespace std;
class BpNode
{
public:
BpNode (int);
BpNode (BayesNode*);
ParamSet getBeliefs (void) const;
double getPiValue (int) const;
void setPiValue (int, double);
double getLambdaValue (int) const;
void setLambdaValue (int, double);
ParamSet& getPiValues (void);
ParamSet& getLambdaValues (void);
double getPiMessageValue (const BayesNode*, int) const;
double getLambdaMessageValue (const BayesNode*, int) const;
const ParamSet& getPiMessage (const BayesNode*) const;
const ParamSet& getLambdaMessage (const BayesNode*) const;
ParamSet& piNextMessageReference (const BayesNode*);
ParamSet& lambdaNextMessageReference (const BayesNode*);
void updatePiMessage (const BayesNode*);
void updateLambdaMessage (const BayesNode*);
double getBeliefChange (void);
void updatePiResidual (const BayesNode*);
void updateLambdaResidual (const BayesNode*);
void clearPiResidual (const BayesNode*);
void clearLambdaResidual (const BayesNode*);
bool hasReceivedChildInfluence (void) const;
// inlines
double getPiResidual (const BayesNode*);
double getLambdaResidual (const BayesNode*);
int getIndex (const BayesNode*) const;
private:
DISALLOW_COPY_AND_ASSIGN (BpNode);
IndexMap indexMap_;
ParamSet piVals_; // pi values
ParamSet ldVals_; // lambda values
vector<ParamSet> currPiMsgs_; // current pi messages
vector<ParamSet> currLdMsgs_; // current lambda messages
vector<ParamSet> nextPiMsgs_;
vector<ParamSet> nextLdMsgs_;
ParamSet oldBeliefs_;
ParamSet piResiduals_;
ParamSet ldResiduals_;
int ds_;
const NodeSet* childs_;
static bool calculateMessageResidual_;
// static const double MAX_CHANGE_ = 10000000.0;
};
inline double
BpNode::getPiResidual (const BayesNode* destination)
{
return piResiduals_[getIndex(destination)];
}
inline double
BpNode::getLambdaResidual (const BayesNode* source)
{
return ldResiduals_[getIndex(source)];
}
inline int
BpNode::getIndex (const BayesNode* node) const
{
assert (node);
//assert (indexMap_.find(node->getVarId()) != indexMap_.end());
//return indexMap_.find(node->getVarId())->second;
for (unsigned i = 0; childs_->size(); i++) {
if ((*childs_)[i]->getVarId() == node->getVarId()) {
return i;
}
}
assert (false);
return -1;
}
#endif

198
packages/CLPBN/clpbn/bp/CountingBP.cpp
View File

@ -1,198 +0,0 @@
#include "CountingBP.h"
CountingBP::~CountingBP (void)
{
delete lfg_;
delete fg_;
for (unsigned i = 0; i < links_.size(); i++) {
delete links_[i];
}
links_.clear();
}
ParamSet
CountingBP::getPosterioriOf (Vid vid) const
{
FgVarNode* var = lfg_->getEquivalentVariable (vid);
ParamSet probs;
if (var->hasEvidence()) {
probs.resize (var->getDomainSize(), 0.0);
probs[var->getEvidence()] = 1.0;
} else {
probs.resize (var->getDomainSize(), 1.0);
CLinkSet links = varsI_[var->getIndex()]->getLinks();
for (unsigned i = 0; i < links.size(); i++) {
ParamSet msg = links[i]->getMessage();
CountingBPLink* l = static_cast<CountingBPLink*> (links[i]);
Util::pow (msg, l->getNumberOfEdges());
for (unsigned j = 0; j < msg.size(); j++) {
probs[j] *= msg[j];
}
}
Util::normalize (probs);
}
return probs;
}
void
CountingBP::initializeSolver (void)
{
lfg_ = new LiftedFG (*fg_);
unsigned nUncVars = fg_->getFgVarNodes().size();
unsigned nUncFactors = fg_->getFactors().size();
CFgVarSet vars = fg_->getFgVarNodes();
unsigned nNeighborLessVars = 0;
for (unsigned i = 0; i < vars.size(); i++) {
CFactorSet factors = vars[i]->getFactors();
if (factors.size() == 1 && factors[0]->getFgVarNodes().size() == 1) {
nNeighborLessVars ++;
}
}
// cout << "UNCOMPRESSED FACTOR GRAPH" << endl;
// fg_->printGraphicalModel();
fg_->exportToDotFormat ("uncompress.dot");
FactorGraph *temp;
temp = fg_;
fg_ = lfg_->getCompressedFactorGraph();
unsigned nCompVars = fg_->getFgVarNodes().size();
unsigned nCompFactors = fg_->getFactors().size();
Statistics::updateCompressingStats (nUncVars,
nUncFactors,
nCompVars,
nCompFactors,
nNeighborLessVars);
cout << "COMPRESSED FACTOR GRAPH" << endl;
fg_->printGraphicalModel();
//fg_->exportToDotFormat ("compress.dot");
SPSolver::initializeSolver();
}
void
CountingBP::createLinks (void)
{
const FactorClusterSet fcs = lfg_->getFactorClusters();
for (unsigned i = 0; i < fcs.size(); i++) {
const VarClusterSet vcs = fcs[i]->getVarClusters();
for (unsigned j = 0; j < vcs.size(); j++) {
unsigned c = lfg_->getGroundEdgeCount (fcs[i], vcs[j]);
links_.push_back (
new CountingBPLink (fcs[i]->getRepresentativeFactor(),
vcs[j]->getRepresentativeVariable(), c));
//cout << (links_.back())->toString() << " edge count =" << c << endl;
}
}
return;
}
void
CountingBP::deleteJunction (Factor* f, FgVarNode*)
{
f->freeDistribution();
}
void
CountingBP::maxResidualSchedule (void)
{
if (nIter_ == 1) {
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->setNextMessage (getFactor2VarMsg (links_[i]));
SortedOrder::iterator it = sortedOrder_.insert (links_[i]);
linkMap_.insert (make_pair (links_[i], it));
if (DL >= 2 && DL < 5) {
cout << "calculating " << links_[i]->toString() << endl;
}
}
return;
}
for (unsigned c = 0; c < links_.size(); c++) {
if (DL >= 2) {
cout << endl << "current residuals:" << endl;
for (SortedOrder::iterator it = sortedOrder_.begin();
it != sortedOrder_.end(); it ++) {
cout << " " << setw (30) << left << (*it)->toString();
cout << "residual = " << (*it)->getResidual() << endl;
}
}
SortedOrder::iterator it = sortedOrder_.begin();
Link* link = *it;
if (DL >= 2) {
cout << "updating " << (*sortedOrder_.begin())->toString() << endl;
}
if (link->getResidual() < SolverOptions::accuracy) {
return;
}
link->updateMessage();
link->clearResidual();
sortedOrder_.erase (it);
linkMap_.find (link)->second = sortedOrder_.insert (link);
// update the messages that depend on message source --> destin
CFactorSet factorNeighbors = link->getVariable()->getFactors();
for (unsigned i = 0; i < factorNeighbors.size(); i++) {
CLinkSet links = factorsI_[factorNeighbors[i]->getIndex()]->getLinks();
for (unsigned j = 0; j < links.size(); j++) {
if (links[j]->getVariable() != link->getVariable()) { //FIXMEFIXME
if (DL >= 2 && DL < 5) {
cout << " calculating " << links[j]->toString() << endl;
}
links[j]->setNextMessage (getFactor2VarMsg (links[j]));
LinkMap::iterator iter = linkMap_.find (links[j]);
sortedOrder_.erase (iter->second);
iter->second = sortedOrder_.insert (links[j]);
}
}
}
}
}
ParamSet
CountingBP::getVar2FactorMsg (const Link* link) const
{
const FgVarNode* src = link->getVariable();
const Factor* dest = link->getFactor();
ParamSet msg;
if (src->hasEvidence()) {
cout << "has evidence" << endl;
msg.resize (src->getDomainSize(), 0.0);
msg[src->getEvidence()] = link->getMessage()[src->getEvidence()];
cout << "-> " << link->getVariable()->getLabel() << " " << link->getFactor()->getLabel() << endl;
cout << "-> p2s " << Util::parametersToString (msg) << endl;
} else {
msg = link->getMessage();
}
const CountingBPLink* l = static_cast<const CountingBPLink*> (link);
Util::pow (msg, l->getNumberOfEdges() - 1);
CLinkSet links = varsI_[src->getIndex()]->getLinks();
for (unsigned i = 0; i < links.size(); i++) {
if (links[i]->getFactor() != dest) {
ParamSet msgFromFactor = links[i]->getMessage();
CountingBPLink* l = static_cast<CountingBPLink*> (links[i]);
Util::pow (msgFromFactor, l->getNumberOfEdges());
for (unsigned j = 0; j < msgFromFactor.size(); j++) {
msg[j] *= msgFromFactor[j];
}
}
}
return msg;
}

45
packages/CLPBN/clpbn/bp/CountingBP.h
View File

@ -1,45 +0,0 @@
#ifndef BP_COUNTING_BP_H
#define BP_COUNTING_BP_H
#include "SPSolver.h"
#include "LiftedFG.h"
class Factor;
class FgVarNode;
class CountingBPLink : public Link
{
public:
CountingBPLink (Factor* f, FgVarNode* v, unsigned c) : Link (f, v)
{
edgeCount_ = c;
}
unsigned getNumberOfEdges (void) const { return edgeCount_; }
private:
unsigned edgeCount_;
};
class CountingBP : public SPSolver
{
public:
CountingBP (FactorGraph& fg) : SPSolver (fg) { }
~CountingBP (void);
ParamSet getPosterioriOf (Vid) const;
private:
void initializeSolver (void);
void createLinks (void);
void deleteJunction (Factor*, FgVarNode*);
void maxResidualSchedule (void);
ParamSet getVar2FactorMsg (const Link*) const;
LiftedFG* lfg_;
};
#endif // BP_COUNTING_BP_H

40
packages/CLPBN/clpbn/bp/Distribution.cpp
View File

@ -1,40 +0,0 @@
#include <vector>
#include <string>
#include <Distribution.h>
Distribution::Distribution (int id,
double* params,
int nParams,
vector<string> domain)
{
this->id = id;
this->params = params;
this->nParams = nParams;
this->domain = domain;
}
Distribution::Distribution (double* params,
int nParams,
vector<string> domain)
{
this->id = -1;
this->params = params;
this->nParams = nParams;
this->domain = domain;
}
/*
Distribution::~Distribution()
{
delete params;
for (unsigned int i = 0; i < cptEntries.size(); i++) {
delete cptEntries[i];
}
}
*/

43
packages/CLPBN/clpbn/bp/FgVarNode.h
View File

@ -1,43 +0,0 @@
#ifndef BP_FG_VAR_NODE_H
#define BP_FG_VAR_NODE_H
#include <vector>
#include "Variable.h"
#include "Shared.h"
using namespace std;
class Factor;
class FgVarNode : public Variable
{
public:
FgVarNode (unsigned vid, unsigned dsize) : Variable (vid, dsize) { }
FgVarNode (const Variable* v) : Variable (v) { }
void addFactor (Factor* f) { factors_.push_back (f); }
CFactorSet getFactors (void) const { return factors_; }
void removeFactor (const Factor* f)
{
if (factors_[factors_.size() -1] == f) {
factors_.pop_back();
} else {
for (unsigned i = 0; i < factors_.size(); i++) {
if (factors_[i] == f) {
factors_.erase (factors_.begin() + i);
return;
}
}
assert (false);
}
}
private:
DISALLOW_COPY_AND_ASSIGN (FgVarNode);
// members
FactorSet factors_;
};
#endif // BP_FG_VAR_NODE_H

278
packages/CLPBN/clpbn/bp/LiftedFG.cpp
View File

@ -1,278 +0,0 @@
#include "LiftedFG.h"
#include "FgVarNode.h"
#include "Factor.h"
#include "Distribution.h"
LiftedFG::LiftedFG (const FactorGraph& fg)
{
groundFg_ = &fg;
freeColor_ = 0;
const FgVarSet& varNodes = fg.getFgVarNodes();
const FactorSet& factors = fg.getFactors();
varColors_.resize (varNodes.size());
factorColors_.resize (factors.size());
for (unsigned i = 0; i < factors.size(); i++) {
factors[i]->setIndex (i);
}
// create the initial variable colors
VarColorMap colorMap;
for (unsigned i = 0; i < varNodes.size(); i++) {
unsigned dsize = varNodes[i]->getDomainSize();
VarColorMap::iterator it = colorMap.find (dsize);
if (it == colorMap.end()) {
it = colorMap.insert (make_pair (
dsize, vector<Color> (dsize + 1,-1))).first;
}
unsigned idx;
if (varNodes[i]->hasEvidence()) {
idx = varNodes[i]->getEvidence();
} else {
idx = dsize;
}
vector<Color>& stateColors = it->second;
if (stateColors[idx] == -1) {
stateColors[idx] = getFreeColor();
}
setColor (varNodes[i], stateColors[idx]);
}
// create the initial factor colors
DistColorMap distColors;
for (unsigned i = 0; i < factors.size(); i++) {
Distribution* dist = factors[i]->getDistribution();
DistColorMap::iterator it = distColors.find (dist);
if (it == distColors.end()) {
it = distColors.insert (make_pair (dist, getFreeColor())).first;
}
setColor (factors[i], it->second);
}
VarSignMap varGroups;
FactorSignMap factorGroups;
bool groupsHaveChanged = true;
unsigned nIter = 0;
while (groupsHaveChanged || nIter == 1) {
nIter ++;
if (Statistics::numCreatedNets == 4) {
cout << "--------------------------------------------" << endl;
cout << "Iteration " << nIter << endl;
cout << "--------------------------------------------" << endl;
}
unsigned prevFactorGroupsSize = factorGroups.size();
factorGroups.clear();
// set a new color to the factors with the same signature
for (unsigned i = 0; i < factors.size(); i++) {
const string& signatureId = getSignatureId (factors[i]);
// cout << factors[i]->getLabel() << " signature: " ;
// cout<< signatureId << endl;
FactorSignMap::iterator it = factorGroups.find (signatureId);
if (it == factorGroups.end()) {
it = factorGroups.insert (make_pair (signatureId, FactorSet())).first;
}
it->second.push_back (factors[i]);
}
if (nIter > 0)
for (FactorSignMap::iterator it = factorGroups.begin();
it != factorGroups.end(); it++) {
Color newColor = getFreeColor();
FactorSet& groupMembers = it->second;
for (unsigned i = 0; i < groupMembers.size(); i++) {
setColor (groupMembers[i], newColor);
}
}
// set a new color to the variables with the same signature
unsigned prevVarGroupsSize = varGroups.size();
varGroups.clear();
for (unsigned i = 0; i < varNodes.size(); i++) {
const string& signatureId = getSignatureId (varNodes[i]);
VarSignMap::iterator it = varGroups.find (signatureId);
// cout << varNodes[i]->getLabel() << " signature: " ;
// cout << signatureId << endl;
if (it == varGroups.end()) {
it = varGroups.insert (make_pair (signatureId, FgVarSet())).first;
}
it->second.push_back (varNodes[i]);
}
if (nIter > 0)
for (VarSignMap::iterator it = varGroups.begin();
it != varGroups.end(); it++) {
Color newColor = getFreeColor();
FgVarSet& groupMembers = it->second;
for (unsigned i = 0; i < groupMembers.size(); i++) {
setColor (groupMembers[i], newColor);
}
}
//if (nIter >= 3) cout << "bigger than three: " << nIter << endl;
groupsHaveChanged = prevVarGroupsSize != varGroups.size()
|| prevFactorGroupsSize != factorGroups.size();
}
printGroups (varGroups, factorGroups);
for (VarSignMap::iterator it = varGroups.begin();
it != varGroups.end(); it++) {
CFgVarSet vars = it->second;
VarCluster* vc = new VarCluster (vars);
for (unsigned i = 0; i < vars.size(); i++) {
vid2VarCluster_.insert (make_pair (vars[i]->getVarId(), vc));
}
varClusters_.push_back (vc);
}
for (FactorSignMap::iterator it = factorGroups.begin();
it != factorGroups.end(); it++) {
VarClusterSet varClusters;
Factor* groundFactor = it->second[0];
FgVarSet groundVars = groundFactor->getFgVarNodes();
for (unsigned i = 0; i < groundVars.size(); i++) {
Vid vid = groundVars[i]->getVarId();
varClusters.push_back (vid2VarCluster_.find (vid)->second);
}
factorClusters_.push_back (new FactorCluster (it->second, varClusters));
}
}
LiftedFG::~LiftedFG (void)
{
for (unsigned i = 0; i < varClusters_.size(); i++) {
delete varClusters_[i];
}
for (unsigned i = 0; i < factorClusters_.size(); i++) {
delete factorClusters_[i];
}
}
string
LiftedFG::getSignatureId (FgVarNode* var) const
{
stringstream ss;
CFactorSet myFactors = var->getFactors();
ss << myFactors.size();
for (unsigned i = 0; i < myFactors.size(); i++) {
ss << "." << getColor (myFactors[i]);
ss << "." << myFactors[i]->getIndexOf(var);
}
ss << "." << getColor (var);
return ss.str();
}
string
LiftedFG::getSignatureId (Factor* factor) const
{
stringstream ss;
CFgVarSet myVars = factor->getFgVarNodes();
ss << myVars.size();
for (unsigned i = 0; i < myVars.size(); i++) {
ss << "." << getColor (myVars[i]);
}
ss << "." << getColor (factor);
return ss.str();
}
FactorGraph*
LiftedFG::getCompressedFactorGraph (void)
{
FactorGraph* fg = new FactorGraph();
for (unsigned i = 0; i < varClusters_.size(); i++) {
FgVarNode* var = varClusters_[i]->getGroundFgVarNodes()[0];
FgVarNode* newVar = new FgVarNode (var);
newVar->setIndex (i);
varClusters_[i]->setRepresentativeVariable (newVar);
fg->addVariable (newVar);
}
for (unsigned i = 0; i < factorClusters_.size(); i++) {
FgVarSet myGroundVars;
const VarClusterSet& myVarClusters = factorClusters_[i]->getVarClusters();
for (unsigned j = 0; j < myVarClusters.size(); j++) {
myGroundVars.push_back (myVarClusters[j]->getRepresentativeVariable());
}
Factor* newFactor = new Factor (myGroundVars,
factorClusters_[i]->getGroundFactors()[0]->getDistribution());
factorClusters_[i]->setRepresentativeFactor (newFactor);
fg->addFactor (newFactor);
}
return fg;
}
unsigned
LiftedFG::getGroundEdgeCount (FactorCluster* fc, VarCluster* vc) const
{
CFactorSet clusterGroundFactors = fc->getGroundFactors();
FgVarNode* var = vc->getGroundFgVarNodes()[0];
unsigned count = 0;
for (unsigned i = 0; i < clusterGroundFactors.size(); i++) {
if (clusterGroundFactors[i]->getIndexOf (var) != -1) {
count ++;
}
}
/*
CFgVarSet vars = vc->getGroundFgVarNodes();
for (unsigned i = 1; i < vars.size(); i++) {
FgVarNode* var = vc->getGroundFgVarNodes()[i];
unsigned count2 = 0;
for (unsigned i = 0; i < clusterGroundFactors.size(); i++) {
if (clusterGroundFactors[i]->getIndexOf (var) != -1) {
count2 ++;
}
}
if (count != count2) { cout << "oops!" << endl; abort(); }
}
*/
return count;
}
void
LiftedFG::printGroups (const VarSignMap& varGroups,
const FactorSignMap& factorGroups) const
{
cout << "variable groups:" << endl;
unsigned count = 0;
for (VarSignMap::const_iterator it = varGroups.begin();
it != varGroups.end(); it++) {
const FgVarSet& groupMembers = it->second;
if (groupMembers.size() > 0) {
cout << ++count << ": " ;
//if (groupMembers.size() > 1) {
for (unsigned i = 0; i < groupMembers.size(); i++) {
cout << groupMembers[i]->getLabel() << " " ;
}
//}
cout << endl;
}
}
cout << endl;
cout << "factor groups:" << endl;
count = 0;
for (FactorSignMap::const_iterator it = factorGroups.begin();
it != factorGroups.end(); it++) {
const FactorSet& groupMembers = it->second;
if (groupMembers.size() > 0) {
cout << ++count << ": " ;
//if (groupMembers.size() > 1) {
for (unsigned i = 0; i < groupMembers.size(); i++) {
cout << groupMembers[i]->getLabel() << " " ;
}
//}
cout << endl;
}
}
}

152
packages/CLPBN/clpbn/bp/LiftedFG.h
View File

@ -1,152 +0,0 @@
#ifndef BP_LIFTED_FG_H
#define BP_LIFTED_FG_H
#include <unordered_map>
#include "FactorGraph.h"
#include "FgVarNode.h"
#include "Factor.h"
#include "Shared.h"
class VarCluster;
class FactorCluster;
class Distribution;
typedef long Color;
typedef vector<Color> Signature;
typedef vector<VarCluster*> VarClusterSet;
typedef vector<FactorCluster*> FactorClusterSet;
typedef map<string, FgVarSet> VarSignMap;
typedef map<string, FactorSet> FactorSignMap;
typedef map<unsigned, vector<Color> > VarColorMap;
typedef map<Distribution*, Color> DistColorMap;
typedef map<Vid, VarCluster*> Vid2VarCluster;
class VarCluster
{
public:
VarCluster (CFgVarSet vs)
{
for (unsigned i = 0; i < vs.size(); i++) {
groundVars_.push_back (vs[i]);
}
}
void addFactorCluster (FactorCluster* fc)
{
factorClusters_.push_back (fc);
}
const FactorClusterSet& getFactorClusters (void) const
{
return factorClusters_;
}
FgVarNode* getRepresentativeVariable (void) const { return representVar_; }
void setRepresentativeVariable (FgVarNode* v) { representVar_ = v; }
CFgVarSet getGroundFgVarNodes (void) const { return groundVars_; }
private:
FgVarSet groundVars_;
FactorClusterSet factorClusters_;
FgVarNode* representVar_;
};
class FactorCluster
{
public:
FactorCluster (CFactorSet groundFactors, const VarClusterSet& vcs)
{
groundFactors_ = groundFactors;
varClusters_ = vcs;
for (unsigned i = 0; i < varClusters_.size(); i++) {
varClusters_[i]->addFactorCluster (this);
}
}
const VarClusterSet& getVarClusters (void) const
{
return varClusters_;
}
bool containsGround (const Factor* f)
{
for (unsigned i = 0; i < groundFactors_.size(); i++) {
if (groundFactors_[i] == f) {
return true;
}
}
return false;
}
Factor* getRepresentativeFactor (void) const { return representFactor_; }
void setRepresentativeFactor (Factor* f) { representFactor_ = f; }
CFactorSet getGroundFactors (void) const { return groundFactors_; }
private:
FactorSet groundFactors_;
VarClusterSet varClusters_;
Factor* representFactor_;
};
class LiftedFG
{
public:
LiftedFG (const FactorGraph&);
~LiftedFG (void);
FactorGraph* getCompressedFactorGraph (void);
unsigned getGroundEdgeCount (FactorCluster*, VarCluster*) const;
void printGroups (const VarSignMap& varGroups,
const FactorSignMap& factorGroups) const;
FgVarNode* getEquivalentVariable (Vid vid)
{
VarCluster* vc = vid2VarCluster_.find (vid)->second;
return vc->getRepresentativeVariable();
}
const VarClusterSet& getVariableClusters (void) { return varClusters_; }
const FactorClusterSet& getFactorClusters (void) { return factorClusters_; }
private:
string getSignatureId (FgVarNode*) const;
string getSignatureId (Factor*) const;
Color getFreeColor (void) { return ++freeColor_ -1; }
Color getColor (FgVarNode* v) const { return varColors_[v->getIndex()]; }
Color getColor (Factor* f) const { return factorColors_[f->getIndex()]; }
void setColor (FgVarNode* v, Color c)
{
varColors_[v->getIndex()] = c;
}
void setColor (Factor* f, Color c)
{
factorColors_[f->getIndex()] = c;
}
VarCluster* getVariableCluster (Vid vid) const
{
return vid2VarCluster_.find (vid)->second;
}
Color freeColor_;
vector<Color> varColors_;
vector<Color> factorColors_;
VarClusterSet varClusters_;
FactorClusterSet factorClusters_;
Vid2VarCluster vid2VarCluster_;
const FactorGraph* groundFg_;
};
#endif // BP_LIFTED_FG_H

470
packages/CLPBN/clpbn/bp/SPSolver.cpp
View File

@ -1,470 +0,0 @@
#include <cassert>
#include <limits>
#include <iostream>
#include "SPSolver.h"
#include "FactorGraph.h"
#include "FgVarNode.h"
#include "Factor.h"
#include "Shared.h"
SPSolver::SPSolver (FactorGraph& fg) : Solver (&fg)
{
fg_ = &fg;
}
SPSolver::~SPSolver (void)
{
for (unsigned i = 0; i < varsI_.size(); i++) {
delete varsI_[i];
}
for (unsigned i = 0; i < factorsI_.size(); i++) {
delete factorsI_[i];
}
for (unsigned i = 0; i < links_.size(); i++) {
delete links_[i];
}
}
void
SPSolver::runTreeSolver (void)
{
CFactorSet factors = fg_->getFactors();
bool finish = false;
while (!finish) {
finish = true;
for (unsigned i = 0; i < factors.size(); i++) {
CLinkSet links = factorsI_[factors[i]->getIndex()]->getLinks();
for (unsigned j = 0; j < links.size(); j++) {
if (!links[j]->messageWasSended()) {
if (readyToSendMessage(links[j])) {
links[j]->setNextMessage (getFactor2VarMsg (links[j]));
links[j]->updateMessage();
}
finish = false;
}
}
}
}
}
bool
SPSolver::readyToSendMessage (const Link* link) const
{
CFgVarSet factorVars = link->getFactor()->getFgVarNodes();
for (unsigned i = 0; i < factorVars.size(); i++) {
if (factorVars[i] != link->getVariable()) {
CLinkSet links = varsI_[factorVars[i]->getIndex()]->getLinks();
for (unsigned j = 0; j < links.size(); j++) {
if (links[j]->getFactor() != link->getFactor() &&
!links[j]->messageWasSended()) {
return false;
}
}
}
}
return true;
}
void
SPSolver::runSolver (void)
{
initializeSolver();
runTreeSolver();
return;
nIter_ = 0;
while (!converged() && nIter_ < SolverOptions::maxIter) {
nIter_ ++;
if (DL >= 2) {
cout << endl;
cout << "****************************************" ;
cout << "****************************************" ;
cout << endl;
cout << " Iteration " << nIter_ << endl;
cout << "****************************************" ;
cout << "****************************************" ;
cout << endl;
}
switch (SolverOptions::schedule) {
case SolverOptions::S_SEQ_RANDOM:
random_shuffle (links_.begin(), links_.end());
// no break
case SolverOptions::S_SEQ_FIXED:
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->setNextMessage (getFactor2VarMsg (links_[i]));
links_[i]->updateMessage();
}
break;
case SolverOptions::S_PARALLEL:
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->setNextMessage (getFactor2VarMsg (links_[i]));
}
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->updateMessage();
}
break;
case SolverOptions::S_MAX_RESIDUAL:
maxResidualSchedule();
break;
}
}
if (DL >= 2) {
cout << endl;
if (nIter_ < SolverOptions::maxIter) {
cout << "Loopy Sum-Product converged in " ;
cout << nIter_ << " iterations" << endl;
} else {
cout << "The maximum number of iterations was hit, terminating..." ;
cout << endl;
}
}
}
ParamSet
SPSolver::getPosterioriOf (Vid vid) const
{
assert (fg_->getFgVarNode (vid));
FgVarNode* var = fg_->getFgVarNode (vid);
ParamSet probs;
if (var->hasEvidence()) {
probs.resize (var->getDomainSize(), 0.0);
probs[var->getEvidence()] = 1.0;
} else {
probs.resize (var->getDomainSize(), 1.0);
CLinkSet links = varsI_[var->getIndex()]->getLinks();
for (unsigned i = 0; i < links.size(); i++) {
CParamSet msg = links[i]->getMessage();
for (unsigned j = 0; j < msg.size(); j++) {
probs[j] *= msg[j];
}
}
Util::normalize (probs);
}
return probs;
}
ParamSet
SPSolver::getJointDistributionOf (const VidSet& jointVids)
{
FgVarSet jointVars;
unsigned dsize = 1;
for (unsigned i = 0; i < jointVids.size(); i++) {
FgVarNode* varNode = fg_->getFgVarNode (jointVids[i]);
dsize *= varNode->getDomainSize();
jointVars.push_back (varNode);
}
unsigned maxVid = std::numeric_limits<unsigned>::max();
FgVarNode* junctionVar = new FgVarNode (maxVid, dsize);
FgVarSet factorVars = { junctionVar };
for (unsigned i = 0; i < jointVars.size(); i++) {
factorVars.push_back (jointVars[i]);
}
unsigned nParams = dsize * dsize;
ParamSet params (nParams);
for (unsigned i = 0; i < nParams; i++) {
unsigned row = i / dsize;
unsigned col = i % dsize;
if (row == col) {
params[i] = 1;
} else {
params[i] = 0;
}
}
Distribution* dist = new Distribution (params, maxVid);
Factor* newFactor = new Factor (factorVars, dist);
fg_->addVariable (junctionVar);
fg_->addFactor (newFactor);
runSolver();
ParamSet results = getPosterioriOf (maxVid);
deleteJunction (newFactor, junctionVar);
return results;
}
void
SPSolver::initializeSolver (void)
{
fg_->setIndexes();
CFgVarSet vars = fg_->getFgVarNodes();
for (unsigned i = 0; i < varsI_.size(); i++) {
delete varsI_[i];
}
varsI_.reserve (vars.size());
for (unsigned i = 0; i < vars.size(); i++) {
varsI_.push_back (new SPNodeInfo());
}
CFactorSet factors = fg_->getFactors();
for (unsigned i = 0; i < factorsI_.size(); i++) {
delete factorsI_[i];
}
factorsI_.reserve (factors.size());
for (unsigned i = 0; i < factors.size(); i++) {
factorsI_.push_back (new SPNodeInfo());
}
for (unsigned i = 0; i < links_.size(); i++) {
delete links_[i];
}
createLinks();
for (unsigned i = 0; i < links_.size(); i++) {
Factor* source = links_[i]->getFactor();
FgVarNode* dest = links_[i]->getVariable();
varsI_[dest->getIndex()]->addLink (links_[i]);
factorsI_[source->getIndex()]->addLink (links_[i]);
}
}
void
SPSolver::createLinks (void)
{
CFactorSet factors = fg_->getFactors();
for (unsigned i = 0; i < factors.size(); i++) {
CFgVarSet neighbors = factors[i]->getFgVarNodes();
for (unsigned j = 0; j < neighbors.size(); j++) {
links_.push_back (new Link (factors[i], neighbors[j]));
}
}
}
void
SPSolver::deleteJunction (Factor* f, FgVarNode* v)
{
fg_->removeFactor (f);
f->freeDistribution();
delete f;
fg_->removeVariable (v);
delete v;
}
bool
SPSolver::converged (void)
{
// this can happen if the graph is fully disconnected
if (links_.size() == 0) {
return true;
}
if (nIter_ == 0 || nIter_ == 1) {
return false;
}
bool converged = true;
if (SolverOptions::schedule == SolverOptions::S_MAX_RESIDUAL) {
Param maxResidual = (*(sortedOrder_.begin()))->getResidual();
if (maxResidual < SolverOptions::accuracy) {
converged = true;
} else {
converged = false;
}
} else {
for (unsigned i = 0; i < links_.size(); i++) {
double residual = links_[i]->getResidual();
if (DL >= 2) {
cout << links_[i]->toString() + " residual = " << residual << endl;
}
if (residual > SolverOptions::accuracy) {
converged = false;
if (DL == 0) break;
}
}
}
return converged;
}
void
SPSolver::maxResidualSchedule (void)
{
if (nIter_ == 1) {
for (unsigned i = 0; i < links_.size(); i++) {
links_[i]->setNextMessage (getFactor2VarMsg (links_[i]));
SortedOrder::iterator it = sortedOrder_.insert (links_[i]);
linkMap_.insert (make_pair (links_[i], it));
if (DL >= 2 && DL < 5) {
cout << "calculating " << links_[i]->toString() << endl;
}
}
return;
}
for (unsigned c = 0; c < links_.size(); c++) {
if (DL >= 2) {
cout << endl << "current residuals:" << endl;
for (SortedOrder::iterator it = sortedOrder_.begin();
it != sortedOrder_.end(); it ++) {
cout << " " << setw (30) << left << (*it)->toString();
cout << "residual = " << (*it)->getResidual() << endl;
}
}
SortedOrder::iterator it = sortedOrder_.begin();
Link* link = *it;
if (DL >= 2) {
cout << "updating " << (*sortedOrder_.begin())->toString() << endl;
}
if (link->getResidual() < SolverOptions::accuracy) {
return;
}
link->updateMessage();
link->clearResidual();
sortedOrder_.erase (it);
linkMap_.find (link)->second = sortedOrder_.insert (link);
// update the messages that depend on message source --> destin
CFactorSet factorNeighbors = link->getVariable()->getFactors();
for (unsigned i = 0; i < factorNeighbors.size(); i++) {
if (factorNeighbors[i] != link->getFactor()) {
CLinkSet links = factorsI_[factorNeighbors[i]->getIndex()]->getLinks();
for (unsigned j = 0; j < links.size(); j++) {
if (links[j]->getVariable() != link->getVariable()) {
if (DL >= 2 && DL < 5) {
cout << " calculating " << links[j]->toString() << endl;
}
links[j]->setNextMessage (getFactor2VarMsg (links[j]));
LinkMap::iterator iter = linkMap_.find (links[j]);
sortedOrder_.erase (iter->second);
iter->second = sortedOrder_.insert (links[j]);
}
}
}
}
}
}
ParamSet
SPSolver::getFactor2VarMsg (const Link* link) const
{
const Factor* src = link->getFactor();
const FgVarNode* dest = link->getVariable();
CFgVarSet neighbors = src->getFgVarNodes();
CLinkSet links = factorsI_[src->getIndex()]->getLinks();
// calculate the product of messages that were sent
// to factor `src', except from var `dest'
Factor result (*src);
Factor temp;
if (DL >= 5) {
cout << "calculating " ;
cout << src->getLabel() << " --> " << dest->getLabel();
cout << endl;
}
for (unsigned i = 0; i < neighbors.size(); i++) {
if (links[i]->getVariable() != dest) {
if (DL >= 5) {
cout << " message from " << links[i]->getVariable()->getLabel();
cout << ": " ;
ParamSet p = getVar2FactorMsg (links[i]);
cout << endl;
Factor temp2 (links[i]->getVariable(), p);
temp.multiplyByFactor (temp2);
temp2.freeDistribution();
} else {
Factor temp2 (links[i]->getVariable(), getVar2FactorMsg (links[i]));
temp.multiplyByFactor (temp2);
temp2.freeDistribution();
}
}
}
if (links.size() >= 2) {
result.multiplyByFactor (temp, &(src->getCptEntries()));
if (DL >= 5) {
cout << " message product: " ;
cout << Util::parametersToString (temp.getParameters()) << endl;
cout << " factor product: " ;
cout << Util::parametersToString (src->getParameters());
cout << " x " ;
cout << Util::parametersToString (temp.getParameters());
cout << " = " ;
cout << Util::parametersToString (result.getParameters()) << endl;
}
temp.freeDistribution();
}
for (unsigned i = 0; i < links.size(); i++) {
if (links[i]->getVariable() != dest) {
result.removeVariable (links[i]->getVariable());
}
}
if (DL >= 5) {
cout << " final message: " ;
cout << Util::parametersToString (result.getParameters()) << endl << endl;
}
ParamSet msg = result.getParameters();
result.freeDistribution();
return msg;
}
ParamSet
SPSolver::getVar2FactorMsg (const Link* link) const
{
const FgVarNode* src = link->getVariable();
const Factor* dest = link->getFactor();
ParamSet msg;
if (src->hasEvidence()) {
msg.resize (src->getDomainSize(), 0.0);
msg[src->getEvidence()] = 1.0;
if (DL >= 5) {
cout << Util::parametersToString (msg);
}
} else {
msg.resize (src->getDomainSize(), 1.0);
}
if (DL >= 5) {
cout << Util::parametersToString (msg);
}
CLinkSet links = varsI_[src->getIndex()]->getLinks();
for (unsigned i = 0; i < links.size(); i++) {
if (links[i]->getFactor() != dest) {
CParamSet msgFromFactor = links[i]->getMessage();
for (unsigned j = 0; j < msgFromFactor.size(); j++) {
msg[j] *= msgFromFactor[j];
}
if (DL >= 5) {
cout << " x " << Util::parametersToString (msgFromFactor);
}
}
}
if (DL >= 5) {
cout << " = " << Util::parametersToString (msg);
}
return msg;
}

130
packages/CLPBN/clpbn/bp/SPSolver.h
View File

@ -1,130 +0,0 @@
#ifndef BP_SP_SOLVER_H
#define BP_SP_SOLVER_H
#include <vector>
#include <set>
#include "Solver.h"
#include "FgVarNode.h"
#include "Factor.h"
using namespace std;
class FactorGraph;
class SPSolver;
class Link
{
public:
Link (Factor* f, FgVarNode* v)
{
factor_ = f;
var_ = v;
currMsg_.resize (v->getDomainSize(), 1);
nextMsg_.resize (v->getDomainSize(), 1);
msgSended_ = false;
residual_ = 0.0;
}
void setMessage (ParamSet msg)
{
Util::normalize (msg);
residual_ = Util::getMaxNorm (currMsg_, msg);
currMsg_ = msg;
}
void setNextMessage (CParamSet msg)
{
nextMsg_ = msg;
Util::normalize (nextMsg_);
residual_ = Util::getMaxNorm (currMsg_, nextMsg_);
}
void updateMessage (void)
{
currMsg_ = nextMsg_;
msgSended_ = true;
}
string toString (void) const
{
stringstream ss;
ss << factor_->getLabel();
ss << " -- " ;
ss << var_->getLabel();
return ss.str();
}
Factor* getFactor (void) const { return factor_; }
FgVarNode* getVariable (void) const { return var_; }
CParamSet getMessage (void) const { return currMsg_; }
bool messageWasSended (void) const { return msgSended_; }
double getResidual (void) const { return residual_; }
void clearResidual (void) { residual_ = 0.0; }
private:
Factor* factor_;
FgVarNode* var_;
ParamSet currMsg_;
ParamSet nextMsg_;
bool msgSended_;
double residual_;
};
class SPNodeInfo
{
public:
void addLink (Link* link) { links_.push_back (link); }
CLinkSet getLinks (void) { return links_; }
private:
LinkSet links_;
};
class SPSolver : public Solver
{
public:
SPSolver (FactorGraph&);
virtual ~SPSolver (void);
void runSolver (void);
virtual ParamSet getPosterioriOf (Vid) const;
ParamSet getJointDistributionOf (CVidSet);
protected:
virtual void initializeSolver (void);
void runTreeSolver (void);
bool readyToSendMessage (const Link*) const;
virtual void createLinks (void);
virtual void deleteJunction (Factor*, FgVarNode*);
bool converged (void);
virtual void maxResidualSchedule (void);
virtual ParamSet getFactor2VarMsg (const Link*) const;
virtual ParamSet getVar2FactorMsg (const Link*) const;
struct CompareResidual {
inline bool operator() (const Link* link1, const Link* link2)
{
return link1->getResidual() > link2->getResidual();
}
};
FactorGraph* fg_;
LinkSet links_;
vector<SPNodeInfo*> varsI_;
vector<SPNodeInfo*> factorsI_;
unsigned nIter_;
typedef multiset<Link*, CompareResidual> SortedOrder;
SortedOrder sortedOrder_;
typedef map<Link*, SortedOrder::iterator> LinkMap;
LinkMap linkMap_;
};
#endif // BP_SP_SOLVER_H

172
packages/CLPBN/clpbn/bp/Variable.h
View File

@ -1,172 +0,0 @@
#ifndef BP_VARIABLE_H
#define BP_VARIABLE_H
#include <algorithm>
#include <sstream>
#include "Shared.h"
using namespace std;
class Variable
{
public:
Variable (const Variable* v)
{
vid_ = v->getVarId();
dsize_ = v->getDomainSize();
if (v->hasDomain()) {
domain_ = v->getDomain();
dsize_ = domain_.size();
} else {
dsize_ = v->getDomainSize();
}
evidence_ = v->getEvidence();
if (v->hasLabel()) {
label_ = new string (v->getLabel());
} else {
label_ = 0;
}
}
Variable (Vid vid)
{
this->vid_ = vid;
this->dsize_ = 0;
this->evidence_ = NO_EVIDENCE;
this->label_ = 0;
}
Variable (Vid vid, unsigned dsize, int evidence = NO_EVIDENCE,
const string& lbl = string())
{
assert (dsize != 0);
assert (evidence < (int)dsize);
this->vid_ = vid;
this->dsize_ = dsize;
this->evidence_ = evidence;
if (!lbl.empty()) {
this->label_ = new string (lbl);
} else {
this->label_ = 0;
}
}
Variable (Vid vid, const Domain& domain, int evidence = NO_EVIDENCE,
const string& lbl = string())
{
assert (!domain.empty());
assert (evidence < (int)domain.size());
this->vid_ = vid;
this->dsize_ = domain.size();
this->domain_ = domain;
this->evidence_ = evidence;
if (!lbl.empty()) {
this->label_ = new string (lbl);
} else {
this->label_ = 0;
}
}
~Variable (void)
{
delete label_;
}
unsigned getVarId (void) const { return vid_; }
unsigned getIndex (void) const { return index_; }
void setIndex (unsigned idx) { index_ = idx; }
unsigned getDomainSize (void) const { return dsize_; }
bool hasEvidence (void) const { return evidence_ != NO_EVIDENCE; }
int getEvidence (void) const { return evidence_; }
bool hasDomain (void) const { return !domain_.empty(); }
bool hasLabel (void) const { return label_ != 0; }
bool isValidStateIndex (int index)
{
return index >= 0 && index < (int)dsize_;
}
bool isValidState (const string& state)
{
return find (domain_.begin(), domain_.end(), state) != domain_.end();
}
Domain getDomain (void) const
{
assert (dsize_ != 0);
if (domain_.size() == 0) {
Domain d;
for (unsigned i = 0; i < dsize_; i++) {
stringstream ss;
ss << "x" << i ;
d.push_back (ss.str());
}
return d;
} else {
return domain_;
}
}
void setDomainSize (unsigned dsize)
{
assert (dsize != 0);
dsize_ = dsize;
}
void setDomain (const Domain& domain)
{
assert (!domain.empty());
domain_ = domain;
dsize_ = domain.size();
}
void setEvidence (int ev)
{
assert (ev < dsize_);
evidence_ = ev;
}
void setEvidence (const string& ev)
{
assert (isValidState (ev));
for (unsigned i = 0; i < domain_.size(); i++) {
if (domain_[i] == ev) {
evidence_ = i;
}
}
}
void setLabel (const string& label)
{
label_ = new string (label);
}
string getLabel (void) const
{
if (label_ == 0) {
stringstream ss;
ss << "v" << vid_;
return ss.str();
} else {
return *label_;
}
}
private:
DISALLOW_COPY_AND_ASSIGN (Variable);
Vid vid_;
unsigned dsize_;
int evidence_;
Domain domain_;
string* label_;
unsigned index_;
};
#endif // BP_VARIABLE_H

147
packages/CLPBN/clpbn/bp/callgrind.h
View File

@ -1,147 +0,0 @@
/*
----------------------------------------------------------------
Notice that the following BSD-style license applies to this one
file (callgrind.h) only. The rest of Valgrind is licensed under the
terms of the GNU General Public License, version 2, unless
otherwise indicated. See the COPYING file in the source
distribution for details.
----------------------------------------------------------------
This file is part of callgrind, a valgrind tool for cache simulation
and call tree tracing.
Copyright (C) 2003-2010 Josef Weidendorfer. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. The origin of this software must not be misrepresented; you must
not claim that you wrote the original software. If you use this
software in a product, an acknowledgment in the product
documentation would be appreciated but is not required.
3. Altered source versions must be plainly marked as such, and must
not be misrepresented as being the original software.
4. The name of the author may not be used to endorse or promote
products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
----------------------------------------------------------------
Notice that the above BSD-style license applies to this one file
(callgrind.h) only. The entire rest of Valgrind is licensed under
the terms of the GNU General Public License, version 2. See the
COPYING file in the source distribution for details.
----------------------------------------------------------------
*/
#ifndef __CALLGRIND_H
#define __CALLGRIND_H
#include "valgrind.h"
/* !! ABIWARNING !! ABIWARNING !! ABIWARNING !! ABIWARNING !!
This enum comprises an ABI exported by Valgrind to programs
which use client requests. DO NOT CHANGE THE ORDER OF THESE
ENTRIES, NOR DELETE ANY -- add new ones at the end.
The identification ('C','T') for Callgrind has historical
reasons: it was called "Calltree" before. Besides, ('C','G') would
clash with cachegrind.
*/
typedef
enum {
VG_USERREQ__DUMP_STATS = VG_USERREQ_TOOL_BASE('C','T'),
VG_USERREQ__ZERO_STATS,
VG_USERREQ__TOGGLE_COLLECT,
VG_USERREQ__DUMP_STATS_AT,
VG_USERREQ__START_INSTRUMENTATION,
VG_USERREQ__STOP_INSTRUMENTATION
} Vg_CallgrindClientRequest;
/* Dump current state of cost centers, and zero them afterwards */
#define CALLGRIND_DUMP_STATS \
{unsigned int _qzz_res; \
VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0, \
VG_USERREQ__DUMP_STATS, \
0, 0, 0, 0, 0); \
}
/* Dump current state of cost centers, and zero them afterwards.
The argument is appended to a string stating the reason which triggered
the dump. This string is written as a description field into the
profile data dump. */
#define CALLGRIND_DUMP_STATS_AT(pos_str) \
{unsigned int _qzz_res; \
VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0, \
VG_USERREQ__DUMP_STATS_AT, \
pos_str, 0, 0, 0, 0); \
}
/* Zero cost centers */
#define CALLGRIND_ZERO_STATS \
{unsigned int _qzz_res; \
VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0, \
VG_USERREQ__ZERO_STATS, \
0, 0, 0, 0, 0); \
}
/* Toggles collection state.
The collection state specifies whether the happening of events
should be noted or if they are to be ignored. Events are noted
by increment of counters in a cost center */
#define CALLGRIND_TOGGLE_COLLECT \
{unsigned int _qzz_res; \
VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0, \
VG_USERREQ__TOGGLE_COLLECT, \
0, 0, 0, 0, 0); \
}
/* Start full callgrind instrumentation if not already switched on.
When cache simulation is done, it will flush the simulated cache;
this will lead to an artifical cache warmup phase afterwards with
cache misses which would not have happened in reality. */
#define CALLGRIND_START_INSTRUMENTATION \
{unsigned int _qzz_res; \
VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0, \
VG_USERREQ__START_INSTRUMENTATION, \
0, 0, 0, 0, 0); \
}
/* Stop full callgrind instrumentation if not already switched off.
This flushes Valgrinds translation cache, and does no additional
instrumentation afterwards, which effectivly will run at the same
speed as the "none" tool (ie. at minimal slowdown).
Use this to bypass Callgrind aggregation for uninteresting code parts.
To start Callgrind in this mode to ignore the setup phase, use
the option "--instr-atstart=no". */
#define CALLGRIND_STOP_INSTRUMENTATION \
{unsigned int _qzz_res; \
VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0, \
VG_USERREQ__STOP_INSTRUMENTATION, \
0, 0, 0, 0, 0); \
}
#endif /* __CALLGRIND_H */

14
packages/CLPBN/clpbn/bp/examples/cbp_example.uai Normal file
View File

@ -0,0 +1,14 @@
MARKOV
3
2 2 2
2
2 0 1
2 2 1
4
1.2 1.4 2.0 0.4
4
1.2 1.4 2.0 0.4

53
packages/CLPBN/clpbn/bp/examples/loop.yap
View File

@ -1,53 +0,0 @@
:- use_module(library(clpbn)).
:- set_clpbn_flag(solver, bp).
%
% A E
% / \ /
% / \ /
% B C
% \ /
% \ /
% D
%
a(A) :-
a_table(ADist),
{ A = a with p([a1, a2], ADist) }.
b(B) :-
a(A),
b_table(BDist),
{ B = b with p([b1, b2], BDist, [A]) }.
c(C) :-
a(A),
c_table(CDist),
{ C = c with p([c1, c2], CDist, [A]) }.
d(D) :-
b(B),
c(C),
d_table(DDist),
{ D = d with p([d1, d2], DDist, [B, C]) }.
e(E) :-
e_table(EDist),
{ E = e with p([e1, e2], EDist) }.
a_table([0.005, 0.995]).
b_table([0.02, 0.97,
0.88, 0.03]).
c_table([0.55, 0.94,
0.45, 0.06]).
d_table([0.192, 0.98, 0.33, 0.013,
0.908, 0.02, 0.77, 0.987]).
e_table([0.055, 0.945]).

60
packages/CLPBN/clpbn/bp/examples/town/run_town_bp_compress.sh Normal file
View File

@ -0,0 +1,60 @@
#!/bin/bash
cp ~/bin/yap ~/bin/town_comp
YAP=~/bin/town_comp
#OUT_FILE_NAME=results`date "+ %H:%M:%S %d-%m-%Y"`.log
OUT_FILE_NAME=bp_compress.log
rm -f $OUT_FILE_NAME
rm -f ignore.$OUT_FILE_NAME
function run_solver
{
if [ $2 = bp ]
then
extra_flag1=clpbn_bp:set_solver_parameter\(run_mode,$4\)
extra_flag2=clpbn_bp:set_solver_parameter\(schedule,$5\)
extra_flag3=clpbn_bp:set_solver_parameter\(always_loopy_solver,$6\)
else
extra_flag1=true
extra_flag2=true
extra_flag3=true
fi
/usr/bin/time -o $OUT_FILE_NAME -a -f "real:%E\tuser:%U\tsys:%S" $YAP << EOF >> $OUT_FILE_NAME 2>> ignore.$OUT_FILE_NAME
[$1].
clpbn:set_clpbn_flag(solver,$2),
clpbn_bp:use_log_space,
$extra_flag1, $extra_flag2, $extra_flag3,
run_query(_R),
open("$OUT_FILE_NAME", 'append',S),
format(S, '$3: ~15+ ',[]),
close(S).
EOF
}
function run_all_graphs
{
echo "*******************************************************************" >> "$OUT_FILE_NAME"
echo "results for solver $2" >> $OUT_FILE_NAME
echo "*******************************************************************" >> "$OUT_FILE_NAME"
run_solver town_1000 $1 town_1000 $3 $4 $5
run_solver town_5000 $1 town_5000 $3 $4 $5
run_solver town_10000 $1 town_10000 $3 $4 $5
run_solver town_50000 $1 town_50000 $3 $4 $5
run_solver town_100000 $1 town_100000 $3 $4 $5
run_solver town_500000 $1 town_500000 $3 $4 $5
run_solver town_1000000 $1 town_1000000 $3 $4 $5
run_solver town_2500000 $1 town_2500000 $3 $4 $5
run_solver town_5000000 $1 town_5000000 $3 $4 $5
run_solver town_7500000 $1 town_7500000 $3 $4 $5
run_solver town_10000000 $1 town_10000000 $3 $4 $5
}
run_solver town_10000 "bp(compress,seq_fixed)" town_10000 compress seq_fixed true
exit
##########
run_all_graphs bp "bp(compress,seq_fixed) " compress seq_fixed true

51
packages/CLPBN/clpbn/bp/examples/town/run_town_bp_convert.sh Normal file
View File

@ -0,0 +1,51 @@
#!/bin/bash
YAP=~/bin/town_conv
#OUT_FILE_NAME=results`date "+ %H:%M:%S %d-%m-%Y"`.log
OUT_FILE_NAME=bp_convert.log
rm -f $OUT_FILE_NAME
rm -f ignore.$OUT_FILE_NAME
function run_solver
{
if [ $2 = bp ]
then
extra_flag1=clpbn_bp:set_solver_parameter\(run_mode,$4\)
extra_flag2=clpbn_bp:set_solver_parameter\(schedule,$5\)
extra_flag3=clpbn_bp:set_solver_parameter\(always_loopy_solver,$6\)
else
extra_flag1=true
extra_flag2=true
extra_flag3=true
fi
/usr/bin/time -o $OUT_FILE_NAME -a -f "real:%E\tuser:%U\tsys:%S" $YAP << EOF >> $OUT_FILE_NAME 2>> ignore.$OUT_FILE_NAME
[$1].
clpbn:set_clpbn_flag(solver,$2),
clpbn_bp:use_log_space,
$extra_flag1, $extra_flag2, $extra_flag3,
run_query(_R),
open("$OUT_FILE_NAME", 'append',S),
format(S, '$3: ~15+ ',[]),
close(S).
EOF
}
function run_all_graphs
{
echo "*******************************************************************" >> "$OUT_FILE_NAME"
echo "results for solver $2" >> $OUT_FILE_NAME
echo "*******************************************************************" >> "$OUT_FILE_NAME"
run_solver town_1000 $1 town_1000 $3 $4 $5
run_solver town_5000 $1 town_5000 $3 $4 $5
run_solver town_10000 $1 town_10000 $3 $4 $5
run_solver town_50000 $1 town_50000 $3 $4 $5
run_solver town_100000 $1 town_100000 $3 $4 $5
run_solver town_500000 $1 town_500000 $3 $4 $5
run_solver town_1000000 $1 town_1000000 $3 $4 $5
}
run_all_graphs bp "bp(convert,seq_fixed) " convert seq_fixed false

50
packages/CLPBN/clpbn/bp/examples/town/run_town_bp_normal.sh Normal file
View File

@ -0,0 +1,50 @@
#!/bin/bash
YAP=~/bin/town_norm
#OUT_FILE_NAME=results`date "+ %H:%M:%S %d-%m-%Y"`.log
OUT_FILE_NAME=bp_normal.log
rm -f $OUT_FILE_NAME
rm -f ignore.$OUT_FILE_NAME
function run_solver
{
if [ $2 = bp ]
then
extra_flag1=clpbn_bp:set_solver_parameter\(run_mode,$4\)
extra_flag2=clpbn_bp:set_solver_parameter\(schedule,$5\)
extra_flag3=clpbn_bp:set_solver_parameter\(always_loopy_solver,$6\)
else
extra_flag1=true
extra_flag2=true
extra_flag3=true
fi
/usr/bin/time -o $OUT_FILE_NAME -a -f "real:%E\tuser:%U\tsys:%S" $YAP << EOF >> $OUT_FILE_NAME 2>> ignore.$OUT_FILE_NAME
[$1].
clpbn:set_clpbn_flag(solver,$2),
clpbn_bp:use_log_space,
$extra_flag1, $extra_flag2, $extra_flag3,
run_query(_R),
open("$OUT_FILE_NAME", 'append',S),
format(S, '$3: ~15+ ',[]),
close(S).
EOF
}
function run_all_graphs
{
echo "*******************************************************************" >> "$OUT_FILE_NAME"
echo "results for solver $2" >> $OUT_FILE_NAME
echo "*******************************************************************" >> "$OUT_FILE_NAME"
run_solver town_1000 $1 town_1000 $3 $4 $5
run_solver town_5000 $1 town_5000 $3 $4 $5
run_solver town_10000 $1 town_10000 $3 $4 $5
run_solver town_50000 $1 town_50000 $3 $4 $5
run_solver town_100000 $1 town_100000 $3 $4 $5
run_solver town_500000 $1 town_500000 $3 $4 $5
run_solver town_1000000 $1 town_1000000 $3 $4 $5
}
run_all_graphs bp "bp(normal,seq_fixed) " normal seq_fixed false

51
packages/CLPBN/clpbn/bp/examples/town/run_town_gibbs_tests.sh Normal file
View File

@ -0,0 +1,51 @@
#!/bin/bash
YAP=~/bin/town_gibbs
#OUT_FILE_NAME=results`date "+ %H:%M:%S %d-%m-%Y"`.log
OUT_FILE_NAME=gibbs.log
rm -f $OUT_FILE_NAME
rm -f ignore.$OUT_FILE_NAME
function run_solver
{
if [ $2 = bp ]
then
extra_flag1=clpbn_bp:set_solver_parameter\(run_mode,$4\)
extra_flag2=clpbn_bp:set_solver_parameter\(schedule,$5\)
extra_flag3=clpbn_bp:set_solver_parameter\(always_loopy_solver,$6\)
else
extra_flag1=true
extra_flag2=true
extra_flag3=true
fi
/usr/bin/time -o $OUT_FILE_NAME -a -f "real:%E\tuser:%U\tsys:%S" $YAP << EOF >> $OUT_FILE_NAME 2>> ignore.$OUT_FILE_NAME
[$1].
clpbn:set_clpbn_flag(solver,$2),
clpbn_bp:use_log_space,
$extra_flag1, $extra_flag2, $extra_flag3,
run_query(_R),
open("$OUT_FILE_NAME", 'append',S),
format(S, '$3: ~15+ ',[]),
close(S).
EOF
}
function run_all_graphs
{
echo "*******************************************************************" >> "$OUT_FILE_NAME"
echo "results for solver $2" >> $OUT_FILE_NAME
echo "*******************************************************************" >> "$OUT_FILE_NAME"
run_solver town_1000 $1 town_1000 $3 $4 $5
run_solver town_5000 $1 town_5000 $3 $4 $5
run_solver town_10000 $1 town_10000 $3 $4 $5
run_solver town_50000 $1 town_50000 $3 $4 $5
run_solver town_100000 $1 town_100000 $3 $4 $5
run_solver town_500000 $1 town_500000 $3 $4 $5
run_solver town_1000000 $1 town_1000000 $3 $4 $5
}
run_all_graphs gibbs "gibbs "

51
packages/CLPBN/clpbn/bp/examples/town/run_town_jt_tests.sh Normal file
View File

@ -0,0 +1,51 @@
#!/bin/bash
YAP=~/bin/town_jt
#OUT_FILE_NAME=results`date "+ %H:%M:%S %d-%m-%Y"`.log
OUT_FILE_NAME=jt.log
rm -f $OUT_FILE_NAME
rm -f ignore.$OUT_FILE_NAME
function run_solver
{
if [ $2 = bp ]
then
extra_flag1=clpbn_bp:set_solver_parameter\(run_mode,$4\)
extra_flag2=clpbn_bp:set_solver_parameter\(schedule,$5\)
extra_flag3=clpbn_bp:set_solver_parameter\(always_loopy_solver,$6\)
else
extra_flag1=true
extra_flag2=true
extra_flag3=true
fi
/usr/bin/time -o $OUT_FILE_NAME -a -f "real:%E\tuser:%U\tsys:%S" $YAP << EOF >> $OUT_FILE_NAME 2>> ignore.$OUT_FILE_NAME
[$1].
clpbn:set_clpbn_flag(solver,$2),
clpbn_bp:use_log_space,
$extra_flag1, $extra_flag2, $extra_flag3,
run_query(_R),
open("$OUT_FILE_NAME", 'append',S),
format(S, '$3: ~15+ ',[]),
close(S).
EOF
}
function run_all_graphs
{
echo "*******************************************************************" >> "$OUT_FILE_NAME"
echo "results for solver $2" >> $OUT_FILE_NAME
echo "*******************************************************************" >> "$OUT_FILE_NAME"
run_solver town_1000 $1 town_1000 $3 $4 $5
run_solver town_5000 $1 town_5000 $3 $4 $5
run_solver town_10000 $1 town_10000 $3 $4 $5
run_solver town_50000 $1 town_50000 $3 $4 $5
run_solver town_100000 $1 town_100000 $3 $4 $5
run_solver town_500000 $1 town_500000 $3 $4 $5
run_solver town_1000000 $1 town_1000000 $3 $4 $5
}
run_all_graphs jt "jt "

51
packages/CLPBN/clpbn/bp/examples/town/run_town_ve_tests.sh Normal file
View File

@ -0,0 +1,51 @@
#!/bin/bash
YAP=~/bin/town_ve
#OUT_FILE_NAME=results`date "+ %H:%M:%S %d-%m-%Y"`.log
OUT_FILE_NAME=ve.log
rm -f $OUT_FILE_NAME
rm -f ignore.$OUT_FILE_NAME
function run_solver
{
if [ $2 = bp ]
then
extra_flag1=clpbn_bp:set_solver_parameter\(run_mode,$4\)
extra_flag2=clpbn_bp:set_solver_parameter\(schedule,$5\)
extra_flag3=clpbn_bp:set_solver_parameter\(always_loopy_solver,$6\)
else
extra_flag1=true
extra_flag2=true
extra_flag3=true
fi
/usr/bin/time -o $OUT_FILE_NAME -a -f "real:%E\tuser:%U\tsys:%S" $YAP << EOF >> $OUT_FILE_NAME 2>> ignore.$OUT_FILE_NAME
[$1].
clpbn:set_clpbn_flag(solver,$2),
clpbn_bp:use_log_space,
$extra_flag1, $extra_flag2, $extra_flag3,
run_query(_R),
open("$OUT_FILE_NAME", 'append',S),
format(S, '$3: ~15+ ',[]),
close(S).
EOF
}
function run_all_graphs
{
echo "*******************************************************************" >> "$OUT_FILE_NAME"
echo "results for solver $2" >> $OUT_FILE_NAME
echo "*******************************************************************" >> "$OUT_FILE_NAME"
run_solver town_1000 $1 town_1000 $3 $4 $5
run_solver town_5000 $1 town_5000 $3 $4 $5
run_solver town_10000 $1 town_10000 $3 $4 $5
run_solver town_50000 $1 town_50000 $3 $4 $5
run_solver town_100000 $1 town_100000 $3 $4 $5
run_solver town_500000 $1 town_500000 $3 $4 $5
run_solver town_1000000 $1 town_1000000 $3 $4 $5
}
run_all_graphs ve "ve "

65
packages/CLPBN/clpbn/bp/examples/town/schema.yap Normal file
View File

@ -0,0 +1,65 @@
conservative_city(City, Cons) :-
cons_table(City, ConsDist),
{ Cons = conservative_city(City) with p([y,n], ConsDist) }.
gender(X, Gender) :-
gender_table(X, GenderDist),
{ Gender = gender(X) with p([m,f], GenderDist) }.
hair_color(X, Color) :-
lives(X, City),
conservative_city(City, Cons),
hair_color_table(X,ColorTable),
{ Color = hair_color(X) with
p([t,f], ColorTable,[Cons]) }.
car_color(X, Color) :-
hair_color(X, HColor),
car_color_table(X,CColorTable),
{ Color = car_color(X) with
p([t,f], CColorTable,[HColor]) }.
height(X, Height) :-
gender(X, Gender),
height_table(X,HeightTable),
{ Height = height(X) with
p([t,f], HeightTable,[Gender]) }.
shoe_size(X, Shoesize) :-
height(X, Height),
shoe_size_table(X,ShoesizeTable),
{ Shoesize = shoe_size(X) with
p([t,f], ShoesizeTable,[Height]) }.
guilty(X, Guilt) :-
guilty_table(X, GuiltDist),
{ Guilt = guilty(X) with p([y,n], GuiltDist) }.
descn(X, Descn) :-
car_color(X, Car),
hair_color(X, Hair),
height(X, Height),
guilty(X, Guilt),
descn_table(X, DescTable),
{ Descn = descn(X) with
p([t,f], DescTable,[Car,Hair,Height,Guilt]) }.
witness(City, Witness) :-
descn(joe, DescnJ),
descn(p2, Descn2),
wit_table(WitTable),
{ Witness = witness(City) with
p([t,f], WitTable,[DescnJ, Descn2]) }.
:- ensure_loaded(tables).

46
packages/CLPBN/clpbn/bp/examples/town/tables.yap Normal file
View File

@ -0,0 +1,46 @@
cons_table(amsterdam, [0.2, 0.8]) :- !.
cons_table(_, [0.8, 0.2]).
gender_table(_, [0.55, 0.44]).
hair_color_table(_,
/* conservative_city */
/* y n */
[ 0.05, 0.1,
0.95, 0.9 ]).
car_color_table(_,
/* t f */
[ 0.9, 0.2,
0.1, 0.8 ]).
height_table(_,
/* m f */
[ 0.6, 0.4,
0.4, 0.6 ]).
shoe_size_table(_,
/* t f */
[ 0.9, 0.1,
0.1, 0.9 ]).
guilty_table(_, [0.23, 0.77]).
descn_table(_,
/* color, hair, height, guilt */
/* ttttt tttf ttft ttff tfttt tftf tfft tfff ttttt fttf ftft ftff ffttt fftf ffft ffff */
[ 0.99, 0.5, 0.23, 0.88, 0.41, 0.3, 0.76, 0.87, 0.44, 0.43, 0.29, 0.72, 0.33, 0.91, 0.95, 0.92,
0.01, 0.5, 0.77, 0.12, 0.59, 0.7, 0.24, 0.13, 0.56, 0.57, 0.61, 0.28, 0.77, 0.09, 0.05, 0.08]).
wit_table([0.2, 0.45, 0.24, 0.34,
0.8, 0.55, 0.76, 0.66]).

59
packages/CLPBN/clpbn/bp/examples/town/town_generator.sh Normal file
View File

@ -0,0 +1,59 @@
#!/home/tgomes/bin/yap -L --
/*
Steps:
1. generate N facts lives(I, nyc), 0 <= I < N.
2. generate evidence on descn for N people, *** except for 1 ***
3. Run query ?- guilty(joe, Guilty), witness(joe, t), descn(2,t), descn(3, f), descn(4, f) ...
*/
:- initialization(main).
main :-
unix(argv([H])),
generate_town(H).
generate_town(N) :-
atomic_concat(['town_', N, '.yap'], FileName),
open(FileName, 'write', S),
write(S, ':- source.\n'),
write(S, ':- style_check(all).\n'),
write(S, ':- yap_flag(unknown,error).\n'),
write(S, ':- yap_flag(write_strings,on).\n'),
write(S, ':- use_module(library(clpbn)).\n'),
write(S, ':- set_clpbn_flag(solver, bp).\n'),
write(S, ':- [-schema].\n\n'),
write(S, 'lives(_joe, nyc).\n'),
atom_number(N, N2),
generate_people(S, N2, 2),
write(S, '\nrun_query(Guilty) :- \n'),
write(S, '\tguilty(joe, Guilty),\n'),
write(S, '\twitness(nyc, t),\n'),
write(S, '\trunall(X, ev(X)).\n\n\n'),
write(S, 'runall(G, Wrapper) :-\n'),
write(S, '\tfindall(G, Wrapper, L),\n'),
write(S, '\texecute_all(L).\n\n\n'),
write(S, 'execute_all([]).\n'),
write(S, 'execute_all(G.L) :-\n'),
write(S, '\tcall(G),\n'),
write(S, '\texecute_all(L).\n\n\n'),
generate_query(S, N2, 2),
close(S).
generate_people(_, N, Counting1) :- !.
generate_people(S, N, Counting) :-
format(S, 'lives(p~w, nyc).~n', [Counting]),
Counting1 is Counting + 1,
generate_people(S, N, Counting1).
generate_query(S, N, Counting) :-
Counting > N, !.
generate_query(S, N, Counting) :- !,
format(S, 'ev(descn(p~w, t)).~n', [Counting]),
Counting1 is Counting + 1,
generate_query(S, N, Counting1).

4536
packages/CLPBN/clpbn/bp/valgrind.h
View File

File diff suppressed because it is too large Load Diff