/** * @file BayesTree.cpp * @brief Bayes Tree is a tree of cliques of a Bayes Chain * @author Frank Dellaert */ #include #include // for operator += using namespace boost::assign; #include "BayesTree.h" #include "inference-inl.h" namespace gtsam { using namespace std; /* ************************************************************************* */ template BayesTree::Clique::Clique(const sharedConditional& conditional) { separator_ = conditional->parents(); this->push_back(conditional); } /* ************************************************************************* */ template Ordering BayesTree::Clique::keys() const { Ordering frontal_keys = this->ordering(), keys = separator_; keys.splice(keys.begin(),frontal_keys); return keys; } /* ************************************************************************* */ template void BayesTree::Clique::print(const string& s) const { cout << s; BOOST_FOREACH(const sharedConditional& conditional, this->conditionals_) cout << " " << conditional->key(); if (!separator_.empty()) { cout << " :"; BOOST_FOREACH(string key, separator_) cout << " " << key; } cout << endl; } /* ************************************************************************* */ template size_t BayesTree::Clique::treeSize() const { size_t size = 1; BOOST_FOREACH(shared_ptr child, children_) size += child->treeSize(); return size; } /* ************************************************************************* */ template void BayesTree::Clique::printTree(const string& indent) const { print(indent); BOOST_FOREACH(shared_ptr child, children_) child->printTree(indent+" "); } /* ************************************************************************* */ // The shortcut density is a conditional P(S|R) of the separator of this // clique on the root. We can compute it recursively from the parent shortcut // P(Sp|R) as \int P(Fp|Sp) P(Sp|R), where Fp are the frontal nodes in p // TODO, why do we actually return a shared pointer, why does eliminate? /* ************************************************************************* */ template template BayesNet BayesTree::Clique::shortcut(shared_ptr R) { // A first base case is when this clique or its parent is the root, // in which case we return an empty Bayes net. if (R.get()==this || parent_==R) { BayesNet empty; return empty; } // The parent clique has a Conditional for each frontal node in Fp // so we can obtain P(Fp|Sp) in factor graph form FactorGraph p_Fp_Sp(*parent_); // If not the base case, obtain the parent shortcut P(Sp|R) as factors FactorGraph p_Sp_R(parent_->shortcut(R)); // now combine P(Cp|R) = P(Fp|Sp) * P(Sp|R) FactorGraph p_Cp_R = combine(p_Fp_Sp, p_Sp_R); // Eliminate into a Bayes net with ordering designed to integrate out // any variables not in *our* separator. Variables to integrate out must be // eliminated first hence the desired ordering is [Cp\S S]. // However, an added wrinkle is that Cp might overlap with the root. // Keys corresponding to the root should not be added to the ordering at all. // Get the key list Cp=Fp+Sp, which will form the basis for the integrands Ordering integrands = parent_->keys(); // Start ordering with the separator Ordering ordering = separator_; // remove any variables in the root, after this integrands = Cp\R, ordering = S\R BOOST_FOREACH(string key, R->ordering()) { integrands.remove(key); ordering.remove(key); } // remove any variables in the separator, after this integrands = Cp\R\S BOOST_FOREACH(string key, separator_) integrands.remove(key); // form the ordering as [Cp\R\S S\R] BOOST_REVERSE_FOREACH(string key, integrands) ordering.push_front(key); // eliminate to get marginal BayesNet p_S_R = eliminate(p_Cp_R,ordering); // remove all integrands BOOST_FOREACH(string key, integrands) p_S_R.pop_front(); // return the parent shortcut P(Sp|R) return p_S_R; } /* ************************************************************************* */ // P(C) = \int_R P(F|S) P(S|R) P(R) // TODO: Maybe we should integrate given parent marginal P(Cp), // \int(Cp\S) P(F|S)P(S|Cp)P(Cp) // Because the root clique could be very big. /* ************************************************************************* */ template template FactorGraph BayesTree::Clique::marginal(shared_ptr R) { // If we are the root, just return this root if (R.get()==this) return *R; // Combine P(F|S), P(S|R), and P(R) BayesNet p_FSR = this->shortcut(R); p_FSR.push_front(*this); p_FSR.push_back(*R); // Find marginal on the keys we are interested in return marginalize(p_FSR,keys()); } /* ************************************************************************* */ // P(C1,C2) = \int_R P(F1|S1) P(S1|R) P(F2|S1) P(S2|R) P(R) /* ************************************************************************* */ template template FactorGraph BayesTree::Clique::joint(shared_ptr C2, shared_ptr R) { // For now, assume neither is the root // Combine P(F1|S1), P(S1|R), P(F2|S2), P(S2|R), and P(R) sharedBayesNet bn(new BayesNet); if (!isRoot()) bn->push_back(*this); // P(F1|S1) if (!isRoot()) bn->push_back(shortcut(R)); // P(S1|R) if (!C2->isRoot()) bn->push_back(*C2); // P(F2|S2) if (!C2->isRoot()) bn->push_back(C2->shortcut(R)); // P(S2|R) bn->push_back(*R); // P(R) // Find the keys of both C1 and C2 Ordering keys12 = keys(); BOOST_FOREACH(string key,C2->keys()) keys12.push_back(key); keys12.unique(); // Calculate the marginal return marginalize(*bn,keys12); } /* ************************************************************************* */ template typename BayesTree::sharedClique BayesTree::addClique (const sharedConditional& conditional, sharedClique parent_clique) { sharedClique new_clique(new Clique(conditional)); nodes_.insert(make_pair(conditional->key(), new_clique)); if (parent_clique != NULL) { new_clique->parent_ = parent_clique; parent_clique->children_.push_back(new_clique); } return new_clique; } /* ************************************************************************* */ template void BayesTree::removeClique(sharedClique clique) { if (clique->isRoot()) root_.reset(); else // detach clique from parent clique->parent_->children_.remove(clique); // orphan my children BOOST_FOREACH(sharedClique child, clique->children_) child->parent_.reset(); BOOST_FOREACH(string key, clique->ordering()) nodes_.erase(key); } /* ************************************************************************* */ template BayesTree::BayesTree() { } /* ************************************************************************* */ template BayesTree::BayesTree(const BayesNet& bayesNet) { typename BayesNet::const_reverse_iterator rit; for ( rit=bayesNet.rbegin(); rit != bayesNet.rend(); ++rit ) insert(*rit); } /* ************************************************************************* */ template void BayesTree::print(const string& s) const { cout << s << ": size == " << size() << endl; if (nodes_.empty()) return; root_->printTree(""); } /* ************************************************************************* */ // binary predicate to test equality of a pair for use in equals template bool check_pair( const pair::sharedClique >& v1, const pair::sharedClique >& v2 ) { return v1.first == v2.first && v1.second->equals(*(v2.second)); } /* ************************************************************************* */ template bool BayesTree::equals(const BayesTree& other, double tol) const { return size()==other.size() && equal(nodes_.begin(),nodes_.end(),other.nodes_.begin(),check_pair); } /* ************************************************************************* */ template void BayesTree::insert(const sharedConditional& conditional) { // get key and parents string key = conditional->key(); list parents = conditional->parents(); // if no parents, start a new root clique if (parents.empty()) { root_ = addClique(conditional); return; } // otherwise, find the parent clique string parent = parents.front(); sharedClique parent_clique = (*this)[parent]; // if the parents and parent clique have the same size, add to parent clique if (parent_clique->size() == parents.size()) { nodes_.insert(make_pair(key, parent_clique)); parent_clique->push_front(conditional); return; } // otherwise, start a new clique and add it to the tree addClique(conditional,parent_clique); } /* ************************************************************************* */ // First finds clique marginal then marginalizes that /* ************************************************************************* */ template template FactorGraph BayesTree::marginal(const string& key) const { // get clique containing key sharedClique clique = (*this)[key]; // calculate or retrieve its marginal FactorGraph cliqueMarginal = clique->marginal(root_); // create an ordering where only the requested key is not eliminated Ordering ord = clique->keys(); ord.remove(key); // partially eliminate, remaining factor graph is requested marginal eliminate(cliqueMarginal,ord); return cliqueMarginal; } /* ************************************************************************* */ template template BayesNet BayesTree::marginalBayesNet(const string& key) const { // calculate marginal as a factor graph FactorGraph fg = this->marginal(key); // eliminate further to Bayes net return eliminate(fg,Ordering(key)); } /* ************************************************************************* */ // Find two cliques, their joint, then marginalizes /* ************************************************************************* */ template template FactorGraph BayesTree::joint(const string& key1, const string& key2) const { // get clique C1 and C2 sharedClique C1 = (*this)[key1], C2 = (*this)[key2]; // calculate joint FactorGraph p_C1C2 = C1->joint(C2,root_); // create an ordering where both requested keys are not eliminated Ordering ord = p_C1C2.keys(); ord.remove(key1); ord.remove(key2); // partially eliminate, remaining factor graph is requested joint // TODO, make eliminate functional eliminate(p_C1C2,ord); return p_C1C2; } /* ************************************************************************* */ template template BayesNet BayesTree::jointBayesNet(const string& key1, const string& key2) const { // calculate marginal as a factor graph FactorGraph fg = this->joint(key1,key2); // eliminate further to Bayes net Ordering ordering; ordering += key1, key2; return eliminate(fg,ordering); } /* ************************************************************************* */ template template pair, list::sharedClique> > BayesTree::removePath(sharedClique clique) { FactorGraph factors; list orphans; // base case is NULL, if so we do nothing and return empties above if (clique!=NULL) { // remove path above me boost::tie(factors,orphans) = removePath(clique->parent_); // add children to list of orphans orphans.insert(orphans.begin(), clique->children_.begin(), clique->children_.end()); // remove me and add my factors removeClique(clique); factors.push_back(*clique); } return make_pair(factors,orphans); } /* ************************************************************************* */ } /// namespace gtsam