feat: global RNG access

[hsloot]

Provide access to dqrng's RNG for expert users:

• Based on daqana/dqrng@access-rng; see https://github.com/hsloot/dqrng/compare/425ae99..f8346b9 for amendments
• Added explanatory text to vignette
• Fixed new function dqrng:: get_rng() (Rcpp::Xptr needs to be constructed without registering a finalizer to avoid accidentally deleting the global RNG); also slightly improve the implementation
• I did not add the tests of my originally linked branch as they do not integrate well into the setup without adjustments

---

Close #41

[rstub]

I am not sure why the CI tests are suddenly failing here. Anyway, right now I am wondering why you advice against wrapping the raw point into a smart pointer. It is clear that one must not do that with a std::unique_ptr. But why not wrapping it inside a std::shared_ptr? What am I missing?

[hsloot]

I am not sure why the CI tests are suddenly failing here.

This appears to happen because the new chunk is evaluated (contrary the the few chunks before). This seems to occur as the vignette building environment in the check does not have the BH dependency available. I will change it so that it is also not executed.

---

Edit: If you actually want to execute those chunks, you might need to add some additional packages to Suggests or Config/Needs/check, and you might want to add % \VignetteDepends{...} with comma-separated dependencies at the top of the vignettes.

[hsloot]

Anyway, right now I am wondering why you advice against wrapping the raw point into a smart pointer. It is clear that one must not do that with a std::unique_ptr. But why not wrapping it inside a std::shared_ptr? What am I missing?

A shared pointer also represents ownership of the underlying object; hence the underlying object is deleted once all ownership claims are gone. In this case, the RNG should be owned only by the dqrng package, and it must not be deleted because some last std::shared_ptr goes out of scope. Using Rcpp::Xptr is also not ideal but the simplest implementation with the interface-option present, I think. IMO, the safest way to make sure nothing wrong happens on accident is for the user to obtain the raw pointer and never do anything except dereferencing it. BTW, std::weak_ptr would also not be good, as it is for temporary ownership.

For background, from https://en.cppreference.com/w/cpp/memory/shared_ptr:

std::shared_ptr is a smart pointer that retains shared ownership of an object through a pointer. Several shared_ptr objects may own the same object. The object is destroyed and its memory deallocated when either of the following happens:

• the last remaining shared_ptr owning the object is destroyed;
• the last remaining shared_ptr owning the object is assigned another pointer via operator= or reset().

The object is destroyed using delete-expression or a custom deleter that is supplied to shared_ptr during construction.

[hsloot]

This answer was for the current implementation. Of course, if the original pointer would be the first std::shared_ptr wrapping the raw pointer, and we somehow manage to expose this shared pointer through the interface, it should be fine. But, this would require a different, more comple implementation without Rcpp::Xptr and the Rcpp interface, I think?

[hsloot]

An easier and more natural way might also be if rng is directly stored as an Rcpp::Xptr which is then also exposed to the used (as I understand it, it is pretty similar to the idea of a shared pointer)?

[rstub]

This has made me think quite a lot about how the RNGs are currently stored (which is good!), and it seems like things can become quite complex.

1. For storage a dqrng::rng64_t == std::shared_ptr<dqrng::random_64bit_generator> is used, which is initialized when the package gets loaded.
2. This shared_ptr is deleted together with the RNG itself when the package is unloaded or (more typically) R exits.
3. Another reason for the RNG to be deleted is dqRNGkind. At the moment this creates a new RNG and wraps it in a new shared_ptr.
4. We can only return Rcpp::Xptr with such a function, since this has to go through R's C API. However instead of a RNG pointer, it should be possible to wrap a pointer that points at a smart pointer, but so far I have not gotten this to work. Have you tried that at some point?
5. Wrapping the raw pointer into a smart pointer can delete the internal RNG, since this new smart pointer would be "disconnected" from the internally used smart pointer. In addition, this smart pointer might give a false sense of security since the raw pointer can be invalidated by calling dqRNGkind.
6. I don't think it would work to use an Rcpp:Xptr internally and hand that out explicitly: The internal Xptr should have a delete-finalizer, while the external one must not as you correctly noticed. Or do you have an idea for that?

In addition one needs to define the use-case properly: With the current "experts only" implementation one can support users that write their own C++ code, understand what is going on and have control over the execution. However, I would like to also support package authors that write code that needs a fast (potentially parallel) RNG. Currently these users have to create their own RNGs and care for properly seeding them. I would like to offer them a mechanism to easily integrate with.

