diff --git a/.Rbuildignore b/.Rbuildignore index 8fa460c..0d639f6 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,2 +1,7 @@ ^wcec\.Rproj$ ^\.Rproj\.user$ +^_pkgdown\.yml$ +^docs$ +^pkgdown$ +^\.github$ +^sandbox$ diff --git a/.github/.gitignore b/.github/.gitignore new file mode 100644 index 0000000..2d19fc7 --- /dev/null +++ b/.github/.gitignore @@ -0,0 +1 @@ +*.html diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml new file mode 100644 index 0000000..d46a617 --- /dev/null +++ b/.github/workflows/R-CMD-check.yaml @@ -0,0 +1,52 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + +name: R-CMD-check.yaml + +permissions: read-all + +jobs: + R-CMD-check: + runs-on: ${{ matrix.config.os }} + + name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + + strategy: + fail-fast: false + matrix: + config: + - {os: macos-latest, r: 'release'} + - {os: windows-latest, r: 'release'} + - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} + - {os: ubuntu-latest, r: 'release'} + - {os: ubuntu-latest, r: 'oldrel-1'} + + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + R_KEEP_PKG_SOURCE: yes + + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + r-version: ${{ matrix.config.r }} + http-user-agent: ${{ matrix.config.http-user-agent }} + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck + needs: check + + - uses: r-lib/actions/check-r-package@v2 + with: + upload-snapshots: true + build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")' diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 0000000..4bbce75 --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,50 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + release: + types: [published] + workflow_dispatch: + +name: pkgdown.yaml + +permissions: read-all + +jobs: + pkgdown: + runs-on: ubuntu-latest + # Only restrict concurrency for non-PR jobs + concurrency: + group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }} + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + permissions: + contents: write + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::pkgdown, local::. + needs: website + + - name: Build site + run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) + shell: Rscript {0} + + - name: Deploy to GitHub pages 🚀 + if: github.event_name != 'pull_request' + uses: JamesIves/github-pages-deploy-action@v4.5.0 + with: + clean: false + branch: gh-pages + folder: docs diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml new file mode 100644 index 0000000..9882260 --- /dev/null +++ b/.github/workflows/test-coverage.yaml @@ -0,0 +1,61 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + +name: test-coverage.yaml + +permissions: read-all + +jobs: + test-coverage: + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::covr, any::xml2 + needs: coverage + + - name: Test coverage + run: | + cov <- covr::package_coverage( + quiet = FALSE, + clean = FALSE, + install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") + ) + covr::to_cobertura(cov) + shell: Rscript {0} + + - uses: codecov/codecov-action@v4 + with: + fail_ci_if_error: ${{ github.event_name != 'pull_request' && true || false }} + file: ./cobertura.xml + plugin: noop + disable_search: true + token: ${{ secrets.CODECOV_TOKEN }} + + - name: Show testthat output + if: always() + run: | + ## -------------------------------------------------------------------- + find '${{ runner.temp }}/package' -name 'testthat.Rout*' -exec cat '{}' \; || true + shell: bash + + - name: Upload test results + if: failure() + uses: actions/upload-artifact@v4 + with: + name: coverage-test-failures + path: ${{ runner.temp }}/package diff --git a/.gitignore b/.gitignore index accfb7b..5bbbb81 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ .httr-oauth .DS_Store .quarto +docs +sandbox diff --git a/DESCRIPTION b/DESCRIPTION index 5e1cb57..117164f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,5 +10,13 @@ License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a Encoding: UTF-8 Roxygen: list(markdown = TRUE) RoxygenNote: 7.3.2 -URL: https://github.com/swarm-lab/wcec +URL: https://github.com/swarm-lab/wcec, https://swarm-lab.github.io/wcec/ BugReports: https://github.com/swarm-lab/wcec/issues +LazyData: true +LinkingTo: + Rcpp, + RcppEigen +Imports: + Rcpp, + RcppEigen, + Rfast diff --git a/NAMESPACE b/NAMESPACE index 6ae9268..42825af 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,2 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(wcec) +import(RcppEigen) +importFrom(Rcpp,evalCpp) +useDynLib(wcec, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..2f49a37 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,33 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#' @name .wcov +#' +#' @title Weighted Covariance Matrix +#' +#' @description \code{wcov} computes the estimates of the weighted covariance +#' matrix and the weighted mean of the data. +#' +#' @param x A matrix with \eqn{m} columns and \eqn{n} rows, where each column +#' represents a different variable and each row a different observation. +#' +#' @param w A non-negative and non-zero vector of weights for each observation. +#' Its length must equal the number of rows of x. +#' +#' @return A list with two components: +#' \itemize{ +#' \item \code{center}: an estimate for the center (mean) of the data. +#' \item \code{cov}: the estimated (weighted) covariance matrix. +#' } +#' +#' @author Simon Garnier, \email{garnier@@njit.edu} +#' +#' @examples +#' m <- matrix(c(rnorm(500, 6), rnorm(500, 11, 3)), ncol = 2) +#' w <- runif(500) +#' gravitree:::.wcov(m, w) +#' +.wcov <- function(x, w) { + .Call('_wcec_wcov', PACKAGE = 'wcec', x, w) +} + diff --git a/R/wcec-package.R b/R/wcec-package.R new file mode 100644 index 0000000..a06bcaa --- /dev/null +++ b/R/wcec-package.R @@ -0,0 +1,57 @@ +#' wcec: ... +#' +#' The \code{wcec} package provides ... +#' +#' @author Simon Garnier +#' @author Jason Graham +#' +#' Maintainer: Simon Garnier +#' +#' @details +#' \tabular{ll}{ +#' Package: \tab wcec\cr +#' Type: \tab Package\cr +#' Version: \tab 0.1\cr +#' Date: \tab 2024-08-07\cr +#' License: \tab GPL-3\cr +#' } +#' +#' @name wcec +#' @import RcppEigen +#' @importFrom Rcpp evalCpp +#' @useDynLib wcec, .registration = TRUE +"_PACKAGE" + + +#' @title Four Gaussian Clusters +#' +#' @description Matrix of 2-dimensional points forming four Gaussian clusters. +#' +#' @name four_Gaussians +#' +#' @docType data +#' +#' @keywords datasets +#' +#' @examples +#' data(four_Gaussians) +#' plot(four_Gaussians, pch = 19) +#' +NULL + + +#' @title Three Gaussian Clusters +#' +#' @description Matrix of 2-dimensional points forming three Gaussian clusters. +#' +#' @name three_Gaussians +#' +#' @docType data +#' +#' @keywords datasets +#' +#' @examples +#' data(three_Gaussians) +#' plot(three_Gaussians, pch = 19) +#' +NULL diff --git a/R/wcec.R b/R/wcec.R new file mode 100644 index 0000000..d97d38b --- /dev/null +++ b/R/wcec.R @@ -0,0 +1,171 @@ +.init_density <- function(x, k, w) { + n_blocks <- k + 2 + block_list <- apply(x, 2, function(y) { + breaks <- seq( + from = min(y, na.rm = TRUE), + to = max(y, na.rm = TRUE), + length.out = n_blocks + ) + cut(y, breaks, labels = FALSE, include.lowest = TRUE, right = FALSE) + }, simplify = FALSE) + fq <- table(block_list) + ord <- order(fq, decreasing = TRUE) + idx <- arrayInd(ord, dim(fq), dimnames(fq)) + init <- matrix(NA, nrow = k, ncol = ncol(x)) + blocks <- simplify2array(block_list) + + for (i in seq_len(k)) { + ix <- Rfast::rowAll(sweep(blocks, 2, idx[1, ], "==")) + init[i, ] <- apply(x[ix, , drop = FALSE], 2, weighted.mean, w = w[ix]) + d <- Rfast::rowMaxs(sweep(idx, 2, idx[1, ])^2, TRUE) + idx <- idx[d > sqrt(n_blocks), , drop = FALSE] + } + + Rfast::rowMins( + Rfast::dista(x, init) + ) +} + +.array_split <- function(x, k) { + nr <- nrow(x) + n <- rep(floor(nr / k), k) + r <- nr - sum(n) + + for (i in seq_len(r)) { + n[i] <- n[i] + 1 + } + + end <- cumsum(n) + start <- c(1, end[1:(k - 1)] + 1) + + l <- lapply(1:k, function(i) { + x[start[i]:end[i], , drop = FALSE] + }) + names(l) <- paste0(start, "-", end) + l +} + +.init_sharding <- function(x, k) { + ord <- order(Rfast::rowsums(x)) + xsplit <- .array_split(x[ord, ], k) + Rfast::rowMins( + Rfast::dista(x, t(sapply(xsplit, Rfast::colmeans))) + ) +} + +.init_clusters <- function(x, k, w, method) { + n <- nrow(x) + + if (n != length(w)) { + stop("The number of elements in `w` is not the same as the number of rows in `x`.") + } + + if (length(k) > 1) { + if (n != length(k)) { + stop("The number of elements in `k` is not the same as the number of rows in `x`.") + } + + if (!all((floor(k) - k) == 0)) { + stop("Not all values in `k` are integers.") + } + + suk <- sort(unique(k)) + lookup <- cbind(suk, as.numeric(as.factor(suk))) + match(k, lookup) + } else if (length(k) == 1) { + if (k > 1) { + if (method == "sharding") { + .init_sharding(x, k) + } else if (method == "density") { + .init_density(x, k, w) + } else if (method == "random") { + sample(1:k, n, TRUE) + } else { + stop("Invalid initialization method.") + } + } else { + rep(1, n) + } + } else { + stop("Invalid 'k'.") + } +} + +.ce_gaussian <- function(n, cov) { + (n * log(2.0 * pi * exp(1)) + log(det(cov))) / 2 +} + +.cost_gaussian <- function(p, n, params) { + sum( + sapply(1:length(p), function(i) { + .ce_gaussian(p[i] * n, params[[i]]$cov) + }) + ) +} + +.dist_gaussian <- function(x, p, mu, sigma) { + -log(p) - log(mvtnorm::dmvnorm(x, mean = mu, sigma = sigma, checkSymmetry = FALSE)) +} + +.wcec <- function(x, k, w, iter_max) { + n <- nrow(x) + out <- list( + clustering = k, + probability = table(k) / n, + params = lapply(unique(k), function(clust) { + idx <- k == clust + wcec:::.wcov( + x[idx, , drop = FALSE], + w[idx] + ) + }), + iter = 0 + ) + d <- matrix(Inf, nrow = n, ncol = length(out$params)) + out$cost <- wcec:::.cost_gaussian(out$probability, n, out$params) + + for (i in 1:iter_max) { + for (j in seq_along(out$params)) { + d[, j] <- wcec:::.dist_gaussian(x, out$probability[j], out$params[[j]]$center, out$params[[j]]$cov) + } + + new_clustering <- Rfast::rowMins(d) + new_clustering <- match(new_clustering, sort(unique(new_clustering))) + test <- all(new_clustering == out$clustering) + out$clustering <- new_clustering + + if (test == TRUE) { + break + } else { + out$params <- lapply(unique(out$clustering), function(clust) { + idx <- out$clustering == clust + wcec:::.wcov( + x[idx, , drop = FALSE], + w[idx] + ) + }) + + d <- d * Inf + out$probability <- table(out$clustering) / n + out$cost <- wcec:::.cost_gaussian(out$probability, n, out$params) + } + } + + out$iter <- i + out +} + +#' @export +wcec <- function( + x, k, w = rep(1 / nrow(x), nrow(x)), iter_max = 10, + init_method = "sharding") { + # CHECKS + + # INITIALIZATION + k <- .init_clusters(x, k, w, method = init_method) + + # FIRST PASS + out <- .wcec(x = x, k = k, w = w, iter_max = iter_max) + + out +} diff --git a/README.md b/README.md index 779d4d7..7061fbe 100644 --- a/README.md +++ b/README.md @@ -1 +1,30 @@ -# wcec \ No newline at end of file + +# wcec + + +[![CRAN status](https://www.r-pkg.org/badges/version/wcec)](https://CRAN.R-project.org/package=wcec) +[![R-CMD-check](https://github.com/swarm-lab/wcec/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/swarm-lab/wcec/actions/workflows/R-CMD-check.yaml) +[![Codecov test coverage](https://codecov.io/gh/swarm-lab/wcec/graph/badge.svg)](https://app.codecov.io/gh/swarm-lab/wcec) +[![test-coverage](https://github.com/swarm-lab/wcec/actions/workflows/test-coverage.yaml/badge.svg)](https://github.com/swarm-lab/wcec/actions/workflows/test-coverage.yaml) + + +The goal of wcec is to ... + +## Installation + +You can install the development version of wcec from [GitHub](https://github.com/) with: + +``` r +# install.packages("pak") +pak::pak("swarm-lab/wcec") +``` + +## Example + +This is a basic example which shows you how to solve a common problem: + +``` r +library(wcec) +## basic example code +``` + diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..60eea50 --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,4 @@ +url: https://swarm-lab.github.io/wcec/ +template: + bootstrap: 5 + diff --git a/data/four_Gaussians.rda b/data/four_Gaussians.rda new file mode 100644 index 0000000..2a95fb6 Binary files /dev/null and b/data/four_Gaussians.rda differ diff --git a/data/three_Gaussians.rda b/data/three_Gaussians.rda new file mode 100644 index 0000000..17b6b5a Binary files /dev/null and b/data/three_Gaussians.rda differ diff --git a/man/dot-wcov.Rd b/man/dot-wcov.Rd new file mode 100644 index 0000000..00eb3c3 --- /dev/null +++ b/man/dot-wcov.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{.wcov} +\alias{.wcov} +\title{Weighted Covariance Matrix} +\usage{ +.wcov(x, w) +} +\arguments{ +\item{x}{A matrix with \eqn{m} columns and \eqn{n} rows, where each column +represents a different variable and each row a different observation.} + +\item{w}{A non-negative and non-zero vector of weights for each observation. +Its length must equal the number of rows of x.} +} +\value{ +A list with two components: +\itemize{ +\item \code{center}: an estimate for the center (mean) of the data. +\item \code{cov}: the estimated (weighted) covariance matrix. +} +} +\description{ +\code{wcov} computes the estimates of the weighted covariance +matrix and the weighted mean of the data. +} +\examples{ +m <- matrix(c(rnorm(500, 6), rnorm(500, 11, 3)), ncol = 2) +w <- runif(500) +gravitree:::.wcov(m, w) + +} +\author{ +Simon Garnier, \email{garnier@njit.edu} +} diff --git a/man/four_Gaussians.Rd b/man/four_Gaussians.Rd new file mode 100644 index 0000000..dc45bf9 --- /dev/null +++ b/man/four_Gaussians.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wcec-package.R +\docType{data} +\name{four_Gaussians} +\alias{four_Gaussians} +\title{Four Gaussian Clusters} +\description{ +Matrix of 2-dimensional points forming four Gaussian clusters. +} +\examples{ +data(four_Gaussians) +plot(four_Gaussians, pch = 19) + +} +\keyword{datasets} diff --git a/man/three_Gaussians.Rd b/man/three_Gaussians.Rd new file mode 100644 index 0000000..ea252fe --- /dev/null +++ b/man/three_Gaussians.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wcec-package.R +\docType{data} +\name{three_Gaussians} +\alias{three_Gaussians} +\title{Three Gaussian Clusters} +\description{ +Matrix of 2-dimensional points forming three Gaussian clusters. +} +\examples{ +data(three_Gaussians) +plot(three_Gaussians, pch = 19) + +} +\keyword{datasets} diff --git a/man/wcec.Rd b/man/wcec.Rd new file mode 100644 index 0000000..e808a7e --- /dev/null +++ b/man/wcec.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wcec-package.R +\docType{package} +\name{wcec} +\alias{wcec-package} +\alias{wcec} +\title{wcec: ...} +\description{ +The \code{wcec} package provides ... +} +\details{ +\tabular{ll}{ +Package: \tab wcec\cr +Type: \tab Package\cr +Version: \tab 0.1\cr +Date: \tab 2024-08-07\cr +License: \tab GPL-3\cr +} +} +\seealso{ +Useful links: +\itemize{ + \item \url{https://github.com/swarm-lab/wcec} + \item \url{https://swarm-lab.github.io/wcec/} + \item Report bugs at \url{https://github.com/swarm-lab/wcec/issues} +} + +} +\author{ +Simon Garnier \href{mailto:garnier@njit.edu}{garnier@njit.edu} + +Jason Graham \href{mailto:jason.graham@scranton.edu}{jason.graham@scranton.edu} + +Maintainer: Simon Garnier \href{mailto:garnier@njit.edu}{garnier@njit.edu} +} diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..22034c4 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,3 @@ +*.o +*.so +*.dll diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..6868846 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,8 @@ +## With Rcpp 0.11.0 and later, we no longer need to set PKG_LIBS as there is +## no user-facing library. The include path to headers is already set by R. +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) + +## With R 3.1.0 or later, you can uncomment the following line to tell R to +## enable compilation with C++11 (or even C++14) where available +#CXX_STD = CXX11 diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..6868846 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,8 @@ +## With Rcpp 0.11.0 and later, we no longer need to set PKG_LIBS as there is +## no user-facing library. The include path to headers is already set by R. +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) + +## With R 3.1.0 or later, you can uncomment the following line to tell R to +## enable compilation with C++11 (or even C++14) where available +#CXX_STD = CXX11 diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..3bc4627 --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,35 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include + +using namespace Rcpp; + +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + +// wcov +Rcpp::List wcov(Eigen::MatrixXd& x, Eigen::VectorXd& w); +RcppExport SEXP _wcec_wcov(SEXP xSEXP, SEXP wSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Eigen::MatrixXd& >::type x(xSEXP); + Rcpp::traits::input_parameter< Eigen::VectorXd& >::type w(wSEXP); + rcpp_result_gen = Rcpp::wrap(wcov(x, w)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_wcec_wcov", (DL_FUNC) &_wcec_wcov, 2}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_wcec(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/fast.cpp b/src/fast.cpp new file mode 100644 index 0000000..55830b6 --- /dev/null +++ b/src/fast.cpp @@ -0,0 +1,47 @@ +#include +// [[Rcpp::depends(RcppEigen)]] + +//' @name .wcov +//' +//' @title Weighted Covariance Matrix +//' +//' @description \code{wcov} computes the estimates of the weighted covariance +//' matrix and the weighted mean of the data. +//' +//' @param x A matrix with \eqn{m} columns and \eqn{n} rows, where each column +//' represents a different variable and each row a different observation. +//' +//' @param w A non-negative and non-zero vector of weights for each observation. +//' Its length must equal the number of rows of x. +//' +//' @return A list with two components: +//' \itemize{ +//' \item \code{center}: an estimate for the center (mean) of the data. +//' \item \code{cov}: the estimated (weighted) covariance matrix. +//' } +//' +//' @author Simon Garnier, \email{garnier@@njit.edu} +//' +//' @examples +//' m <- matrix(c(rnorm(500, 6), rnorm(500, 11, 3)), ncol = 2) +//' w <- runif(500) +//' gravitree:::.wcov(m, w) +//' +// [[Rcpp::export(.wcov)]] +Rcpp::List wcov(Eigen::MatrixXd &x, Eigen::VectorXd &w) { + Eigen::VectorXd center(x.cols()); + double ws = w.sum(); + for (int i = 0; i < x.cols(); i++) { + center(i) = (x.col(i).array() * w.array()).sum() / ws; + } + + int p = x.cols(); + Eigen::VectorXd sqw = (w.array() / ws).cwiseSqrt(); + Eigen::MatrixXd X = x.array().rowwise() - center.transpose().array(); + Eigen::MatrixXd cov = Eigen::MatrixXd(p, p) + .setZero() + .selfadjointView() + .rankUpdate(X.transpose() * sqw.asDiagonal()); + + return Rcpp::List::create(Rcpp::_["center"] = center, Rcpp::_["cov"] = cov); +}