| 1 | /* ---------------------------------------------------------------------------- |
| 2 | |
| 3 | * GTSAM Copyright 2010-2020, Georgia Tech Research Corporation, |
| 4 | * Atlanta, Georgia 30332-0415 |
| 5 | * All Rights Reserved |
| 6 | * Authors: Frank Dellaert, et al. (see THANKS for the full author list) |
| 7 | |
| 8 | * See LICENSE for the license information |
| 9 | |
| 10 | * -------------------------------------------------------------------------- */ |
| 11 | |
| 12 | /** |
| 13 | * @file DiscreteBayesNetExample.cpp |
| 14 | * @brief Hidden Markov Model example, discrete. |
| 15 | * @author Frank Dellaert |
| 16 | * @date July 12, 2020 |
| 17 | */ |
| 18 | |
| 19 | #include <gtsam/discrete/DiscreteFactorGraph.h> |
| 20 | #include <gtsam/discrete/DiscreteMarginals.h> |
| 21 | #include <gtsam/inference/BayesNet.h> |
| 22 | |
| 23 | #include <iomanip> |
| 24 | #include <sstream> |
| 25 | |
| 26 | using namespace std; |
| 27 | using namespace gtsam; |
| 28 | |
| 29 | int main(int argc, char **argv) { |
| 30 | const int nrNodes = 4; |
| 31 | const size_t nrStates = 3; |
| 32 | |
| 33 | // Define variables as well as ordering |
| 34 | Ordering ordering; |
| 35 | vector<DiscreteKey> keys; |
| 36 | for (int k = 0; k < nrNodes; k++) { |
| 37 | DiscreteKey key_i(k, nrStates); |
| 38 | keys.push_back(x: key_i); |
| 39 | ordering.emplace_back(args&: k); |
| 40 | } |
| 41 | |
| 42 | // Create HMM as a DiscreteBayesNet |
| 43 | DiscreteBayesNet hmm; |
| 44 | |
| 45 | // Define backbone |
| 46 | const string transition = "8/1/1 1/8/1 1/1/8" ; |
| 47 | for (int k = 1; k < nrNodes; k++) { |
| 48 | hmm.add(args&: keys[k] | keys[k - 1] = transition); |
| 49 | } |
| 50 | |
| 51 | // Add some measurements, not needed for all time steps! |
| 52 | hmm.add(args: keys[0] % "7/2/1" ); |
| 53 | hmm.add(args: keys[1] % "1/9/0" ); |
| 54 | hmm.add(args: keys.back() % "5/4/1" ); |
| 55 | |
| 56 | // print |
| 57 | hmm.print(s: "HMM" ); |
| 58 | |
| 59 | // Convert to factor graph |
| 60 | DiscreteFactorGraph factorGraph(hmm); |
| 61 | |
| 62 | // Do max-prodcut |
| 63 | auto mpe = factorGraph.optimize(); |
| 64 | GTSAM_PRINT(mpe); |
| 65 | |
| 66 | // Create solver and eliminate |
| 67 | // This will create a DAG ordered with arrow of time reversed |
| 68 | DiscreteBayesNet::shared_ptr chordal = |
| 69 | factorGraph.eliminateSequential(ordering); |
| 70 | chordal->print(s: "Eliminated" ); |
| 71 | |
| 72 | // We can also sample from it |
| 73 | cout << "\n10 samples:" << endl; |
| 74 | for (size_t k = 0; k < 10; k++) { |
| 75 | auto sample = chordal->sample(); |
| 76 | GTSAM_PRINT(sample); |
| 77 | } |
| 78 | |
| 79 | // Or compute the marginals. This re-eliminates the FG into a Bayes tree |
| 80 | cout << "\nComputing Node Marginals .." << endl; |
| 81 | DiscreteMarginals marginals(factorGraph); |
| 82 | for (int k = 0; k < nrNodes; k++) { |
| 83 | Vector margProbs = marginals.marginalProbabilities(key: keys[k]); |
| 84 | stringstream ss; |
| 85 | ss << "marginal " << k; |
| 86 | print(v: margProbs, s: ss.str()); |
| 87 | } |
| 88 | |
| 89 | // TODO(frank): put in the glue to have DiscreteMarginals produce *arbitrary* |
| 90 | // joints efficiently, by the Bayes tree shortcut magic. All the code is there |
| 91 | // but it's not yet connected. |
| 92 | |
| 93 | return 0; |
| 94 | } |
| 95 | |