RcppArmadillo: Pengindeksan negatif dalam loop-for

Question

Saya baru mengenal Rcpp dan mencoba melakukan perhitungan berdasarkan pengindeksan negatif di for()-loop menggunakan RcppArmadillo. Saya sudah menemukan bahwa pengindeksan negatif di RcppArmadillo tidak begitu mudah, tetapi dapat dilakukan melalui vektor elemen yang harus disimpan (seperti yang saya temukan di sini) . Tampaknya agak sulit bagi saya ketika elemen yang akan dihapus adalah indeks perulangan. Saya mencoba menerapkan pendekatan terakhir dalam jawaban ini , tetapi tidak berhasil. Adakah cara mudah untuk menentukan vektor elemen yang tidak termasuk elemen dengan indeks loop?

Jadi, saya mencoba mencari padanan di RcppArmadillo untuk y[-i]di MWE berikut:

# In R:
# Input 
n <- 10
y <- 1:20

# Computation
x <- rep(NA, n)
for(i in 1:n){
x[i] <- sum(y[-i])
}

Kode Rcpp saya sejauh ini:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>   

// [[Rcpp::export]]
arma::vec rcpp_sum (arma::vec y, int n){

arma::vec x(n);
arma::vec ind(n);

for (i=0; i<n; i++){
ind[i] = /*no idea...*/
x[i] = sum(y[ind]);
}

return x;
}

Bantuan apa pun sangat dihargai!

Answer 1

Untuk tugas seperti ini, mungkin lebih baik melewatkan indeks yang dipermasalahkan. Dalam standar C++kami akan memeriksa indeks dan melewatkannya. Sesuatu seperti ini:

// [[Rcpp::export]]
arma::vec rcpp_sum (arma::vec y, int n){
    
    arma::vec x(n);
    
    for (int i = 0; i < n; i++) {
        x[i] = 0; // Initialize value
        
        for (int j = 0; j < y.size(); ++j) {
            if (i != j) {
                x[i] += y[j];
            }
        }
    }
    
    return x;
}

Di atas, kita menjauh dari sintaks gula. IMO tidak apa-apa dalam kasus seperti ini karena alternatifnya tidak terlalu rumit. Sementara kami menyederhanakan, ketergantungan pada RcppArmadillotidak diperlukan karena kami hanya dapat menggunakan pureRcpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector pure_rcpp_sum (NumericVector y, int n){
    
    NumericVector x(n);
    
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < y.size(); ++j) {
            if (i != j) {
                x[i] += y[j];
            }
        }
    }
    
    return x;
}

Memverifikasi keluaran:

all.equal(as.vector(rcpp_sum(y, n)), x)
[1] TRUE

all.equal(pure_rcpp_sum(y, n), x)
[1] TRUE

Memperbarui

Sesuai permintaan OP, kami memiliki pendekatan yang dioptimalkan Runtuk tujuan khusus ini. Di atas menunjukkan bagaimana untuk menyerang masalah yang sangat spesifik hanya menjumlahkan nilai dari vektor yang meninggalkan satu nilai C++. Ini dimaksudkan untuk bersifat pedagogis dan belum tentu cara terbaik untuk tugas khusus ini (seperti yang akan kami tunjukkan di bawah).

Sebelum kami menunjukkan Rkode sederhana , saya ingin menunjukkan bahwa kekhawatiran OP memiliki pernyataan bersyarat sederhana di loop dalam C++tidak perlu ditakuti (seperti kasus di base R). Dari apa yang saya tahu, subset seperti yang ditunjukkan dalam tautan dari OP, adalah O (n) dan memiliki overhead tambahan dari vektor logis ekstra. Solusi yang kami usulkan di atas harus lebih efisien karena pada dasarnya melakukan hal yang sama tanpa objek tambahan.

Sekarang, untuk kode yang diperbarui:

baseR <- function(y, n) {
    mySum <- sum(y)
    vapply(1:n, function(x) mySum - y[x], FUN.VALUE = 1)    
}

## Here is the OP code for reference
OP <- function(y, n) {
    x <- rep(NA, n)
    for(i in 1:n) {x[i] <- sum(y[-i])}
    x
}

Itu dia. Ini secepat kilat juga:

huge_y <- rnorm(1e6)
huge_n <- 1e3

system.time(t1 <- baseR(huge_y, huge_n))
 user  system elapsed 
 0.003   0.000   0.003

system.time(t2 <- pure_rcpp_sum(huge_y, huge_n))
 user  system elapsed 
2.776   0.003   2.779 

system.time(t3 <- OP(huge_y, huge_n))
 user  system elapsed 
9.555   1.248  10.805

all.equal(t1, t2)
[1] TRUE

all.equal(t1, t3)
[1] TRUE