7

I'm working on a program in R to calculate the Gabriel Graph for up to 1000 data points. I used a program I found online first (GabrielGraph based on Bhattacharya et al. 1981 lines 781-830).

Unfortunately it takes quite a bit of time to get the result so I tried reprogramming it using Rcpp. For this I wrote a couple of small programs and a big one called edges which is used to calculate the edges of the Gabriel Graph. I'm also new to programming with Rcpp so I probably did everything more complicated than necessary but I didn't know how to do it any better.

#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;
}

I used these Rcpp programs in this R program:

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)))
}

For a sample data of ten data points it worked fine, for example:

z <- 10
X <- runif(z)*100
Y <- runif(z)*100
GabrielGraphMatrix(X,Y)

returns

> 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

But if I try to put in bigger data sets I get this error message:

Error: Value of SET_STRING_ELT() must be a 'CHARSXP' not a 'builtin'

I would be amazingly grateful if anybody had at least an idea of what I did wrong.

rgum
  • 153
  • 1
  • 1
  • 7
  • 3
    The error means that you are passing a function (not a character vector) to `SET_STRING_ELT`. But that's a lot of code to dig through. Try to narrow it down to a *minimal* example. http://stackoverflow.com/q/5963269/134830 – Richie Cotton Jan 06 '15 at 10:53
  • Thank you for answering so quickly! I know it's a lot of code and I'm sorry but I really don't know where the mistake might be, so I wouldn't know how to minimize the code :/ I'll give it a try though. – rgum Jan 06 '15 at 11:05
  • When I try to run `GabrielGraphMatrix(X,Y)`, I get `Error in Gab(zeros(npts), edges1, edges2, anz) : could not find function "zeros"`. Where is this function? – Richie Cotton Jan 06 '15 at 11:32
  • You are right, that was a function I had programmed somewhere else. It is just supposed to create a matrix filled with zeros. I edited my code and replaced zeros(npts) with matrix(0, npts, npts). Thank you again for your help, it is very much appreciated. – rgum Jan 06 '15 at 12:08
  • I've been trying out some things and now I got the error message `Value of SET_STRING_ELT() must be a 'CHARSXP' not a 'pairlist'` is that a different problem now or the same? @Richie – rgum Jan 08 '15 at 10:01

2 Answers2

4

Just in case anybody has the same problem. Mine was solved quite easily in the end. The mistake was in the function

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

The for-loop was too long. It should have been for(int j=0; j<col; j++)instead of for(int j=0; j<=col; j++).

rgum
  • 153
  • 1
  • 1
  • 7
1

I couldn't reproduce your error, but it threw a variety of similar ones, and often made R crash. Here are a couple of obvious problems.


In your C++ function Gab you have at least two problems:

  1. You don't define the variable anz before you use it.
  2. You are using round rather than square bracket to index Gabriel.

This

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

should be

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

In your R function GabrielGraphMatrix you are growing edges1 and edges2 in a loop. This means that they have to be reallocated in every iteration of the for loop. This will cause problems once you get above trivial loop lengths.

Instead, preallocate them as lists, then call unlist afterwards to get the vector you want.

# 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)
Richie Cotton
  • 118,240
  • 47
  • 247
  • 360
  • Thank you again for your help. I tried changing the function Gab but unfortunately that didn't work for me. The matrix didn't look the way it was supposed to. I did change the edges1/edges2 though, thank you for bringing that to my attention. I still don't know what the problem is but I'll keep trying. Thank you again! – rgum Jan 08 '15 at 12:12