diff --git a/gtsam/base/utilities.h b/gtsam/base/utilities.h index a67c5d1b6..d2039cdea 100644 --- a/gtsam/base/utilities.h +++ b/gtsam/base/utilities.h @@ -1,11 +1,18 @@ #pragma once -#include -#include -#include - #include +#include +#include +#include +#include + +/** + * @brief Global default pseudo-random number generator object. + * In wrappers we can access std::mt19937_64 via gtsam.MT19937 + */ +static std::mt19937_64 kRandomNumberGenerator(42); + namespace gtsam { /** * For Python __str__(). @@ -28,7 +35,7 @@ private: std::streambuf* coutBuffer_; }; -} +} // namespace gtsam namespace gtsam { // Adapted from https://stackoverflow.com/a/32223343/9151520 diff --git a/gtsam/discrete/DiscreteConditional.cpp b/gtsam/discrete/DiscreteConditional.cpp index 959f3f2ce..7577c273b 100644 --- a/gtsam/discrete/DiscreteConditional.cpp +++ b/gtsam/discrete/DiscreteConditional.cpp @@ -18,6 +18,7 @@ #include #include +#include #include #include #include diff --git a/gtsam/discrete/DiscreteConditional.h b/gtsam/discrete/DiscreteConditional.h index 00bfadd14..76f979a23 100644 --- a/gtsam/discrete/DiscreteConditional.h +++ b/gtsam/discrete/DiscreteConditional.h @@ -27,9 +27,6 @@ #include #include -// In wrappers we can access std::mt19937_64 via gtsam.MT19937 -static std::mt19937_64 kRandomNumberGenerator(42); - namespace gtsam { /** diff --git a/gtsam/discrete/TableDistribution.cpp b/gtsam/discrete/TableDistribution.cpp index b5a9f2c14..8c87ceada 100644 --- a/gtsam/discrete/TableDistribution.cpp +++ b/gtsam/discrete/TableDistribution.cpp @@ -17,6 +17,7 @@ #include #include +#include #include #include #include diff --git a/gtsam/hybrid/HybridBayesNet.cpp b/gtsam/hybrid/HybridBayesNet.cpp index 3445ae2da..cf737bbb8 100644 --- a/gtsam/hybrid/HybridBayesNet.cpp +++ b/gtsam/hybrid/HybridBayesNet.cpp @@ -202,16 +202,6 @@ HybridValues HybridBayesNet::sample(std::mt19937_64 *rng) const { return sample(given, rng); } -/* ************************************************************************* */ -HybridValues HybridBayesNet::sample(const HybridValues &given) const { - return sample(given, &kRandomNumberGenerator); -} - -/* ************************************************************************* */ -HybridValues HybridBayesNet::sample() const { - return sample(&kRandomNumberGenerator); -} - /* ************************************************************************* */ AlgebraicDecisionTree HybridBayesNet::errorTree( const VectorValues &continuousValues) const { diff --git a/gtsam/hybrid/HybridBayesNet.h b/gtsam/hybrid/HybridBayesNet.h index 08bab93ec..0058c406c 100644 --- a/gtsam/hybrid/HybridBayesNet.h +++ b/gtsam/hybrid/HybridBayesNet.h @@ -181,10 +181,11 @@ class GTSAM_EXPORT HybridBayesNet : public BayesNet { * auto sample = bn.sample(given, &rng); * * @param given Values of missing variables. - * @param rng The pseudo-random number generator. + * @param rng The optional pseudo-random number generator. * @return HybridValues */ - HybridValues sample(const HybridValues &given, std::mt19937_64 *rng) const; + HybridValues sample(const HybridValues &given, + std::mt19937_64 *rng = nullptr) const; /** * @brief Sample using ancestral sampling. @@ -193,25 +194,10 @@ class GTSAM_EXPORT HybridBayesNet : public BayesNet { * std::mt19937_64 rng(42); * auto sample = bn.sample(&rng); * - * @param rng The pseudo-random number generator. + * @param rng The optional pseudo-random number generator. * @return HybridValues */ - HybridValues sample(std::mt19937_64 *rng) const; - - /** - * @brief Sample from an incomplete BayesNet, use default rng. - * - * @param given Values of missing variables. - * @return HybridValues - */ - HybridValues sample(const HybridValues &given) const; - - /** - * @brief Sample using ancestral sampling, use default rng. - * - * @return HybridValues - */ - HybridValues sample() const; + HybridValues sample(std::mt19937_64 *rng = nullptr) const; /** * @brief Prune the Bayes Net such that we have at most maxNrLeaves leaves. diff --git a/gtsam/hybrid/hybrid.i b/gtsam/hybrid/hybrid.i index bbee89733..dbc9bf5fe 100644 --- a/gtsam/hybrid/hybrid.i +++ b/gtsam/hybrid/hybrid.i @@ -158,10 +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; + gtsam::HybridValues sample(const gtsam::HybridValues& given, std::mt19937_64@ rng = nullptr) const; + gtsam::HybridValues sample(std::mt19937_64@ rng = nullptr) const; void print(string s = "HybridBayesNet\n", const gtsam::KeyFormatter& keyFormatter = diff --git a/gtsam/linear/GaussianBayesNet.cpp b/gtsam/linear/GaussianBayesNet.cpp index f8c4ddc4c..aa11ca911 100644 --- a/gtsam/linear/GaussianBayesNet.cpp +++ b/gtsam/linear/GaussianBayesNet.cpp @@ -26,9 +26,6 @@ using namespace std; using namespace gtsam; -// In Wrappers we have no access to this so have a default ready -static std::mt19937_64 kRandomNumberGenerator(42); - namespace gtsam { // Instantiate base class @@ -76,15 +73,6 @@ namespace gtsam { return result; } - /* ************************************************************************ */ - VectorValues GaussianBayesNet::sample() const { - return sample(&kRandomNumberGenerator); - } - - VectorValues GaussianBayesNet::sample(const VectorValues& given) const { - return sample(given, &kRandomNumberGenerator); - } - /* ************************************************************************ */ VectorValues GaussianBayesNet::optimizeGradientSearch() const { diff --git a/gtsam/linear/GaussianBayesNet.h b/gtsam/linear/GaussianBayesNet.h index 538d1532e..8d22826f8 100644 --- a/gtsam/linear/GaussianBayesNet.h +++ b/gtsam/linear/GaussianBayesNet.h @@ -131,7 +131,7 @@ namespace gtsam { * std::mt19937_64 rng(42); * auto sample = gbn.sample(&rng); */ - VectorValues sample(std::mt19937_64* rng) const; + VectorValues sample(std::mt19937_64* rng = nullptr) const; /** * Sample from an incomplete BayesNet, given missing variables @@ -140,13 +140,7 @@ namespace gtsam { * VectorValues given = ...; * auto sample = gbn.sample(given, &rng); */ - VectorValues sample(const VectorValues& given, std::mt19937_64* rng) const; - - /// Sample using ancestral sampling, use default rng - VectorValues sample() const; - - /// Sample from an incomplete BayesNet, use default rng - VectorValues sample(const VectorValues& given) const; + VectorValues sample(const VectorValues& given, std::mt19937_64* rng = nullptr) const; /** * Return ordering corresponding to a topological sort. diff --git a/gtsam/linear/GaussianConditional.cpp b/gtsam/linear/GaussianConditional.cpp index a77639c00..4430ac0cb 100644 --- a/gtsam/linear/GaussianConditional.cpp +++ b/gtsam/linear/GaussianConditional.cpp @@ -15,11 +15,12 @@ * @author Christian Potthast, Frank Dellaert */ +#include +#include #include #include #include #include -#include #ifdef __GNUC__ #pragma GCC diagnostic push @@ -34,9 +35,6 @@ #include #include -// In wrappers we can access std::mt19937_64 via gtsam.MT19937 -static std::mt19937_64 kRandomNumberGenerator(42); - using namespace std; namespace gtsam { @@ -347,6 +345,10 @@ namespace gtsam { VectorValues solution = solve(parentsValues); Key key = firstFrontalKey(); + + // Check if rng is nullptr, then assign default + rng = (rng == nullptr) ? &kRandomNumberGenerator : rng; + // The vector of sigma values for sampling. // If no model, initialize sigmas to 1, else to model sigmas const Vector& sigmas = (!model_) ? Vector::Ones(rows()) : model_->sigmas(); @@ -359,16 +361,7 @@ namespace gtsam { throw std::invalid_argument( "sample() can only be invoked on no-parent prior"); VectorValues values; - return sample(values); - } - - /* ************************************************************************ */ - VectorValues GaussianConditional::sample() const { - return sample(&kRandomNumberGenerator); - } - - VectorValues GaussianConditional::sample(const VectorValues& given) const { - return sample(given, &kRandomNumberGenerator); + return sample(values, rng); } /* ************************************************************************ */ diff --git a/gtsam/linear/GaussianConditional.h b/gtsam/linear/GaussianConditional.h index f1e2a2684..2a54e9c66 100644 --- a/gtsam/linear/GaussianConditional.h +++ b/gtsam/linear/GaussianConditional.h @@ -217,7 +217,7 @@ namespace gtsam { * std::mt19937_64 rng(42); * auto sample = gc.sample(&rng); */ - VectorValues sample(std::mt19937_64* rng) const; + VectorValues sample(std::mt19937_64* rng = nullptr) const; /** * Sample from conditional, given missing variables @@ -227,13 +227,7 @@ namespace gtsam { * auto sample = gc.sample(given, &rng); */ VectorValues sample(const VectorValues& parentsValues, - std::mt19937_64* rng) const; - - /// Sample, use default rng - VectorValues sample() const; - - /// Sample with given values, use default rng - VectorValues sample(const VectorValues& parentsValues) const; + std::mt19937_64* rng = nullptr) const; /// @} /// @name Linear algebra. diff --git a/gtsam/linear/linear.i b/gtsam/linear/linear.i index 6ccac0902..6f4524b27 100644 --- a/gtsam/linear/linear.i +++ b/gtsam/linear/linear.i @@ -560,10 +560,9 @@ virtual class GaussianConditional : gtsam::JacobianFactor { const gtsam::VectorValues& frontalValues) const; gtsam::JacobianFactor* likelihood(gtsam::Vector frontal) 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; + gtsam::VectorValues sample(std::mt19937_64 @rng = nullptr) const; + gtsam::VectorValues sample(const gtsam::VectorValues& parents, + std::mt19937_64 @rng = nullptr) const; // Advanced Interface gtsam::VectorValues solveOtherRHS(const gtsam::VectorValues& parents, @@ -629,9 +628,10 @@ virtual class GaussianBayesNet { gtsam::VectorValues optimize() const; gtsam::VectorValues optimize(const gtsam::VectorValues& given) const; gtsam::VectorValues optimizeGradientSearch() const; - - gtsam::VectorValues sample(const gtsam::VectorValues& given) const; - gtsam::VectorValues sample() const; + + gtsam::VectorValues sample(const gtsam::VectorValues& given, + std::mt19937_64 @rng = nullptr) const; + gtsam::VectorValues sample(std::mt19937_64 @rng = nullptr) const; gtsam::VectorValues backSubstitute(const gtsam::VectorValues& gx) const; gtsam::VectorValues backSubstituteTranspose(const gtsam::VectorValues& gx) const; diff --git a/gtsam/linear/tests/testGaussianConditional.cpp b/gtsam/linear/tests/testGaussianConditional.cpp index 0252368ea..7133a3437 100644 --- a/gtsam/linear/tests/testGaussianConditional.cpp +++ b/gtsam/linear/tests/testGaussianConditional.cpp @@ -441,7 +441,7 @@ TEST(GaussianConditional, likelihood) { /* ************************************************************************* */ // Test sampling -TEST(GaussianConditional, sample) { +TEST(GaussianConditional, Sample) { Matrix A1 = (Matrix(2, 2) << 1., 2., 3., 4.).finished(); const Vector2 b(20, 40), x1(3, 4); const double sigma = 0.01; @@ -465,8 +465,10 @@ TEST(GaussianConditional, sample) { auto actual3 = conditional.sample(given, &rng); EXPECT_LONGS_EQUAL(1, actual2.size()); // regressions -#if __APPLE__ || _WIN32 - EXPECT(assert_equal(Vector2(31.0111856, 64.9850775), actual2[X(0)], 1e-5)); +#if __APPLE__ +EXPECT(assert_equal(Vector2(31.0111856, 64.9850775), actual2[X(0)], 1e-5)); +#elif _WIN32 + EXPECT(assert_equal(Vector2(30.995317, 64.9943165), actual2[X(0)], 1e-5)); #elif __linux__ EXPECT(assert_equal(Vector2(30.9809331, 64.9927588), actual2[X(0)], 1e-5)); #endif diff --git a/gtsam/sfm/ShonanAveraging.cpp b/gtsam/sfm/ShonanAveraging.cpp index 7d68c7d29..b268fda6f 100644 --- a/gtsam/sfm/ShonanAveraging.cpp +++ b/gtsam/sfm/ShonanAveraging.cpp @@ -42,7 +42,7 @@ namespace gtsam { // In Wrappers we have no access to this so have a default ready -static std::mt19937 kRandomNumberGenerator(42); +static std::mt19937 kPRNG(42); using Sparse = Eigen::SparseMatrix; @@ -869,7 +869,7 @@ Values ShonanAveraging::initializeRandomly(std::mt19937 &rng) const { /* ************************************************************************* */ template Values ShonanAveraging::initializeRandomly() const { - return initializeRandomly(kRandomNumberGenerator); + return initializeRandomly(kPRNG); } /* ************************************************************************* */ @@ -883,7 +883,7 @@ Values ShonanAveraging::initializeRandomlyAt(size_t p, /* ************************************************************************* */ template Values ShonanAveraging::initializeRandomlyAt(size_t p) const { - return initializeRandomlyAt(p, kRandomNumberGenerator); + return initializeRandomlyAt(p, kPRNG); } /* ************************************************************************* */ diff --git a/gtsam/sfm/TranslationRecovery.cpp b/gtsam/sfm/TranslationRecovery.cpp index 02c78133e..828aa6880 100644 --- a/gtsam/sfm/TranslationRecovery.cpp +++ b/gtsam/sfm/TranslationRecovery.cpp @@ -39,7 +39,7 @@ using namespace gtsam; using namespace std; // In Wrappers we have no access to this so have a default ready. -static std::mt19937 kRandomNumberGenerator(42); +static std::mt19937 kPRNG(42); // Some relative translations may be zero. We treat nodes that have a zero // relativeTranslation as a single node. @@ -185,7 +185,7 @@ Values TranslationRecovery::initializeRandomly( const std::vector> &betweenTranslations, const Values &initialValues) const { return initializeRandomly(relativeTranslations, betweenTranslations, - &kRandomNumberGenerator, initialValues); + &kPRNG, initialValues); } Values TranslationRecovery::run( diff --git a/gtsam/sfm/tests/testShonanAveraging.cpp b/gtsam/sfm/tests/testShonanAveraging.cpp index e297678b8..14667269a 100644 --- a/gtsam/sfm/tests/testShonanAveraging.cpp +++ b/gtsam/sfm/tests/testShonanAveraging.cpp @@ -45,7 +45,7 @@ ShonanAveraging3 fromExampleName( static const ShonanAveraging3 kShonan = fromExampleName("toyExample.g2o"); -static std::mt19937 kRandomNumberGenerator(42); +static std::mt19937 kPRNG(42); /* ************************************************************************* */ TEST(ShonanAveraging3, checkConstructor) { @@ -78,7 +78,7 @@ TEST(ShonanAveraging3, buildGraphAt) { /* ************************************************************************* */ TEST(ShonanAveraging3, checkOptimality) { - const Values randomRotations = kShonan.initializeRandomly(kRandomNumberGenerator); + const Values randomRotations = kShonan.initializeRandomly(kPRNG); Values random = ShonanAveraging3::LiftTo(4, randomRotations); // lift to 4! auto Lambda = kShonan.computeLambda(random); EXPECT_LONGS_EQUAL(15, Lambda.rows()); @@ -106,7 +106,7 @@ TEST(ShonanAveraging3, checkSubgraph) { // Create initial random estimation Values initial; - initial = subgraphShonan.initializeRandomly(kRandomNumberGenerator); + initial = subgraphShonan.initializeRandomly(kPRNG); // Run Shonan with SUBGRAPH solver auto result = subgraphShonan.run(initial, 3, 3); @@ -115,7 +115,7 @@ TEST(ShonanAveraging3, checkSubgraph) { /* ************************************************************************* */ TEST(ShonanAveraging3, tryOptimizingAt3) { - const Values randomRotations = kShonan.initializeRandomly(kRandomNumberGenerator); + const Values randomRotations = kShonan.initializeRandomly(kPRNG); Values initial = ShonanAveraging3::LiftTo(3, randomRotations); // convert to SOn EXPECT(!kShonan.checkOptimality(initial)); const Values result = kShonan.tryOptimizingAt(3, initial); @@ -130,7 +130,7 @@ TEST(ShonanAveraging3, tryOptimizingAt3) { /* ************************************************************************* */ TEST(ShonanAveraging3, tryOptimizingAt4) { - const Values randomRotations = kShonan.initializeRandomly(kRandomNumberGenerator); + const Values randomRotations = kShonan.initializeRandomly(kPRNG); Values random = ShonanAveraging3::LiftTo(4, randomRotations); // lift to 4! const Values result = kShonan.tryOptimizingAt(4, random); EXPECT(kShonan.checkOptimality(result)); @@ -228,7 +228,7 @@ TEST(ShonanAveraging3, CheckWithEigen) { /* ************************************************************************* */ TEST(ShonanAveraging3, initializeWithDescent) { - const Values randomRotations = kShonan.initializeRandomly(kRandomNumberGenerator); + const Values randomRotations = kShonan.initializeRandomly(kPRNG); Values random = ShonanAveraging3::LiftTo(3, randomRotations); const Values Qstar3 = kShonan.tryOptimizingAt(3, random); Vector minEigenVector; @@ -240,7 +240,7 @@ TEST(ShonanAveraging3, initializeWithDescent) { /* ************************************************************************* */ TEST(ShonanAveraging3, run) { - auto initial = kShonan.initializeRandomly(kRandomNumberGenerator); + auto initial = kShonan.initializeRandomly(kPRNG); auto result = kShonan.run(initial, 5); EXPECT_DOUBLES_EQUAL(0, kShonan.cost(result.first), 1e-3); EXPECT_DOUBLES_EQUAL(-5.427688831332745e-07, result.second, @@ -295,7 +295,7 @@ TEST(ShonanAveraging3, runKlaus) { EXPECT(assert_equal(R02, wR0.between(wR2), 0.1)); // Run Shonan (with prior on first rotation) - auto initial = shonan.initializeRandomly(kRandomNumberGenerator); + auto initial = shonan.initializeRandomly(kPRNG); auto result = shonan.run(initial, 5); EXPECT_DOUBLES_EQUAL(0, shonan.cost(result.first), 1e-2); EXPECT_DOUBLES_EQUAL(-9.2259161494467889e-05, result.second, @@ -323,7 +323,7 @@ TEST(ShonanAveraging3, runKlausKarcher) { static const ShonanAveraging3 shonan = fromExampleName("Klaus3.g2o"); // Run Shonan (with Karcher mean prior) - auto initial = shonan.initializeRandomly(kRandomNumberGenerator); + auto initial = shonan.initializeRandomly(kPRNG); auto result = shonan.run(initial, 5); EXPECT_DOUBLES_EQUAL(0, shonan.cost(result.first), 1e-2); EXPECT_DOUBLES_EQUAL(-1.361402670507772e-05, result.second, @@ -353,7 +353,7 @@ TEST(ShonanAveraging2, noisyToyGraph) { // Check graph building NonlinearFactorGraph graph = shonan.buildGraphAt(2); EXPECT_LONGS_EQUAL(6, graph.size()); - auto initial = shonan.initializeRandomly(kRandomNumberGenerator); + auto initial = shonan.initializeRandomly(kPRNG); auto result = shonan.run(initial, 2); EXPECT_DOUBLES_EQUAL(0.0008211, shonan.cost(result.first), 1e-6); EXPECT_DOUBLES_EQUAL(0, result.second, 1e-10); // certificate! @@ -391,7 +391,7 @@ TEST(ShonanAveraging2, noisyToyGraphWithHuber) { } // test result - auto initial = shonan.initializeRandomly(kRandomNumberGenerator); + auto initial = shonan.initializeRandomly(kPRNG); auto result = shonan.run(initial, 2,2); EXPECT_DOUBLES_EQUAL(0.0008211, shonan.cost(result.first), 1e-6); EXPECT_DOUBLES_EQUAL(0, result.second, 1e-10); // certificate!