Switch to EliminateDiscrete for max-product

gtsam/hybrid/GaussianMixture.cpp

@@ -51,28 +51,28 @@ GaussianMixture::GaussianMixture(
                       Conditionals(discreteParents, conditionalsList)) {}
 
 /* *******************************************************************************/
-GaussianMixture::Sum GaussianMixture::add(
-    const GaussianMixture::Sum &sum) const {
-  using Y = GaussianMixtureFactor::GraphAndConstant;
+GaussianFactorGraphTree GaussianMixture::add(
+    const GaussianFactorGraphTree &sum) const {
+  using Y = GraphAndConstant;
   auto add = [](const Y &graph1, const Y &graph2) {
     auto result = graph1.graph;
     result.push_back(graph2.graph);
     return Y(result, graph1.constant + graph2.constant);
   };
-  const Sum tree = asGaussianFactorGraphTree();
+  const auto tree = asGaussianFactorGraphTree();
   return sum.empty() ? tree : sum.apply(tree, add);
 }
 
 /* *******************************************************************************/
-GaussianMixture::Sum GaussianMixture::asGaussianFactorGraphTree() const {
+GaussianFactorGraphTree GaussianMixture::asGaussianFactorGraphTree() const {
   auto lambda = [](const GaussianConditional::shared_ptr &conditional) {
     GaussianFactorGraph result;
     result.push_back(conditional);
     if (conditional) {
-      return GaussianMixtureFactor::GraphAndConstant(
+      return GraphAndConstant(
           result, conditional->logNormalizationConstant());
     } else {
-      return GaussianMixtureFactor::GraphAndConstant(result, 0.0);
+      return GraphAndConstant(result, 0.0);
     }
   };
   return {conditionals_, lambda};
@@ -103,7 +103,19 @@ GaussianConditional::shared_ptr GaussianMixture::operator()(
 /* *******************************************************************************/
 bool GaussianMixture::equals(const HybridFactor &lf, double tol) const {
   const This *e = dynamic_cast<const This *>(&lf);
-  return e != nullptr && BaseFactor::equals(*e, tol);
+  if (e == nullptr) return false;
+
+  // This will return false if either conditionals_ is empty or e->conditionals_
+  // is empty, but not if both are empty or both are not empty:
+  if (conditionals_.empty() ^ e->conditionals_.empty()) return false;
+
+  // Check the base and the factors:
+  return BaseFactor::equals(*e, tol) &&
+         conditionals_.equals(e->conditionals_,
+                              [tol](const GaussianConditional::shared_ptr &f1,
+                                    const GaussianConditional::shared_ptr &f2) {
+                                return f1->equals(*(f2), tol);
+                              });
 }
 
 /* *******************************************************************************/

gtsam/hybrid/GaussianMixture.h

@@ -59,9 +59,6 @@ class GTSAM_EXPORT GaussianMixture
   using BaseFactor = HybridFactor;
   using BaseConditional = Conditional<HybridFactor, GaussianMixture>;
 
-  /// Alias for DecisionTree of GaussianFactorGraphs
-  using Sum = DecisionTree<Key, GaussianMixtureFactor::GraphAndConstant>;
-
   /// typedef for Decision Tree of Gaussian Conditionals
   using Conditionals = DecisionTree<Key, GaussianConditional::shared_ptr>;
 
@@ -71,7 +68,7 @@ class GTSAM_EXPORT GaussianMixture
   /**
    * @brief Convert a DecisionTree of factors into a DT of Gaussian FGs.
    */
-  Sum asGaussianFactorGraphTree() const;
+  GaussianFactorGraphTree asGaussianFactorGraphTree() const;
 
   /**
    * @brief Helper function to get the pruner functor.
@@ -172,6 +169,16 @@ class GTSAM_EXPORT GaussianMixture
    */
   double error(const HybridValues &values) const override;
 
+  //   /// Calculate probability density for given values `x`.
+  //   double evaluate(const HybridValues &values) const;
+
+  //   /// Evaluate probability density, sugar.
+  //   double operator()(const HybridValues &values) const { return
+  //   evaluate(values); }
+
+  //   /// Calculate log-density for given values `x`.
+  //   double logDensity(const HybridValues &values) const;
+
   /**
    * @brief Prune the decision tree of Gaussian factors as per the discrete
    * `decisionTree`.
@@ -186,9 +193,9 @@ class GTSAM_EXPORT GaussianMixture
    * maintaining the decision tree structure.
    *
    * @param sum Decision Tree of Gaussian Factor Graphs
-   * @return Sum
+   * @return GaussianFactorGraphTree
    */
-  Sum add(const Sum &sum) const;
+  GaussianFactorGraphTree add(const GaussianFactorGraphTree &sum) const;
   /// @}
 };
 

gtsam/hybrid/GaussianMixtureFactor.cpp

@@ -81,32 +81,36 @@ void GaussianMixtureFactor::print(const std::string &s,
 }
 
 /* *******************************************************************************/
-const GaussianMixtureFactor::Mixture GaussianMixtureFactor::factors() const {
-  return Mixture(factors_, [](const FactorAndConstant &factor_z) {
-    return factor_z.factor;
-  });
+GaussianFactor::shared_ptr GaussianMixtureFactor::factor(
+    const DiscreteValues &assignment) const {
+  return factors_(assignment).factor;
 }
 
 /* *******************************************************************************/
-GaussianMixtureFactor::Sum GaussianMixtureFactor::add(
-    const GaussianMixtureFactor::Sum &sum) const {
-  using Y = GaussianMixtureFactor::GraphAndConstant;
+double GaussianMixtureFactor::constant(const DiscreteValues &assignment) const {
+  return factors_(assignment).constant;
+}
+
+/* *******************************************************************************/
+GaussianFactorGraphTree GaussianMixtureFactor::add(
+    const GaussianFactorGraphTree &sum) const {
+  using Y = GraphAndConstant;
   auto add = [](const Y &graph1, const Y &graph2) {
     auto result = graph1.graph;
     result.push_back(graph2.graph);
     return Y(result, graph1.constant + graph2.constant);
   };
-  const Sum tree = asGaussianFactorGraphTree();
+  const auto tree = asGaussianFactorGraphTree();
   return sum.empty() ? tree : sum.apply(tree, add);
 }
 
 /* *******************************************************************************/
-GaussianMixtureFactor::Sum GaussianMixtureFactor::asGaussianFactorGraphTree()
+GaussianFactorGraphTree GaussianMixtureFactor::asGaussianFactorGraphTree()
     const {
   auto wrap = [](const FactorAndConstant &factor_z) {
     GaussianFactorGraph result;
     result.push_back(factor_z.factor);
-    return GaussianMixtureFactor::GraphAndConstant(result, factor_z.constant);
+    return GraphAndConstant(result, factor_z.constant);
   };
   return {factors_, wrap};
 }

gtsam/hybrid/GaussianMixtureFactor.h

@@ -62,6 +62,7 @@ class GTSAM_EXPORT GaussianMixtureFactor : public HybridFactor {
       // Note: constant is log of normalization constant for probabilities.
       // Errors is the negative log-likelihood,
       // hence we subtract the constant here.
+      if (!factor) return 0.0;  // If nullptr, return 0.0 error
       return factor->error(values) - constant;
     }
 
@@ -71,22 +72,6 @@ class GTSAM_EXPORT GaussianMixtureFactor : public HybridFactor {
     }
   };
 
-  /// Gaussian factor graph and log of normalizing constant.
-  struct GraphAndConstant {
-    GaussianFactorGraph graph;
-    double constant;
-
-    GraphAndConstant(const GaussianFactorGraph &graph, double constant)
-        : graph(graph), constant(constant) {}
-
-    // Check pointer equality.
-    bool operator==(const GraphAndConstant &other) const {
-      return graph == other.graph && constant == other.constant;
-    }
-  };
-
-  using Sum = DecisionTree<Key, GraphAndConstant>;
-
   /// typedef for Decision Tree of Gaussian factors and log-constant.
   using Factors = DecisionTree<Key, FactorAndConstant>;
   using Mixture = DecisionTree<Key, sharedFactor>;
@@ -99,9 +84,9 @@ class GTSAM_EXPORT GaussianMixtureFactor : public HybridFactor {
    * @brief Helper function to return factors and functional to create a
    * DecisionTree of Gaussian Factor Graphs.
    *
-   * @return Sum (DecisionTree<Key, GaussianFactorGraph>)
+   * @return GaussianFactorGraphTree
    */
-  Sum asGaussianFactorGraphTree() const;
+  GaussianFactorGraphTree asGaussianFactorGraphTree() const;
 
  public:
   /// @name Constructors
@@ -151,12 +136,16 @@ class GTSAM_EXPORT GaussianMixtureFactor : public HybridFactor {
   void print(
       const std::string &s = "GaussianMixtureFactor\n",
       const KeyFormatter &formatter = DefaultKeyFormatter) const override;
+
   /// @}
   /// @name Standard API
   /// @{
 
-  /// Getter for the underlying Gaussian Factor Decision Tree.
-  const Mixture factors() const;
+  /// Get factor at a given discrete assignment.
+  sharedFactor factor(const DiscreteValues &assignment) const;
+
+  /// Get constant at a given discrete assignment.
+  double constant(const DiscreteValues &assignment) const;
 
   /**
    * @brief Combine the Gaussian Factor Graphs in `sum` and `this` while
@@ -166,7 +155,7 @@ class GTSAM_EXPORT GaussianMixtureFactor : public HybridFactor {
    * variables.
    * @return Sum
    */
-  Sum add(const Sum &sum) const;
+  GaussianFactorGraphTree add(const GaussianFactorGraphTree &sum) const;
 
   /**
    * @brief Compute error of the GaussianMixtureFactor as a tree.
@@ -184,7 +173,8 @@ class GTSAM_EXPORT GaussianMixtureFactor : public HybridFactor {
   double error(const HybridValues &values) const override;
 
   /// Add MixtureFactor to a Sum, syntactic sugar.
-  friend Sum &operator+=(Sum &sum, const GaussianMixtureFactor &factor) {
+  friend GaussianFactorGraphTree &operator+=(
+      GaussianFactorGraphTree &sum, const GaussianMixtureFactor &factor) {
     sum = factor.add(sum);
     return sum;
   }

gtsam/hybrid/HybridBayesNet.cpp

@@ -26,6 +26,17 @@ static std::mt19937_64 kRandomNumberGenerator(42);
 
 namespace gtsam {
 
+/* ************************************************************************* */
+void HybridBayesNet::print(const std::string &s,
+                           const KeyFormatter &formatter) const {
+  Base::print(s, formatter);
+}
+
+/* ************************************************************************* */
+bool HybridBayesNet::equals(const This &bn, double tol) const {
+  return Base::equals(bn, tol);
+}
+
 /* ************************************************************************* */
 DecisionTreeFactor::shared_ptr HybridBayesNet::discreteConditionals() const {
   AlgebraicDecisionTree<Key> decisionTree;
@@ -271,12 +282,15 @@ double HybridBayesNet::evaluate(const HybridValues &values) const {
 
   // Iterate over each conditional.
   for (auto &&conditional : *this) {
+    // TODO: should be delegated to derived classes.
     if (auto gm = conditional->asMixture()) {
       const auto component = (*gm)(discreteValues);
       logDensity += component->logDensity(continuousValues);
+
     } else if (auto gc = conditional->asGaussian()) {
       // If continuous only, evaluate the probability and multiply.
       logDensity += gc->logDensity(continuousValues);
+
     } else if (auto dc = conditional->asDiscrete()) {
       // Conditional is discrete-only, so return its probability.
       probability *= dc->operator()(discreteValues);

gtsam/hybrid/HybridBayesNet.h

@@ -50,18 +50,14 @@ class GTSAM_EXPORT HybridBayesNet : public BayesNet<HybridConditional> {
   /// @name Testable
   /// @{
 
-  /** Check equality */
-  bool equals(const This &bn, double tol = 1e-9) const {
-    return Base::equals(bn, tol);
-  }
-
-  /// print graph
+  /// GTSAM-style printing
   void print(
       const std::string &s = "",
-      const KeyFormatter &formatter = DefaultKeyFormatter) const override {
-    Base::print(s, formatter);
-  }
+      const KeyFormatter &formatter = DefaultKeyFormatter) const override;
 
+  /// GTSAM-style equals
+  bool equals(const This& fg, double tol = 1e-9) const;
+
   /// @}
   /// @name Standard Interface
   /// @{

gtsam/hybrid/HybridConditional.cpp

@@ -17,6 +17,7 @@
 
 #include <gtsam/hybrid/HybridConditional.h>
 #include <gtsam/hybrid/HybridFactor.h>
+#include <gtsam/hybrid/HybridValues.h>
 #include <gtsam/inference/Conditional-inst.h>
 #include <gtsam/inference/Key.h>
 
@@ -102,7 +103,38 @@ void HybridConditional::print(const std::string &s,
 /* ************************************************************************ */
 bool HybridConditional::equals(const HybridFactor &other, double tol) const {
   const This *e = dynamic_cast<const This *>(&other);
-  return e != nullptr && BaseFactor::equals(*e, tol);
+  if (e == nullptr) return false;
+  if (auto gm = asMixture()) {
+    auto other = e->asMixture();
+    return other != nullptr && gm->equals(*other, tol);
+  }
+  if (auto gc = asGaussian()) {
+    auto other = e->asGaussian();
+    return other != nullptr && gc->equals(*other, tol);
+  }
+  if (auto dc = asDiscrete()) {
+    auto other = e->asDiscrete();
+    return other != nullptr && dc->equals(*other, tol);
+  }
+  return inner_->equals(*(e->inner_), tol);
+
+  return inner_ ? (e->inner_ ? inner_->equals(*(e->inner_), tol) : false)
+                : !(e->inner_);
+}
+
+/* ************************************************************************ */
+double HybridConditional::error(const HybridValues &values) const {
+  if (auto gm = asMixture()) {
+    return gm->error(values);
+  }
+  if (auto gc = asGaussian()) {
+    return gc->error(values.continuous());
+  }
+  if (auto dc = asDiscrete()) {
+    return -log((*dc)(values.discrete()));
+  }
+  throw std::runtime_error(
+      "HybridConditional::error: conditional type not handled");
 }
 
 }  // namespace gtsam

gtsam/hybrid/HybridConditional.h

@@ -176,15 +176,7 @@ class GTSAM_EXPORT HybridConditional
   boost::shared_ptr<Factor> inner() const { return inner_; }
 
   /// Return the error of the underlying conditional.
-  /// Currently only implemented for Gaussian mixture.
-  double error(const HybridValues& values) const override {
-    if (auto gm = asMixture()) {
-      return gm->error(values);
-    } else {
-      throw std::runtime_error(
-          "HybridConditional::error: only implemented for Gaussian mixture");
-    }
-  }
+  double error(const HybridValues& values) const override;
 
   /// @}
 
@@ -195,6 +187,7 @@ class GTSAM_EXPORT HybridConditional
   void serialize(Archive& ar, const unsigned int /*version*/) {
     ar& BOOST_SERIALIZATION_BASE_OBJECT_NVP(BaseFactor);
     ar& BOOST_SERIALIZATION_BASE_OBJECT_NVP(BaseConditional);
+    ar& BOOST_SERIALIZATION_NVP(inner_);
   }
 
 };  // HybridConditional