1
\$\begingroup\$

I am trying to fill missing values with linearly interpolated values in an Rcpp::NumericVector -- leaving in place leading and trailing NA-values.

The desired output is below:

R> x <- c(NA, NA, NA, 1:4, NA, 6:8, NA, 10, NA, NA, NA)
R> naLinInt(x)
R> [1] NA NA NA  1  2  3  4  5  6  7  8  9 10 NA NA NA

I've implemented it via the below, but as i am very new to C++ i'd very much benefit from some help in understanding how to do this better.

I feel that i've really hacked my way through this -- in particular i suspect i've over-used pointers and that there's a better way to do this sort of indexing and offsetting.

#include<Rcpp.h>

//[[Rcpp::export]]
Rcpp::NumericVector naLinInt(Rcpp::NumericVector x) {
    // This function linearly interpolates between NA values

    Rcpp::Rcout << "We begin with: " << std::endl << x << std::endl;

    double * it1 = x.begin();
    double * it2 = x.begin();
    int step = 0;

    it1++;

    while( it1 < x.end() ) {
        if( Rcpp::NumericVector::is_na(*it1)) {
            if( Rcpp::NumericVector::is_na(*(it1 - 1))) {
                it1++;
                continue;
            } else {
            it2 = it1;
            step = 0;
            while(it2 < x.end() && Rcpp::NumericVector::is_na(*it2)) {
                step++;
                it2++;
            }
            if( it2 == x.end()) break;
            // step through missing values and replace
            for(int j = 0; j < step; j++) {
                *it1 = *(it1 - 1) + (*it2 - *(it1 - 1))/(step + 1 - j);
                it1++;
            }
            it1 = it2 + 1;
            }
        } else {
            it1++;
        }
    }

    return x;
}
\$\endgroup\$

1 Answer 1

3
\$\begingroup\$

As you suggest, the code does look very "pointery". I think it might be made clearer with use of some standard algorithms:

  • std::adjacent_find() (with suitable predicate functions) to find a number followed by NA, or a NA followed by a number.
  • std::generate() (with suitable generator function) to populate a series of values from one iterator to another.

With those, you won't need any arithmetic on iterators other than a subtraction to find the number of elements you're interpreting.


The form of such a solution is something like this (untested):

#include <Rcpp.h>
#include <algorithm>

constexpr auto is_na = Rcpp::NumericVector::is_na;

//[[Rcpp::export]]
Rcpp::NumericVector naLinInt(Rcpp::NumericVector x) {
    // This function linearly interpolates to fill sequences of NA
    // values surrounded by valid numbers.
    static auto const detect_start_na = [](auto a, auto b){
        return !is_na(a) && is_na(b);
    };
    static auto const detect_end_na = [](auto a, auto b){
        return is_na(a) && !is_na(b);
    };

    auto start = x.begin();

    while (true) {
        // Find transitions to and from NA values.  If we hit end of
        // vector whilst looking, our work is done.
        auto num_to_na = std::adjacent_find(start, x.end(), detect_start_na);
        auto na_to_num = std::adjacent_find(start, x.end(), detect_end_na);
        if (na_to_num == x.end()) {
            break;
        }

        // At this point, num_to_na points to the last number before
        // an interpolation block, and na_to_num points to the last NA
        // of that block.

        ++na_to_num;            // Now, both iterators point to numbers.
        auto const base = *num_to_na;
        auto const target = *na_to_num;

        // To count rails rather than posts, measure difference before
        // incrementing the start position.
        auto const gaps = std::distance(na_to_num, num_to_na);

        ++num_to_na;
        // Now both iterators point immediately *after* transition.

        auto const make_value = [base, target, gaps, i = std::size_t{0}]()
            mutable { return base + (++i * (target - base) / gaps); };
        std::generate(na_to_num, num_to_na, make_value);

        // Advance onwards
        start = na_to_num;
    }

    return x;
}

There are a couple of iterator increments there (because std::adjacent_find returns an iterator to the first of the matching pair), but there's quite a lot less arithmetic than the original.

Although this is longer than the original, much of the difference is the comments that help the reader understand how the state changes through the loop. That's something that I think is an improvement (though perhaps the terms "rails" and "posts" require a link to a definition of fencepost error for the uninitiated).

\$\endgroup\$
8
  • \$\begingroup\$ Great review. One could also use std::distance, but I guess it will get too verbose. \$\endgroup\$ Commented Nov 27, 2018 at 18:12
  • 1
    \$\begingroup\$ Hi @Incomputable - I did show std::distance in the worked example; just didn't feel the need to pick it out as something that makes the code clearer (because I think that's more arguable than the things from <algorithm>). \$\endgroup\$ Commented Nov 27, 2018 at 18:13
  • \$\begingroup\$ i get the below error when i try to run the code ... naLinInt.cpp:51:16: error: constexpr function's return type '(lambda at naLinInt.cpp:53:12)' is not a literal type constexpr auto na_transition(bool reverse) ^ naLinInt.cpp:53:12: note: '' is not literal because it is not an aggregate and has no constexpr constructors other than copy or move constructors return [reverse](auto a, auto b){ ^ 1 error generated. make: *** [naLinInt.o] Error 1 \$\endgroup\$
    – ricardo
    Commented Nov 28, 2018 at 13:17
  • \$\begingroup\$ I guess you're using an older C++ than me, @ricardo. Probably best to simply inline suitable definitions into the naLinInt() function. \$\endgroup\$ Commented Nov 28, 2018 at 13:22
  • \$\begingroup\$ I've changed the code to a simpler C++14-compliant version that should be easier to work with. \$\endgroup\$ Commented Nov 28, 2018 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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