2

writing a function for columnwise computation of sample skewness in Rcpp I'm in trouble using sqrt()-function. I know that sqrt(x) works for NumericVector types (tested it in a separate file), however in my code (where I'm trying to pass doubles) it does not work.

Here's my code:

#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector colSkew(NumericMatrix x) {
    int nc = x.ncol();
    int nr = x.nrow();
    NumericVector colS(nc);
    for(int i = 0; i < nc; i++){
        double cMean = mean(x( _ ,i));
        double xSq = 0;
        double cSt = 0;
          for(int j = 0; j < nr; j++){
              xSq += std::pow(x(j,i), 2.0);
              cSt += std::pow(x(j,i) - cMean, 3.0);
          }
          double colMsq = nr * std::pow(cMean, 2.0);
          double cTT = sqrt((xSq - colMsq)) / (nr - 1);
          double colVar = cTT / (nr - 1);
          double colNew = nr * std::pow(colVar, 3);
          colS[i] = cSt / colNew;
      }
    return(colS);
}

I tried std::sqrt() instead, giving me the "call to sqrt is ambiguous" error. For that one, neither Getting: error C2668: 'sqrt' : ambiguous call to overloaded function, nor http://dirk.eddelbuettel.com/code/rcpp.examples.html's

inline static double sqrt_double( double x ){ return ::sqrt( x ); }

helped. (The latter one has itself either the "no matching function" or the "ambiguous" problem). The code itself compiles, but does not work properly (I can call the function but results are not fine).

Reproducible example:

# First load the cpp function, however you want to;
Create example data:
A = matrix(rchisq(1000, 5), nrow = 100)
library(timeDate)
skew.1 = apply(A, 2, skewness)
skew.2 = colSkew(A)

# A custom R-function which should mimic the cpp function
colSkew.r = function(x){
  
  nc = ncol(x)
  nr = nrow(x)
  colS = numeric(nc)
  cMean = colMeans(x)
  xSq = colSums(x^2)
  cSt = 0
  for(i in 1:nc){
      cSt = sum((x[,i]-cMean[i])^3)
      colMsq = nr * cMean[i]^2
      cTT = sqrt((xSq[i] - colMsq) / (nr - 1))
      colNew = nr * cTT^3
      colS[i] = cSt / colNew
  }
  return(colS)
}
skew.3 = colSkew.r(A)
Erin Sprünken
  • 400
  • 2
  • 13
  • 1
    Does the compiler list the overloads of `sqrt` that it thinks could match your call, that it can't disambiguate? Does "Rcpp.h" include a `sqrt` function? – 1201ProgramAlarm Dec 26 '20 at 17:02
  • If the problem is not compilation, can you add a small example data to use the function on, as well as your current and desired output for that small example? Thanks! See [mcve] for details. I suspect changing `(nr - 1)` to `(nr - 1.0)` will help, but would need an example to run it on to be sure – duckmayr Dec 26 '20 at 17:03
  • 2
    The only possible answer here is that Rcpp has a version of sqrt that is causing an overload resolution issue. Preprocess the file and look for all declarations of `sqrt` to find out where it's coming from. Removing the `using namespace Rcpp;` may also help. – Tumbleweed53 Dec 26 '20 at 17:05
  • or don't remove it and do `std::sqrt` – Stéphane Laurent Dec 26 '20 at 17:24
  • Agrree with @Tumbleweed53 here and just suggested that in my answer. Doing `using namespace Rcpp;` is almost never a good idea for more complicated code. Ditto for `using namespace std;`. Being explicit is generally preferable (though not required). – Dirk Eddelbuettel Dec 26 '20 at 17:27
  • Thank you all for your fast comments. @duckayr, I will add a MRE and an R-function which does (correctly) what the Rcpp-function should. Switching Namespaces does not solve my problem, unfortunately. As mentioned, std::sqrt only gives me the ambiguous problem and making (nr - 1.0) does not solve the problem; As I said, the function also runs (I don't know why), but results are wrong. Possibly something else gives me the wrong results, but having an error like this makes me skeptical. – Erin Sprünken Dec 26 '20 at 17:34
  • 1
    Yeah guys, thank you so much for your work, the mistake was a goddamn bracket set wrong. :( – Erin Sprünken Dec 26 '20 at 17:41

1 Answers1

4

To be brief, there are few recurrent and common beginner errors we can avoid:

  • do not include unneeded headers (mostly harmless, but Rcpp already brings math headers)

  • do not use a globally flattened namespace Rcpp especially after you get compilation errors on symbol visibility

With that and the required minimal changes it builds and runs (and returns something non-sensical but I leave that for you to tackle :)

There are other ways to sort this out as std:: symbols and Rcpp symbols generally coexist quite happily with distinct signatures.

Code

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector colSkew(Rcpp::NumericMatrix x) {
    int nc = x.ncol();
    int nr = x.nrow();
    Rcpp::NumericVector colS(nc);
    for(int i = 0; i < nc; i++){
      double cMean = Rcpp::mean(x( Rcpp::_ ,i));
        double xSq = 0;
        double cSt = 0;
          for(int j = 0; j < nr; j++){
              xSq += std::pow(x(j,i), 2.0);
              cSt += std::pow(x(j,i) - cMean, 3.0);
          }
          double colMsq = nr * std::pow(cMean, 2.0);
          double cTT = std::sqrt((xSq - colMsq)) / (nr - 1);
          double colVar = cTT / (nr - 1);
          double colNew = nr * std::pow(colVar, 3);
          colS[i] = cSt / colNew;
      }
    return(colS);
}

/*** R
set.seed(42)
colSkew( matrix(rnorm(100), 100, 1) )
*/

Output

> sourceCpp("answer.cpp")

> set.seed(42)

> colSkew( matrix(rnorm(100), 100, 1) )
[1] -448250877
>
Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725