2017-02-15 2 views
0

Обратите внимание, что эта ошибка была взята из более крупного контекста, о котором я не могу полностью отчитываться здесь.Неожиданное поведение в Rcpp

У меня есть следующие функции в файле fun.cpp

#include <RcppArmadilloExtensions/sample.h> 

using namespace Rcpp; 

// [[Rcpp::depends(RcppArmadillo)]] 

arma::vec colMeans(arma::mat data){ 

    int n_0 = data.n_rows; 

    arma::vec xbar(data.n_cols); 
    for(int i = 0; i < data.n_rows; i++){ 
     for(int j = 0; j < data.n_cols; j++){ 
     xbar[j] += data(i,j) /n_0; 
     } 
    } 

    return xbar; 
} 

// [[Rcpp::export]] 
List PosteriorNIW(arma::mat data, arma::vec mu0, double lambda0, 
       double df0, arma::mat V){ 

    // Compute posterior 
    int n = data.n_rows; 

    arma::vec xbar = colMeans(data); 

    double lambdan = lambda0 + n; 

    arma::vec mun = (lambda0 * mu0 + n * xbar)/lambdan; 

    arma::mat S; 
    S.zeros(data.n_cols, data.n_cols); 
    for(int i = 0; i < n; i++){ 
     S += (arma::conv_to<arma::vec>::from(data.row(i)) - xbar) * arma::trans(arma::conv_to<arma::vec>::from(data.row(i)) - xbar); 
    } 

    arma::mat Vn = V + S + ((lambda0*n)/(lambda0 + n)) * (xbar - mu0) * arma::trans(xbar - mu0); 

    return List::create(_["mun"] = mun, 
        _["Vn"] = Vn, 
        _["lambdan"] = lambdan); 
} 

Вызывать Сейчас:

library(Rcpp); library(RcppArmadillo) 
mu0 <- c(3,3) 
V0 <- matrix(c(2.5,0.0,0.0,2.5), nrow = 2) 
sourceCpp("fun.cpp") 

data <- cbind(rep(5,15),rep(0,15)) 
PosteriorNIW(data, mu0, 1, 1, V0) 

дает ожидаемый результат.

$mun 
    [,1] 
[1,] 4.8750 
[2,] 0.1875 

$Vn 
    [,1] [,2] 
[1,] 6.250 -5.6250 
[2,] -5.625 10.9375 

$lambdan 
[1] 16 

Теперь, если добавить в файл fun.cpp следующие функции (опять-таки, они взяты из большего контекста, так что не пытайтесь понять, но просто вставить их) странные вещи происходит:

// [[Rcpp::export]] 
NumericMatrix myFun(arma::mat t_dish, arma::cube data){ 
    int l = 0; 
    for(int j = 0; j < data.n_rows; j++){ 
     l++; 
    } 
    NumericMatrix Dk(l, 2); 
    return Dk; 
} 

// [[Rcpp::export]] 
int myFun2(arma::cube n_cust){ 

    arma::mat temp = n_cust.subcube(arma::span(0), arma::span(), arma::span()); 
    int i; 
    for(i = 0; i < n_cust.n_cols; i++){ 
     arma::rowvec temp2 = temp.row(i); 
    } 

    return i + 1; 
} 

// [[Rcpp::export]] 
arma::vec myFun3(arma::mat k_tables){ 
    arma::vec temp(k_tables.n_cols * k_tables.n_rows); 
    int l = 0; 
    if(!R_IsNA(k_tables(0,0))){ 
     l++; 
    } 

    arma::vec temp2(l); 
    arma::vec tmp3 = sort(temp2); 
    return tmp3; 
} 

double myFun4(arma::vec x, double nu, arma::vec mu, arma::mat Sigma){ 
    arma::vec product = (arma::trans(x - mu) * arma::inv(Sigma) * (x - mu)); 
    double num = pow(1 + (1/nu) * product[0], - (nu + 2)/2); 
    double den = pow(sqrt(M_PI * nu),2) * sqrt(arma::det(Sigma)); 
    return num/den; 
} 

bool myFun5(NumericVector X, double z) { 
    return std::find(X.begin(), X.end(), z)!=X.end(); 
} 

, вызывающий PosteriorNIW(data, mu0, 1, 1, V0), неоднократно начинает давать разные результаты каждый раз. Обратите внимание, что в функциях нет случайности, и очевидно, что эти функции не имеют никакого влияния, поскольку они не вызываются в исходной функции.

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

Я знаю, что удаление этой функции (даже одной из них) устраняет проблему, но ясно, что это не приемлемое решение, когда я работаю с большим количеством функций.

Я хотел бы знать, могут ли другие пользователи реплицировать это поведение, и если да, если есть исправление для него.

Спасибо заранее

EDIT:

Версия R является 3.3.2 и Rtools составляет 3,4. И Rcpp, и RcppArmadillo являются актуальными

+1

Пожалуйста, используйте отступ. Форматирование выше * abysmal * и очень затрудняет чтение вашего кода. – nrussell

+0

Извините, он удалил отступы при вставке. Должно быть проще читать сейчас – adaien

+1

Вы не платите за пробел. Многие из них используют четыре, некоторые бедные души считают, что два лучше. Один по-прежнему по существу не читается. Если вы хотите получить бесплатную помощь от нас, не мешайте нам. –

ответ

1

Вы не обнуляете xbar в своей функции colMeans. Если я это сделать:

arma::vec colMeans(arma::mat data){ 

    int n_0 = data.n_rows; 

    arma::vec xbar; 
    xbar.zeros(data.n_cols); 
    for(int i = 0; i < data.n_rows; i++){ 
     for(int j = 0; j < data.n_cols; j++){ 
     xbar[j] += data(i,j) /n_0; 
     } 
    } 

    return xbar; 
} 

Я получаю каждый раз:

> PosteriorNIW(data, mu0, 1, 1.1, V0) 
$mun 
     [,1] 
[1,] 4.8750 
[2,] 0.1875 

$Vn 
     [,1] [,2] 
[1,] 6.250 -5.6250 
[2,] -5.625 10.9375 

$lambdan 
[1] 16 

Даже когда я добавить свой дополнительный блок кода.

Я не знаю, задокументированы ли эти векторы, чтобы они были инициализированы нолем их конструктором (в этом случае это может быть ошибка там) или нет, и в этом случае это ваша ошибка!

+0

Отлично, спасибо большое! – adaien

+0

Моя техника отладки в этих ситуациях заключается в том, чтобы распечатывать вещи по мере их расчета и выяснять, куда вещи уходят. Это привело меня к функции 'colMeans', а затем я заметил ненулевой вектор. – Spacedman

+0

Я делаю то же самое обычно, но когда я пытался распечатать 'n' в функции' PosteriorNIW', просто добавив печать, чтобы функция работала и удалила ошибку. Cheers – adaien

Смежные вопросы