From 97054723a71d4e54319641f4834b65f3fd1e97fc Mon Sep 17 00:00:00 2001
From: lavakin <nikola@kalabova.eu>
Date: Fri, 18 Aug 2023 11:06:56 +0200
Subject: [PATCH] variance issue fixed

---
 src/rcpp_funcs.cpp | 79 ++++++++++++++++++++++++----------------------
 1 file changed, 42 insertions(+), 37 deletions(-)

diff --git a/src/rcpp_funcs.cpp b/src/rcpp_funcs.cpp
index 55c09f99..922e0247 100644
--- a/src/rcpp_funcs.cpp
+++ b/src/rcpp_funcs.cpp
@@ -39,9 +39,10 @@ int randWrapper(const int& n)
 
 // Initilizing the random number generator outside of the function
 std::random_device rng;
-std::mt19937_64 urng(rng());
+
 std::default_random_engine gn(rng());
 std::minstd_rand mrgn(42);
+std::mt19937_64 urng(rng());
 
 
 void updateProgressBar(int currentProgress, int totalProgress, int barWidth = 40) {
@@ -74,47 +75,51 @@ Eigen::MatrixXd permut_mat(const Eigen::VectorXd& a,const int& permutations) {
   // already added by sourceCpp(), but needed standalone
   Eigen::MatrixXd permutedMat(permutations, a.size());
   Eigen::VectorXd shuffledVec = a;
-
+  
   const int updateFrequency = 200;  // Update progress every x iterations
-  #ifdef _OPENMP  
+#ifdef _OPENMP  
   std::atomic<int> progress(0);
   
-  #pragma omp parallel
-  {
-    int localProgress = 0;
-    #pragma omp for
-      for (int i = 0; i < permutations; i++) {
-        shuffledVec = a;
-        std::shuffle(shuffledVec.data(), shuffledVec.data() + shuffledVec.size(), urng);
-        permutedMat.row(i) = shuffledVec.transpose();
-        localProgress++;
-        
-        
-        if (localProgress % updateFrequency == 0) {
-          progress.fetch_add(updateFrequency, std::memory_order_relaxed);
-          
-          // Display progress
-          
-          int currentProgress = progress.load(std::memory_order_relaxed);
-          
-          updateProgressBar(currentProgress,permutations);
-          }
-      }
-    }
+#pragma omp parallel
+{
+  int localProgress = 0;
+  std::mt19937_64 urngp(rng());
+  Eigen::VectorXd shuff = a;
+  Eigen::VectorXd shuffl = shuff;
+
+#pragma omp for
+  for (int i = 0; i < permutations; i++) {
+    shuffl = shuff;
+    std::shuffle(shuffl.data(), shuffl.data() + shuffl.size(), urngp);
+    permutedMat.row(i) = shuffl.transpose();
+    localProgress++;
     
-  #else
-    for (int i = 0; i < permutations; i++) {
-      shuffledVec = a;
-      std::shuffle(shuffledVec.data(), shuffledVec.data() + shuffledVec.size(), urng);
-      permutedMat.row(i) = shuffledVec.transpose();
-      if (i % updateFrequency == 0){
-        updateProgressBar(i,permutations);
-      }
-        
+    
+    if (localProgress % updateFrequency == 0) {
+      progress.fetch_add(updateFrequency, std::memory_order_relaxed);
+      
+      // Display progress
+      
+      int currentProgress = progress.load(std::memory_order_relaxed);
+      
+      updateProgressBar(currentProgress,permutations);
     }
-  #endif
-    updateProgressBar(permutations,permutations);
-    return permutedMat;
+  }
+}
+
+#else
+for (int i = 0; i < permutations; i++) {
+  shuffledVec = a;
+  std::shuffle(shuffledVec.data(), shuffledVec.data() + shuffledVec.size(), urng);
+  permutedMat.row(i) = shuffledVec.transpose();
+  if (i % updateFrequency == 0){
+    updateProgressBar(i,permutations);
+  }
+  
+}
+#endif
+updateProgressBar(permutations,permutations);
+return permutedMat;
 }
 
 NumericVector permut(const NumericVector& a)