diff options
author | Hind-M <hind.montassif@gmail.com> | 2022-01-18 16:30:33 +0100 |
---|---|---|
committer | Hind-M <hind.montassif@gmail.com> | 2022-01-18 16:30:33 +0100 |
commit | f7c3fd1033e8956777b82d84d231115cb5540bc3 (patch) | |
tree | 1bd6467863b75e7b11247a899c5172f99f14a3b4 /src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h | |
parent | aa600c433e1f756bec4323e29e86786b937d9443 (diff) | |
parent | de5aa9c891ef13c9fc2b2635bcd27ab873b0057b (diff) |
Merge remote-tracking branch 'upstream/master' into fetch_datasets
Diffstat (limited to 'src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h')
-rw-r--r-- | src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h | 72 |
1 files changed, 72 insertions, 0 deletions
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h new file mode 100644 index 00000000..6a896e94 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h @@ -0,0 +1,72 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2019 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_ +#define FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_ + +#include <cstdlib> // for std::size_t +#include <cmath> // for std::cos, std::sin +#include <random> // for std::uniform_real_distribution, std::random_device + +#include <Eigen/Dense> +#include <Eigen/Sparse> +#include <Eigen/SVD> + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> + +#include <boost/math/constants/constants.hpp> // for PI value + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \brief Generates a uniform random orthogonal matrix using the "subgroup algorithm" by + * Diaconis & Shashahani. + * \details Taken from https://en.wikipedia.org/wiki/Rotation_matrix#Uniform_random_rotation_matrices. + * The idea: take a random rotation matrix of dimension d-1, embed it + * as a d*d matrix M with the last column (0,...,0,1). + * Pick a random vector v on a sphere S^d. rotate the matrix M so that its last column is v. + * The determinant of the matrix can be either 1 or -1 + */ +// Note: the householderQR operation at the end seems to take a lot of time at compilation. +// The CGAL headers are another source of long compilation time. +Eigen::MatrixXd random_orthogonal_matrix(std::size_t d) { + typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel; + typedef typename Kernel::Point_d Point_d; + if (d == 1) return Eigen::VectorXd::Constant(1, 1.0); + if (d == 2) { + // 0. < alpha < 2 Pi + std::uniform_real_distribution<double> unif(0., 2 * boost::math::constants::pi<double>()); + std::random_device rand_dev; + std::mt19937 rand_engine(rand_dev()); + double alpha = unif(rand_engine); + + Eigen::Matrix2d rot; + rot << std::cos(alpha), -std::sin(alpha), std::sin(alpha), cos(alpha); + return rot; + } + Eigen::MatrixXd low_dim_rot = random_orthogonal_matrix(d - 1); + Eigen::MatrixXd rot(d, d); + Point_d v = *CGAL::Random_points_on_sphere_d<Point_d>(d, 1); + for (std::size_t i = 0; i < d; ++i) rot(i, 0) = v[i]; + for (std::size_t i = 0; i < d - 1; ++i) + for (std::size_t j = 1; j < d - 1; ++j) rot(i, j) = low_dim_rot(i, j - 1); + for (std::size_t j = 1; j < d; ++j) rot(d - 1, j) = 0; + rot = rot.householderQr() + .householderQ(); // a way to do Gram-Schmidt, see https://forum.kde.org/viewtopic.php?f=74&t=118568#p297246 + return rot; +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif |