From 7c53ca3491c66eb06722bdd921cda66801a3d95f Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 17:04:22 -0400
Subject: [PATCH 01/22] fix docstrings

---
 gtsam/linear/GaussianConditional.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/gtsam/linear/GaussianConditional.h b/gtsam/linear/GaussianConditional.h
index d71119d6a..f1e2a2684 100644
--- a/gtsam/linear/GaussianConditional.h
+++ b/gtsam/linear/GaussianConditional.h
@@ -215,7 +215,7 @@ namespace gtsam {
      * Sample from conditional, zero parent version
      * Example:
      *   std::mt19937_64 rng(42);
-     *   auto sample = gbn.sample(&rng);
+     *   auto sample = gc.sample(&rng);
      */
     VectorValues sample(std::mt19937_64* rng) const;
 
@@ -224,7 +224,7 @@ namespace gtsam {
      * Example:
      *   std::mt19937_64 rng(42);
      *   VectorValues given = ...;
-     *   auto sample = gbn.sample(given, &rng);
+     *   auto sample = gc.sample(given, &rng);
      */
     VectorValues sample(const VectorValues& parentsValues,
                         std::mt19937_64* rng) const;

From 6527437c672e56fbac322e21353e56d65eb32e68 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 17:12:13 -0400
Subject: [PATCH 02/22] create common PRNG object in DiscreteConditional

---
 gtsam/discrete/DiscreteConditional.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/gtsam/discrete/DiscreteConditional.cpp b/gtsam/discrete/DiscreteConditional.cpp
index 9a3b6cf46..27ec69d44 100644
--- a/gtsam/discrete/DiscreteConditional.cpp
+++ b/gtsam/discrete/DiscreteConditional.cpp
@@ -32,6 +32,9 @@
 #include <utility>
 #include <vector>
 
+// In wrappers we can access std::mt19937_64 via gtsam.MT19937
+static std::mt19937_64 kRandomNumberGenerator(2);
+
 using namespace std;
 using std::pair;
 using std::stringstream;
@@ -287,8 +290,6 @@ void DiscreteConditional::sampleInPlace(DiscreteValues* values) const {
 
 /* ************************************************************************** */
 size_t DiscreteConditional::sample(const DiscreteValues& parentsValues) const {
-  static mt19937 rng(2);  // random number generator
-
   // Get the correct conditional distribution
   ADT pFS = choose(parentsValues, true);  // P(F|S=parentsValues)
 
@@ -310,7 +311,7 @@ size_t DiscreteConditional::sample(const DiscreteValues& parentsValues) const {
     }
   }
   std::discrete_distribution<size_t> distribution(p.begin(), p.end());
-  return distribution(rng);
+  return distribution(kRandomNumberGenerator);
 }
 
 /* ************************************************************************** */

From 82190cb7eb409127905404d1226a8c8e37818a76 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 17:15:04 -0400
Subject: [PATCH 03/22] update wrapper comment in GaussianConditional

---
 gtsam/linear/GaussianConditional.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gtsam/linear/GaussianConditional.cpp b/gtsam/linear/GaussianConditional.cpp
index 7735e5129..a77639c00 100644
--- a/gtsam/linear/GaussianConditional.cpp
+++ b/gtsam/linear/GaussianConditional.cpp
@@ -34,7 +34,7 @@
 #include <string>
 #include <cmath>
 
-// In Wrappers we have no access to this so have a default ready
+// In wrappers we can access std::mt19937_64 via gtsam.MT19937
 static std::mt19937_64 kRandomNumberGenerator(42);
 
 using namespace std;

From b58d509b68eaf214a9837a19ff4a3ac45eb403e4 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 17:15:30 -0400
Subject: [PATCH 04/22] add python test for sampling from GaussianConditional
 with a PRNG

---
 python/gtsam/tests/test_GaussianBayesNet.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/python/gtsam/tests/test_GaussianBayesNet.py b/python/gtsam/tests/test_GaussianBayesNet.py
index 05522441b..f06512e93 100644
--- a/python/gtsam/tests/test_GaussianBayesNet.py
+++ b/python/gtsam/tests/test_GaussianBayesNet.py
@@ -23,11 +23,12 @@ _x_ = 11
 _y_ = 22
 _z_ = 33
 
+I_1x1 = np.eye(1, dtype=float)
+
 
 def smallBayesNet():
     """Create a small Bayes Net for testing"""
     bayesNet = GaussianBayesNet()
-    I_1x1 = np.eye(1, dtype=float)
     bayesNet.push_back(GaussianConditional(_x_, [9.0], I_1x1, _y_, I_1x1))
     bayesNet.push_back(GaussianConditional(_y_, [5.0], I_1x1))
     return bayesNet
@@ -51,8 +52,9 @@ class TestGaussianBayesNet(GtsamTestCase):
         values.insert(_x_, np.array([9.0]))
         values.insert(_y_, np.array([5.0]))
         for i in [0, 1]:
-            self.assertAlmostEqual(bayesNet.at(i).logProbability(values),
-                                   np.log(bayesNet.at(i).evaluate(values)))
+            self.assertAlmostEqual(
+                bayesNet.at(i).logProbability(values),
+                np.log(bayesNet.at(i).evaluate(values)))
         self.assertAlmostEqual(bayesNet.logProbability(values),
                                np.log(bayesNet.evaluate(values)))
 
@@ -66,6 +68,16 @@ class TestGaussianBayesNet(GtsamTestCase):
         mean = bayesNet.optimize()
         self.gtsamAssertEquals(sample, mean, tol=3.0)
 
+        # Sample with rng
+        rng = gtsam.MT19937(42)
+        conditional = GaussianConditional(_x_, [9.0], I_1x1)
+        # Sample multiple times and average to gey mean
+        val = 0
+        niters = 10000
+        for _ in range(niters):
+            val += conditional.sample(rng).at(_x_).item()
+        self.assertAlmostEqual(val / niters, 9.0, 2)
+
 
 if __name__ == "__main__":
     unittest.main()

From 84d8c7ed78bdbb36c5bd6b4df7e96d20242c7bb4 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 17:26:22 -0400
Subject: [PATCH 05/22] DiscreteConditional::sample which uses a pseudo RNG

---
 gtsam/discrete/DiscreteConditional.cpp | 23 ++++++++++++++++++++---
 gtsam/discrete/DiscreteConditional.h   | 22 ++++++++++++++++++++++
 2 files changed, 42 insertions(+), 3 deletions(-)

diff --git a/gtsam/discrete/DiscreteConditional.cpp b/gtsam/discrete/DiscreteConditional.cpp
index 27ec69d44..75e090572 100644
--- a/gtsam/discrete/DiscreteConditional.cpp
+++ b/gtsam/discrete/DiscreteConditional.cpp
@@ -290,6 +290,12 @@ void DiscreteConditional::sampleInPlace(DiscreteValues* values) const {
 
 /* ************************************************************************** */
 size_t DiscreteConditional::sample(const DiscreteValues& parentsValues) const {
+  return sample(parentsValues, &kRandomNumberGenerator);
+}
+
+/* ************************************************************************** */
+size_t DiscreteConditional::sample(const DiscreteValues& parentsValues,
+                                   std::mt19937_64* rng) const {
   // Get the correct conditional distribution
   ADT pFS = choose(parentsValues, true);  // P(F|S=parentsValues)
 
@@ -311,27 +317,38 @@ size_t DiscreteConditional::sample(const DiscreteValues& parentsValues) const {
     }
   }
   std::discrete_distribution<size_t> distribution(p.begin(), p.end());
-  return distribution(kRandomNumberGenerator);
+  return distribution(*rng);
 }
 
 /* ************************************************************************** */
 size_t DiscreteConditional::sample(size_t parent_value) const {
+  return sample(parent_value, &kRandomNumberGenerator);
+}
+
+/* ************************************************************************** */
+size_t DiscreteConditional::sample(size_t parent_value,
+                                   std::mt19937_64* rng) const {
   if (nrParents() != 1)
     throw std::invalid_argument(
         "Single value sample() can only be invoked on single-parent "
         "conditional");
   DiscreteValues values;
   values.emplace(keys_.back(), parent_value);
-  return sample(values);
+  return sample(values, rng);
 }
 
 /* ************************************************************************** */
 size_t DiscreteConditional::sample() const {
+  return sample(&kRandomNumberGenerator);
+}
+
+/* ************************************************************************** */
+size_t DiscreteConditional::sample(std::mt19937_64* rng) const {
   if (nrParents() != 0)
     throw std::invalid_argument(
         "sample() can only be invoked on no-parent prior");
   DiscreteValues values;
-  return sample(values);
+  return sample(values, rng);
 }
 
 /* ************************************************************************* */
diff --git a/gtsam/discrete/DiscreteConditional.h b/gtsam/discrete/DiscreteConditional.h
index c22fcdf85..02bd76da1 100644
--- a/gtsam/discrete/DiscreteConditional.h
+++ b/gtsam/discrete/DiscreteConditional.h
@@ -23,6 +23,7 @@
 #include <gtsam/inference/Conditional-inst.h>
 
 #include <memory>
+#include <random>  // for std::mt19937_64
 #include <string>
 #include <vector>
 
@@ -201,12 +202,33 @@ class GTSAM_EXPORT DiscreteConditional
    */
   virtual size_t sample(const DiscreteValues& parentsValues) const;
 
+  /**
+   * Sample from conditional, given missing variables
+   * Example:
+   *   std::mt19937_64 rng(42);
+   *   DiscreteValues given = ...;
+   *   size_t sample = dc.sample(given, &rng);
+   */
+  size_t sample(const DiscreteValues& parentsValues,
+                std::mt19937_64* rng) const;
+
   /// Single parent version.
   size_t sample(size_t parent_value) const;
 
+  /// Single parent version with PRNG
+  size_t sample(size_t parent_value, std::mt19937_64* rng) const;
+
   /// Zero parent version.
   size_t sample() const;
 
+  /**
+   * Sample from conditional, zero parent version
+   * Example:
+   *   std::mt19937_64 rng(42);
+   *   auto sample = dc.sample(&rng);
+   */
+  size_t sample(std::mt19937_64* rng) const;
+
   /**
    * @brief Return assignment for single frontal variable that maximizes value.
    * @param parentsValues Known assignments for the parents.

From 42959035139147dfece24dac96fe874ab9cd5e82 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:08:09 -0400
Subject: [PATCH 06/22] default rng argument to make code DRY

---
 gtsam/discrete/DiscreteBayesNet.cpp    |  7 +++---
 gtsam/discrete/DiscreteBayesNet.h      |  5 ++--
 gtsam/discrete/DiscreteConditional.cpp | 25 +++----------------
 gtsam/discrete/DiscreteConditional.h   | 34 ++++++++++++--------------
 gtsam/discrete/TableDistribution.cpp   |  7 +++---
 gtsam/discrete/TableDistribution.h     |  5 +++-
 6 files changed, 33 insertions(+), 50 deletions(-)

diff --git a/gtsam/discrete/DiscreteBayesNet.cpp b/gtsam/discrete/DiscreteBayesNet.cpp
index 7c6da3dac..ef1286032 100644
--- a/gtsam/discrete/DiscreteBayesNet.cpp
+++ b/gtsam/discrete/DiscreteBayesNet.cpp
@@ -50,12 +50,13 @@ double DiscreteBayesNet::evaluate(const DiscreteValues& values) const {
 }
 
 /* ************************************************************************* */
-DiscreteValues DiscreteBayesNet::sample() const {
+DiscreteValues DiscreteBayesNet::sample(std::mt19937_64* rng) const {
   DiscreteValues result;
   return sample(result);
 }
 
-DiscreteValues DiscreteBayesNet::sample(DiscreteValues result) const {
+DiscreteValues DiscreteBayesNet::sample(DiscreteValues result,
+                                        std::mt19937_64* rng) const {
   // sample each node in turn in topological sort order (parents first)
   for (auto it = std::make_reverse_iterator(end());
        it != std::make_reverse_iterator(begin()); ++it) {
@@ -63,7 +64,7 @@ DiscreteValues DiscreteBayesNet::sample(DiscreteValues result) const {
     // Sample the conditional only if value for j not already in result
     const Key j = conditional->firstFrontalKey();
     if (result.count(j) == 0) {
-      conditional->sampleInPlace(&result);
+      conditional->sampleInPlace(&result, rng);
     }
   }
   return result;
diff --git a/gtsam/discrete/DiscreteBayesNet.h b/gtsam/discrete/DiscreteBayesNet.h
index e15576b37..3d84cd656 100644
--- a/gtsam/discrete/DiscreteBayesNet.h
+++ b/gtsam/discrete/DiscreteBayesNet.h
@@ -112,7 +112,7 @@ class GTSAM_EXPORT DiscreteBayesNet: public BayesNet<DiscreteConditional> {
      *
      * @return a sampled value for all variables.
      */
-    DiscreteValues sample() const;
+    DiscreteValues sample(std::mt19937_64* rng = &kRandomNumberGenerator) const;
 
     /**
      * @brief do ancestral sampling, given certain variables.
@@ -122,7 +122,8 @@ class GTSAM_EXPORT DiscreteBayesNet: public BayesNet<DiscreteConditional> {
      *
      * @return given values extended with sampled value for all other variables.
      */
-    DiscreteValues sample(DiscreteValues given) const;
+    DiscreteValues sample(DiscreteValues given,
+                          std::mt19937_64* rng = &kRandomNumberGenerator) const;
 
     /**
      * @brief Prune the Bayes net
diff --git a/gtsam/discrete/DiscreteConditional.cpp b/gtsam/discrete/DiscreteConditional.cpp
index 75e090572..f5ec20045 100644
--- a/gtsam/discrete/DiscreteConditional.cpp
+++ b/gtsam/discrete/DiscreteConditional.cpp
@@ -32,9 +32,6 @@
 #include <utility>
 #include <vector>
 
-// In wrappers we can access std::mt19937_64 via gtsam.MT19937
-static std::mt19937_64 kRandomNumberGenerator(2);
-
 using namespace std;
 using std::pair;
 using std::stringstream;
@@ -271,7 +268,8 @@ size_t DiscreteConditional::argmax(const DiscreteValues& parentsValues) const {
 }
 
 /* ************************************************************************** */
-void DiscreteConditional::sampleInPlace(DiscreteValues* values) const {
+void DiscreteConditional::sampleInPlace(DiscreteValues* values,
+                                        std::mt19937_64* rng) const {
   // throw if more than one frontal:
   if (nrFrontals() != 1) {
     throw std::invalid_argument(
@@ -284,13 +282,8 @@ void DiscreteConditional::sampleInPlace(DiscreteValues* values) const {
     throw std::invalid_argument(
         "DiscreteConditional::sampleInPlace: values already contains j");
   }
-  size_t sampled = sample(*values);  // Sample variable given parents
-  (*values)[j] = sampled;            // store result in partial solution
-}
-
-/* ************************************************************************** */
-size_t DiscreteConditional::sample(const DiscreteValues& parentsValues) const {
-  return sample(parentsValues, &kRandomNumberGenerator);
+  size_t sampled = sample(*values, rng);  // Sample variable given parents
+  (*values)[j] = sampled;                 // store result in partial solution
 }
 
 /* ************************************************************************** */
@@ -320,11 +313,6 @@ size_t DiscreteConditional::sample(const DiscreteValues& parentsValues,
   return distribution(*rng);
 }
 
-/* ************************************************************************** */
-size_t DiscreteConditional::sample(size_t parent_value) const {
-  return sample(parent_value, &kRandomNumberGenerator);
-}
-
 /* ************************************************************************** */
 size_t DiscreteConditional::sample(size_t parent_value,
                                    std::mt19937_64* rng) const {
@@ -337,11 +325,6 @@ size_t DiscreteConditional::sample(size_t parent_value,
   return sample(values, rng);
 }
 
-/* ************************************************************************** */
-size_t DiscreteConditional::sample() const {
-  return sample(&kRandomNumberGenerator);
-}
-
 /* ************************************************************************** */
 size_t DiscreteConditional::sample(std::mt19937_64* rng) const {
   if (nrParents() != 0)
diff --git a/gtsam/discrete/DiscreteConditional.h b/gtsam/discrete/DiscreteConditional.h
index 02bd76da1..970a0a142 100644
--- a/gtsam/discrete/DiscreteConditional.h
+++ b/gtsam/discrete/DiscreteConditional.h
@@ -27,6 +27,9 @@
 #include <string>
 #include <vector>
 
+// In wrappers we can access std::mt19937_64 via gtsam.MT19937
+static std::mt19937_64 kRandomNumberGenerator(42);
+
 namespace gtsam {
 
 /**
@@ -195,31 +198,23 @@ class GTSAM_EXPORT DiscreteConditional
   /** Single variable version of likelihood. */
   DecisionTreeFactor::shared_ptr likelihood(size_t frontal) const;
 
-  /**
-   * sample
-   * @param parentsValues Known values of the parents
-   * @return sample from conditional
-   */
-  virtual size_t sample(const DiscreteValues& parentsValues) const;
-
   /**
    * Sample from conditional, given missing variables
    * Example:
    *   std::mt19937_64 rng(42);
    *   DiscreteValues given = ...;
    *   size_t sample = dc.sample(given, &rng);
+   *
+   * @param parentsValues Known values of the parents
+   * @param rng Pseudo-Random Number Generator.
+   * @return sample from conditional
    */
-  size_t sample(const DiscreteValues& parentsValues,
-                std::mt19937_64* rng) const;
+  virtual size_t sample(const DiscreteValues& parentsValues,
+                        std::mt19937_64* rng = &kRandomNumberGenerator) const;
 
   /// Single parent version.
-  size_t sample(size_t parent_value) const;
-
-  /// Single parent version with PRNG
-  size_t sample(size_t parent_value, std::mt19937_64* rng) const;
-
-  /// Zero parent version.
-  size_t sample() const;
+  size_t sample(size_t parent_value,
+                std::mt19937_64* rng = &kRandomNumberGenerator) const;
 
   /**
    * Sample from conditional, zero parent version
@@ -227,7 +222,7 @@ class GTSAM_EXPORT DiscreteConditional
    *   std::mt19937_64 rng(42);
    *   auto sample = dc.sample(&rng);
    */
-  size_t sample(std::mt19937_64* rng) const;
+  size_t sample(std::mt19937_64* rng = &kRandomNumberGenerator) const;
 
   /**
    * @brief Return assignment for single frontal variable that maximizes value.
@@ -249,8 +244,9 @@ class GTSAM_EXPORT DiscreteConditional
   /// @name Advanced Interface
   /// @{
 
-  /// sample in place, stores result in partial solution
-  void sampleInPlace(DiscreteValues* parentsValues) const;
+  /// Sample in place with optional PRNG, stores result in partial solution
+  void sampleInPlace(DiscreteValues* parentsValues,
+                     std::mt19937_64* rng = &kRandomNumberGenerator) const;
 
   /// Return all assignments for frontal variables.
   std::vector<DiscreteValues> frontalAssignments() const;
diff --git a/gtsam/discrete/TableDistribution.cpp b/gtsam/discrete/TableDistribution.cpp
index ce0d92bff..614918b74 100644
--- a/gtsam/discrete/TableDistribution.cpp
+++ b/gtsam/discrete/TableDistribution.cpp
@@ -144,9 +144,8 @@ void TableDistribution::prune(size_t maxNrAssignments) {
 }
 
 /* ****************************************************************************/
-size_t TableDistribution::sample(const DiscreteValues& parentsValues) const {
-  static mt19937 rng(2);  // random number generator
-
+size_t TableDistribution::sample(const DiscreteValues& parentsValues,
+                                 std::mt19937_64* rng) const {
   DiscreteKeys parentsKeys;
   for (auto&& [key, _] : parentsValues) {
     parentsKeys.push_back({key, table_.cardinality(key)});
@@ -173,7 +172,7 @@ size_t TableDistribution::sample(const DiscreteValues& parentsValues) const {
     }
   }
   std::discrete_distribution<size_t> distribution(p.begin(), p.end());
-  return distribution(rng);
+  return distribution(*rng);
 }
 
 }  // namespace gtsam
diff --git a/gtsam/discrete/TableDistribution.h b/gtsam/discrete/TableDistribution.h
index 8e28bed5f..a556d6edb 100644
--- a/gtsam/discrete/TableDistribution.h
+++ b/gtsam/discrete/TableDistribution.h
@@ -143,9 +143,12 @@ class GTSAM_EXPORT TableDistribution : public DiscreteConditional {
   /**
    * sample
    * @param parentsValues Known values of the parents
+   * @param rng Pseudo random number generator
    * @return sample from conditional
    */
-  virtual size_t sample(const DiscreteValues& parentsValues) const override;
+  virtual size_t sample(
+      const DiscreteValues& parentsValues,
+      std::mt19937_64* rng = &kRandomNumberGenerator) const override;
 
   /// @}
   /// @name Advanced Interface

From 347f4b2136a6cba76550645c788779a05f6847df Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:08:26 -0400
Subject: [PATCH 07/22] update wrapper

---
 gtsam/discrete/discrete.i | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/gtsam/discrete/discrete.i b/gtsam/discrete/discrete.i
index 27f1fdfa1..2f6687d44 100644
--- a/gtsam/discrete/discrete.i
+++ b/gtsam/discrete/discrete.i
@@ -138,10 +138,13 @@ virtual class DiscreteConditional : gtsam::DecisionTreeFactor {
   gtsam::DecisionTreeFactor* likelihood(
       const gtsam::DiscreteValues& frontalValues) const;
   gtsam::DecisionTreeFactor* likelihood(size_t value) const;
-  size_t sample(const gtsam::DiscreteValues& parentsValues) const;
-  size_t sample(size_t value) const;
-  size_t sample() const;
-  void sampleInPlace(gtsam::DiscreteValues @parentsValues) const;
+  size_t sample(const gtsam::DiscreteValues& parentsValues,
+                std::mt19937_64 @rng = &kRandomNumberGenerator) const;
+  size_t sample(size_t value,
+                std::mt19937_64 @rng = &kRandomNumberGenerator) const;
+  size_t sample(std::mt19937_64 @rng = &kRandomNumberGenerator) const;
+  void sampleInPlace(gtsam::DiscreteValues @parentsValues,
+                     std::mt19937_64 @rng = &kRandomNumberGenerator) const;
   size_t argmax(const gtsam::DiscreteValues& parentsValues) const;
 
   // Markdown and HTML
@@ -233,8 +236,11 @@ class DiscreteBayesNet {
   double evaluate(const gtsam::DiscreteValues& values) const;
   double operator()(const gtsam::DiscreteValues& values) const;
 
-  gtsam::DiscreteValues sample() const;
-  gtsam::DiscreteValues sample(gtsam::DiscreteValues given) const;
+  gtsam::DiscreteValues sample(std::mt19937_64
+                               @rng = &kRandomNumberGenerator) const;
+  gtsam::DiscreteValues sample(gtsam::DiscreteValues given,
+                               std::mt19937_64
+                               @rng = &kRandomNumberGenerator) const;
 
   string dot(
       const gtsam::KeyFormatter& keyFormatter = gtsam::DefaultKeyFormatter,

From 5848b58f1de45f590c700a8c2d981c00b34c7f07 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:10:17 -0400
Subject: [PATCH 08/22] remove PRNG declaration since it is now imported from
 DiscreteConditional

---
 gtsam/hybrid/HybridBayesNet.cpp | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/gtsam/hybrid/HybridBayesNet.cpp b/gtsam/hybrid/HybridBayesNet.cpp
index 5bb0723d2..b8ef7741a 100644
--- a/gtsam/hybrid/HybridBayesNet.cpp
+++ b/gtsam/hybrid/HybridBayesNet.cpp
@@ -25,9 +25,6 @@
 
 #include <memory>
 
-// In Wrappers we have no access to this so have a default ready
-static std::mt19937_64 kRandomNumberGenerator(42);
-
 namespace gtsam {
 
 /* ************************************************************************* */

From 6f8dfdf3f28de0c1d10a4efe7b45249fe7451b22 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:10:33 -0400
Subject: [PATCH 09/22] wrap rng in hybrid

---
 gtsam/hybrid/hybrid.i | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/gtsam/hybrid/hybrid.i b/gtsam/hybrid/hybrid.i
index d308573f1..bbee89733 100644
--- a/gtsam/hybrid/hybrid.i
+++ b/gtsam/hybrid/hybrid.i
@@ -158,6 +158,8 @@ class HybridBayesNet {
   gtsam::HybridValues optimize() const;
   gtsam::VectorValues optimize(const gtsam::DiscreteValues& assignment) const;
 
+  gtsam::HybridValues sample(const gtsam::HybridValues& given, std::mt19937_64@ rng) const;
+  gtsam::HybridValues sample(std::mt19937_64@ rng) const;
   gtsam::HybridValues sample(const gtsam::HybridValues& given) const;
   gtsam::HybridValues sample() const;
 

From eb35ff8df68690c945672cc56d24c1da7e13a919 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:10:54 -0400
Subject: [PATCH 10/22] wrap rng in DiscreteConditional

---
 gtsam/linear/linear.i | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/gtsam/linear/linear.i b/gtsam/linear/linear.i
index 01a28511c..6ccac0902 100644
--- a/gtsam/linear/linear.i
+++ b/gtsam/linear/linear.i
@@ -559,9 +559,12 @@ virtual class GaussianConditional : gtsam::JacobianFactor {
   gtsam::JacobianFactor* likelihood(
       const gtsam::VectorValues& frontalValues) const;
   gtsam::JacobianFactor* likelihood(gtsam::Vector frontal) const;
-  gtsam::VectorValues sample(const gtsam::VectorValues& parents) const;
+
+  gtsam::VectorValues sample(std::mt19937_64@ rng) const;
+  gtsam::VectorValues sample(const gtsam::VectorValues& parents, std::mt19937_64@ rng) const;
   gtsam::VectorValues sample() const;
-  
+  gtsam::VectorValues sample(const gtsam::VectorValues& parents) const;
+
   // Advanced Interface
   gtsam::VectorValues solveOtherRHS(const gtsam::VectorValues& parents,
                                     const gtsam::VectorValues& rhs) const;

From 0d70a475718fa49232692c7b1c6fd12c5904af01 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:11:29 -0400
Subject: [PATCH 11/22] wrap mt19937_64

---
 python/gtsam/specializations/gtsam.h | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/python/gtsam/specializations/gtsam.h b/python/gtsam/specializations/gtsam.h
index da8842eaf..d15a8fa11 100644
--- a/python/gtsam/specializations/gtsam.h
+++ b/python/gtsam/specializations/gtsam.h
@@ -10,3 +10,15 @@
  * with `PYBIND11_MAKE_OPAQUE` this allows the types to be modified with Python,
  * and saves one copy operation.
  */
+
+/*
+ * Custom Pybind11 module for the Mersenne-Twister PRNG object
+ * `std::mt19937_64`.
+ * This can be invoked with `gtsam.MT19937()` and passed
+ * wherever a rng pointer is expected.
+ */
+#include <random>
+py::class_<std::mt19937_64>(m_, "MT19937")
+    .def(py::init<>())
+    .def(py::init<std::mt19937_64::result_type>())
+    .def("__call__", &std::mt19937_64::operator());

From f853b1584b4243073752940a9e9c3aa45cc4b0a0 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:11:40 -0400
Subject: [PATCH 12/22] sampling tests

---
 python/gtsam/tests/test_DiscreteBayesNet.py | 49 ++++++++++++++++++---
 python/gtsam/tests/test_GaussianBayesNet.py |  2 +-
 2 files changed, 43 insertions(+), 8 deletions(-)

diff --git a/python/gtsam/tests/test_DiscreteBayesNet.py b/python/gtsam/tests/test_DiscreteBayesNet.py
index 36bed816d..471d82264 100644
--- a/python/gtsam/tests/test_DiscreteBayesNet.py
+++ b/python/gtsam/tests/test_DiscreteBayesNet.py
@@ -33,6 +33,39 @@ XRay = (2, 2)
 Dyspnea = (1, 2)
 
 
+class TestDiscreteConditional(GtsamTestCase):
+    """Tests for Discrete Conditional"""
+
+    def setUp(self):
+        self.key = (0, 2)
+        self.parent = (1, 2)
+        self.parents = DiscreteKeys()
+        self.parents.push_back(self.parent)
+
+    def test_sample(self):
+        """Tests to check sampling in DiscreteConditionals"""
+        rng = gtsam.MT19937(11)
+        niters = 1000
+
+        # Sample with only 1 variable
+        conditional = DiscreteConditional(self.key, "7/3")
+        p = 0
+        for _ in range(niters):
+            p += conditional.sample(rng)
+
+        self.assertAlmostEqual(p / niters, 0.3, 1)
+
+        # Sample with variable and parent
+        conditional = DiscreteConditional(self.key, self.parents, "7/3 2/8")
+        p = 0
+        parentValues = gtsam.DiscreteValues()
+        parentValues[self.parent[0]] = 1
+        for _ in range(niters):
+            p += conditional.sample(parentValues, rng)
+
+        self.assertAlmostEqual(p / niters, 0.8, 1)
+
+
 class TestDiscreteBayesNet(GtsamTestCase):
     """Tests for Discrete Bayes Nets."""
 
@@ -85,10 +118,12 @@ class TestDiscreteBayesNet(GtsamTestCase):
         # solve
         actualMPE = fg.optimize()
         expectedMPE = DiscreteValues()
-        for key in [Asia, Dyspnea, XRay, Tuberculosis, Smoking, Either, LungCancer, Bronchitis]:
+        for key in [
+                Asia, Dyspnea, XRay, Tuberculosis, Smoking, Either, LungCancer,
+                Bronchitis
+        ]:
             expectedMPE[key[0]] = 0
-        self.assertEqual(list(actualMPE.items()),
-                         list(expectedMPE.items()))
+        self.assertEqual(list(actualMPE.items()), list(expectedMPE.items()))
 
         # Check value for MPE is the same
         self.assertAlmostEqual(asia(actualMPE), fg(actualMPE))
@@ -104,8 +139,7 @@ class TestDiscreteBayesNet(GtsamTestCase):
             expectedMPE2[key[0]] = 0
         for key in [Asia, Dyspnea, Smoking, Bronchitis]:
             expectedMPE2[key[0]] = 1
-        self.assertEqual(list(actualMPE2.items()),
-                         list(expectedMPE2.items()))
+        self.assertEqual(list(actualMPE2.items()), list(expectedMPE2.items()))
 
         # now sample from it
         chordal2 = fg.eliminateSequential(ordering)
@@ -135,8 +169,9 @@ class TestDiscreteBayesNet(GtsamTestCase):
         # self.assertEqual(len(values), 5)
 
         for i in [0, 1, 2]:
-            self.assertAlmostEqual(fragment.at(i).logProbability(values),
-                                   math.log(fragment.at(i).evaluate(values)))
+            self.assertAlmostEqual(
+                fragment.at(i).logProbability(values),
+                math.log(fragment.at(i).evaluate(values)))
         self.assertAlmostEqual(fragment.logProbability(values),
                                math.log(fragment.evaluate(values)))
         actual = fragment.sample(given)
diff --git a/python/gtsam/tests/test_GaussianBayesNet.py b/python/gtsam/tests/test_GaussianBayesNet.py
index f06512e93..468ee8de5 100644
--- a/python/gtsam/tests/test_GaussianBayesNet.py
+++ b/python/gtsam/tests/test_GaussianBayesNet.py
@@ -76,7 +76,7 @@ class TestGaussianBayesNet(GtsamTestCase):
         niters = 10000
         for _ in range(niters):
             val += conditional.sample(rng).at(_x_).item()
-        self.assertAlmostEqual(val / niters, 9.0, 2)
+        self.assertAlmostEqual(val / niters, 9.0, 1)
 
 
 if __name__ == "__main__":

From 2e8f948e6639810d9f6cbb63d3967452d83b79bb Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 18:27:59 -0400
Subject: [PATCH 13/22] fix regressions

---
 gtsam/discrete/tests/testDiscreteBayesNet.cpp | 6 +++---
 gtsam/hybrid/HybridBayesNet.cpp               | 2 +-
 gtsam/hybrid/tests/testHybridBayesNet.cpp     | 8 ++++----
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/gtsam/discrete/tests/testDiscreteBayesNet.cpp b/gtsam/discrete/tests/testDiscreteBayesNet.cpp
index dd5d218f8..61e3a2820 100644
--- a/gtsam/discrete/tests/testDiscreteBayesNet.cpp
+++ b/gtsam/discrete/tests/testDiscreteBayesNet.cpp
@@ -99,9 +99,9 @@ TEST(DiscreteBayesNet, Asia) {
 
   // now sample from it
   DiscreteValues expectedSample{{Asia.first, 1},       {Dyspnea.first, 1},
-                                {XRay.first, 1},       {Tuberculosis.first, 0},
-                                {Smoking.first, 1},    {Either.first, 1},
-                                {LungCancer.first, 1}, {Bronchitis.first, 0}};
+                                {XRay.first, 0},       {Tuberculosis.first, 0},
+                                {Smoking.first, 1},    {Either.first, 0},
+                                {LungCancer.first, 0}, {Bronchitis.first, 1}};
   SETDEBUG("DiscreteConditional::sample", false);
   auto actualSample = chordal2->sample();
   EXPECT(assert_equal(expectedSample, actualSample));
diff --git a/gtsam/hybrid/HybridBayesNet.cpp b/gtsam/hybrid/HybridBayesNet.cpp
index b8ef7741a..3445ae2da 100644
--- a/gtsam/hybrid/HybridBayesNet.cpp
+++ b/gtsam/hybrid/HybridBayesNet.cpp
@@ -188,7 +188,7 @@ HybridValues HybridBayesNet::sample(const HybridValues &given,
     }
   }
   // Sample a discrete assignment.
-  const DiscreteValues assignment = dbn.sample(given.discrete());
+  const DiscreteValues assignment = dbn.sample(given.discrete(), rng);
   // Select the continuous Bayes net corresponding to the assignment.
   GaussianBayesNet gbn = choose(assignment);
   // Sample from the Gaussian Bayes net.
diff --git a/gtsam/hybrid/tests/testHybridBayesNet.cpp b/gtsam/hybrid/tests/testHybridBayesNet.cpp
index 3ddad23ff..2d9a483bf 100644
--- a/gtsam/hybrid/tests/testHybridBayesNet.cpp
+++ b/gtsam/hybrid/tests/testHybridBayesNet.cpp
@@ -90,7 +90,7 @@ TEST(HybridBayesNet, EvaluatePureDiscrete) {
 
   // sample
   std::mt19937_64 rng(42);
-  EXPECT(assert_equal(zero, bayesNet.sample(&rng)));
+  EXPECT(assert_equal(one, bayesNet.sample(&rng)));
   EXPECT(assert_equal(one, bayesNet.sample(one, &rng)));
   EXPECT(assert_equal(zero, bayesNet.sample(zero, &rng)));
 
@@ -616,13 +616,13 @@ TEST(HybridBayesNet, Sampling) {
   double discrete_sum =
       std::accumulate(discrete_samples.begin(), discrete_samples.end(),
                       decltype(discrete_samples)::value_type(0));
-  EXPECT_DOUBLES_EQUAL(0.477, discrete_sum / num_samples, 1e-9);
+  EXPECT_DOUBLES_EQUAL(0.519, discrete_sum / num_samples, 1e-9);
 
   VectorValues expected;
   // regression for specific RNG seed
 #if __APPLE__ || _WIN32
-  expected.insert({X(0), Vector1(-0.0131207162712)});
-  expected.insert({X(1), Vector1(-0.499026377568)});
+  expected.insert({X(0), Vector1(0.0252479903896)});
+  expected.insert({X(1), Vector1(-0.513637101911)});
 #elif __linux__
   expected.insert({X(0), Vector1(-0.00799425182219)});
   expected.insert({X(1), Vector1(-0.526463854268)});

From d01bf5d80b190f2f595c783de85958bac85f6206 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Thu, 15 May 2025 22:59:26 -0400
Subject: [PATCH 14/22] some docs

---
 python/gtsam/tests/test_DiscreteBayesNet.py | 2 ++
 python/gtsam/tests/test_GaussianBayesNet.py | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/python/gtsam/tests/test_DiscreteBayesNet.py b/python/gtsam/tests/test_DiscreteBayesNet.py
index 471d82264..861d6dc49 100644
--- a/python/gtsam/tests/test_DiscreteBayesNet.py
+++ b/python/gtsam/tests/test_DiscreteBayesNet.py
@@ -49,6 +49,7 @@ class TestDiscreteConditional(GtsamTestCase):
 
         # Sample with only 1 variable
         conditional = DiscreteConditional(self.key, "7/3")
+        # Sample multiple times and average to get mean
         p = 0
         for _ in range(niters):
             p += conditional.sample(rng)
@@ -57,6 +58,7 @@ class TestDiscreteConditional(GtsamTestCase):
 
         # Sample with variable and parent
         conditional = DiscreteConditional(self.key, self.parents, "7/3 2/8")
+        # Sample multiple times and average to get mean
         p = 0
         parentValues = gtsam.DiscreteValues()
         parentValues[self.parent[0]] = 1
diff --git a/python/gtsam/tests/test_GaussianBayesNet.py b/python/gtsam/tests/test_GaussianBayesNet.py
index 468ee8de5..31ec1c3bc 100644
--- a/python/gtsam/tests/test_GaussianBayesNet.py
+++ b/python/gtsam/tests/test_GaussianBayesNet.py
@@ -71,7 +71,7 @@ class TestGaussianBayesNet(GtsamTestCase):
         # Sample with rng
         rng = gtsam.MT19937(42)
         conditional = GaussianConditional(_x_, [9.0], I_1x1)
-        # Sample multiple times and average to gey mean
+        # Sample multiple times and average to get mean
         val = 0
         niters = 10000
         for _ in range(niters):

From 22a36b3d7ab08ee89e6fbb75446b98624b7cc2ec Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:07:21 -0400
Subject: [PATCH 15/22] pass in the rng within the DiscreteBayesNet sample call

---
 gtsam/discrete/DiscreteBayesNet.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gtsam/discrete/DiscreteBayesNet.cpp b/gtsam/discrete/DiscreteBayesNet.cpp
index 081177d59..11876c120 100644
--- a/gtsam/discrete/DiscreteBayesNet.cpp
+++ b/gtsam/discrete/DiscreteBayesNet.cpp
@@ -52,7 +52,7 @@ double DiscreteBayesNet::evaluate(const DiscreteValues& values) const {
 /* ************************************************************************* */
 DiscreteValues DiscreteBayesNet::sample(std::mt19937_64* rng) const {
   DiscreteValues result;
-  return sample(result);
+  return sample(result, rng);
 }
 
 DiscreteValues DiscreteBayesNet::sample(DiscreteValues result,

From 8b6ce55151c45fd9d5952129d2cbb5dd6681bcee Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:09:31 -0400
Subject: [PATCH 16/22] set default pointer value to nullptr which is more
 idiomatic

---
 gtsam/discrete/DiscreteBayesNet.h    | 4 ++--
 gtsam/discrete/DiscreteConditional.h | 9 ++++-----
 gtsam/discrete/TableDistribution.h   | 6 ++----
 3 files changed, 8 insertions(+), 11 deletions(-)

diff --git a/gtsam/discrete/DiscreteBayesNet.h b/gtsam/discrete/DiscreteBayesNet.h
index 3d84cd656..01c9c004c 100644
--- a/gtsam/discrete/DiscreteBayesNet.h
+++ b/gtsam/discrete/DiscreteBayesNet.h
@@ -112,7 +112,7 @@ class GTSAM_EXPORT DiscreteBayesNet: public BayesNet<DiscreteConditional> {
      *
      * @return a sampled value for all variables.
      */
-    DiscreteValues sample(std::mt19937_64* rng = &kRandomNumberGenerator) const;
+    DiscreteValues sample(std::mt19937_64* rng = nullptr) const;
 
     /**
      * @brief do ancestral sampling, given certain variables.
@@ -123,7 +123,7 @@ class GTSAM_EXPORT DiscreteBayesNet: public BayesNet<DiscreteConditional> {
      * @return given values extended with sampled value for all other variables.
      */
     DiscreteValues sample(DiscreteValues given,
-                          std::mt19937_64* rng = &kRandomNumberGenerator) const;
+                          std::mt19937_64* rng = nullptr) const;
 
     /**
      * @brief Prune the Bayes net
diff --git a/gtsam/discrete/DiscreteConditional.h b/gtsam/discrete/DiscreteConditional.h
index 970a0a142..00bfadd14 100644
--- a/gtsam/discrete/DiscreteConditional.h
+++ b/gtsam/discrete/DiscreteConditional.h
@@ -210,11 +210,10 @@ class GTSAM_EXPORT DiscreteConditional
    * @return sample from conditional
    */
   virtual size_t sample(const DiscreteValues& parentsValues,
-                        std::mt19937_64* rng = &kRandomNumberGenerator) const;
+                        std::mt19937_64* rng = nullptr) const;
 
   /// Single parent version.
-  size_t sample(size_t parent_value,
-                std::mt19937_64* rng = &kRandomNumberGenerator) const;
+  size_t sample(size_t parent_value, std::mt19937_64* rng = nullptr) const;
 
   /**
    * Sample from conditional, zero parent version
@@ -222,7 +221,7 @@ class GTSAM_EXPORT DiscreteConditional
    *   std::mt19937_64 rng(42);
    *   auto sample = dc.sample(&rng);
    */
-  size_t sample(std::mt19937_64* rng = &kRandomNumberGenerator) const;
+  size_t sample(std::mt19937_64* rng = nullptr) const;
 
   /**
    * @brief Return assignment for single frontal variable that maximizes value.
@@ -246,7 +245,7 @@ class GTSAM_EXPORT DiscreteConditional
 
   /// Sample in place with optional PRNG, stores result in partial solution
   void sampleInPlace(DiscreteValues* parentsValues,
-                     std::mt19937_64* rng = &kRandomNumberGenerator) const;
+                     std::mt19937_64* rng = nullptr) const;
 
   /// Return all assignments for frontal variables.
   std::vector<DiscreteValues> frontalAssignments() const;
diff --git a/gtsam/discrete/TableDistribution.h b/gtsam/discrete/TableDistribution.h
index a556d6edb..bd959d135 100644
--- a/gtsam/discrete/TableDistribution.h
+++ b/gtsam/discrete/TableDistribution.h
@@ -125,7 +125,6 @@ class GTSAM_EXPORT TableDistribution : public DiscreteConditional {
   /// Create new factor by maximizing over all values with the same separator.
   DiscreteFactor::shared_ptr max(const Ordering& keys) const override;
 
-
   /// Multiply by scalar s
   DiscreteFactor::shared_ptr operator*(double s) const override;
 
@@ -146,9 +145,8 @@ class GTSAM_EXPORT TableDistribution : public DiscreteConditional {
    * @param rng Pseudo random number generator
    * @return sample from conditional
    */
-  virtual size_t sample(
-      const DiscreteValues& parentsValues,
-      std::mt19937_64* rng = &kRandomNumberGenerator) const override;
+  virtual size_t sample(const DiscreteValues& parentsValues,
+                        std::mt19937_64* rng = nullptr) const override;
 
   /// @}
   /// @name Advanced Interface

From ab0a4a801e8c5ea1015aba18698181622e8aeedf Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:09:59 -0400
Subject: [PATCH 17/22] nullptr check and assign default

---
 gtsam/discrete/DiscreteConditional.cpp | 4 ++++
 gtsam/discrete/TableDistribution.cpp   | 4 ++++
 2 files changed, 8 insertions(+)

diff --git a/gtsam/discrete/DiscreteConditional.cpp b/gtsam/discrete/DiscreteConditional.cpp
index f5ec20045..959f3f2ce 100644
--- a/gtsam/discrete/DiscreteConditional.cpp
+++ b/gtsam/discrete/DiscreteConditional.cpp
@@ -309,6 +309,10 @@ size_t DiscreteConditional::sample(const DiscreteValues& parentsValues,
       return value;  // shortcut exit
     }
   }
+
+  // Check if rng is nullptr, then assign default
+  rng = (rng == nullptr) ? &kRandomNumberGenerator : rng;
+
   std::discrete_distribution<size_t> distribution(p.begin(), p.end());
   return distribution(*rng);
 }
diff --git a/gtsam/discrete/TableDistribution.cpp b/gtsam/discrete/TableDistribution.cpp
index 614918b74..b5a9f2c14 100644
--- a/gtsam/discrete/TableDistribution.cpp
+++ b/gtsam/discrete/TableDistribution.cpp
@@ -171,6 +171,10 @@ size_t TableDistribution::sample(const DiscreteValues& parentsValues,
       return value;  // shortcut exit
     }
   }
+
+  // Check if rng is nullptr, then assign default
+  rng = (rng == nullptr) ? &kRandomNumberGenerator : rng;
+
   std::discrete_distribution<size_t> distribution(p.begin(), p.end());
   return distribution(*rng);
 }

From 9fea009ac3d59790a9150c60456aab9e8c6bc854 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:10:35 -0400
Subject: [PATCH 18/22] use nullptr for pointer default in discrete wrapping

---
 gtsam/discrete/discrete.i | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/gtsam/discrete/discrete.i b/gtsam/discrete/discrete.i
index 2f6687d44..54d00f82a 100644
--- a/gtsam/discrete/discrete.i
+++ b/gtsam/discrete/discrete.i
@@ -139,12 +139,12 @@ virtual class DiscreteConditional : gtsam::DecisionTreeFactor {
       const gtsam::DiscreteValues& frontalValues) const;
   gtsam::DecisionTreeFactor* likelihood(size_t value) const;
   size_t sample(const gtsam::DiscreteValues& parentsValues,
-                std::mt19937_64 @rng = &kRandomNumberGenerator) const;
+                std::mt19937_64 @rng = nullptr) const;
   size_t sample(size_t value,
-                std::mt19937_64 @rng = &kRandomNumberGenerator) const;
-  size_t sample(std::mt19937_64 @rng = &kRandomNumberGenerator) const;
+                std::mt19937_64 @rng = nullptr) const;
+  size_t sample(std::mt19937_64 @rng = nullptr) const;
   void sampleInPlace(gtsam::DiscreteValues @parentsValues,
-                     std::mt19937_64 @rng = &kRandomNumberGenerator) const;
+                     std::mt19937_64 @rng = nullptr) const;
   size_t argmax(const gtsam::DiscreteValues& parentsValues) const;
 
   // Markdown and HTML
@@ -237,10 +237,10 @@ class DiscreteBayesNet {
   double operator()(const gtsam::DiscreteValues& values) const;
 
   gtsam::DiscreteValues sample(std::mt19937_64
-                               @rng = &kRandomNumberGenerator) const;
+                               @rng = nullptr) const;
   gtsam::DiscreteValues sample(gtsam::DiscreteValues given,
                                std::mt19937_64
-                               @rng = &kRandomNumberGenerator) const;
+                               @rng = nullptr) const;
 
   string dot(
       const gtsam::KeyFormatter& keyFormatter = gtsam::DefaultKeyFormatter,

From 30aeff9c73608cb059b1ca973d3d80feae9abcfe Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:13:31 -0400
Subject: [PATCH 19/22] add print to verify regression values

---
 gtsam/hybrid/tests/testHybridBayesNet.cpp | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/gtsam/hybrid/tests/testHybridBayesNet.cpp b/gtsam/hybrid/tests/testHybridBayesNet.cpp
index 2d9a483bf..d1ec855f4 100644
--- a/gtsam/hybrid/tests/testHybridBayesNet.cpp
+++ b/gtsam/hybrid/tests/testHybridBayesNet.cpp
@@ -628,6 +628,13 @@ TEST(HybridBayesNet, Sampling) {
   expected.insert({X(1), Vector1(-0.526463854268)});
 #endif
 
+  std::cout << std::setprecision(12)
+            << average_continuous.scale(1.0 / num_samples).at(X(0))
+            << std::endl;
+  std::cout << std::setprecision(12)
+            << average_continuous.scale(1.0 / num_samples).at(X(1))
+            << std::endl;
+
   EXPECT(assert_equal(expected, average_continuous.scale(1.0 / num_samples)));
 }
 

From a84ac8f0f2f7a4dcabc1c5d7afb1208cf03f1191 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:58:58 -0400
Subject: [PATCH 20/22] fix regression for linux

---
 gtsam/hybrid/tests/testHybridBayesNet.cpp | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/gtsam/hybrid/tests/testHybridBayesNet.cpp b/gtsam/hybrid/tests/testHybridBayesNet.cpp
index d1ec855f4..40620f13a 100644
--- a/gtsam/hybrid/tests/testHybridBayesNet.cpp
+++ b/gtsam/hybrid/tests/testHybridBayesNet.cpp
@@ -624,17 +624,10 @@ TEST(HybridBayesNet, Sampling) {
   expected.insert({X(0), Vector1(0.0252479903896)});
   expected.insert({X(1), Vector1(-0.513637101911)});
 #elif __linux__
-  expected.insert({X(0), Vector1(-0.00799425182219)});
-  expected.insert({X(1), Vector1(-0.526463854268)});
+  expected.insert({X(0), Vector1(0.0165089744897)});
+  expected.insert({X(1), Vector1(-0.454323399979)});
 #endif
 
-  std::cout << std::setprecision(12)
-            << average_continuous.scale(1.0 / num_samples).at(X(0))
-            << std::endl;
-  std::cout << std::setprecision(12)
-            << average_continuous.scale(1.0 / num_samples).at(X(1))
-            << std::endl;
-
   EXPECT(assert_equal(expected, average_continuous.scale(1.0 / num_samples)));
 }
 

From dcc47d268e12bee810f4a77400c801f2e4e7308b Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 00:59:15 -0400
Subject: [PATCH 21/22] reduce swap file size so CI doesn't fail

---
 .github/workflows/build-special.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/build-special.yml b/.github/workflows/build-special.yml
index c1d321cf9..b4e8fc842 100644
--- a/.github/workflows/build-special.yml
+++ b/.github/workflows/build-special.yml
@@ -201,7 +201,7 @@ jobs:
         if: runner.os == 'Linux'
         uses: pierotofy/set-swap-space@master
         with:
-          swap-size-gb: 12
+          swap-size-gb: 8
 
       - name: Build & Test
         run: |

From c4e1f7ec7fa1b735921a95e07b0556831edde125 Mon Sep 17 00:00:00 2001
From: Varun Agrawal <varagrawal@gmail.com>
Date: Fri, 16 May 2025 01:46:22 -0400
Subject: [PATCH 22/22] fix regression in test_HybridFactorGraph.py

---
 python/gtsam/tests/test_HybridFactorGraph.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/gtsam/tests/test_HybridFactorGraph.py b/python/gtsam/tests/test_HybridFactorGraph.py
index 6edab3449..8650cdcb6 100644
--- a/python/gtsam/tests/test_HybridFactorGraph.py
+++ b/python/gtsam/tests/test_HybridFactorGraph.py
@@ -287,8 +287,8 @@ class TestHybridGaussianFactorGraph(GtsamTestCase):
             print(f"P(mode=1; Z) = {marginals[1]}")
 
         # Check that the estimate is close to the true value.
-        self.assertAlmostEqual(marginals[0], 0.23, delta=0.01)
-        self.assertAlmostEqual(marginals[1], 0.77, delta=0.01)
+        self.assertAlmostEqual(marginals[0], 0.219, delta=0.01)
+        self.assertAlmostEqual(marginals[1], 0.781, delta=0.01)
 
         # Convert to factor graph using measurements.
         fg = bayesNet.toFactorGraph(measurements)