2017-11-10 22 views
2

클러스터링 계수 Kaiser입니다. 원본 코드는 MATLAB에 있으며 저자의 website (마지막 링크, 페이지 하단)에서 사용할 수 있습니다.카이저 Rcpp의 클러스터링 계수

먼저 R으로 원래의 MATLAB 함수를 구현합니다. 아래 코드는 (kaiser_r.R) 붙여 넣습니다.

> A <- matrix(c(0,1,1,0,1,0,1,1,1,1,0,0,0,1,0,0), 4, 4) 
> cc_kaiser(A) 
[1] 0.7777778 

결과는 (내가 MATLAB과 테스트) 올바른 :

# R equivalent of MATLAB find function 
# Find indices of nonzero elements in a vector 
rfind <- function(adj) seq(along = adj)[adj != 0] 

# Main function 
cc_kaiser <- function(adj) { 
    n_count <- nrow(adj) 
    w <- rep(0, n_count) 
    # Number of nodes with at least two neighbors 
    n_neigh = 0 
    for (i in 1:n_count) { 
    n <- rfind(adj[i, ] + t(adj[, i])) 
    n_e <- 0 
    l_n <- length(n) 
    for (j in 1:l_n) { 
     vec <- t(as.matrix(adj[n[j], ])) 
     n_v <- rfind(vec) 
     n_e <- n_e + l_n + length(n_v) - length(union(n, n_v)) 
    } 
    if (l_n > 1) { 
     w[i] <- n_e/(l_n * (l_n - 1)) 
     n_neigh <- n_neigh + 1 
    } 
    } 
    cl <- sum(w)/n_neigh 
    return(cl) 
} 

나는 함께이 기능을 테스트합니다. 그런 다음 Rcpp과 동일한 기능을 구현하려고합니다.

> kaiser(A) 
[1] 0.6666667 

내가 친절하게 도움을 요청 : 나는 kaiser_rcpp.cpp을 실행하면

#include <RcppArmadillo.h> 
using namespace Rcpp; 

// [[Rcpp::depends(RcppArmadillo)]] 
// [[Rcpp::export]] 
double kaiser(arma::mat A) { 
    int n_count = A.n_rows; 
    std::vector<int> w(n_count); 
    int n_neigh = 0; 

    for(int i = 0; i < n_count; i++) { 
    arma::rowvec bla = A.row(i) + A.col(i).t(); 
    arma::uvec n = unique(find(bla > 0)); 
    int n_e = 0; 
    int l_n = n.n_elem; 
    for(int j = 0; j < l_n; j++) { 
     arma::colvec vec = A.row(n(j)).t(); 
     arma::uvec n_v = unique(find(vec > 0)); 
     IntegerVector uni = union_(as<IntegerVector>(wrap(n)), as<IntegerVector>(wrap(n_v))); 
     n_e = n_e + l_n + n_v.n_elem - uni.size(); 
    } 
    if(l_n > 1) { 
     w[i] = n_e/(l_n * (l_n - 1)); 
     n_neigh = n_neigh + 1; 
    } 
    } 
    double s = std::accumulate(w.begin(), w.end(), 0.0); 
    double cl = s/n_neigh; 

    return(cl); 
} 

내가 다른 값을 얻을 : 여기 내 시도 (kaiser_rcpp.cpp)입니다. 내 Rcpp 코드에서 실수가 어디인지 전혀 알 수 없습니다.

답변

2

처음 w은 정수가 아닌 복식의 벡터입니다.

그러면 w[i] = n_e/(l_n * (l_n - 1)); 행이 잘못되었습니다. w[i] = (double)n_e/(l_n * (l_n - 1));으로 교체해야합니다. n_e(l_n * (l_n - 1))은 모두 정수이므로 정수 나누기 (예 : 3/2 = 1)를 수행합니다.

전체 코드 :

#include <RcppArmadillo.h> 
using namespace Rcpp; 

// [[Rcpp::depends(RcppArmadillo)]] 
// [[Rcpp::export]] 
double kaiser(arma::mat A) { 
    int n_count = A.n_rows; 
    std::vector<double> w(n_count);      // CHANGE HERE 
    int n_neigh = 0; 

    for(int i = 0; i < n_count; i++) { 
    arma::rowvec bla = A.row(i) + A.col(i).t(); 
    arma::uvec n = unique(find(bla > 0)); 
    int n_e = 0; 
    int l_n = n.n_elem; 
    for(int j = 0; j < l_n; j++) { 
     arma::colvec vec = A.row(n(j)).t(); 
     arma::uvec n_v = unique(find(vec > 0)); 
     IntegerVector uni = union_(as<IntegerVector>(wrap(n)), as<IntegerVector>(wrap(n_v))); 
     n_e = n_e + l_n + n_v.n_elem - uni.size(); 
     // Rcout << n_e << std::endl; 
    } 
    if(l_n > 1) { 
     w[i] = (double)n_e/(l_n * (l_n - 1));   // CHANGE HERE 
     n_neigh++;          // (CHANGE HERE) 
    } 
    } 
    double s = std::accumulate(w.begin(), w.end(), 0.0); 
    double cl = s/n_neigh; 

    // Rcout << as<NumericVector>(wrap(w)) << std::endl; 

    return(cl); 
}