From 1daa09ab5631620dcfe3e6ecacda3e5899c19849 Mon Sep 17 00:00:00 2001 From: Mauricio 'Pacha' Vargas Sepulveda Date: Mon, 15 Jul 2024 20:48:04 -0400 Subject: [PATCH] some refactors to save memory --- R/cpp11.R | 4 - R/internals.R | 9 +- docs/pkgdown.yml | 2 +- docs/reference/apes.html | 8 +- docs/reference/bias_corr.html | 4 +- docs/reference/feglm.html | 4 +- docs/reference/felm.html | 2 +- docs/reference/fenegbin.html | 2 +- docs/reference/fepoisson.html | 2 +- docs/reference/vcov.feglm.html | 12 +-- docs/reference/vcov.felm.html | 12 +-- man/kendall_cor_test.Rd | 2 +- man/vcov.feglm.Rd | 2 +- man/vcov.felm.Rd | 2 +- src/00_main.h | 16 +--- src/01_center_variables.cpp | 5 +- src/02_get_alpha.cpp | 14 ++- src/03_group_sums.cpp | 88 ++++++++++--------- src/04_linear_algebra.cpp | 43 ++------- src/05_kendall_correlation.cpp | 9 +- src/cpp11.cpp | 12 +-- .../{test-apes.R => test-apes-bias.R} | 0 tests/testthat/test-cpp.R | 28 ------ tests/testthat/test-fepoisson.R | 6 ++ tests/testthat/test-group-sums.R | 32 +++++++ tests/testthat/test-linear_algebra.R | 20 +++++ 26 files changed, 158 insertions(+), 182 deletions(-) rename tests/testthat/{test-apes.R => test-apes-bias.R} (100%) delete mode 100644 tests/testthat/test-cpp.R create mode 100644 tests/testthat/test-group-sums.R diff --git a/R/cpp11.R b/R/cpp11.R index db83d54..8bf0895 100644 --- 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 index bbba560..68279b5 100644 --- 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 index 4397510..82a9933 100644 --- 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 index 6bf4b9d..bec9e49 100644 --- a/docs/reference/apes.html +++ b/docs/reference/apes.html @@ -178,7 +178,7 @@

Examples

summary(mod_ape) #> Estimates: #> Estimate Std. Error z value Pr(>|z|) -#> lang 0.05594 0.01513 3.698 0.000218 *** +#> lang 0.05594 0.01512 3.699 0.000217 *** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 @@ -186,7 +186,7 @@

Examples

mod_bc <- bias_corr(mod) summary(mod_bc) #> Formula: trade ~ lang | year -#> <environment: 0x5feef0344168> +#> <environment: 0x6194fd0f13a0> #> #> Family: Binomial #> @@ -194,7 +194,7 @@

Examples

#> #> | | Estimate | Std. Error | z value | Pr(>|z|) | #> |------|----------|------------|---------|------------| -#> | lang | 0.2393 | 0.0634 | 3.7730 | 0.0002 *** | +#> | lang | 0.2393 | 0.0634 | 3.7724 | 0.0002 *** | #> #> Significance codes: *** 99.9%; ** 99%; * 95%; . 90% #> @@ -207,7 +207,7 @@

Examples

summary(mod_ape_bc) #> Estimates: #> Estimate Std. Error z value Pr(>|z|) -#> lang 0.05594 0.01513 3.698 0.000218 *** +#> lang 0.05593 0.01512 3.698 0.000217 *** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 diff --git a/docs/reference/bias_corr.html b/docs/reference/bias_corr.html index d8028b8..fba8378 100644 --- a/docs/reference/bias_corr.html +++ b/docs/reference/bias_corr.html @@ -146,7 +146,7 @@

Examples

mod_bc <- bias_corr(mod) summary(mod_bc) #> Formula: trade ~ lang | year -#> <environment: 0x5feef0339918> +#> <environment: 0x6194fd0d1c00> #> #> Family: Binomial #> @@ -154,7 +154,7 @@

Examples

#> #> | | Estimate | Std. Error | z value | Pr(>|z|) | #> |------|----------|------------|---------|------------| -#> | lang | 0.2393 | 0.0634 | 3.7730 | 0.0002 *** | +#> | lang | 0.2393 | 0.0634 | 3.7724 | 0.0002 *** | #> #> Significance codes: *** 99.9%; ** 99%; * 95%; . 90% #> diff --git a/docs/reference/feglm.html b/docs/reference/feglm.html index edea9b2..c395427 100644 --- a/docs/reference/feglm.html +++ b/docs/reference/feglm.html @@ -162,7 +162,7 @@

