R - хранить матрицу расстояния в Rcpp

Question

У меня возникли проблемы с определением, как хранить матрицу расстояния в Rcpp.

Давайте представим себе, что я хочу, чтобы хранить следующую функцию в дистанционной матрице n*n лиц (не квадратуру sum, потому что я не уверен, как сделать это в rcpp.

distxy = function(x,y) sum (x - y) 

В этом например, я хочу, чтобы сравнить попарно 3 лиц

 [,1] [,2] [,3] [,4] 
[1,] 24 24 22 20 
[2,] 21 24 30 20 
[3,] 44 34 41 13 

в R я бы петле функции через матрицу, как этот

mat = matrix(0, nrow(d), nrow(d)) 

    len = nrow(d) 
    mat = matrix(0, len, len) 

    for(j in 1:len){ 
    for(i in 1:len){ 
     mat[j,i] = distxy(d[j,], d[i,]) 
    } 
    } 

и получить (я могу квадрат результаты, но это неважно здесь)

 [,1] [,2] [,3] 
[1,] 0 -5 -42 
[2,] 5 0 -37 
[3,] 42 37 0 

У меня возникают проблемы, делая то же самое в rcpp

То, что я достиг до сих пор

// [[Rcpp::export]] 
NumericVector FunCpp(NumericMatrix x) { 
    int nrow = x.nrow(); 
    NumericMatrix out(nrow); 

    for (int i = 0; i < nrow; i++) { 
    for (int j = 0; j < nrow; j++) { 
     out[i,j] = sum(x(i,_) - x(j,_)) ; 
    } 
    } 
    return out; 
} 

Но матрица расстояния неверна. Есть идеи ?

d = rbind(c(24, 24, 22, 20), 
     c(21, 24, 30, 20), 
     c(44, 34, 41, 13)) 

Answer 1

Есть несколько синтаксических ошибок в коде Rcpp:

• Возвращение NumericVector вместо NumericMatrix
• Использование operator[] проиндексировать двумя размерами (out[i,j])

Предлагаемая уборка:

#include <Rcpp.h> 

inline double distxy(Rcpp::NumericVector x, Rcpp::NumericVector y) { 
    return Rcpp::sum(x - y); 
} 

// [[Rcpp::export]] 
Rcpp::NumericMatrix FunCpp(Rcpp::NumericMatrix x) { 
    int nrow = x.nrow(); 
    Rcpp::NumericMatrix out(nrow); 

    for (int i = 0; i < nrow; i++) { 
     for (int j = 0; j < nrow; j++) { 
      out(j, i) = distxy(x.row(j), x.row(i)); 
     } 
    } 

    return out; 
} 

---

Тестирование против вашей функции R,

m <- matrix(
    c(24, 24, 22, 20, 
     21, 24, 30, 20, 
     44, 34, 41, 13), 
    nrow = 3, byrow = TRUE 
) 

all.equal(FunR(m), FunCpp(m)) 
#[1] TRUE 

---

Что касается возведения в квадрат, вы можете использовать std::pow внутри distxy,

return std::pow(Rcpp::sum(x - y), 2); 

или внутри FunCpp во внутреннем петля:

out(j, i) = std::pow(distxy(x.row(j), x.row(i)), 2); 

---

distxy <- function(x,y) sum(x - y) 

FunR <- function(d) { 
    len <- nrow(d) 
    mat <- matrix(0, len, len) 

    for(j in 1:len){ 
     for(i in 1:len){ 
      mat[j,i] <- distxy(d[j,], d[i,]) 
     } 
    } 
    mat 
}