tidy cpp code

Author: pachadotdev

configure

@@ -1,16 +1,106 @@
-#!/bin/sh
-
-# Set the number of cores
-if [ -z "$PKG_NCORES" ]; then
-  PKG_NCORES=$(getconf _NPROCESSORS_ONLN)
-    if [ -z "$PKG_NCORES" ] || [ "$PKG_NCORES" -lt 1 ]; then
-      PKG_NCORES=1
-    elif [ "$PKG_NCORES" -gt 1 ]; then
-      PKG_NCORES=$((PKG_NCORES - 1))
-    fi
+# Anticonf script by Pacha (2025)
+
+PKG_CONFIG_NAME="capybara"
+
+# Parse command line arguments
+for arg in "$@"; do
+  case "$arg" in
+    --enable-optimization)
+      ENABLE_OPTIMIZATION=yes
+      ;;
+  esac
+done
+
+# Find compiler
+CXX=$(${R_HOME}/bin/R CMD config CXX)
+CXXFLAGS=$(${R_HOME}/bin/R CMD config CXXFLAGS)
+CPPFLAGS=$(${R_HOME}/bin/R CMD config CPPFLAGS)
+
+# Set optimization flags if requested
+OPTIMIZATION_FLAGS=""
+if [ "$ENABLE_OPTIMIZATION" = "yes" ]; then
+  OPTIMIZATION_FLAGS="-O3 -funroll-loops"
+  
+  # Check for AVX2 support
+  if grep -q avx2 /proc/cpuinfo 2>/dev/null; then
+    OPTIMIZATION_FLAGS="$OPTIMIZATION_FLAGS -mavx2"
+  fi
+
+  echo "Using optimization flags: $OPTIMIZATION_FLAGS"
+else
+  echo "Not using optimization flags"
+fi
+
+# Check for OpenMP support
+OPENMP_FLAG=""
+for flag in -fopenmp -qopenmp -openmp -xopenmp; do
+  save_CXXFLAGS="$CXXFLAGS"
+  CXXFLAGS="$CXXFLAGS $flag"
+  
+  echo "int main() { return 0; }" | $CXX -xc++ -E $CXXFLAGS - >/dev/null 2>&1
+  if [ $? -eq 0 ]; then
+    OPENMP_FLAG="$flag"
+    break
+  fi
+  CXXFLAGS="$save_CXXFLAGS"
+done
+
+if [ -n "$OPENMP_FLAG" ]; then
+  echo "OpenMP is enabled: $OPENMP_FLAG"
+  PKG_CXXFLAGS="$OPENMP_FLAG $PKG_CXXFLAGS"
+  PKG_LIBS="$OPENMP_FLAG $PKG_LIBS"
+else
+  echo "OpenMP is disabled"
 fi
 
-# Write Makevars
-sed -e "s|@ncores@|$PKG_NCORES|" \
-  src/Makevars.in > src/Makevars
-  
+# Determine number of cores using OpenMP if available
+if [ -n "$CAPYBARA_NCORES" ]; then
+  num_cores="$CAPYBARA_NCORES"
+  echo "using environment value: $num_cores"
+else
+  cat > nthreads.cpp << 'EOL'
+#include <iostream>
+#ifdef _OPENMP
+  #include <omp.h>
+#endif
+
+int main() {
+  int nthreads;
+  #ifdef _OPENMP
+    nthreads = std::max(1, (omp_get_max_threads() + 1) / 2);
+  #else
+    nthreads = 1;
+  #endif
+  std::cout << nthreads << std::endl;
+  return 0;
+}
+EOL
+
+  # Need to compile WITH the OpenMP flag
+  if [ -n "$OPENMP_FLAG" ]; then
+    $CXX nthreads.cpp -o nthreads $OPENMP_FLAG
+  else
+    $CXX nthreads.cpp -o nthreads
+  fi
+  
+  # Check if compilation succeeded
+  if [ $? -eq 0 ] && [ -x ./nthreads ]; then
+    num_cores=$(./nthreads)
+    echo "Number of available cores: $num_cores"
+  else
+    echo "OpenMP test failed, defaulting to 1"
+    num_cores=1
+  fi
+
+  # Clean up
+  rm -f nthreads nthreads.cpp
+fi
+
+# Write to Makevars
+echo "creating src/Makevars"
+sed -e "s|@ncores@|${num_cores}|g" \
+    -e "s|@OPTIMIZATION_FLAGS@|${OPTIMIZATION_FLAGS}|g" \
+    src/Makevars.in > src/Makevars
+
+# Success
+exit 0

