Update sampling_correlation_matrices_test.cpp

test/sampling_correlation_matrices_test.cpp

@@ -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;
@@ -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;
@@ -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;