7

I would write a function (using Rcpp) that removes all the NA values from a R vector.

Before doing so, I did a little test function through Rcpp::cppFunction function.

library(inline)
cppFunction('
  Vector<INTSXP> na_test(const Vector<INTSXP>& x) {
    return setdiff(x, Vector<INTSXP>::create(::traits::get_na<INTSXP>()));
  }
')

That works this way:

na_test(c(1, NA, NA, 1, 2, NA))
# [1] 1 2

After that I tried to generalize this function through the C++ template mechanism.

So, in an external .cpp file (sourced through sourceCpp function), I've written:

template <int RTYPE>
Vector<RTYPE> na_test_template(const Vector<RTYPE>& x) {
    return setdiff(x, Vector<RTYPE>::create(::traits::get_na<RTYPE>()));
}

// [[Rcpp::export(na_test_cpp)]]
SEXP na_test(SEXP x) {
    switch(TYPEOF(x)) {
        case INTSXP:
            return na_test_template<INTSXP>(x);
        case REALSXP:
            return na_test_template<REALSXP>(x);
    }
    return R_NilValue;
}

This code compiles but behaves differently and I cannot explain why.

Infact:

na_test_cpp(c(1, NA, NA, 1, 2, NA))
# [1]  2 NA NA NA  1

Why the same function (apparently) behaves differently? What happens here?

  • 1
    Your example calls na_test_cpp, which your code does not show. – Dirk Eddelbuettel Apr 11 '13 at 16:33
  • 1
    this expression // [[Rcpp::export(na_test_cpp)]] does not need to export the function (with the desired name) ? – leodido Apr 11 '13 at 16:35
  • my bad, overlooked it – Dirk Eddelbuettel Apr 11 '13 at 16:40
7

Following up on your answer, I would use something like this as the template:

template <int RTYPE>
Vector<RTYPE> na_omit_template(const Vector<RTYPE>& x) {
    int n = x.size() ;
    int n_out = n - sum( is_na(x) ) ;

    Vector<RTYPE> out(n_out) ;
    for( int i=0, j=0; i<n; i++){
        if( Vector<RTYPE>::is_na( x[i] ) ) continue ;
        out[j++] = x[i];
    }
    return out ;
}

So the idea is to first calculate the length of the result, and then just use Rcpp vector classes instead of std::vector. This will lead to less copies of the data.


With the development version of Rcpp (svn revision >= 4308), it works for me for all types, and we can then use our RCPP_RETURN_VECTOR dispatch macro instead of writing the switch:

// [[Rcpp::export]]
SEXP na_omit( SEXP x ){
     RCPP_RETURN_VECTOR( na_omit_template, x ) ;   
}

na_omit has been included in Rcpp(svn revision >= 4309), with a few modifications, i.e. it can handle named vectors and arbitrary sugar expressions.

| improve this answer | |
  • Very elegant solution! ;) Any idea to make it to work with STRSXP and CPLXSXP vectors? – leodido Apr 14 '13 at 2:25
  • Works for me for STRSXP. The complex version does not work because there was a bug in Rcpp, which is now fixed in revision 4308. – Romain Francois Apr 14 '13 at 7:51
  • 1
    I have added na_omit to Rcpp, with a few modifications, e.g. it can handle named vectors and arbitrary sugar expressions. – Romain Francois Apr 16 '13 at 9:53
  • very good ! Since I wish to use the na_omit that you added to Rcpp: when the new relase? :) Or I can install it from your SVN? – leodido Apr 16 '13 at 11:08
  • 1
    Now new release is planned at the moment. We don't release for each new feature or bug fix. Just get the source from r-forge and build it yourself: svn checkout svn://svn.r-forge.r-project.org/svnroot/rcpp/pkg/Rcpp – Romain Francois Apr 16 '13 at 11:14
1

I've continued to investigate a solution about the templating problem (i.e. see @Sameer reply).

So I've written another function and now the templating mechanism works.

In an external .cpp file:

#include <Rcpp.h>
template <int RTYPE, class T>
Vector<RTYPE> na_omit_template(const Vector<RTYPE>& x) {
    typedef typename Vector<RTYPE>::iterator rvector_it;
    if (x.size() == 0) {
        return x;
    }
    std::vector<T> out;
    rvector_it it = x.begin();
    for (; it != x.end(); ++it) {
        if (!Vector<RTYPE>::is_na(*it)) {
            out.push_back(*it);
        }
    }
    return wrap(out);
}

// [[Rcpp::export(na_omit_cpp)]]
SEXP na_omit(SEXP x) {
    switch(TYPEOF(x)) {
        case INTSXP:
            return na_omit_template<INTSXP, int>(x);
        case REALSXP:                   
            return na_omit_template<REALSXP, double>(x);
        case LGLSXP:
            return na_omit_template<LGLSXP, bool>(x);
        case CPLXSXP:
            return na_omit_template<CPLXSXP, Rcomplex>(x);
        case RAWSXP:
            return na_omit_template<RAWSXP, Rbyte>(x);
        default:
            stop("unsupported data type");
    }
}

This function removes the NA values, that was my initial purpose.

Unfortunately, at the moment, it does not work for all types of vectors, as the R examples below show.

