CODEWITHSUNDEEP
rinlinercpp
Adrian Albert
Adrian Albert

Reputation: 3

NumericMatrix returns a matrix of all zeros?

I'm trying to implement a simple code in Rcpp that calculates and populates the entries of a distance matrix. The problem is that the Rcpp code (below) returns a matrix D with all elements having a value of zero. This issue does not seem to be addressed anywhere in the forums - I'd appreciate some advice!

  src_d_err_c <- '
  using namespace Rcpp;
  double d_err_c(NumericVector cx, NumericVector csx, NumericVector cy, NumericVector csy) {
    using namespace Rcpp;    
    NumericVector d = (cx - cy)*(cx - cy) / csx;
    double s = std::accumulate(d.begin(), d.end(), 0.0);

    return s;
  }'

  src_d_mat = '   
  using namespace Rcpp;
  // input

  Rcpp::NumericMatrix cX(X);
  Rcpp::NumericMatrix cY(Y);
  Rcpp::NumericMatrix cSX(SX);
  Rcpp::NumericMatrix cSY(SY);  
  int N1 = cX.nrow();
  int N2 = cY.nrow();
  NumericMatrix D(N1, N2);
  NumericVector v(N1);

  for (int x = 0; x++; x<N1){
    v[x] = x;
    for (int y = 0; y++; y<N2) {
      D(x,y) = d_err_c(cX(x,_), cSX(x,_), cY(y,_), cSY(y,_));
    };
  };
  return wrap(v);'  
  fun <- cxxfunction(signature(X = "numeric",  SX = "numeric", 
                               Y = "numeric",  SY = "numeric"),
                     body = src_d_mat, includes = src_d_err_c, 
                     plugin = "Rcpp")

Upvotes: 0

Views: 296

Answers (2)

Dirk is no longer here
Dirk is no longer here

Reputation: 368181

@Vincent correctly pointed out one main error (not actually looping), but there is another major one: you return v when your computations go into D which is what you meant to return. (You also do not need v at all, actually.)

Here is a repaired version, which uses correct loop indexing, returns the correct object, and omits the unused code. It also switched to Rcpp Attributes.

Save this in a file, say "distmat.cpp", and use sourceCpp("distmat.cpp") after which you have a new function d_mat you can call.

#include <Rcpp.h>

using namespace Rcpp;

double d_err_c(NumericVector cx, NumericVector csx, 
               NumericVector cy, NumericVector csy) {
    NumericVector d = (cx - cy)*(cx - cy) / csx;
    return std::accumulate(d.begin(), d.end(), 0.0);
}

// [[Rcpp::export]]
NumericMatrix d_mat(NumericMatrix cX, NumericMatrix cY, 
                    NumericMatrix cSX, NumericMatrix cSY) {
    int N1 = cX.nrow();
    int N2 = cY.nrow();
    NumericMatrix D(N1, N2);

    for (int x = 0; x<N1; x++) {
        for (int y = 0; y<N2; y++) {
            D(x,y) = d_err_c(cX(x,_), cSX(x,_), cY(y,_), cSY(y,_));
        }
    }
    return D;
}

Upvotes: 1

Vincent Zoonekynd
Vincent Zoonekynd

Reputation: 32351

The arguments of your for loops are in the wrong order: the condition should be in the middle and the increment at the end. 

for (int x = 0; x < N1; x++)

Upvotes: 4

Related Questions