Функция Rcpp для добавления элементов вектора

Question

У меня очень длинный вектор параметров (приблизительно 4^10 элементов) и вектор индексов. Моя цель - объединить все значения параметров, индексированных в векторе индексов.

Например, если бы у меня были парас = [1,2,3,4,5,5,5] и индексы = [3,3,1,6], то я бы хотел найти кумулятивную сумму третье значение (3) дважды, первое значение (1) и шестое (5), чтобы получить 12. Существует дополнительно возможность деформирования значений параметров в соответствии с их местоположением.

Я пытаюсь ускорить реализацию R, поскольку я называю это миллионы раз.

Мой текущий код всегда возвращает NA, и я не могу увидеть, где это происходит не так

Вот функция Rcpp:

double dot_prod_c(NumericVector indices, NumericVector paras, 
        NumericVector warp = NA_REAL) { 
int len = indices.size(); 
LogicalVector indices_ok; 
for (int i = 0; i < len; i++){ 
    indices_ok.push_back(R_IsNA(indices[i])); 
} 
if(is_true(any(indices_ok))){ 
    return NA_REAL; 
} 
double counter = 0; 
if(NumericVector::is_na(warp[1])){ 
    for (int i = 0; i < len; i++){ 
     counter += paras[indices[i]]; 
    } 
} else { 
    for (int i = 0; i < len; i++){ 
     counter += paras[indices[i]] * warp[i]; 
    } 
} 
return counter; 
} 

А вот рабочая версия R:

dot_prod <- function(indices, paras, warp = NA){ 
    if(is.na(warp[1])){ 
     return(sum(sapply(indices, function(ind) paras[ind + 1]))) 
    } else { 
     return(sum(sapply(1:length(indices), function(i){ 
      ind <- indices[i] 
      paras[ind + 1] * warp[i] 
     }))) 
    } 
} 

Ниже приведен код для тестирования и бенчмаркинга с использованием пакета микрообъективов:

# testing 
library(Rcpp) 
library(microbenchmark) 

parameters <- list() 
indices <- list() 
indices_trad <- list() 

set.seed(2) 
for (i in 4:12){ 
    size <- 4^i 
    window_size <- 100 
    parameters[[i-3]] <- runif(size) 
    indices[[i-3]] <- floor(runif(window_size)*size) 
    temp <- rep(0, size) 
    for (j in 1:window_size){ 
     temp[indices[[i-3]][j] + 1] <- temp[indices[[i-3]][j] + 1] + 1 
    } 
    indices_trad[[i-3]] <- temp 
} 

microbenchmark(
    x <- sapply(1:9, function(i) dot_prod(indices[[i]], parameters[[i]])), 
    x_c <- sapply(1:9, function(i) dot_prod_c(indices[[i]], parameters[[i]])), 
    x_base <- sapply(1:9, function(i) indices_trad[[i]] %*% parameters[[i]]) 
) 
all.equal(x, x_base) # is true, does work 
all.equal(x_c, x_base) # not true - C++ version returns only NAs 

Answer 1

у меня было мало проблем, пытаясь интерпретировать вашу общую цель через ваш код, так что я просто хочу, чтобы пойти с этим объяснением

Например, если бы я имел пункты = [1,2 , 3,4,5,5,5] и индексы = [3,3,1,6] , тогда я хотел бы найти совокупную сумму третьего значения (3) дважды, первое значение (1) и шестой (5), чтобы получить 12. Существует , а также возможность деформирования значений параметров в соответствии с их местоположением.

так как это было для меня наиболее понятным.

---

Есть некоторые проблемы с кодом на C++. Для начала, вместо этого, используйте Rcpp::Nullable<> (см. Ниже). Это позволит решить несколько проблем:

