From 5344df8b8534403766e8ea28d7b57a0c2ae20583 Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Fri, 19 Jul 2024 13:48:09 +0530 Subject: [PATCH 1/5] modify MT functions --- .../sampling/sample_correlation_matrices.hpp | 36 +++++++++++++++---- 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/include/sampling/sample_correlation_matrices.hpp b/include/sampling/sample_correlation_matrices.hpp index 9a1de2319..e2b0e7e6b 100644 --- a/include/sampling/sample_correlation_matrices.hpp +++ b/include/sampling/sample_correlation_matrices.hpp @@ -41,19 +41,27 @@ template typename WalkTypePolicy, typename PointType, typename RNGType, - typename PointList + typename MT > void uniform_correlation_sampling_MT( const unsigned int &n, - PointList &randPoints, + std::list &randCorMatrices, const unsigned int &walkL, const unsigned int &num_points, unsigned int const& nburns){ + using PointList = std::vector; + PointList randPoints; + CorrelationSpectrahedron_MT P(n); const unsigned int d = P.dimension(); PointType startingPoint(n); RNGType rng(d); uniform_sampling(randPoints, P, rng, walkL, num_points, startingPoint, nburns); + + for(const auto&p : randPoints){ + MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n); + randCorMatrices.push_back(final_cor_mat); + } } template @@ -83,21 +91,29 @@ template typename WalkTypePolicy, typename PointType, typename RNGType, - typename PointList, + typename MT, typename NT > void gaussian_correlation_sampling_MT( const unsigned int &n, - PointList &randPoints, + std::list &randCorMatrices, const unsigned int &walkL, const unsigned int &num_points, const NT &a, unsigned int const& nburns = 0){ + using PointList = std::vector; + PointList randPoints; + CorrelationSpectrahedron_MT P(n); const unsigned int d = P.dimension(); PointType startingPoint(n); RNGType rng(d); gaussian_sampling(randPoints, P, rng, walkL, num_points, a, startingPoint, nburns); + + for(const auto&p : randPoints){ + MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n); + randCorMatrices.push_back(final_cor_mat); + } } template @@ -130,17 +146,20 @@ template typename WalkTypePolicy, typename PointType, typename RNGType, - typename PointList, + typename MT, typename NT, typename VT > void exponential_correlation_sampling_MT( const unsigned int &n, - PointList &randPoints, + std::list &randCorMatrices, const unsigned int &walkL, const unsigned int &num_points, const VT &c, const NT &T, unsigned int const& nburns = 0){ + using PointList = std::vector; + PointList randPoints; + CorrelationSpectrahedron_MT P(n); const unsigned int d = P.dimension(); PointType startingPoint(n); @@ -148,6 +167,11 @@ void exponential_correlation_sampling_MT( const unsigned int &n, PointType _c(c); exponential_sampling(randPoints, P, rng, walkL, num_points, _c, T, startingPoint, nburns); + + for(const auto&p : randPoints){ + MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n); + randCorMatrices.push_back(final_cor_mat); + } } #endif //VOLESTI_SAMPLING_SAMPLE_CORRELATION_MATRICES_HPP From 26a6e278a7a6e809426af19b83bd1410af1273fe Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Fri, 19 Jul 2024 14:03:25 +0530 Subject: [PATCH 2/5] Update sampling_correlation_matrices_test.cpp --- test/sampling_correlation_matrices_test.cpp | 48 +++++++++++++++++++-- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/test/sampling_correlation_matrices_test.cpp b/test/sampling_correlation_matrices_test.cpp index 39bb0e605..8083ad583 100644 --- a/test/sampling_correlation_matrices_test.cpp +++ b/test/sampling_correlation_matrices_test.cpp @@ -36,6 +36,21 @@ MT rebuildMatrix(const VT &xvector, const unsigned int n){ return mat; } +template +Eigen::Matrix getCoefficientsFromMatrix(const MT& mat) { + int n = mat.rows(); + int d = n * (n - 1) / 2; + Eigen::Matrix coeffs(d); + int k = 0; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < i; ++j) { + coeffs(k) = mat(i, j); + ++k; + } + } + return coeffs; +} + template void check_output(PointList &randPoints, int num_points, int n){ int d = n*(n-1)/2, count = 0; @@ -64,6 +79,33 @@ void check_output(PointList &randPoints, int num_points, int n){ CHECK(score.maxCoeff() < 1.1); } +template +void check_output_MT(std::list &randCorMatrices, int num_points, int n){ + int d = n*(n-1)/2, count = 0; + MT A; + Eigen::LDLT mat_ldlt; + for(auto& mat : randCorMatrices){ + mat_ldlt = Eigen::LDLT(mat); + if(mat_ldlt.info() == Eigen::NumericalIssue || !mat_ldlt.isPositive()){ + ++count; + } + } + std::cout << "Fails " << count << " / " << num_points << " samples\n"; + CHECK(count == 0); + + MT samples(d, num_points); + unsigned int jj = 0; + for(const auto& mat : randCorMatrices){ + samples.col(jj) = getCoefficientsFromMatrix(mat); + jj++; + } + + VT score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; + + CHECK(score.maxCoeff() < 1.1); +} + template void test_corre_spectra_classes(unsigned int const n){ typedef Cartesian Kernel; @@ -136,18 +178,18 @@ void test_new_uniform_MT(const unsigned int n, const unsigned int num_points = 1 std::cout << "Test new sampling 2 : "<< num_points << " uniform correlation matrices of size " << n << std::endl; std::chrono::steady_clock::time_point start, end; double time; - std::vector randPoints; + std::list randCorMatrices; unsigned int walkL = 1; start = std::chrono::steady_clock::now(); - uniform_correlation_sampling_MT(n, randPoints, walkL, num_points, 0); + uniform_correlation_sampling_MT(n, randCorMatrices, walkL, num_points, 0); end = std::chrono::steady_clock::now(); time = std::chrono::duration_cast(end - start).count(); std::cout << "Elapsed time : " << time << " (ms)" << std::endl; - check_output(randPoints, num_points, n); + check_output_MT(randCorMatrices, num_points, n); } int n = 3; From 2e4ee750101a7c419ab566f4619de3fb7293686e Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Fri, 19 Jul 2024 15:12:42 +0530 Subject: [PATCH 3/5] Update sampler.cpp --- examples/correlation_matrices/sampler.cpp | 40 ++++++++--------------- 1 file changed, 13 insertions(+), 27 deletions(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 0655767ea..9faaf9553 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -90,28 +90,6 @@ void write_to_file(std::string filename, std::vector const& randPoint std::cout.rdbuf(coutbuf); } -bool is_correlation_matrix(const MT& matrix, const double tol = 1e-8){ - //check if all the diagonal elements are ones - for(int i=0 ; i tol) - { - return false; - } - } - - //check if the matrix is positive semidefinite - using NT = double; - using MatrixType = Eigen::Matrix; - EigenvaluesProblems> solver; - - if(solver.isPositiveSemidefinite(matrix)) - { - return true; - } - return false; -} - template void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){ @@ -138,12 +116,12 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned std::cout << walkname << " samples uniformly "<< num_points << " correlation matrices of size " << n << " with matrix PointType" << std::endl; std::chrono::steady_clock::time_point start, end; double time; - std::vector randPoints; + std::list randCorMatrices; unsigned int walkL = 1; start = std::chrono::steady_clock::now(); - uniform_correlation_sampling_MT(n, randPoints, walkL, num_points, 0); + uniform_correlation_sampling_MT(n, randCorMatrices, walkL, num_points, 0); end = std::chrono::steady_clock::now(); time = std::chrono::duration_cast(end - start).count(); @@ -151,13 +129,21 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned int valid_points = 0; EigenvaluesProblems> solver; - for(const auto& points : randPoints){ - if(solver.is_correlation_matrix(points.mat)){ + for(const auto& matrix : randCorMatrices){ + if(solver.is_correlation_matrix(matrix)){ valid_points++; - } + } } + std::cout << "Number of valid points = " << valid_points << std::endl; + std::vector randPoints; + for(const auto &mat : randCorMatrices){ + PointMT p; + p.mat = mat; + randPoints.push_back(p); + } + write_to_file(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints); } From a703cac576084e18bc146e63dd618c918af231d5 Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Fri, 19 Jul 2024 22:33:36 +0530 Subject: [PATCH 4/5] use list --- include/sampling/sample_correlation_matrices.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/sampling/sample_correlation_matrices.hpp b/include/sampling/sample_correlation_matrices.hpp index e2b0e7e6b..50b7a45da 100644 --- a/include/sampling/sample_correlation_matrices.hpp +++ b/include/sampling/sample_correlation_matrices.hpp @@ -48,7 +48,7 @@ void uniform_correlation_sampling_MT( const unsigned int &n, const unsigned int &walkL, const unsigned int &num_points, unsigned int const& nburns){ - using PointList = std::vector; + using PointList = std::list; PointList randPoints; CorrelationSpectrahedron_MT P(n); @@ -100,7 +100,7 @@ void gaussian_correlation_sampling_MT( const unsigned int &n, const unsigned int &num_points, const NT &a, unsigned int const& nburns = 0){ - using PointList = std::vector; + using PointList = std::list; PointList randPoints; CorrelationSpectrahedron_MT P(n); @@ -157,7 +157,7 @@ void exponential_correlation_sampling_MT( const unsigned int &n, const VT &c, const NT &T, unsigned int const& nburns = 0){ - using PointList = std::vector; + using PointList = std::list; PointList randPoints; CorrelationSpectrahedron_MT P(n); From 66c242cbb6df2435417bc393a9bd1ae86cce7b99 Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Fri, 19 Jul 2024 22:34:59 +0530 Subject: [PATCH 5/5] fix indentation --- test/sampling_correlation_matrices_test.cpp | 34 ++++++++++----------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/test/sampling_correlation_matrices_test.cpp b/test/sampling_correlation_matrices_test.cpp index 8083ad583..a17c36ce6 100644 --- a/test/sampling_correlation_matrices_test.cpp +++ b/test/sampling_correlation_matrices_test.cpp @@ -81,29 +81,29 @@ void check_output(PointList &randPoints, int num_points, int n){ template void check_output_MT(std::list &randCorMatrices, int num_points, int n){ - int d = n*(n-1)/2, count = 0; - MT A; - Eigen::LDLT mat_ldlt; - for(auto& mat : randCorMatrices){ - mat_ldlt = Eigen::LDLT(mat); + int d = n*(n-1)/2, count = 0; + MT A; + Eigen::LDLT mat_ldlt; + for(auto& mat : randCorMatrices){ + mat_ldlt = Eigen::LDLT(mat); if(mat_ldlt.info() == Eigen::NumericalIssue || !mat_ldlt.isPositive()){ - ++count; + ++count; } - } - std::cout << "Fails " << count << " / " << num_points << " samples\n"; - CHECK(count == 0); + } + std::cout << "Fails " << count << " / " << num_points << " samples\n"; + CHECK(count == 0); - MT samples(d, num_points); - unsigned int jj = 0; - for(const auto& mat : randCorMatrices){ - samples.col(jj) = getCoefficientsFromMatrix(mat); + MT samples(d, num_points); + unsigned int jj = 0; + for(const auto& mat : randCorMatrices){ + samples.col(jj) = getCoefficientsFromMatrix(mat); jj++; - } + } - VT score = univariate_psrf(samples); - std::cout << "psrf = " << score.maxCoeff() << std::endl; + VT score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; - CHECK(score.maxCoeff() < 1.1); + CHECK(score.maxCoeff() < 1.1); } template