From ad5b3bfaad3c7976283e36b5b20139d0919f0124 Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Tue, 30 Jul 2024 22:59:39 +0530 Subject: [PATCH] new function to tune the walkL parameter (#328) * function to tune walkL * fix ess calculation * fix indentation * fix calling functions to ess, psrf; print statements * fix indentation * move getCoefficients() to include/sampling/sample_correlation_matrices * reuse getcoefficientsfrommatrix() * reuse getcoefficientsfrommatrix() * fix walkL tuning function * fix tune_parameter_L() * Update sampler.cpp * fix getcoefficientsfrommatrix() * fix tune_parameter_L() * fix L values --- examples/correlation_matrices/sampler.cpp | 70 +++++++++++++++++++ .../sampling/sample_correlation_matrices.hpp | 15 ++++ test/sampling_correlation_matrices_test.cpp | 15 ---- 3 files changed, 85 insertions(+), 15 deletions(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 9faaf9553..22a9b7f3a 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -16,8 +16,11 @@ #include "cartesian_geom/cartesian_kernel.h" #include "convex_bodies/spectrahedra/spectrahedron.h" #include "random_walks/random_walks.hpp" +#include "random_walks/uniform_accelerated_billiard_walk.hpp" #include "sampling/sample_correlation_matrices.hpp" #include "matrix_operations/EigenvaluesProblems.h" +#include "diagnostics/effective_sample_size.hpp" +#include "diagnostics/univariate_psrf.hpp" typedef double NT; typedef Eigen::Matrix MT; @@ -147,6 +150,73 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned write_to_file(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints); } +template +void tune_parameter_L(const int choice, const std::vector& dimensions, const unsigned int num_points, + const unsigned int walkL=1, const unsigned int nburns=0){ + using NT = typename MT::Scalar; + for (unsigned int n : dimensions) { + std::list randCorMatrices; + int d = n*(n-1)/2; + double _L; + switch(choice){ + case 1: + _L = sqrt(NT(d)); + std::cout << "The value of L = sqrt(d) = " << _L << std::endl; + break; + case 2: + _L = 2 * sqrt(NT(d)); + std::cout << "The value of L = 2 * sqrt(d) = " << _L << std::endl; + break; + case 3: + _L = 4 * sqrt(NT(d)); + std::cout << "The value of L = 4 * sqrt(d) = " << _L << std::endl; + break; + case 4: + _L = NT(d) / 4.0; + std::cout << "The value of L = d / 4.0 = " << _L << std::endl; + break; + default: + _L = NT(d) / 10.0; + std::cout << "The value of L = d / 10.0 = " << _L << std::endl; + break; + } + std::cout<< "Testing L = " << _L << std::endl; + WalkTypePolicy walk(_L); + + std::chrono::steady_clock::time_point start, end; + double time; + start = std::chrono::steady_clock::now(); + + uniform_correlation_sampling_MT(n, randCorMatrices, walkL, num_points, nburns); + + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + std::cout << "Elapsed time : " << time << " (ms) for dimension: "<< n << std::endl; + + MT samples(d, num_points); + unsigned int jj = 0; + for(auto& mat : randCorMatrices){ + samples.col(jj) = getCoefficientsFromMatrix(mat); + jj++; + } + + //calculate psrf + VT psrf = univariate_psrf(samples); + double max_psrf = psrf.maxCoeff(); + std::cout << "PSRF = " << max_psrf << std::endl; + + //calculate ess + unsigned int min_ess = 0; + VT ess_vector = effective_sample_size(samples, min_ess); + std::cout << "Effective Sample Size = " << min_ess << std::endl; + std::cout << "Average Effective Sample Size = " << min_ess/(double)num_points << std::endl; + std::cout << std::endl; + + // Clear the matrices for the next iteration + randCorMatrices.clear(); + } +} + int main(int argc, char const *argv[]){ // To enable Intel MKL, change option USE_MKL to ON in CMakeLists.txt diff --git a/include/sampling/sample_correlation_matrices.hpp b/include/sampling/sample_correlation_matrices.hpp index 50b7a45da..b2bd118b3 100644 --- a/include/sampling/sample_correlation_matrices.hpp +++ b/include/sampling/sample_correlation_matrices.hpp @@ -14,6 +14,21 @@ #include +template +MT 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; +} + // New implementations for sampling correlation matrices with CorrelationSpectrahedron and CorrelationSpectrahedron_MT template diff --git a/test/sampling_correlation_matrices_test.cpp b/test/sampling_correlation_matrices_test.cpp index a17c36ce6..f14d89261 100644 --- a/test/sampling_correlation_matrices_test.cpp +++ b/test/sampling_correlation_matrices_test.cpp @@ -36,21 +36,6 @@ 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;