Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add native SIR filter example #85

Merged
merged 5 commits into from
Oct 16, 2024
Merged

Add native SIR filter example #85

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
129705f
Simplify test setup
richfitz Oct 15, 2024
7273af6
Add alternative RNG access
richfitz Oct 15, 2024
1a988c4
Alternative sir implementation
richfitz Oct 15, 2024
0d18d9a
Rename and remove dust dependency
richfitz Oct 15, 2024
4f84a3c
Bump version
richfitz Oct 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 1 addition & 4 deletions DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: monty
Title: Monte Carlo Models
Version: 0.2.18
Version: 0.2.19
Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
email = "[email protected]"),
person("Wes", "Hinsley", role = "aut"),
Expand Down Expand Up @@ -31,7 +31,6 @@ Suggests:
coda,
decor,
knitr,
dust,
mockery,
mvtnorm,
numDeriv,
Expand All @@ -43,5 +42,3 @@ Suggests:
Config/testthat/edition: 3
Language: en-GB
VignetteBuilder: knitr
Remotes:
mrc-ide/dust
12 changes: 12 additions & 0 deletions R/cpp11.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,18 @@ monty_rng_state <- function(ptr, is_float) {
.Call(`_monty_monty_rng_state`, ptr, is_float)
}

cpp_monty_random_real <- function(ptr) {
.Call(`_monty_cpp_monty_random_real`, ptr)
}

cpp_monty_random_binomial <- function(r_size, r_prob, ptr) {
.Call(`_monty_cpp_monty_random_binomial`, r_size, r_prob, ptr)
}

cpp_monty_random_exponential_rate <- function(r_rate, ptr) {
.Call(`_monty_cpp_monty_random_exponential_rate`, r_rate, ptr)
}

monty_rng_pointer_init <- function(n_streams, seed, long_jump, algorithm) {
.Call(`_monty_monty_rng_pointer_init`, n_streams, seed, long_jump, algorithm)
}
Expand Down
23 changes: 23 additions & 0 deletions R/rng2.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
monty_random_alloc <- function(n_streams = 1L, seed = NULL,
deterministic = FALSE) {
float <- FALSE
ptr <- monty_rng_alloc(seed, n_streams, deterministic, float)
ret <- list(ptr = ptr, n_streams = n_streams, deterministic = deterministic)
class(ret) <- "monty_random_state"
ret
}


monty_random_real <- function(state) {
cpp_monty_random_real(state$ptr)
}


monty_random_binomial <- function(size, prob, state) {
cpp_monty_random_binomial(size, prob, state$ptr)
}


