我正在研究R中的程序以計算最多1000個數據點的Gabriel圖。我使用了一個我在網上找到的程序(GabrielGraph based on Bhattacharya et al. 1981第781-830行)。錯誤:SET_STRING_ELT()的值必須是'CHARSXP'而不是'內建'?
不幸的是,得到結果需要相當多的時間,所以我試圖用Rcpp重新編程。爲此,我編寫了幾個小程序和一個叫做Gabriel Graph的邊緣的邊緣圖。我也是使用Rcpp進行編程的新手,所以我可能做了比必要更復雜的任何事情,但我不知道如何做得更好。
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double vecnorm(NumericVector x){
//to calculate the vectornorm sqrt(sum of (vector entries)^2)
double out;
out = sqrt(sum(pow(x,2.0)));
return out;
}
// [[Rcpp::export]]
NumericVector vektorzugriff(NumericMatrix xy,int i){
//to return a row of the Matrix xy
int col = xy.ncol();
NumericVector out(col);
for(int j=0; j<=col; j++){
out[j] = xy(i-1,j);
}
return out;
}
// [[Rcpp::export]]
IntegerVector vergl(NumericVector eins, NumericVector zwei){
//to see if two Vectors have any identical entries
IntegerVector out = match(eins, zwei);
return out;
}
// [[Rcpp::export]]
IntegerVector verglInt(int eins, NumericVector zwei){
NumericVector dummy = NumericVector::create(eins) ;
IntegerVector out = match(dummy, zwei);
return out;
}
// [[Rcpp::export]]
NumericVector toVec(NumericVector excluded, int k){
//to append int k to a Vector excluded
NumericVector dummy = NumericVector::create(k) ;
int len = excluded.size();
int len2 = dummy.size();
int i=0;
NumericVector out(len+len2);
while(i<len+len2){
if(i<len){
out[i]=excluded[i];
i++;
}
else{
out[i]=dummy[i-len];
i++;
}
}
return out;
}
// [[Rcpp::export]]
LogicalVector isNA(IntegerVector x) {
//to see which Vector Entries are NAs
int n = x.size();
LogicalVector out(n);
for (int i = 0; i < n; ++i) {
out[i] = IntegerVector::is_na(x[i]);
}
return out;
}
// [[Rcpp::export]]
NumericMatrix Gab(NumericMatrix Gabriel, NumericVector edges1, NumericVector edges2, int anz){
//to fill a Matrix with the Gabrieledges
for(int i=0; i<anz; i++) {
Gabriel(edges1[i]-1, edges2[i]-1) = 1 ;
Gabriel(edges2[i]-1, edges1[i]-1) = 1 ;
}
return Gabriel;
}
// [[Rcpp::export]]
NumericVector edges(NumericMatrix xy,NumericVector vertices,NumericVector excluded, int i){
//actual function to calculate the edges of the GabrielGraph
int npts = xy.nrow()+1;
double d1;
double d2;
double d3;
for(int r=i+1; r<npts; r++) {
// Skip vertices in excluded
if(!is_true(any(isNA(verglInt(r,excluded))))){
continue;}
d1 = vecnorm(vektorzugriff(xy,i) - vektorzugriff(xy,r));
for(int k=1; k<npts; k++) {
if((k!=r) && (k!=i)){
d2 = vecnorm(vektorzugriff(xy,i) - vektorzugriff(xy,k));
d3 = vecnorm(vektorzugriff(xy,r) - vektorzugriff(xy,k));
//Betrachte vertices, die noch nicht excluded sind
if(!is_true(any(isNA(verglInt(k,vertices[isNA(vergl(vertices,excluded))]))))){
//Wenn d(x,z)^2 > d(x,y)^2+d(y,z)^2 -> Kante gehoert nicht zum GG
if(pow(d2,2.0) > pow(d1,2.0) + pow(d3,2.0)) {
excluded = toVec(excluded,k);
}
}
if(pow(d1,2.0) > pow(d2,2.0) + pow(d3,2.0)){
excluded = toVec(excluded,r);
break;
}
}
}
}
return excluded;
}
我用在該R節目這些RCPP節目:
GabrielGraphMatrix <- function(X,Y,PlotIt=FALSE){
# Heuristic rejection Algorithm for Gabriel Graph Construction (Bhattacharya et al. 1981)
# Algorithm is ~ O(d n^2)
#loading Rcpp functions
library(Rcpp)
sourceCpp("... .cpp")
XY <- cbind(X,Y)
ndim <- ncol(XY)
npts <- nrow(XY)
edges1<- c()
edges2<- c()
for(i in 1:(npts-1)) {
# Candidate set of Gabriel neighbors
vertices <- (i+1):npts
# Initialize list of vertices to be excluded from Ni
excluded <- edges(XY,vertices,vector(),i);
adj <- vertices[which(!match(vertices,excluded,nomatch=F)>0)]
if(length(adj) > 0) {
edges1=c(edges1,rep(i,length(adj)))
edges2=c(edges2,adj)
}
}
anz <- length(edges1)
Gabriel <- Gab(matrix(0, npts, npts),edges1,edges2,anz)
return(list(Gabriel=Gabriel,edges=cbind(edges1,edges2)))
}
對於十個數據的樣本數據點它工作得很好,例如:
z <- 10
X <- runif(z)*100
Y <- runif(z)*100
GabrielGraphMatrix(X,Y)
返回
> GabrielGraphMatrix(X,Y)
$Gabriel
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 1 0 0 0 0 0 0 0 0
[2,] 1 0 0 1 0 0 1 0 0 0
[3,] 0 0 0 1 1 0 0 0 0 1
[4,] 0 1 1 0 0 0 0 0 0 0
[5,] 0 0 1 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 1 0 0
[7,] 0 1 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 1 0 0 1 1
[9,] 0 0 0 0 0 0 0 1 0 1
[10,] 0 0 1 0 0 0 0 1 1 0
$edges
edges1 edges2
[1,] 1 2
[2,] 2 4
[3,] 2 7
[4,] 3 4
[5,] 3 5
[6,] 3 10
[7,] 6 8
[8,] 8 9
[9,] 8 10
[10,] 9 10
但是,如果我嘗試投入更大的數據集,我收到此錯誤信息:
Error: Value of SET_STRING_ELT() must be a 'CHARSXP' not a 'builtin'
,我將不勝感激驚人地如果有人有至少我做錯了什麼的想法。
該錯誤表示您正在將一個函數(而不是一個字符向量)傳遞給'SET_STRING_ELT'。但這是很多需要深入研究的代碼。嘗試縮小到一個*最小*的例子。 http://stackoverflow.com/q/5963269/134830 –
謝謝你回答這麼快!我知道這是很多代碼,我很抱歉,但我真的不知道錯誤的可能性,所以我不知道如何最小化代碼:我會嘗試一下。 – rgum
當我嘗試運行'GabrielGraphMatrix(X,Y)'時,我得到了'錯誤在Gab(零點(npts),邊緣1,邊緣2,anz):找不到函數「零」。這個功能在哪裏? –