时间戳

首页 > 解决方案 > Rcpp 函数抱怨未初始化的变量

r - Rcpp 函数抱怨未初始化的变量

问题描述

在创建可以使用 Rcpp 从 R 调用的 C++ 函数的第一次尝试中，我有一个简单的函数来使用 Prim 算法从距离矩阵计算最小生成树。此函数已从 ANSI C 的旧版本转换为 C++（运行良好）。

这里是：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
DataFrame primlm(const int n, NumericMatrix d)
{
    double const din = 9999999.e0;
    long int i1, nc, nc1;
    double dlarge, dtot;
    NumericVector is, l, lp, dist;
    l(1) = 1;
    is(1) = 1;
    for (int i=2; i <= n; i++) {
        is(i) = 0;
    }
    for (int i=2; i <= n; i++) {
        dlarge = din;
        i1 = i - 1;
        for (int j=1; j <= i1; j++) {
            for (int k=1; k <= n; k++) {
                if (l(j) == k)
                    continue;
                if (d[l(j), k] > dlarge)
                    continue;
                if (is(k) == 1)
                    continue;
                nc = k;
                nc1 = l(j);
                dlarge = d(nc1, nc);
            }
        }
        is(nc) = 1;
        l(i) = nc;
        lp(i) = nc1;
        dist(i) = dlarge;
    }
    dtot = 0.e0;
    for (int i=2; i <= n; i++){
        dtot += dist(i);
    }
    return DataFrame::create(Named("l") = l, 
                          Named("lp") = lp, 
                          Named("dist") = dist, 
                          Named("dtot") = dtot);
}

当我在 RStudio 下使用 Rcpp 编译此函数时，我收到两个警告，抱怨变量 'nc' 和 'nc1' 尚未初始化。坦率地说，我无法理解这一点，因为在我看来，这两个变量都在第三个循环中被初始化。另外，为什么没有关于变量“i1”的类似抱怨？

也许不足为奇的是，当尝试从 R 调用这个函数时，使用下面的代码，我得到的是 R 系统的崩溃！

# Read test data
df <- read.csv("zygo.csv", header=TRUE)
lonlat <- data.frame(df$Longitude, df$Latitude)
colnames(lonlat) <- c("lon", "lat")

# Compute distance matrix using geosphere library
library(geosphere)
d <- distm(lonlat, lonlat, fun=distVincentyEllipsoid)

# Calls Prim minimum spanning tree routine via Rcpp
library(Rcpp)
sourceCpp("Prim.cpp")
n <- nrow(df)
p <- primlm(n, d)

这是我用于测试目的的数据集：

"Scientific name",Locality,Longitude,Latitude Zygodontmys,Bush Bush
Forest,-61.05,10.4 Zygodontmys,Cerro Azul,-79.4333333333,9.15
Zygodontmys,Dividive,-70.6666666667,9.53333333333 Zygodontmys,Hato El
Frio,-63.1166666667,7.91666666667 Zygodontmys,Finca Vuelta
Larga,-63.1166666667,10.55 Zygodontmys,Isla
Cebaco,-81.1833333333,7.51666666667 Zygodontmys,Kayserberg
Airstrip,-56.4833333333,3.1 Zygodontmys,Limao,-60.5,3.93333333333
Zygodontmys,Montijo Bay,-81.0166666667,7.66666666667
Zygodontmys,Parcela 200,-67.4333333333,8.93333333333 Zygodontmys,Rio
Chico,-65.9666666667,10.3166666667 Zygodontmys,San Miguel
Island,-78.9333333333,8.38333333333
Zygodontmys,Tukuko,-72.8666666667,9.83333333333
Zygodontmys,Urama,-68.4,10.6166666667
Zygodontmys,Valledup,-72.9833333333,10.6166666667

谁能给我一个提示？

标签： rrcpp

解决方案

如果三个语句之一为真，则永远不会达到nc和的初始化。您的数据可能无法做到这一点，但编译器无法知道这一点。nc1if

但是，这不是崩溃的原因。当我运行您的代码时，我得到：

索引越界：[index=1; 范围=0]。

这来自这里：

NumericVector is, l, lp, dist;
l(1) = 1;
is(1) = 1;

声明 a 时NumericVector，如果要按索引分配值，则必须告知所需的大小。在你的情况下

NumericVector is(n), l(n), lp(n), dist(n);

可能会奏效。您必须仔细分析 C 代码 wrt 内存分配和数组边界。

或者，您可以按原样使用 C 代码并使用 Rcpp 构建包装函数，例如

#include <array>
#include <Rcpp.h>
using namespace Rcpp;

// One possibility for the function signature ...
double prim(const int n, double *d, double *l, double *lp, double *dist) {
    ....
} 

// [[Rcpp::export]]
List primlm(NumericMatrix d) {
    int n = d.nrow();
    std::array<double, n> lp;    // adjust size as needed!
    std::array<double, n> dist;  // adjust size as needed!
    double dtot = prim(n, d.begin(), l.begin(), lp.begin(), dist.begin()); 

    return List::create(Named("l") = l, 
                        Named("lp") = lp, 
                        Named("dist") = dist, 
                        Named("dtot") = dtot);
}

笔记：

  • 我返回 aList而不是 aDataFrame因为dtot是标量值。
  • 上面的代码是为了说明这个想法。如果不进行调整，它很可能无法正常工作！

推荐阅读