So on the one hand would be clear guidance. In order to protect from the dangers of somebody calling dqRNGkind, it is crucial that users do not save the pointer outside of a function, e.g. foo() is bad while bar() is fine:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng.h>

auto rng1 = dqrng::get_rng();

// [[Rcpp::export]]
void foo() {
  Rcpp::Rcout << (*rng1)() << std::endl;
}

// [[Rcpp::export]]
void bar() {
  auto rng2 = dqrng::get_rng();
  Rcpp::Rcout << (*rng2)() << std::endl;
}

For the second type of problems (do not modify or delete ...), it might make sense to provide a wrapper class that makes this impossible or at least hard to do. And document that as the only supported way to use dqrng's RNG. Something along these lines:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng.h>

namespace dqrng {
namespace external {
class random_64bit_generator : public dqrng::random_64bit_generator {
private:
  Rcpp::XPtr<dqrng::random_64bit_generator> gen;

public:
  random_64bit_generator() : gen(dqrng::get_rng()) {};
  virtual result_type operator() () {return gen->operator()();}
  virtual void seed(result_type seed) {throw std::runtime_error("Seed handling not supported for this class!");}
  virtual void seed(result_type seed, result_type stream) {throw std::runtime_error("Seed handling not supported for this class!");}
};
};
};

// [[Rcpp::export]]
void baz() {
  dqrng::external::random_64bit_generator rng3;
  Rcpp::Rcout << rng3() << std::endl;
}

Notes:

• This does not compile w/o moving the implementations for bit64() and bit32 from random_64bit_wrapper to random_64bit_generator. But that does make sense anyway.
• The naming of classes and namesapces here is very much provisional ...\
• Seeding is turned off as this should be handled via dqset.seed or dqset_seed.

This brings me to the final point: The sampling methods from dqrng_sample.h currently require a reference to a dqrng::rng64_t as input, which one cannot construct with the current pointer. This is something I would like to change. The above wrapper class might help here.

Any comments?

[hsloot]

6. I don't think it would work to use an Rcpp:Xptr internally and hand that out explicitly: The internal Xptr should have a delete-finalizer, while the external one must not as you correctly noticed. Or do you have an idea for that?

I will have to think about this with some time on my hand, but at the top of my head I would try to find out first if 6. is really true. I do not know the implementation well, but I would hope that when the Xptr is copy-constructed, it will not register the finalized again.

---

Edit: Context

[hsloot]

I share your assessments in the beginning, except for no. 6 for which I am not sure (see above). A few comments:

Storing the RNG

The RNG is currently stored as std::shared_ptr<random_64bit_generator>; I wonder, if it would not be more natural to make it a std::unique_ptr<random_64bit_generator> as it is not really shared anywhere.

Supporting parallel RNG usage

However, I would like to also support package authors that write code that needs a fast (potentially parallel) RNG. Currently these users have to create their own RNGs and care for properly seeding them. I would like to offer them a mechanism to easily integrate with.

I am not sure if the global RNG access is really that useful writing parallel sampling algorithms; looking at the sample stream-individual seeding is nothing you could get rid of, right?

Guidance and wrapper

I do not see a usecase for which client-code that would call dqrng::get_rng() (whatever that returns) for anything else than to pass the RNG into distribution objects (see example in details).

• My first idea was to just provide an interface to the raw pointer (Rcpp::Xptr without finalizer is somewhat close to that) and give guidance to use it after seeding and not use it in anyway other than described above. Yes, you could (accidentally) delete the RNG, but for that you would have to do it intentionally and leave the API's intended scope.
• Provide an interface to an external-wrapper class as you sketched it. This is a bit more safe, but the RNG would still implicitly be available (through dqrng::get_rng()). However, it would have the benefit to encapsulate implementation from API.

Idea: Write templated sampling algorithm and inject base-distributions (uniform, exponential, ...) such that it is simple to use it alongside existing frameworks (base R, dqrng, ...).

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng.h>

namespace internal {
  template <typename _RealType, typename _UniformRealDistribution>
  class pareto_distribution {
    public:
    using result_type = _RealType;

    class param_type {
      // implementation
    };

    explicit pareto_distribution(const _RealType alpha,
                                 const _RealType lower_bound)
        : parm_{alpha, lower_bound} {
    // implementation
    }

    // implementation

    template <typename _Engine>
    result_type operator()(_Engine&& engine) {
        return (*this)(std::forward<_Engine>(engine), parm_);
    }

    template <typename _Engine>
    result_type operator()(_Engine&& engine, const param_type& parm) {
        return parm.lower_bound_ /
               std::pow(unit_uniform_real_dist_(std::forward<_Engine>(engine)),
                        1. / parm.alpha_);
    }