Examples

summary(mod) #> Formula: trade ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5feef175de60> +#> <environment: 0x6194fe4f6ed0> #> #> Family: Poisson #> @@ -192,7 +192,7 @@

Examples

summary(mod, type = "clustered") #> Formula: trade ~ log_dist + lang + cntg + clny | exp_year + imp_year | #> pair -#> <environment: 0x5feef175de60> +#> <environment: 0x6194fe4f6ed0> #> #> Family: Poisson #> diff --git a/docs/reference/felm.html b/docs/reference/felm.html index 790db2a..53de5a3 100644 --- a/docs/reference/felm.html +++ b/docs/reference/felm.html @@ -116,7 +116,7 @@

Examples

summary(mod) #> Formula: log(trade) ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5feef1816da8> +#> <environment: 0x6194fe51d928> #> #> Estimates: #> diff --git a/docs/reference/fenegbin.html b/docs/reference/fenegbin.html index 7d51e2a..5170cfd 100644 --- a/docs/reference/fenegbin.html +++ b/docs/reference/fenegbin.html @@ -141,7 +141,7 @@

Examples

summary(mod) #> Formula: trade ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5feef1c5bdd0> +#> <environment: 0x6194fec80d60> #> #> Family: Negative Binomial(1.1839) #> diff --git a/docs/reference/fepoisson.html b/docs/reference/fepoisson.html index 9d50d2c..13dc8c5 100644 --- a/docs/reference/fepoisson.html +++ b/docs/reference/fepoisson.html @@ -128,7 +128,7 @@

Examples

summary(mod) #> Formula: trade ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5feef19487d8> +#> <environment: 0x6194fe9e6d30> #> #> Family: Poisson #> diff --git a/docs/reference/vcov.feglm.html b/docs/reference/vcov.feglm.html index 51a4251..3907120 100644 --- a/docs/reference/vcov.feglm.html +++ b/docs/reference/vcov.feglm.html @@ -120,12 +120,12 @@

Examples

trade_panel ) -vcov(mod, type = "clustered") -#> [,1] [,2] [,3] [,4] -#> [1,] 0.0006847158 0.0000482318 0.0010614958 -0.0006316697 -#> [2,] 0.0000482318 0.0038793033 -0.0011847203 -0.0019330635 -#> [3,] 0.0010614958 -0.0011847203 0.0045202435 -0.0001310871 -#> [4,] -0.0006316697 -0.0019330635 -0.0001310871 0.0085389097 +round(vcov(mod, type = "clustered"), 5) +#> [,1] [,2] [,3] [,4] +#> [1,] 0.00068 0.00005 0.00106 -0.00063 +#> [2,] 0.00005 0.00388 -0.00118 -0.00193 +#> [3,] 0.00106 -0.00118 0.00452 -0.00013 +#> [4,] -0.00063 -0.00193 -0.00013 0.00854 diff --git a/docs/reference/vcov.felm.html b/docs/reference/vcov.felm.html index 62b5795..fb77315 100644 --- a/docs/reference/vcov.felm.html +++ b/docs/reference/vcov.felm.html @@ -120,12 +120,12 @@

Examples

