correct bug with 3 or more FEs

src/01_center_variables.cpp

@@ -1,10 +1,9 @@
 #include "00_main.h"
 
 // Method of alternating projections (Halperin)
-[[cpp11::register]] doubles_matrix<>
-center_variables_(const doubles_matrix<> &V_r, const doubles &v_sum_r,
-                  const doubles &w_r, const list &klist, const double tol,
-                  const int maxiter, bool sum_v) {
+[[cpp11::register]] doubles_matrix<> center_variables_(
+    const doubles_matrix<> &V_r, const doubles &v_sum_r, const doubles &w_r,
+    const list &klist, const double tol, const int maxiter, bool sum_v) {
   // Types conversion
   Mat<double> V = as_Mat(V_r);
   Mat<double> w = as_Mat(w_r);
@@ -32,13 +31,10 @@ center_variables_(const doubles_matrix<> &V_r, const doubles &v_sum_r,
     // Center each variable
     x = V.col(p);
 
-    int interruptCheckCounter = 0;
-
     for (iter = 0; iter < 100000; ++iter) {
       // Check user interrupt
-      if (++interruptCheckCounter == 1000) {
+      if ((iter % 1000) == 0) {
         check_user_interrupt();
-        interruptCheckCounter = 0;
       }
 
       // Store centered vector from the last iteration

src/02_get_alpha.cpp

@@ -18,27 +18,21 @@
   field<Col<double>> Alpha(K);
   for (k = 0; k < K; k++) {
     J = as_cpp<list>(klist[k]).size();
-    Alpha(k) = zeros(J);
+    Alpha(k) = zeros<Col<double>>(J);
   }
 
   // Start alternating between normal equations
   field<Col<double>> Alpha0(size(Alpha));
 
-  // Create vector for alpha
-  Col<double> alpha(J);
-
-  int interruptCheckCounter = 0;
-
   for (iter = 0; iter < 10000; iter++) {
     // Check user interrupt
-    if (++interruptCheckCounter == 1000) {
+    if ((iter % 1000) == 0) {
       check_user_interrupt();
-      interruptCheckCounter = 0;
     }
-    // Store \alpha_{0} of the previous iteration
+
+    // Store alpha_0 of the previous iteration
     Alpha0 = Alpha;
 
-    // Solve normal equations of category k
     for (k = 0; k < K; k++) {
       // Compute adjusted dependent variable
       y = p;
@@ -50,14 +44,15 @@
             integers indexes = as_cpp<list>(klist[l])[j];
             I = indexes.size();
             for (i = 0; i < I; i++) {
-              y(indexes[i]) -= Alpha(j)(j);
+              y(indexes[i]) -= Alpha(l)(j);
             }
           }
         }
       }
 
-      // Compute group mean
       J = as_cpp<list>(klist[k]).size();
+      Col<double> alpha = zeros<Col<double>>(J);
+
       for (j = 0; j < J; j++) {
         // Subset the j-th group of category k
         integers indexes = as_cpp<list>(klist[k])[j];
@@ -73,7 +68,7 @@
         alpha(j) = sum / I;
       }
 
-      // Update \alpha_{k}
+      // Update alpha_k
       Alpha(k) = alpha;
     }
 

src/Makevars

@@ -1,3 +1,2 @@
 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) 
 PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
-

src/Makevars.win

@@ -1,3 +1,2 @@
 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) 
 PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
-