monty_random_exponential_rate <- function(rate, state) {
cpp_monty_random_exponential_rate(rate, state$ptr)
}
70 changes: 47 additions & 23 deletions src/cpp11.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,27 @@ extern "C" SEXP _monty_monty_rng_state(SEXP ptr, SEXP is_float) {
return cpp11::as_sexp(monty_rng_state(cpp11::as_cpp<cpp11::decay_t<SEXP>>(ptr), cpp11::as_cpp<cpp11::decay_t<bool>>(is_float)));
END_CPP11
}
// random2.cpp
cpp11::doubles cpp_monty_random_real(SEXP ptr);
extern "C" SEXP _monty_cpp_monty_random_real(SEXP ptr) {
BEGIN_CPP11
return cpp11::as_sexp(cpp_monty_random_real(cpp11::as_cpp<cpp11::decay_t<SEXP>>(ptr)));
END_CPP11
}
// random2.cpp
cpp11::doubles cpp_monty_random_binomial(cpp11::doubles r_size, cpp11::doubles r_prob, SEXP ptr);
extern "C" SEXP _monty_cpp_monty_random_binomial(SEXP r_size, SEXP r_prob, SEXP ptr) {
BEGIN_CPP11
return cpp11::as_sexp(cpp_monty_random_binomial(cpp11::as_cpp<cpp11::decay_t<cpp11::doubles>>(r_size), cpp11::as_cpp<cpp11::decay_t<cpp11::doubles>>(r_prob), cpp11::as_cpp<cpp11::decay_t<SEXP>>(ptr)));
END_CPP11
}
// random2.cpp
cpp11::doubles cpp_monty_random_exponential_rate(cpp11::doubles r_rate, SEXP ptr);
extern "C" SEXP _monty_cpp_monty_random_exponential_rate(SEXP r_rate, SEXP ptr) {
BEGIN_CPP11
return cpp11::as_sexp(cpp_monty_random_exponential_rate(cpp11::as_cpp<cpp11::decay_t<cpp11::doubles>>(r_rate), cpp11::as_cpp<cpp11::decay_t<SEXP>>(ptr)));
END_CPP11
}
// rng_pointer.cpp
cpp11::sexp monty_rng_pointer_init(int n_streams, cpp11::sexp seed, int long_jump, std::string algorithm);
extern "C" SEXP _monty_monty_rng_pointer_init(SEXP n_streams, SEXP seed, SEXP long_jump, SEXP algorithm) {
Expand Down Expand Up @@ -172,29 +193,32 @@ extern "C" SEXP _monty_test_xoshiro_run(SEXP obj) {

extern "C" {
static const R_CallMethodDef CallEntries[] = {
{"_monty_monty_rng_alloc", (DL_FUNC) &_monty_monty_rng_alloc, 4},
{"_monty_monty_rng_beta", (DL_FUNC) &_monty_monty_rng_beta, 6},
{"_monty_monty_rng_binomial", (DL_FUNC) &_monty_monty_rng_binomial, 6},
{"_monty_monty_rng_cauchy", (DL_FUNC) &_monty_monty_rng_cauchy, 6},
{"_monty_monty_rng_exponential_mean", (DL_FUNC) &_monty_monty_rng_exponential_mean, 5},
{"_monty_monty_rng_exponential_rate", (DL_FUNC) &_monty_monty_rng_exponential_rate, 5},
{"_monty_monty_rng_gamma_rate", (DL_FUNC) &_monty_monty_rng_gamma_rate, 6},
{"_monty_monty_rng_gamma_scale", (DL_FUNC) &_monty_monty_rng_gamma_scale, 6},
{"_monty_monty_rng_hypergeometric", (DL_FUNC) &_monty_monty_rng_hypergeometric, 7},
{"_monty_monty_rng_jump", (DL_FUNC) &_monty_monty_rng_jump, 2},
{"_monty_monty_rng_long_jump", (DL_FUNC) &_monty_monty_rng_long_jump, 2},
{"_monty_monty_rng_multinomial", (DL_FUNC) &_monty_monty_rng_multinomial, 6},
{"_monty_monty_rng_nbinomial", (DL_FUNC) &_monty_monty_rng_nbinomial, 6},
{"_monty_monty_rng_normal", (DL_FUNC) &_monty_monty_rng_normal, 7},
{"_monty_monty_rng_pointer_init", (DL_FUNC) &_monty_monty_rng_pointer_init, 4},
{"_monty_monty_rng_pointer_sync", (DL_FUNC) &_monty_monty_rng_pointer_sync, 2},
{"_monty_monty_rng_poisson", (DL_FUNC) &_monty_monty_rng_poisson, 5},
{"_monty_monty_rng_random_normal", (DL_FUNC) &_monty_monty_rng_random_normal, 5},
{"_monty_monty_rng_random_real", (DL_FUNC) &_monty_monty_rng_random_real, 4},
{"_monty_monty_rng_state", (DL_FUNC) &_monty_monty_rng_state, 2},
{"_monty_monty_rng_uniform", (DL_FUNC) &_monty_monty_rng_uniform, 6},
{"_monty_test_rng_pointer_get", (DL_FUNC) &_monty_test_rng_pointer_get, 2},
{"_monty_test_xoshiro_run", (DL_FUNC) &_monty_test_xoshiro_run, 1},
{"_monty_cpp_monty_random_binomial", (DL_FUNC) &_monty_cpp_monty_random_binomial, 3},
{"_monty_cpp_monty_random_exponential_rate", (DL_FUNC) &_monty_cpp_monty_random_exponential_rate, 2},
{"_monty_cpp_monty_random_real", (DL_FUNC) &_monty_cpp_monty_random_real, 1},
{"_monty_monty_rng_alloc", (DL_FUNC) &_monty_monty_rng_alloc, 4},
{"_monty_monty_rng_beta", (DL_FUNC) &_monty_monty_rng_beta, 6},
{"_monty_monty_rng_binomial", (DL_FUNC) &_monty_monty_rng_binomial, 6},
{"_monty_monty_rng_cauchy", (DL_FUNC) &_monty_monty_rng_cauchy, 6},
{"_monty_monty_rng_exponential_mean", (DL_FUNC) &_monty_monty_rng_exponential_mean, 5},
{"_monty_monty_rng_exponential_rate", (DL_FUNC) &_monty_monty_rng_exponential_rate, 5},
{"_monty_monty_rng_gamma_rate", (DL_FUNC) &_monty_monty_rng_gamma_rate, 6},
{"_monty_monty_rng_gamma_scale", (DL_FUNC) &_monty_monty_rng_gamma_scale, 6},
{"_monty_monty_rng_hypergeometric", (DL_FUNC) &_monty_monty_rng_hypergeometric, 7},
{"_monty_monty_rng_jump", (DL_FUNC) &_monty_monty_rng_jump, 2},
{"_monty_monty_rng_long_jump", (DL_FUNC) &_monty_monty_rng_long_jump, 2},
{"_monty_monty_rng_multinomial", (DL_FUNC) &_monty_monty_rng_multinomial, 6},
{"_monty_monty_rng_nbinomial", (DL_FUNC) &_monty_monty_rng_nbinomial, 6},
{"_monty_monty_rng_normal", (DL_FUNC) &_monty_monty_rng_normal, 7},
{"_monty_monty_rng_pointer_init", (DL_FUNC) &_monty_monty_rng_pointer_init, 4},
{"_monty_monty_rng_pointer_sync", (DL_FUNC) &_monty_monty_rng_pointer_sync, 2},
{"_monty_monty_rng_poisson", (DL_FUNC) &_monty_monty_rng_poisson, 5},
{"_monty_monty_rng_random_normal", (DL_FUNC) &_monty_monty_rng_random_normal, 5},
{"_monty_monty_rng_random_real", (DL_FUNC) &_monty_monty_rng_random_real, 4},
{"_monty_monty_rng_state", (DL_FUNC) &_monty_monty_rng_state, 2},
{"_monty_monty_rng_uniform", (DL_FUNC) &_monty_monty_rng_uniform, 6},
{"_monty_test_rng_pointer_get", (DL_FUNC) &_monty_test_rng_pointer_get, 2},
{"_monty_test_xoshiro_run", (DL_FUNC) &_monty_test_xoshiro_run, 1},
{NULL, NULL, 0}
};
}
Expand Down
73 changes: 73 additions & 0 deletions src/random2.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
#include <cpp11/doubles.hpp>
#include <cpp11/external_pointer.hpp>
#include <cpp11/integers.hpp>
#include <cpp11/raws.hpp>

#include <monty/r/random.hpp>
#include <monty/random/random.hpp>
#include <monty/utils.hpp>

using default_rng64 = monty::random::prng<monty::random::generator<double>>;

template <typename T>
void check_length(T x, size_t len, const char * name) {
const size_t len_given = Rf_length(x);
if (len_given != len) {
cpp11::stop("Expected '%s' to have length %d, not %d",

Check warning on line 16 in src/random2.cpp

View check run for this annotation

Codecov / codecov/patch

src/random2.cpp#L16 

Added line #L16 was not covered by tests
name, len, len_given);
}
}


[[cpp11::register]]
cpp11::doubles cpp_monty_random_real(SEXP ptr) {
default_rng64 *rng = cpp11::as_cpp<cpp11::external_pointer<default_rng64>>(ptr).get();
const int n = rng->size();
cpp11::writable::doubles r_y = cpp11::writable::doubles(n);
double *y = REAL(r_y);
for (int i = 0; i < n; ++i) {
y[i] = monty::random::random_real<double>(rng->state(i));
}
return r_y;
}


[[cpp11::register]]
cpp11::doubles cpp_monty_random_binomial(cpp11::doubles r_size,
cpp11::doubles r_prob,
SEXP ptr) {
default_rng64 *rng = cpp11::as_cpp<cpp11::external_pointer<default_rng64>>(ptr).get();
const int n = rng->size();
check_length(r_size, n, "size");
check_length(r_prob, n, "size");
const double * size = REAL(r_size);
const double * prob = REAL(r_prob);

cpp11::writable::doubles r_y = cpp11::writable::doubles(n);
double *y = REAL(r_y);
for (int i = 0; i < n; ++i) {
auto &state = rng->state(i);
y[i] = monty::random::binomial<double>(state, size[i], prob[i]);
}

return r_y;
}


[[cpp11::register]]
cpp11::doubles cpp_monty_random_exponential_rate(cpp11::doubles r_rate,
SEXP ptr) {
default_rng64 *rng = cpp11::as_cpp<cpp11::external_pointer<default_rng64>>(ptr).get();
const int n = rng->size();
check_length(r_rate, n, "size");
const double * rate = REAL(r_rate);

cpp11::writable::doubles r_y = cpp11::writable::doubles(n);
double *y = REAL(r_y);
for (int i = 0; i < n; ++i) {
auto &state = rng->state(i);
y[i] = monty::random::exponential_rate<double>(state, rate[i]);
}

return r_y;
}
Loading
Loading