Reduce the number of individual PRNG instances Use same instance of a PRNG throughout by passing it to methods and functions as an argument to generate random numbers without breaking the sequence. Change-Id: Ib024bbc1ea2d14e4b9afb71857856a5fb77b1667
diff --git a/internal/ceres/autodiff_test.cc b/internal/ceres/autodiff_test.cc index 629eea8..0916a45 100644 --- a/internal/ceres/autodiff_test.cc +++ b/internal/ceres/autodiff_test.cc
@@ -30,6 +30,8 @@ #include "ceres/internal/autodiff.h" +#include <algorithm> +#include <iterator> #include <random> #include "gtest/gtest.h" @@ -173,9 +175,9 @@ // Make random P and X, in a single vector. double PX[12 + 4]; - for (double& PX_i : PX) { - PX_i = uniform01(prng); - } + std::generate(std::begin(PX), std::end(PX), [&prng, &uniform01] { + return uniform01(prng); + }); // Handy names for the P and X parts. double* P = PX + 0; @@ -294,7 +296,10 @@ double qcX[4 + 3 + 3]; std::mt19937 prng; std::uniform_real_distribution<double> uniform01(0.0, 1.0); - for (double& qcX_i : qcX) qcX_i = uniform01(prng); + + std::generate(std::begin(qcX), std::end(qcX), [&prng, &uniform01] { + return uniform01(prng); + }); // Handy names. double* q = qcX;
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index ae6bd3a..90bd926 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc
@@ -375,7 +375,7 @@ } std::unique_ptr<BlockSparseMatrix> BlockSparseMatrix::CreateRandomMatrix( - const BlockSparseMatrix::RandomMatrixOptions& options) { + const BlockSparseMatrix::RandomMatrixOptions& options, std::mt19937& prng) { CHECK_GT(options.num_row_blocks, 0); CHECK_GT(options.min_row_block_size, 0); CHECK_GT(options.max_row_block_size, 0); @@ -383,7 +383,6 @@ CHECK_GT(options.block_density, 0.0); CHECK_LE(options.block_density, 1.0); - std::mt19937 prng; std::uniform_int_distribution<int> col_distribution( options.min_col_block_size, options.max_col_block_size); std::uniform_int_distribution<int> row_distribution( @@ -436,9 +435,10 @@ auto matrix = std::make_unique<BlockSparseMatrix>(bs); double* values = matrix->mutable_values(); std::normal_distribution<double> standard_normal_distribution; - for (int i = 0; i < matrix->num_nonzeros(); ++i) { - values[i] = standard_normal_distribution(prng); - } + std::generate_n( + values, matrix->num_nonzeros(), [&standard_normal_distribution, &prng] { + return standard_normal_distribution(prng); + }); return matrix; }
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h index 7cef18d..44cb5e5 100644 --- a/internal/ceres/block_sparse_matrix.h +++ b/internal/ceres/block_sparse_matrix.h
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -35,6 +35,7 @@ #define CERES_INTERNAL_BLOCK_SPARSE_MATRIX_H_ #include <memory> +#include <random> #include "ceres/block_structure.h" #include "ceres/crs_matrix.h" @@ -123,7 +124,7 @@ // distributed and whose structure is determined by // RandomMatrixOptions. static std::unique_ptr<BlockSparseMatrix> CreateRandomMatrix( - const RandomMatrixOptions& options); + const RandomMatrixOptions& options, std::mt19937& prng); private: int num_rows_;
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc index 574b0c0..ce587da 100644 --- a/internal/ceres/compressed_row_sparse_matrix.cc +++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -119,11 +119,12 @@ transpose_rows[0] = 0; } +template <class RandomNormalFunctor> void AddRandomBlock(const int num_rows, const int num_cols, const int row_block_begin, const int col_block_begin, - std::function<double()> dist, + RandomNormalFunctor&& randn, std::vector<int>* rows, std::vector<int>* cols, std::vector<double>* values) { @@ -131,20 +132,21 @@ for (int c = 0; c < num_cols; ++c) { rows->push_back(row_block_begin + r); cols->push_back(col_block_begin + c); - values->push_back(dist()); + values->push_back(randn()); } } } +template <class RandomNormalFunctor> void AddSymmetricRandomBlock(const int num_rows, const int row_block_begin, - std::function<double()> dist, + RandomNormalFunctor&& randn, std::vector<int>* rows, std::vector<int>* cols, std::vector<double>* values) { for (int r = 0; r < num_rows; ++r) { for (int c = r; c < num_rows; ++c) { - const double v = dist(); + const double v = randn(); rows->push_back(row_block_begin + r); cols->push_back(row_block_begin + c); values->push_back(v); @@ -625,7 +627,8 @@ std::unique_ptr<CompressedRowSparseMatrix> CompressedRowSparseMatrix::CreateRandomMatrix( - CompressedRowSparseMatrix::RandomMatrixOptions options) { + CompressedRowSparseMatrix::RandomMatrixOptions options, + std::mt19937& prng) { CHECK_GT(options.num_row_blocks, 0); CHECK_GT(options.min_row_block_size, 0); CHECK_GT(options.max_row_block_size, 0); @@ -651,7 +654,6 @@ vector<int> col_blocks; col_blocks.reserve(options.num_col_blocks); - std::mt19937 prng; std::uniform_int_distribution<int> col_distribution( options.min_col_block_size, options.max_col_block_size); std::uniform_int_distribution<int> row_distribution( @@ -708,6 +710,9 @@ // Randomly determine if this block is present or not. if (uniform01(prng) <= options.block_density) { + auto randn = [&standard_normal_distribution, &prng] { + return standard_normal_distribution(prng); + }; // If the matrix is symmetric, then we take care to generate // symmetric diagonal blocks. if (options.storage_type == StorageType::UNSYMMETRIC || r != c) { @@ -715,14 +720,14 @@ col_blocks[c], row_block_begin, col_block_begin, - values_dist, + randn, &tsm_rows, &tsm_cols, &tsm_values); } else { AddSymmetricRandomBlock(row_blocks[r], row_block_begin, - values_dist, + randn, &tsm_rows, &tsm_cols, &tsm_values);
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h index 2580045..d164bef 100644 --- a/internal/ceres/compressed_row_sparse_matrix.h +++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -32,6 +32,7 @@ #define CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_ #include <memory> +#include <random> #include <vector> #include "ceres/internal/disable_warnings.h" @@ -195,7 +196,7 @@ // normally distributed and whose structure is determined by // RandomMatrixOptions. static std::unique_ptr<CompressedRowSparseMatrix> CreateRandomMatrix( - RandomMatrixOptions options); + RandomMatrixOptions options, std::mt19937& prng); private: static std::unique_ptr<CompressedRowSparseMatrix> FromTripletSparseMatrix(
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc index 42f5498..718a659 100644 --- a/internal/ceres/compressed_row_sparse_matrix_test.cc +++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -328,6 +328,7 @@ } TEST(CompressedRowSparseMatrix, FromTripletSparseMatrix) { + std::mt19937 prng; TripletSparseMatrix::RandomMatrixOptions options; options.num_rows = 5; options.num_cols = 7; @@ -336,7 +337,7 @@ const int kNumTrials = 10; for (int i = 0; i < kNumTrials; ++i) { std::unique_ptr<TripletSparseMatrix> tsm = - TripletSparseMatrix::CreateRandomMatrix(options); + TripletSparseMatrix::CreateRandomMatrix(options, prng); std::unique_ptr<CompressedRowSparseMatrix> crsm = CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm); @@ -354,6 +355,7 @@ } TEST(CompressedRowSparseMatrix, FromTripletSparseMatrixTransposed) { + std::mt19937 prng; TripletSparseMatrix::RandomMatrixOptions options; options.num_rows = 5; options.num_cols = 7; @@ -362,7 +364,7 @@ const int kNumTrials = 10; for (int i = 0; i < kNumTrials; ++i) { std::unique_ptr<TripletSparseMatrix> tsm = - TripletSparseMatrix::CreateRandomMatrix(options); + TripletSparseMatrix::CreateRandomMatrix(options, prng); std::unique_ptr<CompressedRowSparseMatrix> crsm = CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm); @@ -421,7 +423,7 @@ options.block_density = uniform(prng); options.storage_type = ::testing::get<0>(param); std::unique_ptr<CompressedRowSparseMatrix> matrix = - CompressedRowSparseMatrix::CreateRandomMatrix(options); + CompressedRowSparseMatrix::CreateRandomMatrix(options, prng); const int num_rows = matrix->num_rows(); const int num_cols = matrix->num_cols(); @@ -491,7 +493,7 @@ options.block_density = uniform(prng); options.storage_type = ::testing::get<0>(param); std::unique_ptr<CompressedRowSparseMatrix> matrix = - CompressedRowSparseMatrix::CreateRandomMatrix(options); + CompressedRowSparseMatrix::CreateRandomMatrix(options, prng); const int num_rows = matrix->num_rows(); const int num_cols = matrix->num_cols(); @@ -561,7 +563,7 @@ options.block_density = uniform(prng); options.storage_type = ::testing::get<0>(param); std::unique_ptr<CompressedRowSparseMatrix> matrix = - CompressedRowSparseMatrix::CreateRandomMatrix(options); + CompressedRowSparseMatrix::CreateRandomMatrix(options, prng); const int num_cols = matrix->num_cols(); Vector actual(num_cols);
diff --git a/internal/ceres/corrector_test.cc b/internal/ceres/corrector_test.cc index 7f58efa..a5dd29c 100644 --- a/internal/ceres/corrector_test.cc +++ b/internal/ceres/corrector_test.cc
@@ -32,7 +32,6 @@ #include <algorithm> #include <cmath> -#include <cstdlib> #include <cstring> #include <random>
diff --git a/internal/ceres/gradient_checker_test.cc b/internal/ceres/gradient_checker_test.cc index 3c5f5da..8d7d10f 100644 --- a/internal/ceres/gradient_checker_test.cc +++ b/internal/ceres/gradient_checker_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2021 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -33,7 +33,6 @@ #include "ceres/gradient_checker.h" #include <cmath> -#include <cstdlib> #include <random> #include <utility> #include <vector> @@ -59,16 +58,16 @@ // version, they are both block vectors, of course. class GoodTestTerm : public CostFunction { public: - GoodTestTerm(int arity, int const* dim) : arity_(arity), return_value_(true) { - std::mt19937 prng; - std::uniform_real_distribution<double> distribution(-1.0, 1.0); - + template <class UniformRandomFunctor> + GoodTestTerm(int arity, int const* dim, UniformRandomFunctor&& randu) + : arity_(arity), return_value_(true) { + std::uniform_real_distribution distribution(-1.0, 1.0); // Make 'arity' random vectors. a_.resize(arity_); for (int j = 0; j < arity_; ++j) { a_[j].resize(dim[j]); for (int u = 0; u < dim[j]; ++u) { - a_[j][u] = distribution(prng); + a_[j][u] = randu(); } } @@ -121,15 +120,15 @@ class BadTestTerm : public CostFunction { public: - BadTestTerm(int arity, int const* dim) : arity_(arity) { - std::mt19937 prng; - std::uniform_real_distribution<double> distribution(-1.0, 1.0); + template <class UniformRandomFunctor> + BadTestTerm(int arity, int const* dim, UniformRandomFunctor&& randu) + : arity_(arity) { // Make 'arity' random vectors. a_.resize(arity_); for (int j = 0; j < arity_; ++j) { a_[j].resize(dim[j]); for (int u = 0; u < dim[j]; ++u) { - a_[j][u] = distribution(prng); + a_[j][u] = randu(); } } @@ -210,10 +209,11 @@ FixedArray<double*> parameters(num_parameters); std::mt19937 prng; std::uniform_real_distribution<double> distribution(-1.0, 1.0); + auto randu = [&prng, &distribution] { return distribution(prng); }; for (int j = 0; j < num_parameters; ++j) { parameters[j] = new double[parameter_sizes[j]]; for (int u = 0; u < parameter_sizes[j]; ++u) { - parameters[j][u] = distribution(prng); + parameters[j][u] = randu(); } } @@ -221,7 +221,7 @@ GradientChecker::ProbeResults results; // Test that Probe returns true for correct Jacobians. - GoodTestTerm good_term(num_parameters, parameter_sizes.data()); + GoodTestTerm good_term(num_parameters, parameter_sizes.data(), randu); std::vector<const Manifold*>* manifolds = nullptr; GradientChecker good_gradient_checker( &good_term, manifolds, numeric_diff_options); @@ -258,7 +258,7 @@ EXPECT_FALSE(results.error_log.empty()); // Test that Probe returns false for incorrect Jacobians. - BadTestTerm bad_term(num_parameters, parameter_sizes.data()); + BadTestTerm bad_term(num_parameters, parameter_sizes.data(), randu); GradientChecker bad_gradient_checker( &bad_term, manifolds, numeric_diff_options); EXPECT_FALSE(
diff --git a/internal/ceres/gradient_checking_cost_function_test.cc b/internal/ceres/gradient_checking_cost_function_test.cc index caf4f25..87c3278 100644 --- a/internal/ceres/gradient_checking_cost_function_test.cc +++ b/internal/ceres/gradient_checking_cost_function_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -69,15 +69,15 @@ public: // The constructor of this function needs to know the number // of blocks desired, and the size of each block. - TestTerm(int arity, int const* dim) : arity_(arity) { - std::mt19937 prng; - std::uniform_real_distribution<double> uniform01(-1.0, 1.0); + template <class UniformRandomFunctor> + TestTerm(int arity, int const* dim, UniformRandomFunctor&& randu) + : arity_(arity) { // Make 'arity' random vectors. a_.resize(arity_); for (int j = 0; j < arity_; ++j) { a_[j].resize(dim[j]); for (int u = 0; u < dim[j]; ++u) { - a_[j][u] = uniform01(prng); + a_[j][u] = randu(); } } @@ -139,11 +139,12 @@ // Make a random set of blocks. vector<double*> parameters(arity); std::mt19937 prng; - std::uniform_real_distribution<double> uniform01(-1.0, 1.0); + std::uniform_real_distribution<double> distribution(-1.0, 1.0); + auto randu = [&prng, &distribution] { return distribution(prng); }; for (int j = 0; j < arity; ++j) { parameters[j] = new double[dim[j]]; for (int u = 0; u < dim[j]; ++u) { - parameters[j][u] = uniform01(prng); + parameters[j][u] = randu(); } } @@ -162,7 +163,7 @@ const double kRelativeStepSize = 1e-6; const double kRelativePrecision = 1e-4; - TestTerm<-1, -1> term(arity, dim); + TestTerm<-1, -1> term(arity, dim, randu); GradientCheckingIterationCallback callback; auto gradient_checking_cost_function = CreateGradientCheckingCostFunction(&term, @@ -196,11 +197,12 @@ // Make a random set of blocks. vector<double*> parameters(arity); std::mt19937 prng; - std::uniform_real_distribution<double> uniform01(-1.0, 1.0); + std::uniform_real_distribution<double> distribution(-1.0, 1.0); + auto randu = [&prng, &distribution] { return distribution(prng); }; for (int j = 0; j < arity; ++j) { parameters[j] = new double[dim[j]]; for (int u = 0; u < dim[j]; ++u) { - parameters[j][u] = uniform01(prng); + parameters[j][u] = randu(); } } @@ -218,7 +220,7 @@ // Should have one term that's bad, causing everything to get dumped. LOG(INFO) << "Bad gradient"; { - TestTerm<1, 2> term(arity, dim); + TestTerm<1, 2> term(arity, dim, randu); GradientCheckingIterationCallback callback; auto gradient_checking_cost_function = CreateGradientCheckingCostFunction(&term, @@ -238,7 +240,7 @@ // The gradient is correct, so no errors are reported. LOG(INFO) << "Good gradient"; { - TestTerm<-1, -1> term(arity, dim); + TestTerm<-1, -1> term(arity, dim, randu); GradientCheckingIterationCallback callback; auto gradient_checking_cost_function = CreateGradientCheckingCostFunction(&term,
diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc index 5bd4c50..4573b6b 100644 --- a/internal/ceres/inner_product_computer_test.cc +++ b/internal/ceres/inner_product_computer_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2017 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -109,7 +109,7 @@ VLOG(2) << "block density: " << options.block_density; std::unique_ptr<BlockSparseMatrix> random_matrix( - BlockSparseMatrix::CreateRandomMatrix(options)); + BlockSparseMatrix::CreateRandomMatrix(options, prng)); TripletSparseMatrix tsm(random_matrix->num_rows(), random_matrix->num_cols(), @@ -169,7 +169,7 @@ VLOG(2) << "block density: " << options.block_density; std::unique_ptr<BlockSparseMatrix> random_matrix( - BlockSparseMatrix::CreateRandomMatrix(options)); + BlockSparseMatrix::CreateRandomMatrix(options, prng)); const std::vector<CompressedRow>& row_blocks = random_matrix->block_structure()->rows;
diff --git a/internal/ceres/normal_prior_test.cc b/internal/ceres/normal_prior_test.cc index 25e6322..881038f 100644 --- a/internal/ceres/normal_prior_test.cc +++ b/internal/ceres/normal_prior_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -42,6 +42,7 @@ TEST(NormalPriorTest, ResidualAtRandomPosition) { std::mt19937 prng; std::uniform_real_distribution<double> distribution(-1.0, 1.0); + auto randu = [&distribution, &prng] { return distribution(prng); }; for (int num_rows = 1; num_rows < 5; ++num_rows) { for (int num_cols = 1; num_cols < 5; ++num_cols) { Vector b(num_cols); @@ -50,7 +51,7 @@ A.setRandom(); auto* x = new double[num_cols]; - for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng); + std::generate_n(x, num_cols, randu); auto* jacobian = new double[num_rows * num_cols]; Vector residuals(num_rows); @@ -77,6 +78,7 @@ TEST(NormalPriorTest, ResidualAtRandomPositionNullJacobians) { std::mt19937 prng; std::uniform_real_distribution<double> distribution(-1.0, 1.0); + auto randu = [&distribution, &prng] { return distribution(prng); }; for (int num_rows = 1; num_rows < 5; ++num_rows) { for (int num_cols = 1; num_cols < 5; ++num_cols) { Vector b(num_cols); @@ -85,7 +87,7 @@ A.setRandom(); auto* x = new double[num_cols]; - for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng); + std::generate_n(x, num_cols, randu); double* jacobians[1]; jacobians[0] = nullptr;
diff --git a/internal/ceres/numeric_diff_cost_function_test.cc b/internal/ceres/numeric_diff_cost_function_test.cc index b9a0beb..dcffcff 100644 --- a/internal/ceres/numeric_diff_cost_function_test.cc +++ b/internal/ceres/numeric_diff_cost_function_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -35,6 +35,7 @@ #include <array> #include <cmath> #include <memory> +#include <random> #include <string> #include <vector> @@ -310,6 +311,7 @@ } TEST(NumericDiffCostFunction, RandomizedFunctorRidders) { + std::mt19937 prng; NumericDiffOptions options; // Larger initial step size is chosen to produce robust results in the // presence of random noise. @@ -321,15 +323,16 @@ 1, // number of residuals 1 // size of x1 >>( - new RandomizedFunctor(kNoiseFactor, kRandomSeed), + new RandomizedFunctor(kNoiseFactor, prng), TAKE_OWNERSHIP, 1, options); - RandomizedFunctor functor(kNoiseFactor, kRandomSeed); + RandomizedFunctor functor(kNoiseFactor, prng); functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function); } TEST(NumericDiffCostFunction, RandomizedCostFunctionRidders) { + std::mt19937 prng; NumericDiffOptions options; // Larger initial step size is chosen to produce robust results in the // presence of random noise. @@ -341,12 +344,12 @@ 1, // number of residuals 1 // size of x1 >>( - new RandomizedCostFunction(kNoiseFactor, kRandomSeed), + new RandomizedCostFunction(kNoiseFactor, prng), TAKE_OWNERSHIP, 1, options); - RandomizedFunctor functor(kNoiseFactor, kRandomSeed); + RandomizedFunctor functor(kNoiseFactor, prng); functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function); }
diff --git a/internal/ceres/numeric_diff_test_utils.cc b/internal/ceres/numeric_diff_test_utils.cc index 3d78e36..9da514e 100644 --- a/internal/ceres/numeric_diff_test_utils.cc +++ b/internal/ceres/numeric_diff_test_utils.cc
@@ -226,14 +226,7 @@ } bool RandomizedFunctor::operator()(const double* x1, double* residuals) const { - double random_value = - static_cast<double>(rand()) / static_cast<double>(RAND_MAX); - - // Normalize noise to [-factor, factor]. - random_value *= 2.0; - random_value -= 1.0; - random_value *= noise_factor_; - + double random_value = uniform_distribution_(*prng_); residuals[0] = x1[0] * x1[0] + random_value; return true; } @@ -244,9 +237,6 @@ const double kTolerance = 2e-4; - // Initialize random number generator with given seed. - srand(random_seed_); - for (double& test : kTests) { double* parameters[] = {&test}; double dydx;
diff --git a/internal/ceres/numeric_diff_test_utils.h b/internal/ceres/numeric_diff_test_utils.h index 68ae5ed..8d3f14b 100644 --- a/internal/ceres/numeric_diff_test_utils.h +++ b/internal/ceres/numeric_diff_test_utils.h
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -31,6 +31,8 @@ #ifndef CERES_INTERNAL_NUMERIC_DIFF_TEST_UTILS_H_ #define CERES_INTERNAL_NUMERIC_DIFF_TEST_UTILS_H_ +#include <random> + #include "ceres/cost_function.h" #include "ceres/internal/export.h" #include "ceres/sized_cost_function.h" @@ -116,8 +118,10 @@ // y = x^2 + [random noise], dy/dx ~ 2x class CERES_NO_EXPORT RandomizedFunctor { public: - RandomizedFunctor(double noise_factor, unsigned int random_seed) - : noise_factor_(noise_factor), random_seed_(random_seed) {} + RandomizedFunctor(double noise_factor, std::mt19937& prng) + : noise_factor_(noise_factor), + prng_(&prng), + uniform_distribution_{-noise_factor, noise_factor} {} bool operator()(const double* x1, double* residuals) const; void ExpectCostFunctionEvaluationIsNearlyCorrect( @@ -125,13 +129,16 @@ private: double noise_factor_; - unsigned int random_seed_; + // Store the generator as a pointer to be able to modify the instance the + // pointer is pointing to. + std::mt19937* prng_; + mutable std::uniform_real_distribution<> uniform_distribution_; }; class CERES_NO_EXPORT RandomizedCostFunction : public SizedCostFunction<1, 1> { public: - RandomizedCostFunction(double noise_factor, unsigned int random_seed) - : functor_(noise_factor, random_seed) {} + RandomizedCostFunction(double noise_factor, std::mt19937& prng) + : functor_(noise_factor, prng) {} bool Evaluate(double const* const* parameters, double* residuals,
diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc index cb3dd14..e4484dd 100644 --- a/internal/ceres/partitioned_matrix_view_test.cc +++ b/internal/ceres/partitioned_matrix_view_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without
diff --git a/internal/ceres/reorder_program_test.cc b/internal/ceres/reorder_program_test.cc index 872f381..d6a6577 100644 --- a/internal/ceres/reorder_program_test.cc +++ b/internal/ceres/reorder_program_test.cc
@@ -265,7 +265,7 @@ problem.GetResidualBlocks(&residual_block_ids); std::vector<ResidualBlock*> residual_blocks = problem.program().residual_blocks(); - auto rng = std::default_random_engine{}; + auto rng = std::mt19937{}; for (int i = 1; i < 6; ++i) { std::shuffle( std::begin(residual_block_ids), std::end(residual_block_ids), rng);
diff --git a/internal/ceres/rotation_test.cc b/internal/ceres/rotation_test.cc index d6de650..97ea514 100644 --- a/internal/ceres/rotation_test.cc +++ b/internal/ceres/rotation_test.cc
@@ -32,6 +32,7 @@ #include <cmath> #include <limits> +#include <random> #include <string> #include "ceres/internal/eigen.h" @@ -56,11 +57,6 @@ const double kPi = 3.14159265358979323846; const double kHalfSqrt2 = 0.707106781186547524401; -static double RandDouble() { - double r = rand(); - return r / RAND_MAX; -} - // A tolerance value for floating-point comparisons. static double const kTolerance = numeric_limits<double>::epsilon() * 10; @@ -322,20 +318,22 @@ // Takes a bunch of random axis/angle values, converts them to quaternions, // and back again. TEST(Rotation, AngleAxisToQuaterionAndBack) { - srand(5); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; for (int i = 0; i < kNumTrials; i++) { double axis_angle[3]; // Make an axis by choosing three random numbers in [-1, 1) and // normalizing. double norm = 0; for (double& coeff : axis_angle) { - coeff = RandDouble() * 2 - 1; + coeff = uniform_distribution(prng); norm += coeff * coeff; } norm = sqrt(norm); // Angle in [-pi, pi). - double theta = kPi * 2 * RandDouble() - kPi; + double theta = uniform_distribution( + prng, std::uniform_real_distribution<double>::param_type{-kPi, kPi}); for (double& coeff : axis_angle) { coeff = coeff * theta / norm; } @@ -355,13 +353,14 @@ // Takes a bunch of random quaternions, converts them to axis/angle, // and back again. TEST(Rotation, QuaterionToAngleAxisAndBack) { - srand(5); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; for (int i = 0; i < kNumTrials; i++) { double quaternion[4]; // Choose four random numbers in [-1, 1) and normalize. double norm = 0; for (double& coeff : quaternion) { - coeff = RandDouble() * 2 - 1; + coeff = uniform_distribution(prng); norm += coeff * coeff; } norm = sqrt(norm); @@ -432,20 +431,24 @@ double matrix[9]; double out_axis_angle[3]; - srand(5); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; for (int i = 0; i < kNumTrials; i++) { // Make an axis by choosing three random numbers in [-1, 1) and // normalizing. double norm = 0; for (double& coeff : in_axis_angle) { - coeff = RandDouble() * 2 - 1; + coeff = uniform_distribution(prng); norm += coeff * coeff; } norm = sqrt(norm); // Angle in [pi - kMaxSmallAngle, pi). - const double kMaxSmallAngle = 1e-8; - double theta = kPi - kMaxSmallAngle * RandDouble(); + constexpr double kMaxSmallAngle = 1e-8; + double theta = + uniform_distribution(prng, + std::uniform_real_distribution<double>::param_type{ + kPi - kMaxSmallAngle, kPi}); for (double& coeff : in_axis_angle) { coeff *= (theta / norm); @@ -531,20 +534,22 @@ // Takes a bunch of random axis/angle values, converts them to rotation // matrices, and back again. TEST(Rotation, AngleAxisToRotationMatrixAndBack) { - srand(5); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; for (int i = 0; i < kNumTrials; i++) { double axis_angle[3]; // Make an axis by choosing three random numbers in [-1, 1) and // normalizing. double norm = 0; for (double& i : axis_angle) { - i = RandDouble() * 2 - 1; + i = uniform_distribution(prng); norm += i * i; } norm = sqrt(norm); // Angle in [-pi, pi). - double theta = kPi * 2 * RandDouble() - kPi; + double theta = uniform_distribution( + prng, std::uniform_real_distribution<double>::param_type{-kPi, kPi}); for (double& i : axis_angle) { i = i * theta / norm; } @@ -564,20 +569,25 @@ // Takes a bunch of random axis/angle values near zero, converts them // to rotation matrices, and back again. TEST(Rotation, AngleAxisToRotationMatrixAndBackNearZero) { - srand(5); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; for (int i = 0; i < kNumTrials; i++) { double axis_angle[3]; // Make an axis by choosing three random numbers in [-1, 1) and // normalizing. double norm = 0; for (double& i : axis_angle) { - i = RandDouble() * 2 - 1; + i = uniform_distribution(prng); norm += i * i; } norm = sqrt(norm); // Tiny theta. - double theta = 1e-16 * (kPi * 2 * RandDouble() - kPi); + constexpr double kScale = 1e-16; + double theta = + uniform_distribution(prng, + std::uniform_real_distribution<double>::param_type{ + -kScale * kPi, kScale * kPi}); for (double& i : axis_angle) { i = i * theta / norm; } @@ -647,11 +657,12 @@ // Test that a random rotation produces an orthonormal rotation // matrix. TEST(EulerAnglesToRotationMatrix, IsOrthonormal) { - srand(5); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-180.0, 180.0}; for (int trial = 0; trial < kNumTrials; ++trial) { double euler_angles_degrees[3]; for (double& euler_angles_degree : euler_angles_degrees) { - euler_angles_degree = RandDouble() * 360.0 - 180.0; + euler_angles_degree = uniform_distribution(prng); } double rotation_matrix[9]; EulerAnglesToRotationMatrix(euler_angles_degrees, 3, rotation_matrix); @@ -906,13 +917,15 @@ // Verify that (a * b) * c == a * (b * c). TEST(Quaternion, MultiplicationIsAssociative) { + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; double a[4]; double b[4]; double c[4]; for (int i = 0; i < 4; ++i) { - a[i] = 2 * RandDouble() - 1; - b[i] = 2 * RandDouble() - 1; - c[i] = 2 * RandDouble() - 1; + a[i] = uniform_distribution(prng); + b[i] = uniform_distribution(prng); + c[i] = uniform_distribution(prng); } double ab[4]; @@ -932,6 +945,8 @@ } TEST(AngleAxis, RotatePointGivesSameAnswerAsRotationMatrix) { + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; double angle_axis[3]; double R[9]; double p[3]; @@ -943,8 +958,8 @@ for (int j = 0; j < 50; ++j) { double norm2 = 0.0; for (int k = 0; k < 3; ++k) { - angle_axis[k] = 2.0 * RandDouble() - 1.0; - p[k] = 2.0 * RandDouble() - 1.0; + angle_axis[k] = uniform_distribution(prng); + p[k] = uniform_distribution(prng); norm2 = angle_axis[k] * angle_axis[k]; } @@ -976,6 +991,8 @@ } TEST(AngleAxis, NearZeroRotatePointGivesSameAnswerAsRotationMatrix) { + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0}; double angle_axis[3]; double R[9]; double p[3]; @@ -985,8 +1002,8 @@ for (int i = 0; i < 10000; ++i) { double norm2 = 0.0; for (int k = 0; k < 3; ++k) { - angle_axis[k] = 2.0 * RandDouble() - 1.0; - p[k] = 2.0 * RandDouble() - 1.0; + angle_axis[k] = uniform_distribution(prng); + p[k] = uniform_distribution(prng); norm2 = angle_axis[k] * angle_axis[k]; } @@ -1111,7 +1128,12 @@ } TEST(RotationMatrixToAngleAxis, ExhaustiveRoundTrip) { - const double kMaxSmallAngle = 1e-8; + constexpr double kMaxSmallAngle = 1e-8; + std::mt19937 prng; + std::uniform_real_distribution<double> uniform_distribution1{ + kPi - kMaxSmallAngle, kPi}; + std::uniform_real_distribution<double> uniform_distribution2{ + -1.0, 2.0 * kMaxSmallAngle - 1.0}; const int kNumSteps = 1000; for (int i = 0; i < kNumSteps; ++i) { const double theta = static_cast<double>(i) / kNumSteps * 2.0 * kPi; @@ -1121,10 +1143,10 @@ CheckRotationMatrixToAngleAxisRoundTrip(theta, phi, kPi); // Rotation of angle approximately Pi. CheckRotationMatrixToAngleAxisRoundTrip( - theta, phi, kPi - kMaxSmallAngle * RandDouble()); + theta, phi, uniform_distribution1(prng)); // Rotations of angle approximately zero. CheckRotationMatrixToAngleAxisRoundTrip( - theta, phi, kMaxSmallAngle * 2.0 * RandDouble() - 1.0); + theta, phi, uniform_distribution2(prng)); } } }
diff --git a/internal/ceres/schur_eliminator_benchmark.cc b/internal/ceres/schur_eliminator_benchmark.cc index 0973e73..c6425ce 100644 --- a/internal/ceres/schur_eliminator_benchmark.cc +++ b/internal/ceres/schur_eliminator_benchmark.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2019 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -91,11 +91,9 @@ matrix_ = std::make_unique<BlockSparseMatrix>(bs); double* values = matrix_->mutable_values(); - std::mt19937 prng; - std::normal_distribution<double> standard_normal; - for (int i = 0; i < matrix_->num_nonzeros(); ++i) { - values[i] = standard_normal(prng); - } + std::generate_n(values, matrix_->num_nonzeros(), [this] { + return standard_normal(prng); + }); b_.resize(matrix_->num_rows()); b_.setRandom(); @@ -128,6 +126,8 @@ Vector diagonal_; Vector z_; Vector y_; + std::mt19937 prng; + std::normal_distribution<> standard_normal; }; static void BM_SchurEliminatorEliminate(benchmark::State& state) {
diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc index da27019..17ce67d 100644 --- a/internal/ceres/schur_eliminator_test.cc +++ b/internal/ceres/schur_eliminator_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2019 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -284,10 +284,10 @@ BlockSparseMatrix matrix(bs); double* values = matrix.mutable_values(); std::mt19937 prng; - std::normal_distribution<double> standard_normal; - for (int i = 0; i < matrix.num_nonzeros(); ++i) { - values[i] = standard_normal(prng); - } + std::normal_distribution<> standard_normal; + std::generate_n(values, matrix.num_nonzeros(), [&prng, &standard_normal] { + return standard_normal(prng); + }); Vector b(matrix.num_rows()); b.setRandom();
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc index 73e1455..2397fdd 100644 --- a/internal/ceres/sparse_cholesky_test.cc +++ b/internal/ceres/sparse_cholesky_test.cc
@@ -55,7 +55,8 @@ const int num_col_blocks, const int min_col_block_size, const int max_col_block_size, - const double block_density) { + const double block_density, + std::mt19937& prng) { // Create a random matrix BlockSparseMatrix::RandomMatrixOptions options; options.num_col_blocks = num_col_blocks; @@ -66,7 +67,7 @@ options.min_row_block_size = 1; options.max_row_block_size = max_col_block_size; options.block_density = block_density; - auto random_matrix = BlockSparseMatrix::CreateRandomMatrix(options); + auto random_matrix = BlockSparseMatrix::CreateRandomMatrix(options, prng); // Add a diagonal block sparse matrix to make it full rank. Vector diagonal = Vector::Ones(random_matrix->num_cols()); @@ -110,7 +111,8 @@ const int num_blocks, const int min_block_size, const int max_block_size, - const double block_density) { + const double block_density, + std::mt19937& prng) { LinearSolver::Options sparse_cholesky_options; sparse_cholesky_options.sparse_linear_algebra_library_type = sparse_linear_algebra_library_type; @@ -120,7 +122,7 @@ sparse_cholesky->StorageType(); auto m = CreateRandomFullRankMatrix( - num_blocks, min_block_size, max_block_size, block_density); + num_blocks, min_block_size, max_block_size, block_density, prng); auto inner_product_computer = InnerProductComputer::Create(*m, storage_type); inner_product_computer->Compute(); CompressedRowSparseMatrix* lhs = inner_product_computer->mutable_result(); @@ -184,7 +186,8 @@ num_blocks, kMinBlockSize, kMaxBlockSize, - block_density); + block_density, + prng); } } }
diff --git a/internal/ceres/subset_preconditioner_test.cc b/internal/ceres/subset_preconditioner_test.cc index cd86695..5fb7f26 100644 --- a/internal/ceres/subset_preconditioner_test.cc +++ b/internal/ceres/subset_preconditioner_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2017 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -31,6 +31,7 @@ #include "ceres/subset_preconditioner.h" #include <memory> +#include <random> #include "Eigen/Dense" #include "Eigen/SparseCore" @@ -100,13 +101,13 @@ options.max_row_block_size = 4; options.block_density = 0.9; - m_ = BlockSparseMatrix::CreateRandomMatrix(options); + m_ = BlockSparseMatrix::CreateRandomMatrix(options, prng); start_row_block_ = m_->block_structure()->rows.size(); // Ensure that the bottom part of the matrix has the same column // block structure. options.col_blocks = m_->block_structure()->cols; - b_ = BlockSparseMatrix::CreateRandomMatrix(options); + b_ = BlockSparseMatrix::CreateRandomMatrix(options, prng); m_->AppendRows(*b_); // Create a Identity block diagonal matrix with the same column @@ -132,6 +133,7 @@ std::unique_ptr<Preconditioner> preconditioner_; Vector diagonal_; int start_row_block_; + std::mt19937 prng; }; TEST_P(SubsetPreconditionerTest, foo) {
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc index 49d367a..37a779f 100644 --- a/internal/ceres/triplet_sparse_matrix.cc +++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -310,7 +310,8 @@ } std::unique_ptr<TripletSparseMatrix> TripletSparseMatrix::CreateRandomMatrix( - const TripletSparseMatrix::RandomMatrixOptions& options) { + const TripletSparseMatrix::RandomMatrixOptions& options, + std::mt19937& prng) { CHECK_GT(options.num_rows, 0); CHECK_GT(options.num_cols, 0); CHECK_GT(options.density, 0.0); @@ -319,7 +320,6 @@ std::vector<int> rows; std::vector<int> cols; std::vector<double> values; - std::mt19937 prng; std::uniform_real_distribution<double> uniform01(0.0, 1.0); std::normal_distribution<double> standard_normal; while (rows.empty()) {
diff --git a/internal/ceres/triplet_sparse_matrix.h b/internal/ceres/triplet_sparse_matrix.h index c962455..549966a 100644 --- a/internal/ceres/triplet_sparse_matrix.h +++ b/internal/ceres/triplet_sparse_matrix.h
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2022 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -32,6 +32,7 @@ #define CERES_INTERNAL_TRIPLET_SPARSE_MATRIX_H_ #include <memory> +#include <random> #include <vector> #include "ceres/crs_matrix.h" @@ -135,7 +136,8 @@ // normally distributed and whose structure is determined by // RandomMatrixOptions. static std::unique_ptr<TripletSparseMatrix> CreateRandomMatrix( - const TripletSparseMatrix::RandomMatrixOptions& options); + const TripletSparseMatrix::RandomMatrixOptions& options, + std::mt19937& prng); // Load a triplet sparse matrix from a text file. static std::unique_ptr<TripletSparseMatrix> CreateFromTextFile(FILE* file);