c++ - 尝试在 SET_VECTOR_ELT 中设置索引 0/0

Question

亲爱的。在花了几个月的时间学习 C++ 以编写可与 R 软件一起使用的代码并成功编写代码之后，我感到沮丧，因为我意识到我需要使用 Rcpp 命名法重写我的代码。在阅读各种 Rcpp 资料和讨论论坛数周之后，我能够使 sourcecpp 功能正常工作。但是，使用该函数出现了一个错误，我不知道它的来源。

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS:   /usr/lib/libblas/libblas.so.3.7.0
LAPACK: /usr/lib/lapack/liblapack.so.3.7.0

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] Rcpp_1.0.1           lba_2.4.4            rgl_0.100.26        
 [4] scatterplot3d_0.3-41 plotrix_3.7-6        alabama_2015.3-1    
 [7] numDeriv_2016.8-1.1  MASS_7.3-51.4        nvimcom_0.9-82      
[10] colorout_1.2-1      

loaded via a namespace (and not attached):
 [1] knitr_1.23              magrittr_1.5            xtable_1.8-4           
 [4] R6_2.4.0                tools_3.6.1             webshot_0.5.1          
 [7] xfun_0.8                miniUI_0.1.1.1          htmltools_0.3.6        
[10] crosstalk_1.0.0         digest_0.6.20           shiny_1.3.2            
[13] later_0.8.0             htmlwidgets_1.3         promises_1.0.1         
[16] manipulateWidget_0.10.0 mime_0.7                compiler_3.6.1         
[19] jsonlite_1.6            httpuv_1.5.1            

saa.cpp 是:

#include <Rcpp.h>
#include <iostream>
#include <cmath>
#include <vector>
#include <cstdlib>

using namespace std;
using namespace Rcpp;

// [[Rcpp::plugins(cpp11)]]

NumericMatrix Tfunc(NumericMatrix Told, double amm, Function funT) {
  NumericMatrix TT = funT(Told, amm);
  return TT;
}

/*
 * Criando uma função para fazer o produto entre duas matrizes
 */
NumericMatrix multiplica_mat(NumericMatrix A, NumericMatrix B)
{
  int i = 0, j = 0, p = 0;
  int I = A.nrow();
  int J = B.ncol();
  int K = B.nrow();
  double sum = 0;
  NumericMatrix C(I,J);

  for(i = 0; i < I; i++)
  {
    for(j = 0; j < J; j++)
    {
      for(p = 0; p < K; p++)
      {
        sum += A(i,p) * B(p,j);
      }
      C(i,j) = sum;
      sum = 0;
    }
  }

  return C;
}

/*
 * Criando função para calcular inversa
 */
NumericMatrix inversa(NumericMatrix mat){
  int K = mat.ncol();
  NumericMatrix matinversa(K,K);
  double determinante = 0;

  for(int i = 0; i < K; i++){
    determinante = determinante + (mat(0,i) * (mat(1,(i+1)%3)* mat(2,(i+2)%3) - mat(1,(i+2)%3) * mat(2,(i+1)%3)));
  }

  for(int i = 0; i < K; i++){
    for(int j = 0; j < K; j++){
      matinversa(i,j) =((mat((j+1)%3,(i+1)%3) * mat((j+2)%3,(i+2)%3)) - (mat((j+1)%3,(i+2)%3) * mat((j+2)%3,(i+1)%3)))/determinante;
    }
  }
  return matinversa;
}

// [[Rcpp::export]]
List saacpp(Function funT,
    NumericMatrix A,
    NumericMatrix T,
    double am,
    double Da)

{
  NumericMatrix Tn = T; 
  List Am;

  for (int j = 0; j < 10; j++) {
    Am[j] = multiplica_mat(A,inversa(Tn));
    Tn = Tfunc(Tn, am, funT);
    am = am*Da;
  }

  List result;
  result["Am"] = Am;

  return result;
} 

R代码是:

Tfunc <- function(Told,amm){
  k <- ncol(Told)
  id1 <- sample(k, size = 1, replace = TRUE)
  id2 <- sample(k, size = 1, replace = TRUE)
  id3 <- sample(c(1:k)[-id2], size = 1, replace = TRUE)
  Told[id1, id2] <- Told[id1, id2] + amm
  Told[id1, id3] <- Told[id1, id3] - amm

  return(Told)
}

A <- matrix(c(0.079,0.141,0.231,0.338,0.437,0.616,0.774,0.803,0.152,0.399,0.32,0.044,0.146,0.056,0.617,0.263,0.244,0.34),ncol=3) 
T0 <- diag(3)

sourceCpp("saa_teste.cpp")
saacpp(Tfunc,
       A,
       T0,
       am=40,
       Da=0.95)

谢谢!

Answer 1

问题是您创建了一个空的 Rcpp::List“Am”，然后尝试设置 Am 的元素 j。您需要更改行

List Am;

到

List Am(10);

或者换行

Am[j] = multiplica_mat(A,inversa(Tn));

到

Am.push_back(multiplica_mat(A,inversa(Tn)));

我会选择第一个解决方案(然后再一次，我会遵循 coatless 的评论并移至 RcppArmadillo 以解决任何涉及矩阵计算的问题)。

在这些修复之后，您将获得正确的输出

$Am
$Am[[1]]
      [,1]  [,2]  [,3]
[1,] 0.079 0.774 0.146
[2,] 0.141 0.803 0.056
[3,] 0.231 0.152 0.617
[4,] 0.338 0.399 0.263
[5,] 0.437 0.320 0.244
[6,] 0.616 0.044 0.340

# [Some output not shown here]

$Am[[10]]
         [,1]      [,2]       [,3]
[1,] 1.343188 0.2328353 -0.5770233
[2,] 1.347821 0.2307821 -0.5786032
[3,] 1.328719 0.2510093 -0.5797280
[4,] 1.341242 0.2403128 -0.5815551
[5,] 1.343531 0.2411428 -0.5836740
[6,] 1.339382 0.2464418 -0.5858239