Rcpp：循环中的计算因错误停止 "Not a matrix"

Question

在 R 脚本中，我获取了一个 cpp 文件来进行一些计算。在该 R 脚本中，调用了 cpp 文件中定义的函数，并提供了一个矩阵和一个整数。在循环几轮之后，它会给出错误“不是矩阵”（在代码行 resid = (x(_,j) - x(_,i))*(x(_,j) - x(_,i)); 中），即使在它工作之前的轮次也是如此。 R 脚本：

## all together
# rm(list=ls())
library(RcppArmadillo)
library(Rcpp)

sourceCpp("~/test.cpp",verbose = FALSE)

cat("start loop")
for(n in c(45:46)){
  cat("\n", n, "\n")
  p_m <- matrix(data=rnorm(n^2,1,1),nrow = n, ncol=n)
  print(class(p_m))
  print(some_function(p_m,nosamples=10))
}
cat("finished")

我通过命令行启动这个 R 脚本。 R 版本 R-4.1.0。在 R-Studio 中，它因致命错误而崩溃。

cpp 文件：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector some_function(NumericMatrix x,int nosamples) {
  int ncol = x.ncol();
  NumericVector out2(nosamples);
  int loops;
  int loops2;
  double result=0;
  NumericVector::iterator it;
  double acc = 0;
  NumericVector resid(ncol);
  NumericVector out(ncol*(ncol-1)/2);
  loops2=0;
  std::cout << nosamples << std::endl;
  std::cout << (ncol-1) << std::endl;
  std::cout << ncol*(ncol-1)/2 << std::endl;
  while(loops2 < (nosamples)){
    std::cout << "loops2:" << std::endl;
    std::cout << loops2 << std::endl;
    loops=0;
    int i;
    int j;
    for(j=0;j<(ncol-1);++j){
      std::cout << " j: " << j << std::endl;
      for (i = (j+1); i < (ncol); ++i) {
        std::cout << " i: " << i << std::endl;
          resid = (x(_,j) - x(_,i))*(x(_,j) - x(_,i)); //here it stops
          std::cout << " i: " << i << std::endl; 
          for(int ii=0; ii<ncol;++ii){
            acc += resid[i];
          }
         result=sqrt(acc);
          loops += 1;
          out[loops] = result;
          std::cout << " i: " << i << std::endl;
      }
    }
    std::cout << "loops:" << std::endl;
    std::cout << loops << std::endl;
    out = out[out > 0];
    it = std::min_element(out.begin(), out.end());
    out= *it;
    std::cout << out << std::endl;
    loops2 += 1;
    out2[loops2]=out[0];
  }
  std::cout << "cpp finished" << std::endl;
  return(out2);
}

谁能解释一下问题出在哪里？

谢谢和亲切的问候

Edit

我修改了cpp文件中的一些东西（如下所示），错误消失了。首先我想，一切都很好。但是当我增加循环次数时，另一个问题出现了：函数中断，但没有显示错误。它在循环编号 543（“循环 2：543”）之后中断。至少它在每一轮 while 循环中使用相同的数据做同样的事情。

我修改了 R 脚本和 ccp 文件，使这个问题（至少在我的机器上）可以重现。

我知道这段代码似乎毫无意义，但它是一个更大程序的一部分，我想在这里举一个最简单的例子。

R 脚本：

## all together
# rm(list=ls())
library(RcppArmadillo)
library(Rcpp)

sourceCpp("~/test.cpp",verbose = FALSE)

cat("start loop")
for(n in c(100:101)){
  cat("\n", n, "\n")
  p_m <- matrix(data=rnorm(n^2,1,1),nrow = n, ncol=n)
  print(class(p_m))
  print(some_function(p_m,nosamples=800))
}
cat("finished")

cpp 文件：

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
using namespace Rcpp;
using  Eigen::Map;
using  Eigen::VectorXd;
typedef  Map<VectorXd>  MapVecd;


// [[Rcpp::export]]
NumericVector some_function(NumericMatrix x,int nosamples) {
  
  int ncol = x.ncol();
  NumericVector out(ncol*(ncol-1)/2);
  NumericVector out2(nosamples);
  NumericVector out3(ncol*(ncol-1)/2);
  NumericVector resid(ncol);
  int loops;
  int loops2;
  double result=0;
  double acc = 0;
  
  int show_cout=0;
  
  loops2=0;
  
  std::cout << nosamples << std::endl;
  std::cout << (ncol-1) << std::endl;
  std::cout << ncol*(ncol-1)/2 << std::endl;
  
  while(loops2 < (nosamples)){
    
    std::cout << "loops2:" << loops2 << std::endl;
    
    loops=0;
    int i;
    int j;
    
    for(j=0;j<(ncol-1);++j){
      // std::cout << " j: " << j << std::endl;
      for (i = (j+1); i < (ncol); ++i) {
        
        if(show_cout==1){
          std::cout << " i: " << i << std::endl;
        }
        
        resid = (x(_,j) - x(_,i))*(x(_,j) - x(_,i));
        
        if(show_cout==1){
          std::cout << " i: " << i << std::endl;
        }
        
        for(int ii=0; ii<ncol;++ii){
          acc += resid[ii];
          
        }
        
        result=sqrt(acc);
        
        loops += 1;
        
        out[loops] = result;
        
        if(show_cout==1){
          std::cout << " i: " << i << std::endl;
        }
      }
      
    }
    
    // std::cout << "loops:" << loops << std::endl;
    
    // 
    out = out[out > 0];
    
    const MapVecd xy(as<MapVecd>(out));
    
    out3=xy.minCoeff();
    
    out2[loops2]=out3[0];
    
    loops2 += 1;
  }
  std::cout << "cpp finished" << std::endl;
  return(out2);
  
}

Answer 1

这里有两件事：

使用 out[loops++] = result; 而不是 loops += 1; out[loops] = result; 因为您从 1 开始，并且可能访问了该向量范围之外的最后一个元素。
使用 for(int ii=0; ii<ncol;++ii){ double eps = x(ii, j) - x(ii, i); acc += eps * eps; } 而不是依赖这个 resid 向量。

Rcpp：循环中的计算因错误停止 "Not a matrix"

Rcpp: Calculation in loop stops with error "Not a matrix"

loops

r

matrix

rcpp