1. Это более читаемо. Если вы не знакомы с классом Nullable, это в значительной степени то, на что это похоже - объект, который может быть или не быть нулевым.
2. Вам не придется выполнять какие-либо неудобные инициализации, такие как NumericVector warp = NA_REAL. Честно говоря, я был удивлен, что компилятор принял это.
3. Вам не придется беспокоиться о том, чтобы случайно забыть, что C++ использует индексирование с нулевым индексом, в отличие от R, как в этой строке: if(NumericVector::is_na(warp[1])){. Это неопределенное поведение, написанное во всем этом.

---

Вот пересмотренная версия, идя от вашего цитируемого описания проблемы выше:

#include <Rcpp.h> 

typedef Rcpp::Nullable<Rcpp::NumericVector> nullable_t; 
// [[Rcpp::export]] 
double DotProd(Rcpp::NumericVector indices, Rcpp::NumericVector params, nullable_t warp_ = R_NilValue) { 
    R_xlen_t i = 0, n = indices.size(); 
    double result = 0.0; 

    if (warp_.isNull()) { 
    for (; i < n; i++) { 
     result += params[indices[i]]; 
    }  
    } else { 
    Rcpp::NumericVector warp(warp_); 
    for (; i < n; i++) { 
     result += params[indices[i]] * warp[i]; 
    } 
    } 

    return result; 
} 

---

Вы имели некоторый сложный код для генерации данных выборки. Я не нашел времени, чтобы пройти через это, потому что это было необязательно, и не было бенчмаркинга. Вы заявили, что версия C++ не дает правильных результатов. Ваш первый приоритет должен состоять в том, чтобы заставить ваш код работать с простыми данными. Затем загрузите несколько более сложных данных. Тогда бенчмарк. Пересмотренный вариант выше работает на простых данных:

---

args <- list(
    indices = c(3, 3, 1, 6), 
    params = c(1, 2, 3, 4, 5, 5, 5), 
    warp = c(.25, .75, 1.25, 1.75) 
) 

all.equal(
    DotProd(args[[1]], args[[2]]), 
    dot_prod(args[[1]], args[[2]])) 
#[1] TRUE 

all.equal(
    DotProd(args[[1]], args[[2]], args[[3]]), 
    dot_prod(args[[1]], args[[2]], args[[3]])) 
#[1] TRUE 

---

Это также быстрее, чем версия R на этом образце данных. У меня нет оснований полагать, что это не было бы для более крупных и более сложных данных - нет ничего волшебного или особенно эффективного в отношении функций * apply; они просто более идиоматические/читаемого Р.

---

microbenchmark::microbenchmark(
    "Rcpp" = DotProd(args[[1]], args[[2]]), 
    "R" = dot_prod(args[[1]], args[[2]])) 
#Unit: microseconds 
#expr min  lq  mean median  uq max neval 
#Rcpp 2.463 2.8815 3.52907 3.3265 3.8445 18.823 100 
#R 18.869 20.0285 21.60490 20.4400 21.0745 66.531 100 
# 
microbenchmark::microbenchmark(
    "Rcpp" = DotProd(args[[1]], args[[2]], args[[3]]), 
    "R" = dot_prod(args[[1]], args[[2]], args[[3]])) 
#Unit: microseconds 
#expr min  lq  mean median  uq max neval 
#Rcpp 2.680 3.0430 3.84796 3.701 4.1360 12.304 100 
#R 21.587 22.6855 23.79194 23.342 23.8565 68.473 100 

---

Я опустил NA чек из приведенного выше примера, но это тоже может быть пересмотрены в нечто более идиоматическое, используя немного Rcpp сахара. Раньше вы делали это:

LogicalVector indices_ok; 
for (int i = 0; i < len; i++){ 
    indices_ok.push_back(R_IsNA(indices[i])); 
} 
if(is_true(any(indices_ok))){ 
    return NA_REAL; 
} 

Это немного агрессивны - вы проверяете целый вектор значений (с R_IsNA), а затем применяя is_true(any(indices_ok)) - когда вы могли бы просто разорвать преждевременно и вернуться NA_REAL по первой инстанции от R_IsNA(indices[i]), что приводит к true. Кроме того, использование push_back замедлит вашу функцию совсем немного - вам было бы лучше инициализировать indices_ok до известного размера и заполнить его доступом по индексу в вашем цикле. Тем не менее, вот один из способов конденсироваться операции:

if (Rcpp::na_omit(indices).size() != indices.size()) return NA_REAL; 

---

Для полноты картины, вот полностью сахарная роскопию версия, которая позволяет избежать петель полностью:

#include <Rcpp.h> 

typedef Rcpp::Nullable<Rcpp::NumericVector> nullable_t; 
// [[Rcpp::export]] 
double DotProd3(Rcpp::NumericVector indices, Rcpp::NumericVector params, nullable_t warp_ = R_NilValue) { 
    if (Rcpp::na_omit(indices).size() != indices.size()) return NA_REAL; 

    if (warp_.isNull()) { 
    Rcpp::NumericVector tmp = params[indices]; 
    return Rcpp::sum(tmp);  
    } else { 
    Rcpp::NumericVector warp(warp_), tmp = params[indices]; 
    return Rcpp::sum(tmp * warp); 
    } 
} 

/*** R 

all.equal(
    DotProd3(args[[1]], args[[2]]), 
    dot_prod(args[[1]], args[[2]])) 
#[1] TRUE 

all.equal(
    DotProd3(args[[1]], args[[2]], args[[3]]), 
    dot_prod(args[[1]], args[[2]], args[[3]])) 
#[1] TRUE 

*/