library(Rcpp)
sourceCpp('file.cpp')
na_omit_cpp(as.integer(c(1, NA, NA, 1, 2, NA))) # OK
# [1] 1 1 2
na_omit_cpp(as.numeric(c(1, NA, NA, 1, 2, NA)))
# [1] 1 1 2
na_omit_cpp(c(NA, 1L, NA, 3L, NA)) # OK
# [1] 1 3
na_omit_cpp(c(NA, 2L, 1, NA)) # OK
# [1] 2 1
na_omit_cpp(c(1.0, 1.1, 2.2, NA, 3, NA, 4)) # OK
# [1] 1.0 1.1 2.2 3.0 4.0
na_omit_cpp(c(1L, NaN, NaN, 0, NA)) # OK
# [1]   1 NaN NaN   0
na_omit_cpp(c(NA, NaN, 1.0, 0.0, 2.2, NA, 3.3, NA, 4.4)) # OK
# [1] NaN 1.0 0.0 2.2 3.3 4.4
na_omit_cpp(as.logical(c(1, 0, 1, NA))) # OK
# [1]  TRUE FALSE  TRUE
na_omit_cpp(as.logical(c(TRUE, FALSE, NA, TRUE, NA))) # OK
# [1]  TRUE FALSE  TRUE
# empty vectors ?
na_omit_cpp(c(NA)) # OK
# logical(0)
na_omit_cpp(numeric(0)) # OK
# numeric(0)
na_omit_cpp(logical(0)) # OK
# logical(0)
na_omit_cpp(raw(0)) # OK
# raw(0)
na_omit_cpp(as.raw(c(40,16,NA,0,2))) # NO! (R converts it to 00)
# [1] 28 10 00 00 02
# Warning message ...
na_omit_cpp(as.complex(c(-1, 2, 1, NA, 0, NA, -1))) # NO! 
# [1] -1+0i  2+0i  1+0i    NA  0+0i    NA -1+0i

So, this function works in almost all cases except for raw vectors and complex vectors.

The current open issues are:

  1. I'm not sure why this bug and I'd like to find out why. Any idea?
  2. The previous templated function, as showed @Sameer, have a strange behaviour.
  3. How to do to accept character vectors?

I thought clearly to case STRSXP: return na_omit_template<STRSXP, ?>(x); but this statement doesn't work substituting std::string, Rcpp:String to ?.

| improve this answer | |
  • I will make an actual reply later. There is no such thing as NA in the case of raw vectors. – Romain Francois Apr 12 '13 at 16:19
  • @RomainFrancois you're right, I should have know (dirk.eddelbuettel.com/code/rcpp/html/get__na_8h.html). – leodido Apr 12 '13 at 21:37
  • Anyway I'm still trying to make this function to work for character vectors. At the moment: no progress. – leodido Apr 12 '13 at 21:38
  • case STRSXP: return na_omit_template<STRSXP, String>(x); works. – Artem Klevtsov May 20 '16 at 0:30
0

The templating mechanism seems to work fine.

> na_test_cpp(as.numeric(c(1, NA, NA, 1, 2, NA)))
[1]  2 NA NA NA  1

> na_test_cpp(as.integer(c(1, NA, NA, 1, 2, NA)))
[1] 1 2

It is this code that is working for INTSXP but not for REALSXP

Vector<REALSXP> na_test_real(const Vector<REALSXP>& x) {
  return setdiff(x, Vector<REALSXP>::create(::traits::get_na<REALSXP>()));
}
| improve this answer | |
  • Good point! I'm start to think also that can be a problem related to setdiff function, but I'm not sure. That's why I wrote here, hoping some guru help me :) – leodido Apr 11 '13 at 18:41
  • That is because > evalCpp( "NA_REAL == NA_REAL" ) gives [1] FALSE. you cannot test for equality with numeric NA, you must use is_na. – Romain Francois Apr 13 '13 at 6:15
0

Some implementations:

// naive
template <int RTYPE>
Vector<RTYPE> na_omit_impl(const Vector<RTYPE>& x) {
    std::size_t n = x.size();
    // Estimate out length
    std::size_t n_out = 0;
    for(std::size_t i = 0; i < n; ++i) {
        if (Vector<RTYPE>::is_na(x[i])) continue;
        ++n_out;
    }
    // exit if no NAs
    if (n_out == n) return x;
    // allocate vector without filling
    Vector<RTYPE> res = no_init(n_out);
    // fill result vector
    for(std::size_t i = 0, j = 0; i < n; ++i) {
        if (Vector<RTYPE>::is_na(x[i])) continue;
        res[j] = x[i];
        ++j;
    }
    return res;
}

// STL algorithms
template <int RTYPE>
struct not_na {
    typedef typename Vector<RTYPE>::stored_type type;
    bool operator() (const type& i) {
        return !Vector<RTYPE>::is_na(i);
    }
};

template <int RTYPE>
Vector<RTYPE> na_omit_impl(const Vector<RTYPE>& x) {
    // Estimate out length
    std::size_t n_out = std::count_if(x.begin(), x.end(), not_na<RTYPE>());
    // exit if no NAs
    if (n_out == x.size()) return x;
    // allocate vector without filling
    Vector<RTYPE> res = no_init(n_out);
    // fill result vector
    std::copy_if(x.begin(), x.end(), res.begin(), not_na<RTYPE>());
    return res;
}

// Rcpp sugar
template <class T>
T na_omit_impl(const T& x) {
    return x[!is_na(x)];
}

// Rcpp sugar
template <class T>
T na_omit_impl(const T& x) {
    return Rcpp::na_omit(x);
}

All implementations works with the RCPP_RETURN_VECTOR macro:

// [[Rcpp::export]]
RObject na_omit(RObject x){
    RCPP_RETURN_VECTOR(na_omit_impl, x);
}
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.