    private:
    param_type parm_{};
    _UniformRealDistribution unit_uniform_real_dist_{_RealType{0}, _RealType{1}};

    // implementation
  };

  class r_uniform_real_distribution {
    // implementation to sample via R with an std-like interface 
  };

  class r_engine {
    // implementation for an std-like interface for R engine
  };
};

// [[Rcpp::export]]
Rcpp::NumericVector rpareto_r(const double alpha, const double lower_bound) {
  using pareto_distribution = internal:: pareto_distribution<double, internal::r_uniform_real_distribution>;
  std::unique_ptr rng = std::make_unique(engine{});
  
  // wrapper-implementation
}

// [[Rcpp::export]]
Rcpp::NumericVector rpareto_dqrng(const double alpha, const double lower_bound) {
  using pareto_distribution = internal:: pareto_distribution<double, dqrng::r_uniform_distribution>;
  Rcpp::Xptr rng = dqrng::get_rng();
  
  // wrapper-implementation
}

Problem: Currently, this would require duplicating all of dqrng's sampling algorithms and seeding interface to obtain the functionality. In addition, this could not be done by different packages at the same time without copying each other. Getting a global access would solve this. How this is done is irrelevant; we only need to pass the underlying generator into the operator()(const Engine& ) method. Obviously, this can be simplified a lot if you only target a dqrng variant.

Call to sampling routines in dqrng_sample.h

I think (if you would be comfortable with changing this) the easiest thing would be to not pass a reference to a pointer into these methods, but a reference to the object, e.g.,

// in dqrng_sampe.h
namespace sample {

inline Rcpp::Vector<RTYPE> replacement(dqrng::rng64_t::element_type &rng, INT m, INT n, int offset) {
  using storage_t = typename Rcpp::traits::storage_type<RTYPE>::type;
  Rcpp::Vector<RTYPE> result(Rcpp::no_init(n));
  std::generate(result.begin(), result.end(),
                [m, offset, rng] () {return static_cast<storage_t>(offset + rng(m));});
  return result;
}

// other implementations

template<int RTYPE, typename INT>
inline Rcpp::Vector<RTYPE> sample(dqrng::rng64_t::element_type &rng, INT m, INT n, bool replace, int offset = 0) {
  if (replace || n <= 1) {
    return dqrng::sample::replacement<RTYPE, INT>(rng, m, n, offset);
  } else {
    if (!(m >= n))
      Rcpp::stop("Argument requirements not fulfilled: m >= n");
    if (m < 2 * n) {
      return dqrng::sample::no_replacement_shuffle<RTYPE, INT>(rng, m, n, offset);
    } else if (m < 1000 * n) {
      return dqrng::sample::no_replacement_set<RTYPE, INT, dqrng::minimal_bit_set>(rng, m, n, offset);
    } else {
      return dqrng::sample::no_replacement_set<RTYPE, INT, dqrng::minimal_hash_set<INT>>(rng, m, n, offset);
    }
  }
}

}; // sample

// in dqrng.cpp
// [[Rcpp::export(rng = false)]]
Rcpp::IntegerVector dqsample_int(int m,
                                 int n,
                                 bool replace = false,
                                 Rcpp::Nullable<Rcpp::NumericVector> probs = R_NilValue,
                                 int offset = 0) {
    if (!(m > 0 && n >= 0))
        Rcpp::stop("Argument requirements not fulfilled: m > 0 && n >= 0");
    return dqrng::sample::sample<INTSXP, uint32_t>(*rng, uint32_t(m), uint32_t(n), replace, offset);
}

Where does this leave us?

I think the next steps would be (in that order):

• Add some tests?
• Move bit32 and.bit64 implementations to dqrng:: random_64bit_generator.
• Solve the problem with dqrng_sample.h
• Create the wrapper class for external use (do you have any preference for naming or in which namespace to put it?).
• Add documentation.

If that is fine for you, I could start with that.

[hsloot]

Another thought: As your proposed wrapper encapsulates implementation, how to store the RNG internally and how to improve the user experience for setting up and seeding the RNG in parallel usage are sort of independent of exposing the RNG for external use. Maybe this should be discussed in separate issues?

[rstub]

Thanks @hsloot! You are right that these other things should be separate issues. So what do we need here:

• One question: Why are you using dqrng::rng64_t::element_type instead of dqrng::random_64bit_generator?
• Tests would be good. What adaptions would be needed for the tests in your original branch?