src/00_main.h

@@ -1,3 +1,4 @@
+#include <armadillo.hpp>
 #include <cpp11.hpp>
 #include <cpp11armadillo.hpp>
 #include <regex>
@@ -8,44 +9,26 @@ using namespace cpp11;
 
 // used across the scripts
 
-#ifdef _OPENMP
-const size_t n_threads = omp_get_max_threads();
-#endif
-
 void center_variables_(mat &V, const vec &w, const list &klist,
-                       const double &tol, const size_t &maxiter,
-                       const size_t &interrupt_iter);
+                       const double &tol, const int &maxiter);
 
 vec solve_beta_(mat MX, const mat &MNU, const vec &w);
 
 vec solve_eta_(const mat &MX, const mat &MNU, const vec &nu, const vec &beta);
 
 mat crossprod_(const mat &X, const vec &w);
 
-// Enum for GLM family types
-enum FamilyType {
-  GAUSSIAN,
-  POISSON,
-  BINOMIAL,
-  GAMMA,
-  INV_GAUSSIAN,
-  NEG_BIN,
-  UNKNOWN
-};
-
 std::string tidy_family_(const std::string &family);
 
-FamilyType get_family_type(const std::string &fam);
-
-vec link_inv_(const vec &eta, const FamilyType &fam);
+vec link_inv_(const vec &eta, const std::string &fam);
 
 double dev_resids_(const vec &y, const vec &mu, const double &theta,
-                   const vec &wt, const FamilyType &fam);
+                   const vec &wt, const std::string &fam);
 
-vec mu_eta_(const vec &eta, const FamilyType &fam);
+vec mu_eta_(const vec &eta, const std::string &fam);
 
-vec variance_(const vec &mu, const double &theta, const FamilyType &fam);
+vec variance_(const vec &mu, const double &theta, const std::string &fam);
 
-bool valid_eta_(const vec &eta, const FamilyType &fam);
+bool valid_eta_(const vec &eta, const std::string &fam);
 
-bool valid_mu_(const vec &mu, const FamilyType &fam);
+bool valid_mu_(const vec &mu, const std::string &fam);

src/01_linear_algebra.cpp

@@ -6,7 +6,8 @@ mat crossprod_(const mat &X, const vec &w) {
   return Y.t() * Y;
 }
 
-// vec solve_beta_(mat MX, const mat &MNU, const vec &w) {
+// vec solve_beta_(mat MX, const mat &MNU,
+//                         const vec &w) {
 //   const vec sqrt_w = sqrt(w);
 //   MX.each_col() %= sqrt_w;
 
@@ -15,27 +16,29 @@ mat crossprod_(const mat &X, const vec &w) {
 //     stop("QR decomposition failed");
 //   }
 
-//   return solve(trimatu(R), Q.t() * (MNU.each_col() % sqrt_w), solve_opts::fast);
+//   return solve(trimatu(R), Q.t() * (MNU.each_col() % sqrt_w),
+//   solve_opts::fast);
 // }
 
+// Cholesky decomposition
 vec solve_beta_(mat MX, const mat &MNU, const vec &w) {
   const vec sqrt_w = sqrt(w);
+
   MX.each_col() %= sqrt_w;
+  mat WMNU = MNU.each_col() % sqrt_w;
+
   mat XtX = MX.t() * MX;
+  vec XtY = MX.t() * (MNU.each_col() % sqrt_w);
 
-  // Cholesky decomposition: XtX = L * L.t()
+  // XtX = L * L.t()
   mat L;
   if (!chol(L, XtX, "lower")) {
     stop("Cholesky decomposition failed.");
   }
 
-  vec Xty = MX.t() * (MNU.each_col() % sqrt_w);
-
   // Solve L * z = Xty
-  vec z = solve(trimatl(L), Xty, solve_opts::fast);
-
-  // Solve L.t() * beta = z
-  vec beta = solve(trimatu(L.t()), z, solve_opts::fast);
+  vec z = solve(trimatl(L), XtY, solve_opts::fast);
 
-  return beta;
+  // Solve Lt * beta = z
+  return solve(trimatu(L.t()), z, solve_opts::fast);
 }