Remove ceres/internal/random.h in favor of <random> Fixes https://github.com/ceres-solver/ceres-solver/issues/854 Change-Id: Id30b8dc2221f9afe4eb83f3a9304b9b2bc7e05d4
diff --git a/examples/BUILD b/examples/BUILD index 7bacf22..00b799e 100644 --- a/examples/BUILD +++ b/examples/BUILD
@@ -47,7 +47,6 @@ "bal_problem.cc", "bal_problem.h", "bundle_adjuster.cc", - "random.h", "snavely_reprojection_error.h", ], copts = EXAMPLE_COPTS, @@ -69,7 +68,6 @@ cc_binary( name = "robot_pose_mle", srcs = [ - "random.h", "robot_pose_mle.cc", ], copts = EXAMPLE_COPTS,
diff --git a/examples/bal_problem.cc b/examples/bal_problem.cc index 99bc075..2423fb6 100644 --- a/examples/bal_problem.cc +++ b/examples/bal_problem.cc
@@ -32,15 +32,15 @@ #include <algorithm> #include <cstdio> -#include <cstdlib> #include <fstream> +#include <functional> +#include <random> #include <string> #include <vector> #include "Eigen/Core" #include "ceres/rotation.h" #include "glog/logging.h" -#include "random.h" namespace ceres::examples { namespace { @@ -55,9 +55,9 @@ } } -void PerturbPoint3(const double sigma, double* point) { +void PerturbPoint3(std::function<double()> dist, double* point) { for (int i = 0; i < 3; ++i) { - point[i] += RandNormal() * sigma; + point[i] += dist(); } } @@ -297,14 +297,19 @@ CHECK_GE(point_sigma, 0.0); CHECK_GE(rotation_sigma, 0.0); CHECK_GE(translation_sigma, 0.0); - + std::mt19937 prng; + std::normal_distribution point_noise_distribution(0.0, point_sigma); double* points = mutable_points(); if (point_sigma > 0) { for (int i = 0; i < num_points_; ++i) { - PerturbPoint3(point_sigma, points + 3 * i); + PerturbPoint3(std::bind(point_noise_distribution, std::ref(prng)), + points + 3 * i); } } + std::normal_distribution rotation_noise_distribution(0.0, point_sigma); + std::normal_distribution translation_noise_distribution(0.0, + translation_sigma); for (int i = 0; i < num_cameras_; ++i) { double* camera = mutable_cameras() + camera_block_size() * i; @@ -314,12 +319,14 @@ // representation. CameraToAngleAxisAndCenter(camera, angle_axis, center); if (rotation_sigma > 0.0) { - PerturbPoint3(rotation_sigma, angle_axis); + PerturbPoint3(std::bind(rotation_noise_distribution, std::ref(prng)), + angle_axis); } AngleAxisAndCenterToCamera(angle_axis, center, camera); if (translation_sigma > 0.0) { - PerturbPoint3(translation_sigma, camera + camera_block_size() - 6); + PerturbPoint3(std::bind(translation_noise_distribution, std::ref(prng)), + camera + camera_block_size() - 6); } } }
diff --git a/examples/robot_pose_mle.cc b/examples/robot_pose_mle.cc index c199b89..f924bc3 100644 --- a/examples/robot_pose_mle.cc +++ b/examples/robot_pose_mle.cc
@@ -127,13 +127,13 @@ #include <cmath> #include <cstdio> +#include <random> #include <vector> #include "ceres/ceres.h" #include "ceres/dynamic_autodiff_cost_function.h" #include "gflags/gflags.h" #include "glog/logging.h" -#include "random.h" using ceres::AutoDiffCostFunction; using ceres::CauchyLoss; @@ -143,7 +143,6 @@ using ceres::Problem; using ceres::Solve; using ceres::Solver; -using ceres::examples::RandNormal; using std::min; using std::vector; @@ -247,6 +246,10 @@ const int num_steps = static_cast<int>(ceil(CERES_GET_FLAG(FLAGS_corridor_length) / CERES_GET_FLAG(FLAGS_pose_separation))); + std::mt19937 prng; + std::normal_distribution odometry_noise( + 0.0, CERES_GET_FLAG(FLAGS_odometry_stddev)); + std::normal_distribution range_noise(0.0, CERES_GET_FLAG(FLAGS_range_stddev)); // The robot starts out at the origin. double robot_location = 0.0; @@ -258,10 +261,8 @@ const double actual_range = CERES_GET_FLAG(FLAGS_corridor_length) - robot_location; const double observed_odometry = - RandNormal() * CERES_GET_FLAG(FLAGS_odometry_stddev) + - actual_odometry_value; - const double observed_range = - RandNormal() * CERES_GET_FLAG(FLAGS_range_stddev) + actual_range; + actual_odometry_value + odometry_noise(prng); + const double observed_range = actual_range + range_noise(prng); odometry_values->push_back(observed_odometry); range_readings->push_back(observed_range); }
diff --git a/internal/ceres/autodiff_test.cc b/internal/ceres/autodiff_test.cc index 5f769a0..629eea8 100644 --- a/internal/ceres/autodiff_test.cc +++ b/internal/ceres/autodiff_test.cc
@@ -30,7 +30,8 @@ #include "ceres/internal/autodiff.h" -#include "ceres/random.h" +#include <random> + #include "gtest/gtest.h" namespace ceres::internal { @@ -162,17 +163,18 @@ // Test projective camera model projector. TEST(AutoDiff, ProjectiveCameraModel) { - srand(5); double const tol = 1e-10; // floating-point tolerance. double const del = 1e-4; // finite-difference step. double const err = 1e-6; // finite-difference tolerance. Projective b; + std::mt19937 prng; + std::uniform_real_distribution<double> uniform01(0.0, 1.0); // Make random P and X, in a single vector. double PX[12 + 4]; for (double& PX_i : PX) { - PX_i = RandDouble(); + PX_i = uniform01(prng); } // Handy names for the P and X parts. @@ -282,16 +284,17 @@ // This test is similar in structure to the previous one. TEST(AutoDiff, Metric) { - srand(5); double const tol = 1e-10; // floating-point tolerance. double const del = 1e-4; // finite-difference step. - double const err = 1e-5; // finite-difference tolerance. + double const err = 2e-5; // finite-difference tolerance. Metric b; // Make random parameter vector. double qcX[4 + 3 + 3]; - for (double& qcX_i : qcX) qcX_i = RandDouble(); + std::mt19937 prng; + std::uniform_real_distribution<double> uniform01(0.0, 1.0); + for (double& qcX_i : qcX) qcX_i = uniform01(prng); // Handy names. double* q = qcX;
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index 22224a3..9ae4384 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc
@@ -33,12 +33,12 @@ #include <algorithm> #include <cstddef> #include <memory> +#include <random> #include <vector> #include "ceres/block_structure.h" #include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" -#include "ceres/random.h" #include "ceres/small_blas.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" @@ -382,6 +382,11 @@ CHECK_GT(options.block_density, 0.0); CHECK_LE(options.block_density, 1.0); + std::mt19937 prng; + std::uniform_int_distribution col_distribution(options.min_col_block_size, + options.max_col_block_size); + std::uniform_int_distribution row_distribution(options.min_row_block_size, + options.max_row_block_size); auto* bs = new CompressedRowBlockStructure(); if (options.col_blocks.empty()) { CHECK_GT(options.num_col_blocks, 0); @@ -392,10 +397,7 @@ // Generate the col block structure. int col_block_position = 0; for (int i = 0; i < options.num_col_blocks; ++i) { - // Generate a random integer in [min_col_block_size, max_col_block_size] - const int delta_block_size = - Uniform(options.max_col_block_size - options.min_col_block_size); - const int col_block_size = options.min_col_block_size + delta_block_size; + const int col_block_size = col_distribution(prng); bs->cols.emplace_back(col_block_size, col_block_position); col_block_position += col_block_size; } @@ -404,22 +406,21 @@ } bool matrix_has_blocks = false; + std::uniform_real_distribution uniform01(0.0, 1.0); while (!matrix_has_blocks) { VLOG(1) << "Clearing"; bs->rows.clear(); int row_block_position = 0; int value_position = 0; for (int r = 0; r < options.num_row_blocks; ++r) { - const int delta_block_size = - Uniform(options.max_row_block_size - options.min_row_block_size); - const int row_block_size = options.min_row_block_size + delta_block_size; + const int row_block_size = row_distribution(prng); bs->rows.emplace_back(); CompressedRow& row = bs->rows.back(); row.block.size = row_block_size; row.block.position = row_block_position; row_block_position += row_block_size; for (int c = 0; c < bs->cols.size(); ++c) { - if (RandDouble() > options.block_density) continue; + if (uniform01(prng) > options.block_density) continue; row.cells.emplace_back(); Cell& cell = row.cells.back(); @@ -433,8 +434,9 @@ auto matrix = std::make_unique<BlockSparseMatrix>(bs); double* values = matrix->mutable_values(); + std::normal_distribution standard_normal_distribution; for (int i = 0; i < matrix->num_nonzeros(); ++i) { - values[i] = RandNormal(); + values[i] = standard_normal_distribution(prng); } return matrix;
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc index d663f4a..bd94fdc 100644 --- a/internal/ceres/compressed_row_sparse_matrix.cc +++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -31,13 +31,14 @@ #include "ceres/compressed_row_sparse_matrix.h" #include <algorithm> +#include <functional> #include <memory> #include <numeric> +#include <random> #include <vector> #include "ceres/crs_matrix.h" #include "ceres/internal/export.h" -#include "ceres/random.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" @@ -122,6 +123,7 @@ const int num_cols, const int row_block_begin, const int col_block_begin, + std::function<double()> dist, std::vector<int>* rows, std::vector<int>* cols, std::vector<double>* values) { @@ -129,19 +131,20 @@ 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(RandNormal()); + values->push_back(dist()); } } } void AddSymmetricRandomBlock(const int num_rows, const int row_block_begin, + std::function<double()> dist, 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 = RandNormal(); + const double v = dist(); rows->push_back(row_block_begin + r); cols->push_back(row_block_begin + c); values->push_back(v); @@ -644,21 +647,26 @@ vector<int> row_blocks; vector<int> col_blocks; + std::mt19937 prng; + std::uniform_int_distribution col_distribution(options.min_col_block_size, + options.max_col_block_size); + std::uniform_int_distribution row_distribution(options.min_row_block_size, + options.max_row_block_size); + std::uniform_real_distribution uniform01(0.0, 1.0); + std::normal_distribution standard_normal_distribution; + auto values_dist = std::bind(standard_normal_distribution, std::ref(prng)); + // Generate the row block structure. for (int i = 0; i < options.num_row_blocks; ++i) { // Generate a random integer in [min_row_block_size, max_row_block_size] - const int delta_block_size = - Uniform(options.max_row_block_size - options.min_row_block_size); - row_blocks.push_back(options.min_row_block_size + delta_block_size); + row_blocks.push_back(row_distribution(prng)); } if (options.storage_type == StorageType::UNSYMMETRIC) { // Generate the col block structure. for (int i = 0; i < options.num_col_blocks; ++i) { // Generate a random integer in [min_col_block_size, max_col_block_size] - const int delta_block_size = - Uniform(options.max_col_block_size - options.min_col_block_size); - col_blocks.push_back(options.min_col_block_size + delta_block_size); + col_blocks.push_back(col_distribution(prng)); } } else { // Symmetric matrices (LOWER_TRIANGULAR or UPPER_TRIANGULAR); @@ -695,7 +703,7 @@ } // Randomly determine if this block is present or not. - if (RandDouble() <= options.block_density) { + if (uniform01(prng) <= options.block_density) { // If the matrix is symmetric, then we take care to generate // symmetric diagonal blocks. if (options.storage_type == StorageType::UNSYMMETRIC || r != c) { @@ -703,12 +711,14 @@ col_blocks[c], row_block_begin, col_block_begin, + values_dist, &tsm_rows, &tsm_cols, &tsm_values); } else { AddSymmetricRandomBlock(row_blocks[r], row_block_begin, + values_dist, &tsm_rows, &tsm_cols, &tsm_values);
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc index 4151845..db90d8b 100644 --- a/internal/ceres/compressed_row_sparse_matrix_test.cc +++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -33,6 +33,7 @@ #include <algorithm> #include <memory> #include <numeric> +#include <random> #include <string> #include "Eigen/SparseCore" @@ -40,7 +41,6 @@ #include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/linear_least_squares_problems.h" -#include "ceres/random.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" #include "gtest/gtest.h" @@ -404,7 +404,8 @@ const int kMinBlockSize = 1; const int kMaxBlockSize = 5; const int kNumTrials = 10; - + std::mt19937 prng; + std::uniform_real_distribution uniform(0.5, 1.0); for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks; ++num_blocks) { for (int trial = 0; trial < kNumTrials; ++trial) { @@ -416,7 +417,7 @@ options.num_row_blocks = 2 * num_blocks; options.min_row_block_size = kMinBlockSize; options.max_row_block_size = kMaxBlockSize; - options.block_density = std::max(0.5, RandDouble()); + options.block_density = uniform(prng); options.storage_type = ::testing::get<0>(param); std::unique_ptr<CompressedRowSparseMatrix> matrix = CompressedRowSparseMatrix::CreateRandomMatrix(options); @@ -473,7 +474,8 @@ const int kMinBlockSize = 1; const int kMaxBlockSize = 5; const int kNumTrials = 10; - + std::mt19937 prng; + std::uniform_real_distribution uniform(0.5, 1.0); for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks; ++num_blocks) { for (int trial = 0; trial < kNumTrials; ++trial) { @@ -485,7 +487,7 @@ options.num_row_blocks = 2 * num_blocks; options.min_row_block_size = kMinBlockSize; options.max_row_block_size = kMaxBlockSize; - options.block_density = std::max(0.5, RandDouble()); + options.block_density = uniform(prng); options.storage_type = ::testing::get<0>(param); std::unique_ptr<CompressedRowSparseMatrix> matrix = CompressedRowSparseMatrix::CreateRandomMatrix(options); @@ -542,7 +544,8 @@ const int kMinBlockSize = 1; const int kMaxBlockSize = 5; const int kNumTrials = 10; - + std::mt19937 prng; + std::uniform_real_distribution uniform(0.5, 1.0); for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks; ++num_blocks) { for (int trial = 0; trial < kNumTrials; ++trial) { @@ -554,7 +557,7 @@ options.num_row_blocks = 2 * num_blocks; options.min_row_block_size = kMinBlockSize; options.max_row_block_size = kMaxBlockSize; - options.block_density = std::max(0.5, RandDouble()); + options.block_density = uniform(prng); options.storage_type = ::testing::get<0>(param); std::unique_ptr<CompressedRowSparseMatrix> matrix = CompressedRowSparseMatrix::CreateRandomMatrix(options);
diff --git a/internal/ceres/corrector_test.cc b/internal/ceres/corrector_test.cc index d2f7c84..b1e9bee 100644 --- a/internal/ceres/corrector_test.cc +++ b/internal/ceres/corrector_test.cc
@@ -34,9 +34,9 @@ #include <cmath> #include <cstdlib> #include <cstring> +#include <random> #include "ceres/internal/eigen.h" -#include "ceres/random.h" #include "gtest/gtest.h" namespace ceres { @@ -160,18 +160,19 @@ // and hessians. Matrix c_hess(2, 2); Vector c_grad(2); - - srand(5); + std::mt19937 prng; + std::uniform_real_distribution uniform01(0.0, 1.0); for (int iter = 0; iter < 10000; ++iter) { // Initialize the jacobian and residual. - for (double& jacobian_entry : jacobian) jacobian_entry = RandDouble(); - for (double& residual : residuals) residual = RandDouble(); + for (double& jacobian_entry : jacobian) jacobian_entry = uniform01(prng); + for (double& residual : residuals) residual = uniform01(prng); const double sq_norm = res.dot(res); rho[0] = sq_norm; - rho[1] = RandDouble(); - rho[2] = 2.0 * RandDouble() - 1.0; + rho[1] = uniform01(prng); + rho[2] = uniform01( + prng, std::uniform_real_distribution<double>::param_type(-1, 1)); // If rho[2] > 0, then the curvature correction to the correction // and the gauss newton approximation will match. Otherwise, we @@ -227,10 +228,11 @@ Matrix c_hess(2, 2); Vector c_grad(2); - srand(5); + std::mt19937 prng; + std::uniform_real_distribution uniform01(0.0, 1.0); for (int iter = 0; iter < 10000; ++iter) { // Initialize the jacobian. - for (double& jacobian_entry : jacobian) jacobian_entry = RandDouble(); + for (double& jacobian_entry : jacobian) jacobian_entry = uniform01(prng); // Zero residuals res.setZero(); @@ -238,8 +240,9 @@ const double sq_norm = res.dot(res); rho[0] = sq_norm; - rho[1] = RandDouble(); - rho[2] = 2 * RandDouble() - 1.0; + rho[1] = uniform01(prng); + rho[2] = uniform01( + prng, std::uniform_real_distribution<double>::param_type(-1, 1)); // Ground truth values. g_res = sqrt(rho[1]) * res;
diff --git a/internal/ceres/dense_linear_solver_benchmark.cc b/internal/ceres/dense_linear_solver_benchmark.cc index 53628ad..a01e0d6 100644 --- a/internal/ceres/dense_linear_solver_benchmark.cc +++ b/internal/ceres/dense_linear_solver_benchmark.cc
@@ -34,7 +34,6 @@ #include "ceres/dense_sparse_matrix.h" #include "ceres/internal/config.h" #include "ceres/linear_solver.h" -#include "ceres/random.h" namespace ceres::internal {
diff --git a/internal/ceres/gradient_checker_test.cc b/internal/ceres/gradient_checker_test.cc index aedb0b9..5de4778 100644 --- a/internal/ceres/gradient_checker_test.cc +++ b/internal/ceres/gradient_checker_test.cc
@@ -34,12 +34,12 @@ #include <cmath> #include <cstdlib> +#include <random> #include <utility> #include <vector> #include "ceres/cost_function.h" #include "ceres/problem.h" -#include "ceres/random.h" #include "ceres/solver.h" #include "ceres/test_util.h" #include "glog/logging.h" @@ -60,12 +60,15 @@ class GoodTestTerm : public CostFunction { public: GoodTestTerm(int arity, int const* dim) : arity_(arity), return_value_(true) { + std::mt19937 prng; + 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] = 2.0 * RandDouble() - 1.0; + a_[j][u] = distribution(prng); } } @@ -119,12 +122,14 @@ class BadTestTerm : public CostFunction { public: BadTestTerm(int arity, int const* dim) : arity_(arity) { + std::mt19937 prng; + 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] = 2.0 * RandDouble() - 1.0; + a_[j][u] = distribution(prng); } } @@ -194,8 +199,6 @@ } TEST(GradientChecker, SmokeTest) { - srand(5); - // Test with 3 blocks of size 2, 3 and 4. int const num_parameters = 3; std::vector<int> parameter_sizes(3); @@ -205,10 +208,12 @@ // Make a random set of blocks. FixedArray<double*> parameters(num_parameters); + std::mt19937 prng; + std::uniform_real_distribution distribution(-1.0, 1.0); 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] = 2.0 * RandDouble() - 1.0; + parameters[j][u] = distribution(prng); } }
diff --git a/internal/ceres/gradient_checking_cost_function_test.cc b/internal/ceres/gradient_checking_cost_function_test.cc index 5c6046f..99e0727 100644 --- a/internal/ceres/gradient_checking_cost_function_test.cc +++ b/internal/ceres/gradient_checking_cost_function_test.cc
@@ -33,6 +33,7 @@ #include <cmath> #include <cstdint> #include <memory> +#include <random> #include <vector> #include "ceres/cost_function.h" @@ -41,7 +42,6 @@ #include "ceres/parameter_block.h" #include "ceres/problem_impl.h" #include "ceres/program.h" -#include "ceres/random.h" #include "ceres/residual_block.h" #include "ceres/sized_cost_function.h" #include "ceres/types.h" @@ -70,12 +70,14 @@ // 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 uniform01(-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] = 2.0 * RandDouble() - 1.0; + a_[j][u] = uniform01(prng); } } @@ -130,18 +132,18 @@ }; TEST(GradientCheckingCostFunction, ResidualsAndJacobiansArePreservedTest) { - srand(5); - // Test with 3 blocks of size 2, 3 and 4. int const arity = 3; int const dim[arity] = {2, 3, 4}; // Make a random set of blocks. vector<double*> parameters(arity); + std::mt19937 prng; + std::uniform_real_distribution uniform01(-1.0, 1.0); for (int j = 0; j < arity; ++j) { parameters[j] = new double[dim[j]]; for (int u = 0; u < dim[j]; ++u) { - parameters[j][u] = 2.0 * RandDouble() - 1.0; + parameters[j][u] = uniform01(prng); } } @@ -187,18 +189,18 @@ } TEST(GradientCheckingCostFunction, SmokeTest) { - srand(5); - // Test with 3 blocks of size 2, 3 and 4. int const arity = 3; int const dim[arity] = {2, 3, 4}; // Make a random set of blocks. vector<double*> parameters(arity); + std::mt19937 prng; + std::uniform_real_distribution uniform01(-1.0, 1.0); for (int j = 0; j < arity; ++j) { parameters[j] = new double[dim[j]]; for (int u = 0; u < dim[j]; ++u) { - parameters[j][u] = 2.0 * RandDouble() - 1.0; + parameters[j][u] = uniform01(prng); } }
diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc index a1f5c49..9c0382b 100644 --- a/internal/ceres/inner_product_computer_test.cc +++ b/internal/ceres/inner_product_computer_test.cc
@@ -32,11 +32,11 @@ #include <memory> #include <numeric> +#include <random> #include "Eigen/SparseCore" #include "ceres/block_sparse_matrix.h" #include "ceres/internal/eigen.h" -#include "ceres/random.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" #include "gtest/gtest.h" @@ -77,11 +77,11 @@ } TEST(InnerProductComputer, NormalOperation) { - // "Randomly generated seed." - SetRandomState(29823); const int kMaxNumRowBlocks = 10; const int kMaxNumColBlocks = 10; const int kNumTrials = 10; + std::mt19937 prng; + std::uniform_real_distribution distribution(0.01, 1.0); // Create a random matrix, compute its outer product using Eigen and // ComputeOuterProduct. Convert both matrices to dense matrices and @@ -98,7 +98,8 @@ options.max_row_block_size = 5; options.min_col_block_size = 1; options.max_col_block_size = 10; - options.block_density = std::max(0.1, RandDouble()); + options.block_density = distribution(prng); + VLOG(2) << "num row blocks: " << options.num_row_blocks; VLOG(2) << "num col blocks: " << options.num_col_blocks; VLOG(2) << "min row block size: " << options.min_row_block_size; @@ -143,6 +144,8 @@ const int kNumRowBlocks = 10; const int kNumColBlocks = 20; const int kNumTrials = 5; + std::mt19937 prng; + std::uniform_real_distribution distribution(0.01, 1.0); // Create a random matrix, compute its outer product using Eigen and // ComputeInnerProductComputer. Convert both matrices to dense matrices and @@ -155,7 +158,7 @@ options.max_row_block_size = 5; options.min_col_block_size = 1; options.max_col_block_size = 10; - options.block_density = std::min(0.01, RandDouble()); + options.block_density = distribution(prng); VLOG(2) << "num row blocks: " << options.num_row_blocks; VLOG(2) << "num col blocks: " << options.num_col_blocks;
diff --git a/internal/ceres/normal_prior_test.cc b/internal/ceres/normal_prior_test.cc index c66f9d6..9c84555 100644 --- a/internal/ceres/normal_prior_test.cc +++ b/internal/ceres/normal_prior_test.cc
@@ -31,43 +31,26 @@ #include "ceres/normal_prior.h" #include <cstddef> +#include <random> #include "ceres/internal/eigen.h" -#include "ceres/random.h" #include "gtest/gtest.h" namespace ceres { namespace internal { -namespace { - -void RandomVector(Vector* v) { - for (int r = 0; r < v->rows(); ++r) (*v)[r] = 2 * RandDouble() - 1; -} - -void RandomMatrix(Matrix* m) { - for (int r = 0; r < m->rows(); ++r) { - for (int c = 0; c < m->cols(); ++c) { - (*m)(r, c) = 2 * RandDouble() - 1; - } - } -} - -} // namespace - TEST(NormalPriorTest, ResidualAtRandomPosition) { - srand(5); - + std::mt19937 prng; + std::uniform_real_distribution distribution(-1.0, 1.0); for (int num_rows = 1; num_rows < 5; ++num_rows) { for (int num_cols = 1; num_cols < 5; ++num_cols) { Vector b(num_cols); - RandomVector(&b); - + b.setRandom(); Matrix A(num_rows, num_cols); - RandomMatrix(&A); + A.setRandom(); auto* x = new double[num_cols]; - for (int i = 0; i < num_cols; ++i) x[i] = 2 * RandDouble() - 1; + for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng); auto* jacobian = new double[num_rows * num_cols]; Vector residuals(num_rows); @@ -92,18 +75,17 @@ } TEST(NormalPriorTest, ResidualAtRandomPositionNullJacobians) { - srand(5); - + std::mt19937 prng; + std::uniform_real_distribution distribution(-1.0, 1.0); for (int num_rows = 1; num_rows < 5; ++num_rows) { for (int num_cols = 1; num_cols < 5; ++num_cols) { Vector b(num_cols); - RandomVector(&b); - + b.setRandom(); Matrix A(num_rows, num_cols); - RandomMatrix(&A); + A.setRandom(); auto* x = new double[num_cols]; - for (int i = 0; i < num_cols; ++i) x[i] = 2 * RandDouble() - 1; + for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng); double* jacobians[1]; jacobians[0] = nullptr;
diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc index d01dd3b..6e9dd92 100644 --- a/internal/ceres/partitioned_matrix_view_test.cc +++ b/internal/ceres/partitioned_matrix_view_test.cc
@@ -31,13 +31,13 @@ #include "ceres/partitioned_matrix_view.h" #include <memory> +#include <random> #include <vector> #include "ceres/block_structure.h" #include "ceres/casts.h" #include "ceres/internal/eigen.h" #include "ceres/linear_least_squares_problems.h" -#include "ceres/random.h" #include "ceres/sparse_matrix.h" #include "glog/logging.h" #include "gtest/gtest.h" @@ -50,7 +50,6 @@ class PartitionedMatrixViewTest : public ::testing::Test { protected: void SetUp() final { - srand(5); std::unique_ptr<LinearLeastSquaresProblem> problem = CreateLinearLeastSquaresProblemFromId(2); CHECK(problem != nullptr); @@ -65,11 +64,16 @@ options, *down_cast<BlockSparseMatrix*>(A_.get())); } + double RandDouble() { return distribution_(prng_); } + int num_rows_; int num_cols_; int num_eliminate_blocks_; std::unique_ptr<SparseMatrix> A_; std::unique_ptr<PartitionedMatrixViewBase> pmv_; + std::mt19937 prng_; + std::uniform_real_distribution<double> distribution_ = + std::uniform_real_distribution(0.0, 1.0); }; TEST_F(PartitionedMatrixViewTest, DimensionsTest) {
diff --git a/internal/ceres/random.h b/internal/ceres/random.h deleted file mode 100644 index 0495d67..0000000 --- a/internal/ceres/random.h +++ /dev/null
@@ -1,73 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. -// http://ceres-solver.org/ -// -// Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions are met: -// -// * Redistributions of source code must retain the above copyright notice, -// this list of conditions and the following disclaimer. -// * Redistributions in binary form must reproduce the above copyright notice, -// this list of conditions and the following disclaimer in the documentation -// and/or other materials provided with the distribution. -// * Neither the name of Google Inc. nor the names of its contributors may be -// used to endorse or promote products derived from this software without -// specific prior written permission. -// -// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE -// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. -// -// Author: keir@google.com (Keir Mierle) -// sameeragarwal@google.com (Sameer Agarwal) - -#ifndef CERES_INTERNAL_RANDOM_H_ -#define CERES_INTERNAL_RANDOM_H_ - -#include <cmath> -#include <cstdlib> - -#include "ceres/internal/export.h" - -namespace ceres { - -inline void SetRandomState(int state) { srand(state); } - -inline int Uniform(int n) { - if (n) { - return rand() % n; - } else { - return 0; - } -} - -inline double RandDouble() { - auto r = static_cast<double>(rand()); - return r / RAND_MAX; -} - -// Box-Muller algorithm for normal random number generation. -// http://en.wikipedia.org/wiki/Box-Muller_transform -inline double RandNormal() { - double x1, x2, w; - do { - x1 = 2.0 * RandDouble() - 1.0; - x2 = 2.0 * RandDouble() - 1.0; - w = x1 * x1 + x2 * x2; - } while (w >= 1.0 || w == 0.0); - - w = sqrt((-2.0 * log(w)) / w); - return x1 * w; -} - -} // namespace ceres - -#endif // CERES_INTERNAL_RANDOM_H_
diff --git a/internal/ceres/schur_eliminator_benchmark.cc b/internal/ceres/schur_eliminator_benchmark.cc index 2bb18ca..0973e73 100644 --- a/internal/ceres/schur_eliminator_benchmark.cc +++ b/internal/ceres/schur_eliminator_benchmark.cc
@@ -29,13 +29,13 @@ // Authors: sameeragarwal@google.com (Sameer Agarwal) #include <memory> +#include <random> #include "Eigen/Dense" #include "benchmark/benchmark.h" #include "ceres/block_random_access_dense_matrix.h" #include "ceres/block_sparse_matrix.h" #include "ceres/block_structure.h" -#include "ceres/random.h" #include "ceres/schur_eliminator.h" namespace ceres::internal { @@ -91,8 +91,10 @@ 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] = RandNormal(); + values[i] = standard_normal(prng); } b_.resize(matrix_->num_rows());
diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc index a8c5f71..da27019 100644 --- a/internal/ceres/schur_eliminator_test.cc +++ b/internal/ceres/schur_eliminator_test.cc
@@ -31,6 +31,7 @@ #include "ceres/schur_eliminator.h" #include <memory> +#include <random> #include "Eigen/Dense" #include "ceres/block_random_access_dense_matrix.h" @@ -41,7 +42,6 @@ #include "ceres/detect_structure.h" #include "ceres/internal/eigen.h" #include "ceres/linear_least_squares_problems.h" -#include "ceres/random.h" #include "ceres/test_util.h" #include "ceres/triplet_sparse_matrix.h" #include "ceres/types.h" @@ -283,8 +283,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] = RandNormal(); + values[i] = standard_normal(prng); } Vector b(matrix.num_rows());
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc index 0522232..e99d59b 100644 --- a/internal/ceres/sparse_cholesky_test.cc +++ b/internal/ceres/sparse_cholesky_test.cc
@@ -32,6 +32,7 @@ #include <memory> #include <numeric> +#include <random> #include <vector> #include "Eigen/Dense" @@ -42,7 +43,6 @@ #include "ceres/internal/config.h" #include "ceres/internal/eigen.h" #include "ceres/iterative_refiner.h" -#include "ceres/random.h" #include "glog/logging.h" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -165,17 +165,18 @@ class SparseCholeskyTest : public ::testing::TestWithParam<Param> {}; TEST_P(SparseCholeskyTest, FactorAndSolve) { - SetRandomState(2982); const int kMinNumBlocks = 1; const int kMaxNumBlocks = 10; const int kNumTrials = 10; const int kMinBlockSize = 1; const int kMaxBlockSize = 5; + std::mt19937 prng; + std::uniform_real_distribution distribution(0.1, 1.0); for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks; ++num_blocks) { for (int trial = 0; trial < kNumTrials; ++trial) { - const double block_density = std::max(0.1, RandDouble()); + const double block_density = distribution(prng); Param param = GetParam(); SparseCholeskySolverUnitTest(::testing::get<0>(param), ::testing::get<1>(param),
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc index 86a444e..c5ff2e5 100644 --- a/internal/ceres/triplet_sparse_matrix.cc +++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -32,12 +32,12 @@ #include <algorithm> #include <memory> +#include <random> #include "ceres/compressed_row_sparse_matrix.h" #include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/internal/export.h" -#include "ceres/random.h" #include "ceres/types.h" #include "glog/logging.h" @@ -282,8 +282,8 @@ } } -std::unique_ptr<TripletSparseMatrix> -TripletSparseMatrix::CreateFromTextFile(FILE* file) { +std::unique_ptr<TripletSparseMatrix> TripletSparseMatrix::CreateFromTextFile( + FILE* file) { CHECK(file != nullptr); int num_rows = 0; int num_cols = 0; @@ -317,16 +317,19 @@ std::vector<int> rows; std::vector<int> cols; std::vector<double> values; + std::mt19937 prng; + std::uniform_real_distribution uniform01(0.0, 1.0); + std::uniform_real_distribution standard_normal; while (rows.empty()) { rows.clear(); cols.clear(); values.clear(); for (int r = 0; r < options.num_rows; ++r) { for (int c = 0; c < options.num_cols; ++c) { - if (RandDouble() <= options.density) { + if (uniform01(prng) <= options.density) { rows.push_back(r); cols.push_back(c); - values.push_back(RandNormal()); + values.push_back(standard_normal(prng)); } } }