trade_panel ) -vcov(mod, type = "clustered") -#> [,1] [,2] [,3] [,4] -#> [1,] 0.0006847158 0.0000482318 0.0010614958 -0.0006316697 -#> [2,] 0.0000482318 0.0038793033 -0.0011847203 -0.0019330635 -#> [3,] 0.0010614958 -0.0011847203 0.0045202435 -0.0001310871 -#> [4,] -0.0006316697 -0.0019330635 -0.0001310871 0.0085389097 +round(vcov(mod, type = "clustered"), 5) +#> [,1] [,2] [,3] [,4] +#> [1,] 0.00068 0.00005 0.00106 -0.00063 +#> [2,] 0.00005 0.00388 -0.00118 -0.00193 +#> [3,] 0.00106 -0.00118 0.00452 -0.00013 +#> [4,] -0.00063 -0.00193 -0.00013 0.00854 diff --git a/man/kendall_cor_test.Rd b/man/kendall_cor_test.Rd index 0bec3a2..b817470 100644 --- a/man/kendall_cor_test.Rd +++ b/man/kendall_cor_test.Rd @@ -4,7 +4,7 @@ \alias{kendall_cor_test} \title{Kendall Correlation Test} \usage{ -kendall_cor_test(x, y, alternative = "two.sided") +kendall_cor_test(x, y, alternative = c("two.sided", "greater", "less")) } \arguments{ \item{x}{a numeric vector.} diff --git a/man/vcov.feglm.Rd b/man/vcov.feglm.Rd index c14d92f..d81d66f 100644 --- a/man/vcov.feglm.Rd +++ b/man/vcov.feglm.Rd @@ -38,7 +38,7 @@ mod <- fepoisson( trade_panel ) -vcov(mod, type = "clustered") +round(vcov(mod, type = "clustered"), 5) } \references{ diff --git a/man/vcov.felm.Rd b/man/vcov.felm.Rd index 22239f7..e97c1de 100644 --- a/man/vcov.felm.Rd +++ b/man/vcov.felm.Rd @@ -38,7 +38,7 @@ mod <- fepoisson( trade_panel ) -vcov(mod, type = "clustered") +round(vcov(mod, type = "clustered"), 5) } \references{ diff --git a/src/00_main.h b/src/00_main.h index ac702ea..cf902ac 100644 --- a/src/00_main.h +++ b/src/00_main.h @@ -1,21 +1,13 @@ -// #include - #include #include #include #include #include -// #include +#include #include +#include "Rmath.h" + +// #include 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 index 9554f74..66250d1 100644 --- 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 C(N, P); Mat x(N, 1); Mat 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 index 560ca63..d95736c 100644 --- a/src/02_get_alpha.cpp +++ b/src/02_get_alpha.cpp @@ -27,7 +27,7 @@ } // Start alternating between normal equations - field> Alpha0(size(Alpha)); + field> Alpha0(K); for (iter = 0; iter < 10000; ++iter) { if ((iter % 1000) == 0) { @@ -42,7 +42,7 @@ 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(klist_l[j])); y(indexes) -= Alpha(l)(j); @@ -50,8 +50,8 @@ } } - list klist_k = as_cpp(klist[k]); - Mat alpha(list_sizes[k], 1); + const list &klist_k = as_cpp(klist[k]); + Mat &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 diff = Alpha(k) - Alpha0(k); + const Mat &diff = Alpha(k) - Alpha0(k); num += accu(diff % diff); denom += accu(Alpha0(k) % Alpha0(k)); - Alpha0(k) = Alpha(k); } crit = sqrt(num / denom); diff --git a/src/03_group_sums.cpp b/src/03_group_sums.cpp index 664d548..8600533 100644 --- a/src/03_group_sums.cpp +++ b/src/03_group_sums.cpp @@ -12,31 +12,29 @@ const int P = M.n_cols; // Auxiliary variables (storage) - int i, j, p; - Mat num(P, 1, fill::zeros); - integers indexes; + int i, j, I; + Mat b(P, 1, fill::zeros); + Mat num(P, 1); + uvec indexes; + double denom; // Compute sum of weighted group sums - double denom = 0.0; - for (j = 0; j < J; j++) { - uvec arma_indexes = as_uvec(as_cpp(jlist[j])); - int I = arma_indexes.size(); + for (j = 0; j < J; ++j) { + denom = 0.0; + + indexes = as_uvec(as_cpp(jlist[j])); + I = indexes.size(); num.zeros(); - for (p = 0; p < P; ++p) { - for (i = 0; i < I; i++) { - num(p, 0) += M(arma_indexes[i], p); - } + for (i = 0; i < I; ++i) { + num += M.row(indexes[i]).t(); + denom += w(indexes[i]); } - for (i = 0; i < I; i++) { - denom += w(arma_indexes[i]); - } + b += num / denom; } - num = num / denom; - - return as_doubles_matrix(num); + return as_doubles_matrix(b); } [[cpp11::register]] doubles_matrix<> @@ -53,25 +51,33 @@ group_sums_spectral_(const doubles_matrix<> &M_r, const doubles_matrix<> &v_r, const int P = M.n_cols; // Auxiliary variables (storage) - int j; - Mat num(P, 1, fill::zeros); + int i, j, k, I; + Mat b(P, 1, fill::zeros); + Mat num(P, 1); + double denom; // Compute sum of weighted group sums - double denom = 0.0; for (j = 0; j < J; j++) { - uvec arma_indexes = as_uvec(as_cpp(jlist[j])); - // arma_indexes -= 1; + uvec indexes = as_uvec(as_cpp(jlist[j])); + I = indexes.size(); + + num.zeros(); + denom = 0.0; - Mat M_sub = M.rows(arma_indexes); - Mat w_sub = w.rows(arma_indexes); + for (i = 1; i < I; ++i) { + for (k = 1; k <= K; ++k) { + num += M.row(indexes[i]) * v(indexes[i - k], 0) * I / (I - 1); + } + } - num += sum(M_sub.each_col() % w_sub, 0).t(); - denom += accu(w_sub); - } + for (i = 0; i < I; ++i) { + denom += w(indexes[i]); + } - num = num / denom; + b += num / denom; + } - return as_doubles_matrix(num); + return as_doubles_matrix(b); } [[cpp11::register]] doubles_matrix<> @@ -85,15 +91,15 @@ group_sums_var_(const doubles_matrix<> &M_r, const list &jlist) { // Auxiliary variables (storage) int j; + Mat v(P, 1); Mat V(P, P, fill::zeros); // Compute covariance matrix for (j = 0; j < J; ++j) { - uvec arma_indexes = as_uvec(as_cpp(jlist[j])); - // arma_indexes -= 1; + uvec indexes = as_uvec(as_cpp(jlist[j])); - Mat M_sub = M.rows(arma_indexes); - vec v = sum(M_sub, 0).t(); + Mat M_sub = M.rows(indexes); + v = sum(M_sub, 0).t(); V += v * v.t(); } @@ -112,19 +118,19 @@ group_sums_cov_(const doubles_matrix<> &M_r, const doubles_matrix<> &N_r, const int P = M.n_cols; // Auxiliary variables (storage) - int j, p, q; + int j; + size_t i, k, I; + uvec indexes; Mat V(P, P, fill::zeros); // Compute covariance matrix for (j = 0; j < J; ++j) { - uvec arma_indexes = as_uvec(as_cpp(jlist[j])); - - Mat M_sub = M.rows(arma_indexes); - Mat N_sub = N.rows(arma_indexes); + indexes = as_uvec(as_cpp(jlist[j])); + I = indexes.n_elem; - for (p = 0; p < P; p++) { - for (q = 0; q < P; q++) { - V(q, p) += accu(M_sub.col(q) * N_sub.col(p).t()); + for (i = 0; i < I; ++i) { + for (k = i + 1; k < I; ++k) { + V += M.row(indexes[i]).t() * N.row(indexes[k]); } } } diff --git a/src/04_linear_algebra.cpp b/src/04_linear_algebra.cpp index c44cc6e..07d1619 100644 --- a/src/04_linear_algebra.cpp +++ b/src/04_linear_algebra.cpp @@ -148,6 +148,7 @@ update_beta_eta_(const doubles &old, const doubles &upd, const double ¶m) { // Weight the X and Y matrices if (weighted) { Mat w = as_Mat(wtilde); + w = sqrt(w); X = X.each_col() % w; // element-wise multiplication Y = Y.each_col() % w; } @@ -183,40 +184,12 @@ update_beta_eta_(const doubles &old, const doubles &upd, const double ¶m) { // eta.upd <- yadj - as.vector(Myadj) + offset - eta -[[cpp11::register]] doubles solve_eta2_(const SEXP &yadj, const SEXP &myadj, - const SEXP &offset, const SEXP &eta) { - int N = Rf_length(yadj); - writable::doubles res(N); - - double *Yadj_data = REAL(yadj); - double *Myadj_data = REAL(myadj); - double *Offset_data = REAL(offset); - double *Eta_data = REAL(eta); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) -#endif - for (int n = 0; n < N; ++n) { - res[n] = Yadj_data[n] - Myadj_data[n] + Offset_data[n] - Eta_data[n]; - } - - return res; -} +[[cpp11::register]] doubles solve_eta2_(const doubles &yadj, const doubles_matrix<> &myadj, + const doubles &offset, const doubles &eta) { + Mat Yadj = as_Mat(yadj); + Mat Myadj = as_Mat(myadj); + Mat Offset = as_Mat(offset); + Mat Eta = as_Mat(eta); -// w <- sqrt(w) - -[[cpp11::register]] doubles sqrt_(const SEXP &w) { - int n = Rf_length(w); - writable::doubles res(n); - - double *w_data = REAL(w); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) -#endif - for (int i = 0; i < n; ++i) { - res[i] = sqrt(w_data[i]); - } - - return res; + return as_doubles(Yadj - Myadj + Offset - Eta); } diff --git a/src/05_kendall_correlation.cpp b/src/05_kendall_correlation.cpp index b22b7fc..12b1c9d 100644 --- a/src/05_kendall_correlation.cpp +++ b/src/05_kendall_correlation.cpp @@ -6,14 +6,7 @@ // note: the len < 2 conditions are commented out because the R function checks // for this condition before calling the C++ functions -#include -#include -#include -#include -#include -#include "Rmath.h" - -using namespace cpp11; +#include "00_main.h" uint64_t insertion_sort_(double *arr, size_t len) { // if (len < 2) { diff --git a/src/cpp11.cpp b/src/cpp11.cpp index db176ce..2c74811 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -125,17 +125,10 @@ extern "C" SEXP _capybara_solve_eta_(SEXP mx, SEXP mnu, SEXP nu, SEXP beta) { END_CPP11 } // 04_linear_algebra.cpp -doubles solve_eta2_(const SEXP & yadj, const SEXP & myadj, const SEXP & offset, const SEXP & eta); +doubles solve_eta2_(const doubles & yadj, const doubles_matrix<> & myadj, const doubles & offset, const doubles & eta); extern "C" SEXP _capybara_solve_eta2_(SEXP yadj, SEXP myadj, SEXP offset, SEXP eta) { BEGIN_CPP11 - return cpp11::as_sexp(solve_eta2_(cpp11::as_cpp>(yadj), cpp11::as_cpp>(myadj), cpp11::as_cpp>(offset), cpp11::as_cpp>(eta))); - END_CPP11 -} -// 04_linear_algebra.cpp -doubles sqrt_(const SEXP & w); -extern "C" SEXP _capybara_sqrt_(SEXP w) { - BEGIN_CPP11 - return cpp11::as_sexp(sqrt_(cpp11::as_cpp>(w))); + return cpp11::as_sexp(solve_eta2_(cpp11::as_cpp>(yadj), cpp11::as_cpp &>>(myadj), cpp11::as_cpp>(offset), cpp11::as_cpp>(eta))); END_CPP11 } // 05_kendall_correlation.cpp @@ -173,7 +166,6 @@ static const R_CallMethodDef CallEntries[] = { {"_capybara_solve_eta2_", (DL_FUNC) &_capybara_solve_eta2_, 4}, {"_capybara_solve_eta_", (DL_FUNC) &_capybara_solve_eta_, 4}, {"_capybara_solve_y_", (DL_FUNC) &_capybara_solve_y_, 2}, - {"_capybara_sqrt_", (DL_FUNC) &_capybara_sqrt_, 1}, {"_capybara_update_beta_eta_", (DL_FUNC) &_capybara_update_beta_eta_, 3}, {"_capybara_update_nu_", (DL_FUNC) &_capybara_update_nu_, 3}, {NULL, NULL, 0} diff --git a/tests/testthat/test-apes.R b/tests/testthat/test-apes-bias.R similarity index 100% rename from tests/testthat/test-apes.R rename to tests/testthat/test-apes-bias.R diff --git a/tests/testthat/test-cpp.R b/tests/testthat/test-cpp.R deleted file mode 100644 index 2b53a70..0000000 --- a/tests/testthat/test-cpp.R +++ /dev/null @@ -1,28 +0,0 @@ -test_that("crossprod works", { - A <- matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2) - expect_equal(crossprod(A), crossprod_(A, NA_real_, FALSE, FALSE)) - - b <- c(1, 2) - expect_equal(crossprod(A * b), crossprod_(A, b, TRUE, FALSE)) - expect_equal(crossprod(A * sqrt(b)), crossprod_(A, b, TRUE, TRUE)) -}) - -test_that("solve_bias_ works", { - A <- matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2) - x <- c(2, 2) - expect_equal(as.vector(A %*% x), solve_y_(A, x)) - expect_equal(x - solve(A, x), solve_bias_(x, A, 1, x)) -}) - -test_that("inv_ works", { - A <- matrix(c(1, 0, 0, 1, 1, 0, 0, 1, 1), nrow = 3, ncol = 3, byrow = TRUE) - expect_equal(solve(A), inv_(A)) - - # non-invertible matrix - A <- matrix(c(1, 0, 0, 1, 0, 0, 0, 1, 1), nrow = 3, ncol = 3, byrow = TRUE) - expect_error(inv_(A)) - - # non-square matrix - A <- matrix(c(1, 0, 0, 1, 1, 0), nrow = 2, ncol = 3, byrow = TRUE) - expect_error(inv_(A)) -}) diff --git a/tests/testthat/test-fepoisson.R b/tests/testthat/test-fepoisson.R index 2bd9fb2..c4500e7 100644 --- a/tests/testthat/test-fepoisson.R +++ b/tests/testthat/test-fepoisson.R @@ -40,4 +40,10 @@ test_that("fepoisson is similar to fixest", { expect_output(summary_r2_(smod, 3)) expect_output(summary_nobs_(smod)) expect_output(summary_fisher_(smod)) + + fe <- fixed_effects(mod) + + expect_equal(length(fe), 2) + expect_equal(round(fe$exp_year[1:3], 3), c(10.195, 11.081, 11.260)) + expect_equal(round(fe$imp_year[1:3], 3), c(0.226, -0.254, 1.115)) }) diff --git a/tests/testthat/test-group-sums.R b/tests/testthat/test-group-sums.R new file mode 100644 index 0000000..86b1411 --- /dev/null +++ b/tests/testthat/test-group-sums.R @@ -0,0 +1,32 @@ +test_that("group_sums_* works", { + set.seed(123) + M <- matrix(rnorm(9), ncol = 3, nrow = 3) + w <- matrix(rnorm(3), ncol = 1, nrow = 3) + v <- matrix(rnorm(3), ncol = 1, nrow = 3) + K <- 2L + jlist <- list(1L, 2L) + + expect_equal( + round(group_sums_(M, w, jlist), 3), + # round(alpaca:::groupSums(M, w, jlist), 3) + matrix(c(4.144, 4.872, -2.942), ncol = 1, nrow = 3) + ) + + expect_equal( + round(group_sums_spectral_(M, w, v, K, jlist), 3), + # round(alpaca:::groupSumsSpectral(M, w, v, K, jlist), 3) + matrix(c(0, 0, 0), ncol = 1, nrow = 3) + ) + + expect_equal( + round(group_sums_var_(M, jlist)[, 1], 3), + # round(alpaca:::groupSumsVar(M, jlist)[, 1], 3) + c(2.483, 2.644, -0.779) + ) + + expect_equal( + round(group_sums_cov_(M, M, jlist)[, 1], 3), + # round(alpaca:::groupSumsCov(M, M, jlist)[, 1], 3) + c(0, 0, 0) + ) +}) diff --git a/tests/testthat/test-linear_algebra.R b/tests/testthat/test-linear_algebra.R index 3a4a036..8a2b453 100644 --- a/tests/testthat/test-linear_algebra.R +++ b/tests/testthat/test-linear_algebra.R @@ -1,3 +1,23 @@ +test_that("solve_bias_ works", { + A <- matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2) + x <- c(2, 2) + expect_equal(as.vector(A %*% x), solve_y_(A, x)) + expect_equal(x - solve(A, x), solve_bias_(x, A, 1, x)) +}) + +test_that("inv_ works", { + A <- matrix(c(1, 0, 0, 1, 1, 0, 0, 1, 1), nrow = 3, ncol = 3, byrow = TRUE) + expect_equal(solve(A), inv_(A)) + + # non-invertible matrix + A <- matrix(c(1, 0, 0, 1, 0, 0, 0, 1, 1), nrow = 3, ncol = 3, byrow = TRUE) + expect_error(inv_(A)) + + # non-square matrix + A <- matrix(c(1, 0, 0, 1, 1, 0), nrow = 2, ncol = 3, byrow = TRUE) + expect_error(inv_(A)) +}) + test_that("crossprod_ works", { set.seed(123) A <- matrix(rnorm(4), nrow = 2)