BTW, I had another look at 6. from above and it looks like you might be correct. When an XPtr, which has been copy-constructed, goes out of scope, the original XPtr is unaffected in my tests:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng_generator.h>

Rcpp::XPtr<dqrng::random_64bit_generator> rng(new dqrng::random_64bit_wrapper<>());

Rcpp::XPtr<dqrng::random_64bit_generator> get_ptr() {
return rng;
}

// [[Rcpp::export]]
void foo() {
  Rcpp::XPtr<dqrng::random_64bit_generator> lrng(get_ptr());
  Rcpp::Rcout << (*lrng.checked_get())() << std::endl;
}

// [[Rcpp::export]]
void bar() {
  Rcpp::Rcout << (*rng.checked_get())() << std::endl;
}

// [[Rcpp::export]]
void baz(bool finalizer) {
  Rcpp::XPtr<dqrng::random_64bit_generator> lrng(get_ptr().get(), finalizer);
  Rcpp::Rcout << (*lrng.checked_get())() << std::endl;
}

Calling foo() or baz(FALSE) does not cause any problems. Calling baz(TRUE) many times eventually invalidates rng. I probably should ask on the rcpp-devel ML to get a definitive answer.

So storing internally in a XPtr might indeed simplify things. But as you said, this belongs into a different issue.

[hsloot]

One question: Why are you using dqrng::rng64_t::element_type instead of dqrng::random_64bit_generator?

As you have to know that they are the same, I find it more clear to only use one type throughout so that this dependence is clear.

Tests would be good. What adaptions would be needed for the tests in your original branch?

More or less only renaming a few things. But I would try to simplify it and remove this distribution caller approach and just duplicate a bit of code.

[hsloot]

The proposed changes to the dqrng-sample algorithm are not really part of any documented public API, right? Can I just change them or do we have to be careful with revdeps?

[rstub]

The sampling methods are used in https://daqana.github.io/dqrng/articles/parallel.html#pcg-multiple-streams-with-rcppparallel. So at least this would need to be updated. But I think we should leave those changes to the "usability enhancements". For now it would be great if you could include some tests and add yourself as contributor to DESCRIPTION.

[hsloot]

The new commits contain:

• DESCRIPTION update
• Added tests (requires adding BH to suggested packages)
• Add class for external RNG and use in tests

Still pending:

• Documentation update
• The last commit was required to make the test work, however, I have no clue why it was not working before. I want to investigate if there is some underlying error; maybe you can take a look, too.

[hsloot]

The reason why this weird fix is needed is a template vs. polymorphism issue; see #65.

[rstub]

Great find w.r.t. the boost specialization. Will have to look at this in more detail tomorrow.

[rstub]

This looks great. Thanks a lot!
I am contemplating a slightly different ordering for the header files. For example, it would be nice if the accessor class were directly available when one loads dqrng.h. But at the same time, it would be nice if the accessor class would be usable w/o modifying the distributions from boost. But I need to try out a few things first.

[hsloot]

Maybe consider putting all boost modifications in a single header and introduce a flag to activate/deactivate them?

[rstub]

What I have now:

• Move random_64bit_generator and random_64bit_accessor w/o the constructor definition into dqrng_types.h.
• This gets automatically included in dqrng_RcppExports.h.
• Put the constructor definition for random_64bit_accessor into dqrng.h.
• Put the boost modifications into dqrng_distributions.h

With these modifications one can use

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng.h>
#include <boost/random/exponential_distribution.hpp>

// [[Rcpp::export(rng = false)]]
Rcpp::NumericVector dqrexp_boost(const std::size_t n, const double rate = 1.0) {
  using dist_t = boost::random::exponential_distribution<double>;
  using parm_t = typename dist_t::param_type;

  const auto parm = parm_t{rate};
  auto dist = dist_t{};;
  auto out = Rcpp::NumericVector(Rcpp::no_init(n));

  auto engine = dqrng::random_64bit_accessor{};
  std::generate(out.begin(), out.end(), [&dist, &parm, &engine]() {
    return dist(engine, parm);
  });
  return out;
}

And get the unmodified exponential distribution from boost.

I feel uneasy about putting the the constructor into dqrng.h like this:

#ifndef dqrng_H
#define dqrng_H

#include "dqrng_RcppExports.h"

namespace dqrng {
random_64bit_accessor::random_64bit_accessor() : gen(dqrng::get_rng().get()) {}
} // namespace dqrng
#endif // dqrng_H

But this was the only way I found to get out of the circular dependency w/o introducing even more header files. I can add these directly here in the PR if that's ok with you.

[rstub]

Thanks again @hsloot!