[Request]: Sugar for pow(int/double/IntegerVector/NumericVector x, IntegerVector/NumericVector z)

[kendonB]

Currently, it seems that only int z and double z work.

These, for example, would work:

// [[Rcpp::export]]
NumericVector pow_fails_double_vector(double x, NumericVector z) {
  return pow(x, z);
}

// This would do element-wise pow, as in R
// [[Rcpp::export]]
NumericVector pow_fails_vector_vector(NumericVector x, NumericVector z) {
  return pow(x, z);
}

Is making sugar functions for more classes hard? Is there a resource to help with creating them?

[eddelbuettel]

It's not hard if you have a wee bit of patience. Many of the functions are being constructed with macros. Do you want to give it a try?

Also, you can of course also just define them locally as helper functions in your package(s).

Now, before we get to carried away, the first one make no sense. Just wrap a rep(...) around the double to make it a vector, and the would have vector, vector. In general you do NOT want too many combinations of these things. I could possibly see a use for that case, but then we also have a templates to worry about here. You probably don't want too many combinations here...

[kendonB]

Could you point me to a file with vector-vector implemented that I can try and replicate?

This? https://github.com/RcppCore/Rcpp/blob/master/inst/include/Rcpp/sugar/functions/pmin.h

[thirdwing]

What about the Plus_Vector_Vector from https://github.com/RcppCore/Rcpp/blob/master/inst/include/Rcpp/sugar/operators/plus.h ?

pow can be a kind of operator between two vector.

[nathan-russell]

Not thoroughly tested or anything, but something like this might work as well:

#include <Rcpp.h>

namespace impl {

template <int RTYPE, bool LHS_NA, typename LHS_T, bool RHS_NA, typename RHS_T>
class Pow : public Rcpp::VectorBase<RTYPE, true, Pow<RTYPE, LHS_NA, LHS_T, RHS_NA, RHS_T> > {
public:
    typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
    typedef Rcpp::VectorBase<RTYPE, LHS_NA, LHS_T> LHS_Type;
    typedef Rcpp::VectorBase<RTYPE, RHS_NA, RHS_T> RHS_Type;

    Pow(const LHS_Type& lhs_, const RHS_Type& rhs_)
        : lhs(lhs_), rhs(rhs_),
          lhs_sz(lhs.size()),
          rhs_sz(rhs.size()),
          sz(lhs_sz >= rhs_sz ? lhs_sz : rhs_sz)
    {}

    stored_type operator[](const R_xlen_t i) const {
        if (Rcpp::traits::is_na<RTYPE>(lhs[i % lhs_sz])) return lhs[i % lhs_sz];
        if (Rcpp::traits::is_na<RTYPE>(rhs[i % rhs_sz])) return rhs[i % rhs_sz];
        return std::pow(lhs[i % lhs_sz], rhs[i % rhs_sz]);
    }

    R_xlen_t size() const {
        return sz;
    }

private:
    const LHS_Type& lhs;
    const RHS_Type& rhs;

    R_xlen_t lhs_sz;
    R_xlen_t rhs_sz;
    R_xlen_t sz;
};

} // impl

template <int RTYPE, bool LHS_NA, typename LHS_T, bool RHS_NA, typename RHS_T>
inline impl::Pow<RTYPE, LHS_NA, LHS_T, RHS_NA, RHS_T>
Pow(const Rcpp::VectorBase<RTYPE, LHS_NA, LHS_T>& lhs,
    const Rcpp::VectorBase<RTYPE, RHS_NA, RHS_T>& rhs)
{
    return impl::Pow<RTYPE, LHS_NA, LHS_T, RHS_NA, RHS_T>(lhs, rhs);
}

// [[Rcpp::export]]
Rcpp::NumericVector vpow(Rcpp::NumericVector lhs, Rcpp::NumericVector rhs) {
    return Pow(lhs, rhs);
}

/*** R

all.equal(2 ^ (1:5), vpow(2, 1:5))
# [1] TRUE

all.equal((1:5) ^ (1:5), vpow(1:5, 1:5))
# [1] TRUE

all.equal((1:5) ^ 2, vpow(1:5, 2))
# [1] TRUE

*/

Maybe just add some specializations for when LHS_NA and / or RHS_NA are false.

[eddelbuettel]

Thanks!

Now the bigger question for Team Rcpp: Do we want / need this?

[nathan-russell]

@eddelbuettel I have no strong feelings either way (as usual); adding this feature would involve replacing the current Rcpp::pow, which means that in addition to testing the added functionality, the new implementation would have to support the current behavior.

[eddelbuettel]

Right. I am marginally against it because of this maintenance burden, and the general lack of use, the case of @kendonB notwithstanding. Or can anyone of you recall ever wanting this more than once?

Maybe just make it an example in, say, the Rcpp Gallery?

Thinking of this, though: using mapply() for this may even be a decent use case to sweep such a simple operation across two vectors.

[thirdwing]

Personally I prefer adding this as an example.

It can be in Rcpp Gallery. Or maybe we can even add a paragraph into the sugar vignette, like "How to write your own Rcpp sugar".

[eddelbuettel]

Oh, I like that idea too! It may be too advanced for the Gallery, but it is pretty spot-on for the vignette.

[coatless]

Just a heads up, there does exist some documentation already on implementing sugar operators within the sugar vignette. Specifically, the document provides an overview on the curiously recurring template pattern (CRTP), VectorBase class, and an example of how sapply was created (c.f http://dirk.eddelbuettel.com/code/rcpp/Rcpp-sugar.pdf#page=8) with subsequent step by step breakdown of each line of code. So, this would serve as a second, more advanced example if added to the sugar vignette.

[dcdillon]

FWIW, and I apologize for the shameless plug, https://github.com/dcdillon/RcppHoney would give you this function.

#include <RcppHoney.hpp>

// [[Rcpp::depends(RcppHoney)]]

// [[Rcpp::export]]
Rcpp::NumericVector pow_fails_double_vector(double x, Rcpp::NumericVector z) {
  return Rcpp::wrap(RcppHoney::pow(x, z));
}

// This would do element-wise pow, as in R
// [[Rcpp::export]]
Rcpp::NumericVector pow_fails_vector_vector(Rcpp::NumericVector x, Rcpp::NumericVector z) {
  return Rcpp::wrap(RcppHoney::pow(x, z));
}

> pow_fails_double_vector(2.0, c(1, 2, 3))
[1] 2 4 8
> pow_fails_vector_vector(c(2.0, 3.0, 4.0), c(1, 2, 3))
[1]  2  9 64
> 

[eddelbuettel]

Working example or the answer didn't exist.

[coatless]

So... Add this as a second example to the Sugar vignette or a gallery post?

[eddelbuettel]

Both?