Skip to content

Commit

Permalink
Feature/symmetric correlation matrix (#321)
Browse files Browse the repository at this point in the history 
* modify MT functions

* Update sampling_correlation_matrices_test.cpp

* Update sampler.cpp

* use list

* fix indentation
  • Loading branch information
@atrayees
atrayees authored Jul 22, 2024
1 parent 3325d6a commit c2504a5
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 36 deletions.
40 changes: 13 additions & 27 deletions examples/correlation_matrices/sampler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,28 +90,6 @@ void write_to_file(std::string filename, std::vector<PointType> 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<matrix.rows() ; i++)
{
if(std::abs(matrix(i, i)-1.0) > tol)
{
return false;
}
}

//check if the matrix is positive semidefinite
using NT = double;
using MatrixType = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>;
EigenvaluesProblems<NT, MatrixType, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;

if(solver.isPositiveSemidefinite(matrix))
{
return true;
}
return false;
}

template<typename WalkType>
void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){

Expand All @@ -138,26 +116,34 @@ 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<PointMT> randPoints;
std::list<MT> randCorMatrices;
unsigned int walkL = 1;

start = std::chrono::steady_clock::now();

uniform_correlation_sampling_MT<WalkType, PointMT, RNGType>(n, randPoints, walkL, num_points, 0);
uniform_correlation_sampling_MT<WalkType, PointMT, RNGType>(n, randCorMatrices, walkL, num_points, 0);

end = std::chrono::steady_clock::now();
time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << "Elapsed time : " << time << " (ms)" << std::endl;

int valid_points = 0;
EigenvaluesProblems<NT, MT, Eigen::Matrix<NT, Eigen::Dynamic, 1>> 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<PointMT> randPoints;
for(const auto &mat : randCorMatrices){
PointMT p;
p.mat = mat;
randPoints.push_back(p);
}

write_to_file<PointMT>(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints);
}

Expand Down
36 changes: 30 additions & 6 deletions include/sampling/sample_correlation_matrices.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<MT> &randCorMatrices,
const unsigned int &walkL,
const unsigned int &num_points,
unsigned int const& nburns){
using PointList = std::list<PointType>;
PointList randPoints;

CorrelationSpectrahedron_MT<PointType> P(n);
const unsigned int d = P.dimension();
PointType startingPoint(n);
RNGType rng(d);

uniform_sampling<WalkTypePolicy>(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
Expand Down Expand Up @@ -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<MT> &randCorMatrices,
const unsigned int &walkL,
const unsigned int &num_points,
const NT &a,
unsigned int const& nburns = 0){
using PointList = std::list<PointType>;
PointList randPoints;

CorrelationSpectrahedron_MT<PointType> P(n);
const unsigned int d = P.dimension();
PointType startingPoint(n);
RNGType rng(d);

gaussian_sampling<WalkTypePolicy>(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
Expand Down Expand Up @@ -130,24 +146,32 @@ 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<MT> &randCorMatrices,
const unsigned int &walkL,
const unsigned int &num_points,
const VT &c,
const NT &T,
unsigned int const& nburns = 0){
using PointList = std::list<PointType>;
PointList randPoints;

CorrelationSpectrahedron_MT<PointType> P(n);
const unsigned int d = P.dimension();
PointType startingPoint(n);
RNGType rng(d);
PointType _c(c);

exponential_sampling<WalkTypePolicy>(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
48 changes: 45 additions & 3 deletions test/sampling_correlation_matrices_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,21 @@ MT rebuildMatrix(const VT &xvector, const unsigned int n){
return mat;
}

template<typename NT, typename MT>
Eigen::Matrix<NT, Eigen::Dynamic, 1> getCoefficientsFromMatrix(const MT& mat) {
int n = mat.rows();
int d = n * (n - 1) / 2;
Eigen::Matrix<NT, Eigen::Dynamic, 1> 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<typename NT, typename VT, typename MT, typename PointList>
void check_output(PointList &randPoints, int num_points, int n){
int d = n*(n-1)/2, count = 0;
Expand Down Expand Up @@ -64,6 +79,33 @@ void check_output(PointList &randPoints, int num_points, int n){
CHECK(score.maxCoeff() < 1.1);
}

template<typename NT, typename VT, typename MT>
void check_output_MT(std::list<MT> &randCorMatrices, int num_points, int n){
int d = n*(n-1)/2, count = 0;
MT A;
Eigen::LDLT<MT> mat_ldlt;
for(auto& mat : randCorMatrices){
mat_ldlt = Eigen::LDLT<MT>(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<NT, MT>(mat);
jj++;
}

VT score = univariate_psrf<NT, VT>(samples);
std::cout << "psrf = " << score.maxCoeff() << std::endl;

CHECK(score.maxCoeff() < 1.1);
}

template <typename NT>
void test_corre_spectra_classes(unsigned int const n){
typedef Cartesian<NT> Kernel;
Expand Down Expand Up @@ -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<Point> randPoints;
std::list<MT> randCorMatrices;
unsigned int walkL = 1;

start = std::chrono::steady_clock::now();

uniform_correlation_sampling_MT<WalkType, Point, RNGType>(n, randPoints, walkL, num_points, 0);
uniform_correlation_sampling_MT<WalkType, Point, RNGType>(n, randCorMatrices, walkL, num_points, 0);

end = std::chrono::steady_clock::now();
time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << "Elapsed time : " << time << " (ms)" << std::endl;

check_output<NT, VT, MT>(randPoints, num_points, n);
check_output_MT<NT, VT, MT>(randCorMatrices, num_points, n);
}

int n = 3;
Expand Down

0 comments on commit c2504a5

Please sign in to comment.