some refactors to save memory

pachadotdev · Jul 16, 2024 · 1daa09a · 1daa09a
1 parent ccf2773
commit 1daa09a
Show file tree

Hide file tree

Showing 26 changed files with 158 additions and 182 deletions.
diff --git a/R/cpp11.R b/R/cpp11.R
@@ -72,10 +72,6 @@ solve_eta2_ <- function(yadj, myadj, offset, eta) {
   .Call(`_capybara_solve_eta2_`, yadj, myadj, offset, eta)
 }
 
-sqrt_ <- function(w) {
-  .Call(`_capybara_sqrt_`, w)
-}
-
 kendall_cor_ <- function(m) {
   .Call(`_capybara_kendall_cor_`, m)
 }

diff --git a/R/internals.R b/R/internals.R
@@ -71,7 +71,6 @@ feglm_fit_ <- function(beta, eta, y, X, wt, k.list, family, control) {
     # Compute weights and dependent variable
     mu.eta <- family[["mu.eta"]](eta)
     w <- (wt * mu.eta^2) / family[["variance"]](mu)
-    w.tilde <- sqrt_(w)
     nu <- (y - mu) / mu.eta
 
     # Centering variables
@@ -81,7 +80,7 @@ feglm_fit_ <- function(beta, eta, y, X, wt, k.list, family, control) {
     # Compute update step and update eta
     # beta.upd <- as.vector(qr.solve(MX * w.tilde, Mnu * w.tilde, epsilon))
     # eta.upd <- nu - as.vector(Mnu - MX %*% beta.upd)
-    beta.upd <- solve_beta_(MX, Mnu, w.tilde, TRUE)
+    beta.upd <- solve_beta_(MX, Mnu, w, TRUE)
     eta.upd <- solve_eta_(MX, Mnu, nu, beta.upd)
 
     # Step-halving with three checks
@@ -101,7 +100,7 @@ feglm_fit_ <- function(beta, eta, y, X, wt, k.list, family, control) {
       val.crit <- family[["valideta"]](eta) && family[["validmu"]](mu)
       imp.crit <- (dev - dev.old) / (0.1 + abs(dev)) <= -dev.tol
       if (dev.crit && val.crit && imp.crit) break
-      rho <- rho / 2.0
+      rho <- rho * 0.5
     }
 
     # Check if step-halving failed (deviance and invalid \eta or \mu)
@@ -226,7 +225,7 @@ feglm_offset_ <- function(object, offset) {
 
     # Centering dependent variable and compute \eta update
     Myadj <- center_variables_(Myadj, yadj, w, k.list, center.tol, 10000L, TRUE)
-    # eta.upd <- yadj - as.vector(Myadj) + offset - eta
+    # eta.upd <- yadj - drop(Myadj) + offset - eta
     eta.upd <- solve_eta2_(yadj, Myadj, offset, eta)
 
     # Step-halving with three checks
@@ -273,7 +272,7 @@ get_index_list_ <- function(k.vars, data) {
 
 # Compute score matrix ----
 
-getScoreMatrix <- function(object) {
+get_score_matrix_ <- function(object) {
   # Extract required quantities from result list
   control <- object[["control"]]
   data <- object[["data"]]

diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
@@ -3,5 +3,5 @@ pkgdown: 2.0.7
 pkgdown_sha: ~
 articles:
   intro: intro.html
-last_built: 2024-07-04T20:13Z
+last_built: 2024-07-16T00:21Z
 
diff --git a/docs/reference/apes.html b/docs/reference/apes.html
diff --git a/docs/reference/bias_corr.html b/docs/reference/bias_corr.html
diff --git a/docs/reference/feglm.html b/docs/reference/feglm.html
diff --git a/docs/reference/felm.html b/docs/reference/felm.html
diff --git a/docs/reference/fenegbin.html b/docs/reference/fenegbin.html
diff --git a/docs/reference/fepoisson.html b/docs/reference/fepoisson.html
diff --git a/docs/reference/vcov.feglm.html b/docs/reference/vcov.feglm.html
diff --git a/docs/reference/vcov.felm.html b/docs/reference/vcov.felm.html
diff --git a/man/kendall_cor_test.Rd b/man/kendall_cor_test.Rd
diff --git a/man/vcov.feglm.Rd b/man/vcov.feglm.Rd
diff --git a/man/vcov.felm.Rd b/man/vcov.felm.Rd
diff --git a/src/00_main.h b/src/00_main.h
@@ -1,21 +1,13 @@
-// #include <omp.h>
-
 #include <algorithm>
 #include <armadillo.hpp>
 #include <cmath>
 #include <cpp11.hpp>
 #include <cpp11armadillo.hpp>
-// #include <iostream>
+#include <numeric>
 #include <vector>
+#include "Rmath.h"
+
+// #include <iostream>
 
 using namespace arma;
 using namespace cpp11;
-
-// helpers used across scripts
-
-#ifndef HELPERS_H
-#define HELPERS_H
-
-uvec as_uvec(const cpp11::integers &x);
-
-#endif  // HELPERS_H
diff --git a/src/01_center_variables.cpp b/src/01_center_variables.cpp
@@ -24,7 +24,6 @@ center_variables_(const doubles_matrix<> &V_r, const doubles &v_sum_r,
   // Auxiliary variables (storage)
   int iter, j, k, p, J;
   double delta, meanj;
-  Mat<double> C(N, P);
   Mat<double> x(N, 1);
   Mat<double> x0(N, 1);
 
@@ -79,9 +78,9 @@ center_variables_(const doubles_matrix<> &V_r, const doubles &v_sum_r,
         break;
       }
     }
-    C.col(p) = x;
+    V.col(p) = x;
   }
 
   // Return matrix with centered variables
-  return as_doubles_matrix(C);
+  return as_doubles_matrix(V);
 }
diff --git a/src/02_get_alpha.cpp b/src/02_get_alpha.cpp
@@ -27,7 +27,7 @@
   }
 
   // Start alternating between normal equations
-  field<Mat<double>> Alpha0(size(Alpha));
+  field<Mat<double>> Alpha0(K);
 
   for (iter = 0; iter < 10000; ++iter) {
     if ((iter % 1000) == 0) {
@@ -42,16 +42,16 @@
 
       for (l = 0; l < K; ++l) {
         if (l != k) {
-          list klist_l = klist[l];
+          const list &klist_l = klist[l];
           for (int j = 0; j < list_sizes[l]; ++j) {
             uvec indexes = as_uvec(as_cpp<integers>(klist_l[j]));
             y(indexes) -= Alpha(l)(j);
           }
         }
       }
 
-      list klist_k = as_cpp<list>(klist[k]);
-      Mat<double> alpha(list_sizes[k], 1);
+      const list &klist_k = as_cpp<list>(klist[k]);
+      Mat<double> &alpha = Alpha(k);
 
       for (int j = 0; j < list_sizes[k]; ++j) {
         // Subset the j-th group of category k
@@ -60,19 +60,15 @@
         // Store group mean
         alpha(j) = mean(y(indexes));
       }
-
-      // Update alpha_k
-      Alpha(k) = alpha;
     }
 
     // Compute termination criterion and check convergence
     num = 0.0;
     denom = 0.0;
     for (k = 0; k < K; ++k) {
-      Mat<double> diff = Alpha(k) - Alpha0(k);
+      const Mat<double> &diff = Alpha(k) - Alpha0(k);
       num += accu(diff % diff);
       denom += accu(Alpha0(k) % Alpha0(k));
-      Alpha0(k) = Alpha(k);
     }
 
     crit = sqrt(num / denom);