Блог

Функция C , которую я написал с использованием Rcpp, выдает разные выходные данные в зависимости от того, есть ли у меня оператор Rcout or Rprintf в коде. Приведенный ниже код возвращает 1 правильное значение для функции с оператором печати H_sigma_1() . Однако для H_sigma_2() функции без оператора печати функция возвращает 2. Я тестировал это на Ubuntu 16.04.1, а также CentOS 6.8. Хотя я не смог воспроизвести эту ошибку в Windows 10. Следовательно, это, похоже, проблема Linux.

Код Rcpp

 library(Rcpp)

cppFunction ( 
  "double H_sigma_1(IntegerVector sigma, NumericMatrix J, NumericVector h)
  {
    double first_sum, second_sum = 0;
    int n = sigma.size();

    for(int i = 0; i < n; i  )
    {
      for(int j = 0; j < n; j  )
      {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum  = J(i, j) * sigma[i] * sigma[j];
        Rcout << first_sum << std::endl;
      }
      second_sum  = h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
  }"
)

cppFunction (
  "double H_sigma_2(IntegerVector sigma, NumericMatrix J, NumericVector h)
  {
    double first_sum, second_sum = 0;
    int n = sigma.size();

    for(int i = 0; i < n; i  )
    {
      for(int j = 0; j < n; j  )
      {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum  = J(i, j) * sigma[i] * sigma[j];
        // Rcout << first_sum << std::endl;
      }
      second_sum  = h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
  }"
)
  

Пример вызова

Настройка:

 n = 2
params = rep(1, n) 
h = rep(params[1], n)
J = toeplitz(c(0, params[2], rep(0, n - 2)))
  

Тест 1:

 H_sigma_1(c(-1, -1), J, h)
  

Вывод:

 1
[1] 1
  

Тест 2

 H_sigma_2(c(-1, -1), J, h)
  

Вывод

 [1] 2
  

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

 double first_sum, second_sum = 0;
  

не совпадает с:

 double first_sum = 0, second_sum = 0;
  

Смотрите Флаги предупреждения компилятора:

 file11d36f564b15.cpp:17:5: warning: variable 'first_sum' is uninitialized when used here [-Wuninitialized]
    first_sum  = J(i, j) * sigma[i] * sigma[j];
    ^~~~~~~~~
file11d36f564b15.cpp:8:21: note: initialize the variable 'first_sum' to silence this warning
    double first_sum, second_sum = 0;
  

Вы сразу же используете:

 first_sum  = J(i, j) * sigma[i] * sigma[j];
  

без установки first_sum = ...;

Кроме того, другая проблема заключается в:

 second_sum = 0;
  

инициализирует double целочисленным значением. Хотя эта проблема незначительна по объему, чтобы исправить это, все, что нужно сделать, это использовать 0.0 вместо 0 .

 second_sum = 0.0;
  

Это также относится к first_sum .

Код с вышеуказанными исправлениями:

 library(Rcpp)

cppFunction ( 
    "double H_sigma_1(IntegerVector sigma, NumericMatrix J, NumericVector h)
    {
    double first_sum = 0.0, second_sum = 0.0;
    int n = sigma.size();

    for(int i = 0; i < n; i  ) {
      for(int j = 0; j < n; j  ) {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum  = J(i, j) * sigma[i] * sigma[j];
        Rcout << first_sum << std::endl;
      }
      second_sum  = h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
    }"
)

cppFunction ( 
    "double H_sigma_2(IntegerVector sigma, NumericMatrix J, NumericVector h)
    {
    double first_sum = 0.0, second_sum = 0.0;
    int n = sigma.size();

    for(int i = 0; i < n; i  ) {
      for(int j = 0; j < n; j  ) {
        // skip inside loop if i >= j to stop double counting
        if(i >= j) {continue;}
        first_sum  = J(i, j) * sigma[i] * sigma[j];
      }
      second_sum  = h[i] * sigma[i];
    }
    return(-1.0 * first_sum - second_sum);
    }"
)
  

Тест:

 n = 2
params = rep(1, n) 
h = rep(params[1], n)
J = toeplitz(c(0, params[2], rep(0, n - 2)))

H_sigma_1(c(-1, -1), J, h)
  

Вывод:

 1
[1] 1
  

Тест 2:

 H_sigma_2(c(-1, -1), J, h)
  

Вывод:

 [1] 1