Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Feature/symmetric correlation matrix #321

Merged
merged 5 commits into from
Jul 22, 2024
Merged

Feature/symmetric correlation matrix #321

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5344df8
modify MT functions
atrayees Jul 19, 2024
26a6e27
Update sampling_correlation_matrices_test.cpp
atrayees Jul 19, 2024
2e4ee75
Update sampler.cpp
atrayees Jul 19, 2024
a703cac
use list
atrayees Jul 19, 2024
66c242c
fix indentation
atrayees Jul 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 13 additions & 27 deletions examples/correlation_matrices/sampler.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why do you need this function? Notice that the class CorreMatrix includes this representation. Look at the member function getCoefficients().

However, it's not clear why this function is needed. Can you comment on this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here in the uniform_correlation_sampling_MT() function

for(const auto&p : randPoints){
        MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n);
    	randCorMatrices.push_back(final_cor_mat);
    }

p.mat, etc. return type MT, so I'm actually storing MT and not CorreMatrix in std::list<MT> randCorMatrices.

but getCoefficients() is a member function of the CorreMatrix class, so in the testing function in the loop

for(const auto& mat : randCorMatrices){
    	    samples.col(jj) = getCoefficientsFromMatrix<NT, MT>(mat);
    	    jj++;
	}

since every member of randCorMatrices is MT, I couldn't call the getCoefficients() function which is why I wrote the separate function to get the coefficients

I thought we need this loop for populating the samples variable

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ic, then let's keep it in the test folder

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){
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

plz fix indentation in this function

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

really sorry for overlooking this. i fixed it

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
Loading