2017-12-31 52 views
0

Ich habe mit dem folgenden C++ - Code für die Integration mit R-Code, den ich geschrieben habe (zu viel, um hier aufzunehmen), aber immer einen Fehler bekommen, dass der Cube: : Operator() Index ist außerhalb der Grenzen und ich bin unsicher, warum dies auftritt. Mein Verdacht ist, dass der 3D-Array ist nicht richtig, wie beschrieben inRcppAramadillo Cube :: operator(): index außerhalb der Grenzen

gefüllt wird

making 3d array with arma::cube in Rcpp shows cube error

aber ich bin mir nicht sicher, wie man richtig das Problem zu lösen.

Unten ist mein voller C++ Code:

// [[Rcpp::depends(RcppArmadillo)]] 
#define ARMA_DONT_PRINT_OPENMP_WARNING 
#include <RcppArmadillo.h> 
#include <RcppArmadilloExtensions/sample.h> 
#include <set> 
using namespace Rcpp; 


int sample_one(int n) { 
    return n * unif_rand(); 
} 

int sample_n_distinct(const IntegerVector& x, 
        int k, 
        const int * pop_ptr) { 

    IntegerVector ind_index = RcppArmadillo::sample(x, k, false); 
    std::set<int> distinct_container; 

    for (int i = 0; i < k; i++) { 
    distinct_container.insert(pop_ptr[ind_index[i]]); 
    } 

    return distinct_container.size(); 
} 

// [[Rcpp::export]] 
arma::Cube<int> fillCube(const arma::Cube<int>& pop, 
        const IntegerVector& specs, 
        int perms, 
        int K) { 

int num_specs = specs.size(); 
arma::Cube<int> res(perms, num_specs, K); 

IntegerVector specs_C = specs - 1; 
const int * pop_ptr; 
int i, j, k; 

for (i = 0; i < K; i++) { 
    for (k = 0; k < num_specs; k++) { 
     for (j = 0; j < perms; j++) { 
      pop_ptr = &(pop(0, sample_one(perms), sample_one(K))); 
      res(j, k, i) = sample_n_distinct(specs_C, k + 1, pop_ptr); 
     } 
    } 
} 

return res; 
} 

Hat jemand eine Idee, wie haben, was kann den besagten Fehler produzieren?

Unten ist der R-Code mit einem Aufruf der C++ - Funktion (einschließlich einer auskommentierten Triple-geschachtelte 'für' Schleife, die der C++ - Code reproduziert).

## Set up container(s) to hold the identity of each individual from each permutation ## 

num.specs <- ceiling(N/K) 

## Create an ID for each haplotype ## 

haps <- 1:Hstar 

## Assign individuals (N) to each subpopulation (K) ## 

specs <- 1:num.specs 

## Generate permutations, assume each permutation has N individuals, and sample those individuals' haplotypes from the probabilities ## 

gen.perms <- function() { 
    sample(haps, size = num.specs, replace = TRUE, prob = probs) 
} 

pop <- array(dim = c(perms, num.specs, K)) 

for (i in 1:K) { 
    pop[,, i] <- replicate(perms, gen.perms()) 
} 

## Make a matrix to hold individuals from each permutation ## 

# HAC.mat <- array(dim = c(perms, num.specs, K)) 

## Perform haplotype accumulation ## 

# for (k in specs) { 
    # for (j in 1:perms) { 
     # for (i in 1:K) { 
      # select.perm <- sample(1:nrow(pop), size = 1, replace = TRUE) # randomly sample a permutation 
      # ind.index <- sample(specs, size = k, replace = FALSE) # randomly sample individuals 
      # select.subpop <- sample(i, size = 1, replace = TRUE) # randomly sample a subpopulation 
      # hap.plot <- pop[select.perm, ind.index, select.subpop] # extract data 
      # HAC.mat[j, k, i] <- length(unique(hap.plot)) # how many haplotypes are recovered 
     # } 
    # } 
# } 

HAC.mat <- fillCube(pop, specs, perms, K) 
+2

Könnten Sie Beispiel _R_ Code hinzufügen, der den _C++ _ Code aufruft, um diesen Fehler auszulösen? – coatless

+1

Sie haben fast sicher einen logischen Fehler in Ihrer Indizierung. Ich würde alles einfacher machen - beginnen Sie einfach mit einem Würfel von 5 x 6 x 7 und stellen Sie sicher, dass das, was Sie über die Indizes denken, tatsächlich gilt. –

+0

@coatless Siehe die aktualisierte Frage für meinen R-Code. –

Antwort

1

Dies ist ein Fehler außerhalb der Grenzen. Der Kern des Problems ist der Aufruf

pop_ptr = &(pop(0, sample_one(perms), sample_one(K))); 

seit

sample_one(perms) 

als Zugriffsindex platziert wird, wo die maximale Länge num_specs ist. Dies zeigt sich durch, wie res definiert:

arma::Cube<int> res(perms, num_specs, K); 

So heraus perms aus num_specs Ort bewegen sollte das Problem beheben.

// [[Rcpp::export]] 
arma::Cube<int> fillCube(const arma::Cube<int>& pop, 
         const IntegerVector& specs, 
         int perms, 
         int K) { 

    int num_specs = specs.size(); 
    arma::Cube<int> res(perms, num_specs, K); 

    IntegerVector specs_C = specs - 1; 
    const int * pop_ptr; 
    int i, j, k; 

    for (i = 0; i < K; i++) { 
    for (k = 0; k < num_specs; k++) { 
     for (j = 0; j < perms; j++) { 

     // swapped location 
     pop_ptr = &(pop(sample_one(perms), 0, sample_one(K))); 
     // should the middle index be 0? 
     res(j, k, i) = sample_n_distinct(specs_C, k + 1, pop_ptr); 
     } 
    } 
    } 

    return res; 
} 
+0

Nizza. Und genau wie wir es vermutet haben. –

+0

Vielen Dank !! Klappt wunderbar! Ich bin völlig beeindruckt, wie viel Rcpp/RcppArmadillo meinen Algorithmus beschleunigt hat! Es ist wirklich ein tolles Paket! –

Verwandte Themen