Is this correct in Rcpp?-CodePudding

I want to compare each columns, and return all the results after calculating. I try to write the codes, but the outcome was not resonable. Because if there are 5 columns in a matrix, the number of result will will be 5*4/2=10 rather than 5. I think the problem is the m in codes. I don't know whether it is correct. Thanks.

library(Rcpp)
sourceCpp(code='
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h> 

double KS(arma::colvec x, arma::colvec y) {
  int n = x.n_rows;
  arma::colvec w = join_cols(x, y);
  arma::uvec z = arma::sort_index(w);
  w.fill(-1); w.elem( find(z <= n-1) ).ones();
  return max(abs(cumsum(w)))/n;
}
// [[Rcpp::export]]
Rcpp::NumericVector K_S(arma::mat mt) {
  int n = mt.n_cols;
  int m = 1;
  Rcpp::NumericVector results(n);
  for (int i = 0; i < n-1; i  ) {
    for (int j = i 1; j < n; j  ){
      arma::colvec x=mt.col(i);
      arma::colvec y=mt.col(j);
      results[m] = KS(x, y);
      m   ;
    }
  }
  return results;
}
')

set.seed(1)
mt <- matrix(rnorm(400*5), ncol=5)
result <- K_S(t(mt))

> result
[1] 0.0000 0.1050 0.0675 0.0475 0.0650

CodePudding user response：

You had a couple of small errors. In fixing it, an intermediate version I had just filled a similar n by n matrix -- that made indexing errors obvious. Returning an arma::rowvec also helps with possible out-of-bounds index errors (it errors by default) but lastly you (in this case !!) can actually just grow a std::vector instead.

Code

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

double KS(arma::colvec x, arma::colvec y) {
  int n = x.n_rows;
  arma::colvec w = join_cols(x, y);
  arma::uvec z = arma::sort_index(w);
  w.fill(-1); w.elem( find(z <= n-1) ).ones();
  return max(abs(cumsum(w)))/n;
}
// [[Rcpp::export]]
std::vector<double> K_S(arma::mat mt) {
  int n = mt.n_cols;
  std::vector<double> res;
  for (int i = 0; i < n; i  ) {
    for (int j = i 1; j < n; j  ){
      arma::colvec x=mt.col(i);
      arma::colvec y=mt.col(j);
      res.push_back(KS(x, y));
    }
  }
  return res;
}

/*** R
set.seed(1)
mt <- matrix(rnorm(400*5), ncol=5)
result <- K_S(mt)
result
*/

Output

> Rcpp::sourceCpp("~/git/stackoverflow/73916783/answer.cpp")
> set.seed(1)
> mt <- matrix(rnorm(400*5), ncol=5)
> result <- K_S(mt)
> result
 [1] 0.1050 0.0675 0.0475 0.0650 0.0500 0.0775 0.0575 0.0500 0.0475 0.0600
>