Search code examples
rcpp

Ignore missing values while finding min in Rcpp


I am trying to find the minimum of 2 values from 2 vectors in Rcpp. But the following does not compile:

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



// [[Rcpp::export]]
NumericVector timesTwo(int time_length, double BXadd, 
                       NumericVector vn_complete, NumericVector vn1_complete) {
  
  // Empty vectors
  NumericVector BX (time_length);
  
  for(int t = 0; t < time_length; t++) {

    BX[t] = BXadd * sqrt(std::min(na_omit(vn_complete[t], vn1_complete[t])));
    
  }

  return BX;
  
  // return vn_complete[0];
}

  Error 1 occurred building shared library.

It works if I don't use na_omit.

R code for running the function:

Rcpp::sourceCpp("test.cpp")

timesTwo(5, 2, 5:9, 1:5)

Solution

  • What follows below is a slightly non-sensical answer (as it only works with vectors not containing NA) but it has all the component your code has, and it compiles.

    #include <Rcpp.h>
    using namespace Rcpp;
    
    // [[Rcpp::export]]
    NumericVector foo(int time_length, double BXadd,
                      NumericVector vn_complete, NumericVector vn1_complete) {
    
        // Empty vectors
        NumericVector BX (time_length);
        vn_complete = na_omit(vn_complete);
        vn1_complete = na_omit(vn1_complete);
        for(int t = 0; t < time_length; t++) {
            double a = vn_complete[t];
            double b = vn1_complete[t];
            BX[t] = BXadd * std::sqrt(std::min(a,b));
        }
        return BX;
    }
    
    // Edited version with new function added 
    // [[Rcpp::export]]
    NumericVector foo2(double BXadd, NumericVector vn_complete, 
                       NumericVector vn1_complete) {
        return BXadd * sqrt(pmin(vn_complete, vn1_complete));
    }
    
    /*** R
    foo(5, 2, 5:9, 1:5)
    foo2(5, 2, 5:9, 1:5)
    */
    

    For a real solution you will have to think harder about what the removal of NA is supposed to do as it will alter the lenth of your vectors too. So this still needs work.

    Lastly, your whole function can be written in R as

    2 * sqrt(pmin(5:9, 1:5))
    

    and I think you could write that same expression using Rcpp sugar too as we have pmin() and sqrt():

    // [[Rcpp::export]]
    NumericVector foo2(double BXadd, NumericVector vn_complete, 
                       NumericVector vn1_complete) {
        return BXadd * sqrt(pmin(vn_complete, vn1_complete));
    }