1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
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
|
