diff --git a/DESCRIPTION b/DESCRIPTION index f0f53f3..633f384 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: capybara Type: Package Title: Fast and Memory Efficient Fitting of Linear Models With High-Dimensional Fixed Effects -Version: 0.3.5 +Version: 0.4 Authors@R: c( person( given = "Mauricio", diff --git a/NEWS.md b/NEWS.md index 6475074..aa2eb56 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# capybara 0.4 + +* Uses R's C API efficiently to add a bit more of memory optimizations + # capybara 0.3.5 * Uses Mat consistently for all matrix operations (avoids vectors) diff --git a/R/cpp11.R b/R/cpp11.R index 4c95ebd..fc3c509 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -4,22 +4,10 @@ center_variables_ <- function(V_r, v_sum_r, w_r, klist, tol, maxiter, sum_v) { .Call(`_capybara_center_variables_`, V_r, v_sum_r, w_r, klist, tol, maxiter, sum_v) } -solve_beta_ <- function(mx, mnu, wtilde, epsilon, weighted) { - .Call(`_capybara_solve_beta_`, mx, mnu, wtilde, epsilon, weighted) -} - get_alpha_ <- function(p_r, klist, tol) { .Call(`_capybara_get_alpha_`, p_r, klist, tol) } -solve_eta_ <- function(mx, mnu, nu, beta) { - .Call(`_capybara_solve_eta_`, mx, mnu, nu, beta) -} - -solve_eta2_ <- function(yadj, myadj, offset, eta) { - .Call(`_capybara_solve_eta2_`, yadj, myadj, offset, eta) -} - group_sums_ <- function(M_r, w_r, jlist) { .Call(`_capybara_group_sums_`, M_r, w_r, jlist) } @@ -72,6 +60,30 @@ sandwich_ <- function(a, b) { .Call(`_capybara_sandwich_`, a, b) } +update_beta_eta_ <- function(old, upd, param) { + .Call(`_capybara_update_beta_eta_`, old, upd, param) +} + +update_nu_ <- function(y, mu, mu_eta) { + .Call(`_capybara_update_nu_`, y, mu, mu_eta) +} + +solve_beta_ <- function(mx, mnu, wtilde, epsilon, weighted) { + .Call(`_capybara_solve_beta_`, mx, mnu, wtilde, epsilon, weighted) +} + +solve_eta_ <- function(mx, mnu, nu, beta) { + .Call(`_capybara_solve_eta_`, mx, mnu, nu, beta) +} + +solve_eta2_ <- function(yadj, myadj, offset, eta) { + .Call(`_capybara_solve_eta2_`, yadj, myadj, offset, eta) +} + +sqrt_ <- function(w) { + .Call(`_capybara_sqrt_`, w) +} + pairwise_cor_ <- function(y, yhat) { .Call(`_capybara_pairwise_cor_`, y, yhat) } diff --git a/R/internals.R b/R/internals.R index b611300..8b3efa4 100644 --- a/R/internals.R +++ b/R/internals.R @@ -24,6 +24,7 @@ felm_fit_ <- function(y, X, wt, k.list, control) { MX <- center_variables_(MX, NA_real_, wt, k.list, center.tol, 10000L, FALSE) # Compute the OLS estimate + # beta <- as.vector(qr.solve(MX, y, epsilon)) beta <- solve_beta_(MX, y, NA_real_, epsilon, FALSE) # Generate result list @@ -70,7 +71,7 @@ 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) + w.tilde <- sqrt_(w) nu <- (y - mu) / mu.eta # Centering variables @@ -78,8 +79,9 @@ feglm_fit_ <- function(beta, eta, y, X, wt, k.list, family, control) { MX <- center_variables_(MX, NA_real_, w, k.list, center.tol, 10000L, FALSE) # 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, epsilon, TRUE) - eta.upd <- solve_eta_(MX, Mnu, nu, beta.upd) # Step-halving with three checks @@ -89,8 +91,10 @@ feglm_fit_ <- function(beta, eta, y, X, wt, k.list, family, control) { rho <- 1.0 for (inner.iter in seq.int(50L)) { - eta <- eta.old + rho * eta.upd - beta <- beta.old + rho * beta.upd + # eta <- eta.old + rho * eta.upd + # beta <- beta.old + rho * beta.upd + eta <- update_beta_eta_(eta.old, eta.upd, rho) + beta <- update_beta_eta_(beta.old, beta.upd, rho) mu <- family[["linkinv"]](eta) dev <- sum(family[["dev.resids"]](y, mu, wt)) dev.crit <- is.finite(dev) @@ -231,7 +235,8 @@ feglm_offset_ <- function(object, offset) { # 3. improvement as in glm2 rho <- 1.0 for (inner.iter in seq.int(50L)) { - eta <- eta.old + rho * eta.upd + # eta <- eta.old + rho * eta.upd + eta <- update_beta_eta_(eta.old, eta.upd, rho) mu <- family[["linkinv"]](eta) dev <- sum(family[["dev.resids"]](y, mu, wt)) dev.crit <- is.finite(dev) @@ -281,7 +286,8 @@ getScoreMatrix <- function(object) { mu <- family[["linkinv"]](eta) mu.eta <- family[["mu.eta"]](eta) w <- (wt * mu.eta^2) / family[["variance"]](mu) - nu <- (y - mu) / mu.eta + # nu <- (y - mu) / mu.eta + nu <- update_nu_(y, mu, mu.eta) # Center regressor matrix (if required) if (control[["keep.mx"]]) { diff --git a/README.Rmd b/README.Rmd index b51b4cb..d112205 100644 --- a/README.Rmd +++ b/README.Rmd @@ -124,7 +124,7 @@ Median time for the different models in the book |:------------|-------:|---------------:|-----------:|-----------------:|--------------------------:|-------------:| |Alpaca | 213.4ms| 2.3s| 1.35s| 1.86s| 2.59s| 4.96s| |Base R | 1.5m | 1.53m| 23.43m| 23.52m| 23.16m| 24.85m| -|**Capybara** | 263ms | 2.98s| 1.43s| 1.78s| 2.58s| 4.51s| +|**Capybara** | 323ms | 3.24s| 1.47s| 1.84s| 2.44s| 4.56s| |Fixest | 67.4ms| 477.08ms| 95.88ms| 136.21ms| 206.12ms| 415.31ms| Memory allocation for the same models @@ -133,7 +133,7 @@ Memory allocation for the same models |:------------|--------:|---------------:|-----------:|-----------------:|--------------------------:|-------------:| |Alpaca | 304.8MB | 339.8MB| 306.3MB| 335.61MB| 393.86MB| 539.49MB| |Base R | 2.73GB| 2.6GB| 11.9GB| 11.94GB| 11.95GB| 11.97GB| -|**Capybara** | 207MB | 231MB| 237MB| 244MB| 258MB| 293MB| +|**Capybara** | 204MB | 231MB| 237MB| 244MB| 258MB| 293MB| |Fixest | 44.59MB| 36.59MB| 28.1MB| 32.43MB| 41.12MB| 62.87MB| # Debugging diff --git a/README.md b/README.md index bdf9017..1778c4f 100644 --- a/README.md +++ b/README.md @@ -118,7 +118,7 @@ Analysis](https://www.wto.org/english/res_e/publications_e/advancedguide2016_e.h | :----------- | ------: | --------------: | ----------: | ----------------: | -------------------------: | ------------: | | Alpaca | 213.4ms | 2.3s | 1.35s | 1.86s | 2.59s | 4.96s | | Base R | 1.5m | 1.53m | 23.43m | 23.52m | 23.16m | 24.85m | -| **Capybara** | 263ms | 2.98s | 1.43s | 1.78s | 2.58s | 4.51s | +| **Capybara** | 323ms | 3.24s | 1.47s | 1.84s | 2.44s | 4.56s | | Fixest | 67.4ms | 477.08ms | 95.88ms | 136.21ms | 206.12ms | 415.31ms | Memory allocation for the same models @@ -127,7 +127,7 @@ Memory allocation for the same models | :----------- | ------: | --------------: | ----------: | ----------------: | -------------------------: | ------------: | | Alpaca | 304.8MB | 339.8MB | 306.3MB | 335.61MB | 393.86MB | 539.49MB | | Base R | 2.73GB | 2.6GB | 11.9GB | 11.94GB | 11.95GB | 11.97GB | -| **Capybara** | 207MB | 231MB | 237MB | 244MB | 258MB | 293MB | +| **Capybara** | 204MB | 231MB | 237MB | 244MB | 258MB | 293MB | | Fixest | 44.59MB | 36.59MB | 28.1MB | 32.43MB | 41.12MB | 62.87MB | # Debugging diff --git a/dev/benchmarks_tests_agtpa_capybara_only.R b/dev/benchmarks_tests_agtpa_capybara_only.R index d52e47b..236574f 100644 --- a/dev/benchmarks_tests_agtpa_capybara_only.R +++ b/dev/benchmarks_tests_agtpa_capybara_only.R @@ -157,7 +157,7 @@ bench_phasing <- readRDS("dev/bench_phasing.rds") bench_globalization <- readRDS("dev/bench_globalization.rds") -bench_ppml %>% +t1 <- bench_ppml %>% mutate(package = "**Capybara**") %>% mutate(model = "PPML") %>% select(model, package, median) %>% @@ -200,7 +200,7 @@ bench_ppml %>% arrange(package) %>% kable() -bench_ppml %>% +t2 <- bench_ppml %>% mutate(package = "**Capybara**") %>% mutate(model = "PPML") %>% select(model, package, mem_alloc) %>% @@ -242,3 +242,6 @@ bench_ppml %>% ) %>% arrange(package) %>% kable() + +t1 +t2 diff --git a/dev/check_bottlenecks.R b/dev/check_bottlenecks.R new file mode 100644 index 0000000..80b7c59 --- /dev/null +++ b/dev/check_bottlenecks.R @@ -0,0 +1,32 @@ +load_all() + +# d <- trade_panel +# d$trade_100 <- ifelse(d$trade > 100, 1L, 0L) + +# simulate a data frame with 1,000,000 rows with: +# trade_100: 0/1 +# lang: 0/1 +# clny: 0/1 +# rta: 0/1 +# year: 2000-2010 +set.seed(200100) +d <- data.frame( + trade_100 = sample(0:1, 1e6, replace = TRUE), + lang = sample(0:1, 1e6, replace = TRUE), + clny = sample(0:1, 1e6, replace = TRUE), + rta = sample(0:1, 1e6, replace = TRUE), + year = sample(2000:2010, 1e6, replace = TRUE) +) + +unique(d$trade_100) +unique(d$lang) +unique(d$clny) +unique(d$rta) +unique(d$year) + +# Fit 'feglm()' +load_all() +profvis::profvis(feglm(trade_100 ~ lang + clny + rta | year, d, family = binomial())) + +# Compute average partial effects +# bench::mark(apes(mod)) diff --git a/docs/index.html b/docs/index.html index dde47fa..fa71a2e 100644 --- a/docs/index.html +++ b/docs/index.html @@ -178,12 +178,12 @@

Benchmarks Capybara -263ms -2.98s -1.43s -1.78s -2.58s -4.51s +323ms +3.24s +1.47s +1.84s +2.44s +4.56s Fixest @@ -237,7 +237,7 @@

Benchmarks Capybara -207MB +204MB 231MB 237MB 244MB diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index b3ef834..17fa80f 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-03-03T21:12Z +last_built: 2024-03-04T01:22Z diff --git a/docs/reference/apes.html b/docs/reference/apes.html index 2d9df06..6f0e2f0 100644 --- a/docs/reference/apes.html +++ b/docs/reference/apes.html @@ -187,7 +187,7 @@

Examples

mod_bc <- bias_corr(mod) summary(mod_bc) #> Formula: trade ~ lang | year -#> <environment: 0x5669d9680358> +#> <environment: 0x6307ba4ace40> #> #> Family: Binomial #> diff --git a/docs/reference/bias_corr.html b/docs/reference/bias_corr.html index e1c5e3e..58049ae 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: 0x5669d8df0f50> +#> <environment: 0x6307b9d69e00> #> #> Family: Binomial #> diff --git a/docs/reference/feglm.html b/docs/reference/feglm.html index 8471b3f..330dba9 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: 0x5669d9d18e40> +#> <environment: 0x6307bac30440> #> #> Family: Poisson #> diff --git a/docs/reference/felm.html b/docs/reference/felm.html index 1cbb6b9..b1e3283 100644 --- a/docs/reference/felm.html +++ b/docs/reference/felm.html @@ -117,7 +117,7 @@

Examples

summary(mod) #> Formula: log(trade) ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5669da368040> +#> <environment: 0x6307bb219c28> #> #> Estimates: #> diff --git a/docs/reference/fenegbin.html b/docs/reference/fenegbin.html index 7888176..9fb3c9e 100644 --- a/docs/reference/fenegbin.html +++ b/docs/reference/fenegbin.html @@ -135,7 +135,7 @@

Examples

summary(mod) #> Formula: trade ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5669da988dd8> +#> <environment: 0x6307bba9e790> #> #> Family: Negative Binomial(1.1839) #> diff --git a/docs/reference/fepoisson.html b/docs/reference/fepoisson.html index 5568620..98edb37 100644 --- a/docs/reference/fepoisson.html +++ b/docs/reference/fepoisson.html @@ -122,7 +122,7 @@

Examples

summary(mod) #> Formula: trade ~ log_dist + lang + cntg + clny | exp_year + imp_year -#> <environment: 0x5669d99dfb98> +#> <environment: 0x6307bbf01968> #> #> Family: Poisson #> diff --git a/src/03_get_alpha.cpp b/src/02_get_alpha.cpp similarity index 100% rename from src/03_get_alpha.cpp rename to src/02_get_alpha.cpp diff --git a/src/02_solve_beta.cpp b/src/02_solve_beta.cpp deleted file mode 100644 index 254efe7..0000000 --- a/src/02_solve_beta.cpp +++ /dev/null @@ -1,38 +0,0 @@ -#include "00_main.h" - -[[cpp11::register]] doubles solve_beta_(const doubles_matrix<> &mx, - const doubles_matrix<> &mnu, - const doubles wtilde, double epsilon, - bool weighted) { - // Types conversion - Mat X = as_Mat(mx); - Mat Y = as_Mat(mnu); - - // Weight the X and Y matrices - if (weighted) { - // Additional type conversion - Mat W = as_Mat(wtilde); - - // Multiply each column of X by W pair-wise - X = X.each_col() % W; - - // Multiply each column of Y by W pair-wise - Y = Y.each_col() % W; - } - - // Now we need to solve the system X * beta = Y - // We proceed with the Economic QR - - Mat Q; - Mat R; - - bool computable = qr_econ(Q, R, X); - - if (!computable) { - stop("QR decomposition failed"); - } else { - // backsolve - Mat beta = solve(R, Q.t() * Y); - return as_doubles(beta); - } -} diff --git a/src/04_group_sums.cpp b/src/03_group_sums.cpp similarity index 100% rename from src/04_group_sums.cpp rename to src/03_group_sums.cpp diff --git a/src/03_solve_eta.cpp b/src/03_solve_eta.cpp deleted file mode 100644 index 4e7a290..0000000 --- a/src/03_solve_eta.cpp +++ /dev/null @@ -1,34 +0,0 @@ -// eta.upd <- nu - as.vector(Mnu - MX %*% beta.upd) - -#include "00_main.h" - -[[cpp11::register]] doubles solve_eta_(const doubles_matrix<> &mx, - const doubles_matrix<> &mnu, - const doubles &nu, const doubles &beta) { - // Types conversion - Mat MX = as_Mat(mx); - Mat Mnu = as_Mat(mnu); - Mat Nu = as_Mat(nu); - Mat Beta = as_Mat(beta); - - Mat res = Nu - (Mnu - (MX * Beta)); - - return as_doubles(res); -} - -// eta.upd <- yadj - as.vector(Myadj) + offset - eta - -[[cpp11::register]] doubles solve_eta2_(const doubles &yadj, - const doubles_matrix<> &myadj, - const doubles &offset, - const doubles &eta) { - // Types conversion - Mat Yadj = as_Mat(yadj); - Mat Myadj = as_Mat(myadj); - Mat Offset = as_Mat(offset); - Mat Eta = as_Mat(eta); - - Mat res = Yadj - Myadj + Offset - Eta; - - return as_doubles(res); -} diff --git a/src/04_linear_algebra.cpp b/src/04_linear_algebra.cpp new file mode 100644 index 0000000..14526e9 --- /dev/null +++ b/src/04_linear_algebra.cpp @@ -0,0 +1,259 @@ +#include "00_main.h" + +// Y <- crossprod(X) +// Y <- t(X) %*% X + +[[cpp11::register]] doubles_matrix<> crossprod_(const doubles_matrix<> &x, + const doubles &w, bool weighted, + bool root_weights) { + // Types conversion + Mat X = as_Mat(x); + + if (weighted) { + // Additional type conversion + Mat W = as_Mat(w); + + if (root_weights) { + W = sqrt(W); + } + + // Multiply each column of X by W pair-wise + X = X.each_col() % W; + } + + Mat res = X.t() * X; + return as_doubles_matrix(res); +} + +// WinvJ < -solve(object[["Hessian"]] / nt.full, J) +// Gamma < -(MX %*% WinvJ - PPsi) * v / nt.full +// V < -crossprod(Gamma) + +[[cpp11::register]] doubles_matrix<> gamma_(const doubles_matrix<> &mx, + const doubles_matrix<> &hessian, + const doubles_matrix<> j, + const doubles_matrix<> &ppsi, + const doubles &v, + const SEXP &nt_full) { + // Types conversion + Mat MX = as_Mat(mx); + Mat H = as_Mat(hessian); + Mat J = as_Mat(j); + Mat PPsi = as_Mat(ppsi); + Mat V = as_Mat(v); + + double N = as_cpp(nt_full); + + Mat res = (MX * solve(H / N, J)) - PPsi; + res = (res.each_col() % V) / N; + return as_doubles_matrix(res); +} + +// chol(crossprod(X)) + +[[cpp11::register]] doubles_matrix<> chol_crossprod_( + const doubles_matrix<> &x) { + // Types conversion + Mat X = as_Mat(x); + + // Crossprod + Cholesky decomposition + Mat res = chol(X.t() * X); + return as_doubles_matrix(res); +} + +// chol2inv(X) +// r comes from a Cholesky decomposition in the R code +// no need to check upper triangularity + +[[cpp11::register]] doubles_matrix<> chol2inv_(const doubles_matrix<> &r) { + // Types conversion + Mat R = as_Mat(r); + + // (X'X)^(-1) from the R part of the Cholesky decomposition + Mat res = inv(R.t() * R); + return as_doubles_matrix(res); +} + +// chol(X) + +[[cpp11::register]] doubles_matrix<> chol_(const doubles_matrix<> &x) { + // Types conversion + Mat X = as_Mat(x); + + // Cholesky decomposition + + Mat res = chol(X); + return as_doubles_matrix(res); +} + +// qr(X)$rank + +[[cpp11::register]] int qr_rank_(const doubles_matrix<> &x) { + Mat X = as_Mat(x); + + Mat Q; + Mat R; + + bool computable = qr_econ(Q, R, X); + + if (!computable) { + stop("QR decomposition failed"); + } else { + // rank = non-zero diagonal elements + int rank = sum(R.diag() != 0.0); + return rank; + } +} + +// Beta_uncorr - solve(H / nt, B) + +[[cpp11::register]] doubles solve_bias_(const doubles &beta_uncorr, + const doubles_matrix<> &hessian, + const double &nt, const doubles &b) { + // Types conversion + Mat Beta_uncorr = as_Mat(beta_uncorr); + Mat H = as_Mat(hessian); + Mat B = as_Mat(b); + + // Solve + return as_doubles(Beta_uncorr - solve(H / nt, B)); +} + +// A %*% x + +[[cpp11::register]] doubles solve_y_(const doubles_matrix<> &a, + const doubles &x) { + // Types conversion + Mat A = as_Mat(a); + Mat X = as_Mat(x); + + // Solve + return as_doubles(A * X); +} + +// A %*% B %*% A + +[[cpp11::register]] doubles_matrix<> sandwich_(const doubles_matrix<> &a, + const doubles_matrix<> &b) { + // Types conversion + Mat A = as_Mat(a); + Mat B = as_Mat(b); + + // Sandwich + + Mat res = A * B * A; + return as_doubles_matrix(res); +} + +// eta <- eta.old + rho * eta.upd + +[[cpp11::register]] doubles update_beta_eta_(const doubles &old, + const doubles &upd, + const double ¶m) { + int n = old.size(); + writable::doubles res(n); + + double *old_data = REAL(old); + double *upd_data = REAL(upd); + + for (int i = 0; i < n; i++) { + res[i] = old_data[i] + (upd_data[i] * param); + } + + return res; +} + +// nu <- (y - mu) / mu.eta + +[[cpp11::register]] doubles update_nu_(const SEXP &y, const SEXP &mu, + const SEXP &mu_eta) { + int n = Rf_length(y); + writable::doubles res(n); + + double *y_data = REAL(y); + double *mu_data = REAL(mu); + double *mu_eta_data = REAL(mu_eta); + + for (int i = 0; i < n; i++) { + res[i] = (y_data[i] - mu_data[i]) / mu_eta_data[i]; + } + + return res; +} + +[[cpp11::register]] doubles solve_beta_(const doubles_matrix<> &mx, + const doubles_matrix<> &mnu, + const doubles wtilde, double epsilon, + bool weighted) { + // Types conversion + Mat X = as_Mat(mx); + Mat Y = as_Mat(mnu); + + // Weight the X and Y matrices + if (weighted) { + // Additional type conversion + Mat W = as_Mat(wtilde); + + // Multiply each column of X by W pair-wise + X = X.each_col() % W; + + // Multiply each column of Y by W pair-wise + Y = Y.each_col() % W; + } + + // Now we need to solve the system X * beta = Y + // We proceed with the Economic QR + + Mat Q; + Mat R; + + bool computable = qr_econ(Q, R, X); + + if (!computable) { + stop("QR decomposition failed"); + } else { + // backsolve + return as_doubles(solve(R, Q.t() * Y)); + } +} + +// eta.upd <- nu - as.vector(Mnu - MX %*% beta.upd) + +[[cpp11::register]] doubles solve_eta_(const doubles_matrix<> &mx, + const doubles_matrix<> &mnu, + const doubles &nu, const doubles &beta) { + // Types conversion + Mat MX = as_Mat(mx); + Mat Mnu = as_Mat(mnu); + Mat Nu = as_Mat(nu); + Mat Beta = as_Mat(beta); + + return as_doubles(Nu - (Mnu - (MX * Beta))); +} + +// eta.upd <- yadj - as.vector(Myadj) + offset - eta + +[[cpp11::register]] doubles solve_eta2_(const doubles &yadj, + const doubles_matrix<> &myadj, + const doubles &offset, + const doubles &eta) { + // Types conversion + int n = yadj.size(); + writable::doubles res(n); + for (int i = 0; i < n; i++) { + res[i] = yadj[i] - myadj(i, 0) + offset[i] - eta[i]; + } + return res; +} + +// w <- sqrt(w) + +[[cpp11::register]] doubles sqrt_(const doubles &w) { + // Types conversion + int n = w.size(); + writable::doubles res(n); + for (int i = 0; i < n; i++) { + res[i] = sqrt(w[i]); + } + return res; +} diff --git a/src/05_linear_algebra.cpp b/src/05_linear_algebra.cpp deleted file mode 100644 index acb634f..0000000 --- a/src/05_linear_algebra.cpp +++ /dev/null @@ -1,144 +0,0 @@ -// Y <- crossprod(X) -// Y <- t(X) %*% X - -#include "00_main.h" - -[[cpp11::register]] doubles_matrix<> crossprod_(const doubles_matrix<> &x, - const doubles &w, bool weighted, - bool root_weights) { - // Types conversion - Mat X = as_Mat(x); - - if (weighted) { - // Additional type conversion - Mat W = as_Mat(w); - - if (root_weights) { - W = sqrt(W); - } - - // Multiply each column of X by W pair-wise - X = X.each_col() % W; - } - - Mat Y = X.t() * X; - - return as_doubles_matrix(Y); -} - -// WinvJ < -solve(object[["Hessian"]] / nt.full, J) -// Gamma < -(MX % * % WinvJ - PPsi) * v / nt.full -// V < -crossprod(Gamma) - -[[cpp11::register]] doubles_matrix<> -gamma_(const doubles_matrix<> &mx, const doubles_matrix<> &hessian, - const doubles_matrix<> j, const doubles_matrix<> &ppsi, const doubles &v, - const int &nt_full) { - // Types conversion - Mat MX = as_Mat(mx); - Mat H = as_Mat(hessian); - Mat J = as_Mat(j); - Mat PPsi = as_Mat(ppsi); - Mat V = as_Mat(v); - - double N = static_cast(nt_full); - Mat res = (MX * solve(H / N, J) - PPsi); - res = (res.each_col() % V) / N; - - return as_doubles_matrix(res); -} - -[[cpp11::register]] doubles_matrix<> -chol_crossprod_(const doubles_matrix<> &x) { - // Types conversion - Mat X = as_Mat(x); - - // Crossprod - Mat Y = X.t() * X; - - // Cholesky decomposition - Mat res = chol(Y); - - return as_doubles_matrix(res); -} - -// chol2inv(X) - -// r comes from a Cholesky decomposition in the R code -// no need to check upper triangularity -[[cpp11::register]] doubles_matrix<> chol2inv_(const doubles_matrix<> &r) { - // Types conversion - Mat R = as_Mat(r); - - // (X'X)^(-1) from the R part of the Cholesky decomposition - Mat res = inv(R.t() * R); - - return as_doubles_matrix(res); -} - -// chol(X) - -[[cpp11::register]] doubles_matrix<> chol_(const doubles_matrix<> &x) { - // Types conversion - Mat X = as_Mat(x); - - // Cholesky decomposition - Mat res = chol(X); - - return as_doubles_matrix(res); -} - -[[cpp11::register]] int qr_rank_(const doubles_matrix<> &x) { - Mat X = as_Mat(x); - - Mat Q; - Mat R; - - bool computable = qr_econ(Q, R, X); - - if (!computable) { - stop("QR decomposition failed"); - } else { - // rank = non-zero diagonal elements - int rank = sum(R.diag() != 0.0); - return rank; - } -} - -[[cpp11::register]] doubles solve_bias_(const doubles &beta_uncorr, - const doubles_matrix<> &hessian, - const double &nt, const doubles &b) { - // Types conversion - Mat Beta_uncorr = as_Mat(beta_uncorr); - Mat H = as_Mat(hessian); - Mat B = as_Mat(b); - - // Solve - Mat res = Beta_uncorr - solve(H / nt, B); - - return as_doubles(res); -} - -[[cpp11::register]] doubles solve_y_(const doubles_matrix<> &a, - const doubles &x) { - // Types conversion - Mat A = as_Mat(a); - Mat X = as_Mat(x); - - // Solve - Mat res = A * X; - - return as_doubles(res); -} - -[[cpp11::register]] doubles_matrix<> sandwich_(const doubles_matrix<> &a, - const doubles_matrix<> &b) { - // Types conversion - Mat A = as_Mat(a); - Mat B = as_Mat(b); - - // Sandwich - Mat res = (A * B) * A; - - return as_doubles_matrix(res); -} diff --git a/src/06_pairwise_correlation.cpp b/src/05_pairwise_correlation.cpp similarity index 100% rename from src/06_pairwise_correlation.cpp rename to src/05_pairwise_correlation.cpp diff --git a/src/cpp11.cpp b/src/cpp11.cpp index e4de51b..1d9d68c 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -12,126 +12,147 @@ extern "C" SEXP _capybara_center_variables_(SEXP V_r, SEXP v_sum_r, SEXP w_r, SE return cpp11::as_sexp(center_variables_(cpp11::as_cpp &>>(V_r), cpp11::as_cpp>(v_sum_r), cpp11::as_cpp>(w_r), cpp11::as_cpp>(klist), cpp11::as_cpp>(tol), cpp11::as_cpp>(maxiter), cpp11::as_cpp>(sum_v))); END_CPP11 } -// 02_solve_beta.cpp -doubles solve_beta_(const doubles_matrix<> & mx, const doubles_matrix<> & mnu, const doubles wtilde, double epsilon, bool weighted); -extern "C" SEXP _capybara_solve_beta_(SEXP mx, SEXP mnu, SEXP wtilde, SEXP epsilon, SEXP weighted) { - BEGIN_CPP11 - return cpp11::as_sexp(solve_beta_(cpp11::as_cpp &>>(mx), cpp11::as_cpp &>>(mnu), cpp11::as_cpp>(wtilde), cpp11::as_cpp>(epsilon), cpp11::as_cpp>(weighted))); - END_CPP11 -} -// 03_get_alpha.cpp +// 02_get_alpha.cpp list get_alpha_(const doubles_matrix<> & p_r, const list & klist, const double tol); extern "C" SEXP _capybara_get_alpha_(SEXP p_r, SEXP klist, SEXP tol) { BEGIN_CPP11 return cpp11::as_sexp(get_alpha_(cpp11::as_cpp &>>(p_r), cpp11::as_cpp>(klist), cpp11::as_cpp>(tol))); END_CPP11 } -// 03_solve_eta.cpp -doubles solve_eta_(const doubles_matrix<> & mx, const doubles_matrix<> & mnu, const doubles & nu, const doubles & beta); -extern "C" SEXP _capybara_solve_eta_(SEXP mx, SEXP mnu, SEXP nu, SEXP beta) { - BEGIN_CPP11 - return cpp11::as_sexp(solve_eta_(cpp11::as_cpp &>>(mx), cpp11::as_cpp &>>(mnu), cpp11::as_cpp>(nu), cpp11::as_cpp>(beta))); - END_CPP11 -} -// 03_solve_eta.cpp -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_group_sums.cpp +// 03_group_sums.cpp doubles group_sums_(const doubles_matrix<> & M_r, const doubles_matrix<> & w_r, const list & jlist); extern "C" SEXP _capybara_group_sums_(SEXP M_r, SEXP w_r, SEXP jlist) { BEGIN_CPP11 return cpp11::as_sexp(group_sums_(cpp11::as_cpp &>>(M_r), cpp11::as_cpp &>>(w_r), cpp11::as_cpp>(jlist))); END_CPP11 } -// 04_group_sums.cpp +// 03_group_sums.cpp doubles group_sums_spectral_(const doubles_matrix<> & M_r, const doubles_matrix<> & v_r, const doubles_matrix<> & w_r, const int K, const list & jlist); extern "C" SEXP _capybara_group_sums_spectral_(SEXP M_r, SEXP v_r, SEXP w_r, SEXP K, SEXP jlist) { BEGIN_CPP11 return cpp11::as_sexp(group_sums_spectral_(cpp11::as_cpp &>>(M_r), cpp11::as_cpp &>>(v_r), cpp11::as_cpp &>>(w_r), cpp11::as_cpp>(K), cpp11::as_cpp>(jlist))); END_CPP11 } -// 04_group_sums.cpp +// 03_group_sums.cpp doubles_matrix<> group_sums_var_(const doubles_matrix<> & M_r, const list & jlist); extern "C" SEXP _capybara_group_sums_var_(SEXP M_r, SEXP jlist) { BEGIN_CPP11 return cpp11::as_sexp(group_sums_var_(cpp11::as_cpp &>>(M_r), cpp11::as_cpp>(jlist))); END_CPP11 } -// 04_group_sums.cpp +// 03_group_sums.cpp doubles_matrix<> group_sums_cov_(const doubles_matrix<> & M_r, const doubles_matrix<> & N_r, const list & jlist); extern "C" SEXP _capybara_group_sums_cov_(SEXP M_r, SEXP N_r, SEXP jlist) { BEGIN_CPP11 return cpp11::as_sexp(group_sums_cov_(cpp11::as_cpp &>>(M_r), cpp11::as_cpp &>>(N_r), cpp11::as_cpp>(jlist))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles_matrix<> crossprod_(const doubles_matrix<> & x, const doubles & w, bool weighted, bool root_weights); extern "C" SEXP _capybara_crossprod_(SEXP x, SEXP w, SEXP weighted, SEXP root_weights) { BEGIN_CPP11 return cpp11::as_sexp(crossprod_(cpp11::as_cpp &>>(x), cpp11::as_cpp>(w), cpp11::as_cpp>(weighted), cpp11::as_cpp>(root_weights))); END_CPP11 } -// 05_linear_algebra.cpp -doubles_matrix<> gamma_(const doubles_matrix<> & mx, const doubles_matrix<> & hessian, const doubles_matrix<> j, const doubles_matrix<> & ppsi, const doubles & v, const int & nt_full); +// 04_linear_algebra.cpp +doubles_matrix<> gamma_(const doubles_matrix<> & mx, const doubles_matrix<> & hessian, const doubles_matrix<> j, const doubles_matrix<> & ppsi, const doubles & v, const SEXP & nt_full); extern "C" SEXP _capybara_gamma_(SEXP mx, SEXP hessian, SEXP j, SEXP ppsi, SEXP v, SEXP nt_full) { BEGIN_CPP11 - return cpp11::as_sexp(gamma_(cpp11::as_cpp &>>(mx), cpp11::as_cpp &>>(hessian), cpp11::as_cpp>>(j), cpp11::as_cpp &>>(ppsi), cpp11::as_cpp>(v), cpp11::as_cpp>(nt_full))); + return cpp11::as_sexp(gamma_(cpp11::as_cpp &>>(mx), cpp11::as_cpp &>>(hessian), cpp11::as_cpp>>(j), cpp11::as_cpp &>>(ppsi), cpp11::as_cpp>(v), cpp11::as_cpp>(nt_full))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles_matrix<> chol_crossprod_(const doubles_matrix<> & x); extern "C" SEXP _capybara_chol_crossprod_(SEXP x) { BEGIN_CPP11 return cpp11::as_sexp(chol_crossprod_(cpp11::as_cpp &>>(x))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles_matrix<> chol2inv_(const doubles_matrix<> & r); extern "C" SEXP _capybara_chol2inv_(SEXP r) { BEGIN_CPP11 return cpp11::as_sexp(chol2inv_(cpp11::as_cpp &>>(r))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles_matrix<> chol_(const doubles_matrix<> & x); extern "C" SEXP _capybara_chol_(SEXP x) { BEGIN_CPP11 return cpp11::as_sexp(chol_(cpp11::as_cpp &>>(x))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp int qr_rank_(const doubles_matrix<> & x); extern "C" SEXP _capybara_qr_rank_(SEXP x) { BEGIN_CPP11 return cpp11::as_sexp(qr_rank_(cpp11::as_cpp &>>(x))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles solve_bias_(const doubles & beta_uncorr, const doubles_matrix<> & hessian, const double & nt, const doubles & b); extern "C" SEXP _capybara_solve_bias_(SEXP beta_uncorr, SEXP hessian, SEXP nt, SEXP b) { BEGIN_CPP11 return cpp11::as_sexp(solve_bias_(cpp11::as_cpp>(beta_uncorr), cpp11::as_cpp &>>(hessian), cpp11::as_cpp>(nt), cpp11::as_cpp>(b))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles solve_y_(const doubles_matrix<> & a, const doubles & x); extern "C" SEXP _capybara_solve_y_(SEXP a, SEXP x) { BEGIN_CPP11 return cpp11::as_sexp(solve_y_(cpp11::as_cpp &>>(a), cpp11::as_cpp>(x))); END_CPP11 } -// 05_linear_algebra.cpp +// 04_linear_algebra.cpp doubles_matrix<> sandwich_(const doubles_matrix<> & a, const doubles_matrix<> & b); extern "C" SEXP _capybara_sandwich_(SEXP a, SEXP b) { BEGIN_CPP11 return cpp11::as_sexp(sandwich_(cpp11::as_cpp &>>(a), cpp11::as_cpp &>>(b))); END_CPP11 } -// 06_pairwise_correlation.cpp +// 04_linear_algebra.cpp +doubles update_beta_eta_(const doubles & old, const doubles & upd, const double & param); +extern "C" SEXP _capybara_update_beta_eta_(SEXP old, SEXP upd, SEXP param) { + BEGIN_CPP11 + return cpp11::as_sexp(update_beta_eta_(cpp11::as_cpp>(old), cpp11::as_cpp>(upd), cpp11::as_cpp>(param))); + END_CPP11 +} +// 04_linear_algebra.cpp +doubles update_nu_(const SEXP & y, const SEXP & mu, const SEXP & mu_eta); +extern "C" SEXP _capybara_update_nu_(SEXP y, SEXP mu, SEXP mu_eta) { + BEGIN_CPP11 + return cpp11::as_sexp(update_nu_(cpp11::as_cpp>(y), cpp11::as_cpp>(mu), cpp11::as_cpp>(mu_eta))); + END_CPP11 +} +// 04_linear_algebra.cpp +doubles solve_beta_(const doubles_matrix<> & mx, const doubles_matrix<> & mnu, const doubles wtilde, double epsilon, bool weighted); +extern "C" SEXP _capybara_solve_beta_(SEXP mx, SEXP mnu, SEXP wtilde, SEXP epsilon, SEXP weighted) { + BEGIN_CPP11 + return cpp11::as_sexp(solve_beta_(cpp11::as_cpp &>>(mx), cpp11::as_cpp &>>(mnu), cpp11::as_cpp>(wtilde), cpp11::as_cpp>(epsilon), cpp11::as_cpp>(weighted))); + END_CPP11 +} +// 04_linear_algebra.cpp +doubles solve_eta_(const doubles_matrix<> & mx, const doubles_matrix<> & mnu, const doubles & nu, const doubles & beta); +extern "C" SEXP _capybara_solve_eta_(SEXP mx, SEXP mnu, SEXP nu, SEXP beta) { + BEGIN_CPP11 + return cpp11::as_sexp(solve_eta_(cpp11::as_cpp &>>(mx), cpp11::as_cpp &>>(mnu), cpp11::as_cpp>(nu), cpp11::as_cpp>(beta))); + END_CPP11 +} +// 04_linear_algebra.cpp +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 doubles & w); +extern "C" SEXP _capybara_sqrt_(SEXP w) { + BEGIN_CPP11 + return cpp11::as_sexp(sqrt_(cpp11::as_cpp>(w))); + END_CPP11 +} +// 05_pairwise_correlation.cpp double pairwise_cor_(const doubles & y, const doubles & yhat); extern "C" SEXP _capybara_pairwise_cor_(SEXP y, SEXP yhat) { BEGIN_CPP11 @@ -160,6 +181,9 @@ 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} }; }