random_orthogonal_matrix.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): Siargey Kachanovich
4  *
5  * Copyright (C) 2019 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  */
10 
11 #ifndef FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_
12 #define FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_
13 
14 #include <cstdlib> // for std::size_t
15 #include <cmath> // for std::cos, std::sin
16 #include <random> // for std::uniform_real_distribution, std::random_device
17 
18 #include <Eigen/Dense>
19 #include <Eigen/Sparse>
20 #include <Eigen/SVD>
21 
22 #include <CGAL/Epick_d.h>
23 #include <CGAL/point_generators_d.h>
24 
25 #include <boost/math/constants/constants.hpp> // for PI value
26 
27 namespace Gudhi {
28 
29 namespace coxeter_triangulation {
30 
39 // Note: the householderQR operation at the end seems to take a lot of time at compilation.
40 // The CGAL headers are another source of long compilation time.
41 Eigen::MatrixXd random_orthogonal_matrix(std::size_t d) {
42  typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel;
43  typedef typename Kernel::Point_d Point_d;
44  if (d == 1) return Eigen::VectorXd::Constant(1, 1.0);
45  if (d == 2) {
46  // 0. < alpha < 2 Pi
47  std::uniform_real_distribution<double> unif(0., 2 * boost::math::constants::pi<double>());
48  std::random_device rand_dev;
49  std::mt19937 rand_engine(rand_dev());
50  double alpha = unif(rand_engine);
51 
52  Eigen::Matrix2d rot;
53  rot << std::cos(alpha), -std::sin(alpha), std::sin(alpha), cos(alpha);
54  return rot;
55  }
56  Eigen::MatrixXd low_dim_rot = random_orthogonal_matrix(d - 1);
57  Eigen::MatrixXd rot(d, d);
58  Point_d v = *CGAL::Random_points_on_sphere_d<Point_d>(d, 1);
59  for (std::size_t i = 0; i < d; ++i) rot(i, 0) = v[i];
60  for (std::size_t i = 0; i < d - 1; ++i)
61  for (std::size_t j = 1; j < d - 1; ++j) rot(i, j) = low_dim_rot(i, j - 1);
62  for (std::size_t j = 1; j < d; ++j) rot(d - 1, j) = 0;
63  rot = rot.householderQr()
64  .householderQ(); // a way to do Gram-Schmidt, see https://forum.kde.org/viewtopic.php?f=74&t=118568#p297246
65  return rot;
66 }
67 
68 } // namespace coxeter_triangulation
69 
70 } // namespace Gudhi
71 
72 #endif
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Fri Jan 14 2022 18:28:42 for GUDHIdev by Doxygen 1.9.1