Функция Rcpp, жалующаяся на неинициализированные переменные

Question

В самой первой попытке создать функцию C ++, которую можно вызывать из R с использованием Rcpp, у меня есть простая функция для вычисления минимального остовного дерева из матрицы расстояний с использованием алгоритма Прима.Эта функция была преобразована в C ++ из предыдущей версии в ANSI 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);
}

Когда я компилирую эту функцию, используя Rcpp под RStudio, яполучить два предупреждения, жалуясь, что переменные 'nc' и 'nc1' не были инициализированы.Честно говоря, я не мог этого понять, поскольку мне кажется, что обе переменные инициализируются внутри третьего цикла.Кроме того, почему нет аналогичных жалоб на переменную 'i1'?

Возможно, неудивительно, что при попытке вызвать эту функцию из R с использованием приведенного ниже кода, что яget - это сбой системы 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

Может кто-нибудь дать мне подсказку?

Answer 1

Инициализация nc и nc1 никогда не будет достигнута, если один из трех операторов if верен.Возможно, это невозможно с вашими данными, но компилятор не знает, как это сделать.

Однако это не является причиной сбоя.Когда я запускаю ваш код, я получаю:

Индекс за пределами: [index = 1;экстент = 0].

Это происходит здесь:

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

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

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

может работать.Вы должны тщательно проанализировать код C по отношению к выделению памяти и границам массивов.

В качестве альтернативы вы можете использовать код 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);
}

Примечания:

• Я возвращаю List вместо DataFrame, поскольку dtot является скалярным значением.
• Приведенный выше код предназначен для иллюстрации идеи.Скорее всего он не будет работать без настроек!