2015-01-06 360 views
3

我正在研究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' 

,我將不勝感激驚人地如果有人有至少我做錯了什麼的想法。

+3

該錯誤表示您正在將一個函數(而不是一個字符向量)傳遞給'SET_STRING_ELT'。但這是很多需要深入研究的代碼。嘗試縮小到一個*最小*的例子。 http://stackoverflow.com/q/5963269/134830 –

+0

謝謝你回答這麼快!我知道這是很多代碼,我很抱歉,但我真的不知道錯誤的可能性,所以我不知道如何最小化代碼:我會嘗試一下。 – rgum

+0

當我嘗試運行'GabrielGraphMatrix(X,Y)'時,我得到了'錯誤在Gab(零點(npts),邊緣1,邊緣2,anz):找不到函數「零」。這個功能在哪裏? –

回答

1

我無法重現您的錯誤,但它拋出了各種類似的錯誤,並且經常使R崩潰。這裏有幾個明顯的問題。


在你的C++函數Gab至少有兩個問題:

  1. 您使用它之前,您不定義變量anz
  2. 您正在使用圓括號而非方括號來索引Gabriel

Gabriel(edges1[i]-1, edges2[i]-1) 

應該

Gabriel[edges1[i]-1, edges2[i]-1] 

在你一個R函數GabrielGraphMatrix你是在一個循環中不斷增長的edges1edges2。這意味着它們必須在for循環的每次迭代中重新分配。這會導致問題,一旦你得到以上瑣碎的循環長度。

取而代之的是,將它們預先分配爲列表,然後調用unlist以獲得所需的矢量。

# before the loop 
edges1 <- vector("list", npts - 1) 
edges2 <- vector("list", npts - 1) 

# in the loop 
if(length(adj) > 0) { 
    edges1[[i]] <- rep(i,length(adj)) 
    edges2[[i]] <- adj 
} 

# after the loop 
edges1 <- unlist(edges1) 
edges2 <- unlist(edges2) 
+0

再次感謝您的幫助。我試圖改變功能Gab,但不幸的是,這並不適合我。矩陣看起來並不像它應該的樣子。 雖然我確實改變了邊緣1 /邊緣2,但感謝您將它引入我的視線。 我仍然不知道問題是什麼,但我會繼續嘗試。再次感謝你! – rgum

2

只是爲了防止任何人有同樣的問題。最終我的礦難很容易解決。功能上的錯誤是

// [[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; 
} 

for循環太長。它應該是for(int j=0; j<col; j++)而不是for(int j=0; j<=col; j++)

相關問題