## -----------------------------------------------------------------------------
#| include: false
library(vws)
library(tidyverse)

set.seed(1234)


## -----------------------------------------------------------------------------
#| eval: false
#| prompt: true
# file.path(path.package("vws"), "doc", "examples")


## // [[Rcpp::depends(vws, fntl)]]                 // <1>
## #include "vws.h"                                // <2>
## #include "MyRegion.h"
## 
## // [[Rcpp::export]]                             // <3>
## Rcpp::List sample(unsigned int n)
## {
##     MyRegion supp( ... );                       // <4>
##     vws::FMMProposal<double, MyRegion> h(supp); // <5>
##     h.refine(N - 1, 0.01);                      // <6>
##     auto out = vws::rejection(h, n);            // <7>
## 
##     return Rcpp::List::create(                  // <8>
##         Rcpp::Named("draws") = out.draws,
##         Rcpp::Named("rejects") = out.rejects
##     );
## }

## -----------------------------------------------------------------------------
#| eval: false
#| prompt: true
# Rcpp::sourceCpp("sample.cpp")
# out = sample(n = 100)
# head(out$draws)
# head(out$rejects)


## double z = 3;
## typedef function<double(double)> dfdd;
## dfdd f1 = [&](double x, double y) -> double { return x*y*z; };
## dfdd f2 = [&](double x, double y) { return x*y*z; };
## dfdd f3 = [=](double x, double y) { return x*y*z; };
## dfdd f4 = [](double x, double y) { return x*y*3; };

## double out = f1(1, 2);

## typedef std::function<double(double x, bool log)> dfdb;

## typedef std::function<double(
## 	const dfdb& w,  // <1>
## 	double lo,      // <2>
## 	double hi,      // <3>
## 	bool log        // <4>
## )> optimizer;

## typedef std::function<double(double a, double b)> midpoint;

## template <class T, class R>
## class FMMProposal { ... };

## FMMProposal(const std::vector<R>& regions);

## FMMProposal(const R& region);

## std::vector<T> r(unsigned int n = 1) const;                           // <1>
## std::pair<std::vector<T>, std::vector<unsigned int>>                  // <2>
## 	r_ext(unsigned int n = 1) const;
## double d(const T& x, bool normalize = true, bool log = false) const;  // <3>
## double w_major(const T& x, bool log = true) const;                    // <4>
## double d_target_unnorm(const T& x, bool log = true) const;            // <5>

## Rcpp::NumericVector xi_upper(bool log = true) const;        // <1>
## Rcpp::NumericVector xi_lower(bool log = true) const;        // <2>
## Rcpp::LogicalVector bifurcatable() const;                   // <3>
## Rcpp::NumericVector pi(bool log = false) const;             // <4>
## Rcpp::NumericVector bound_contrib(bool log = false) const;  // <5>
## double bound(bool log = false) const;                       // <6>
## double nc(bool log = false) const;                          // <7>
## unsigned int size() const;                                  // <8>

## std::set<R>::const_iterator regions_begin() const;               // <1>
## std::set<R>::const_iterator regions_end() const;
## Rcpp::NumericVector::const_iterator log_xi_upper_begin() const;  // <2>
## Rcpp::NumericVector::const_iterator log_xi_upper_end() const;
## Rcpp::NumericVector::const_iterator log_xi_lower_begin() const;  // <3>
## Rcpp::NumericVector::const_iterator log_xi_lower_end() const;
## Rcpp::LogicalVector::const_iterator bifurcatable_begin() const;  // <4>
## Rcpp::LogicalVector::const_iterator bifurcatable_end() const;

## Rcpp::NumericVector refine(
## 	const std::vector<T>& knots, // <1>
## 	bool log = true              // <2>
## );

## Rcpp::NumericVector refine(
## 	unsigned int N,                        // <1>
## 	double tol = 0,                        // <2>
## 	bool greedy = false,                   // <3>
## 	unsigned int report = fntl::uint_max,  // <4>
## 	bool log = true                        // <5>
## );

## std::vector<double> seq(double lo, double hi, unsigned int N,
## 	bool endpoints = false);

## Rcpp::DataFrame summary() const;       // <1>
## void print(unsigned int n = 5) const;  // <2>

## template <class T>
## class Region { ... };

## virtual double d_base(const T& x, bool log = false) const = 0; // <1>
## virtual std::vector<T> r(unsigned int n) const = 0;            // <2>
## virtual bool s(const T& x) const = 0;                          // <3>
## virtual double w(const T& x, bool log = true) const = 0;       // <4>
## virtual double w_major(const T& x, bool log = true) const = 0; // <5>
## virtual bool is_bifurcatable() const = 0;                      // <6>
## virtual double xi_upper(bool log = true) const = 0;            // <7>
## virtual double xi_lower(bool log = true) const = 0;            // <8>
## virtual std::string description() const = 0;                   // <9>

## class MyRegion : public Region<type>
## {
## public:
##     std::pair<MyRegion,MyRegion> bifurcate() const;              // <1>
##     std::pair<MyRegion,MyRegion> bifurcate(const type& x) const; // <2>
##     MyRegion singleton(const type& x) const;                     // <3>
##     bool operator<(const MyRegion& x) const;                     // <4>
##     bool operator==(const MyRegion& x) const;                    // <5>
##     const MyRegion& operator=(const MyRegion& x);                // <6>
## ...
## };

## RealConstRegion(
## 	double a,                                 // <1>
## 	double b,                                 // <2>
## 	const dfdb& w,                            // <3>
## 	const UnivariateHelper& helper,           // <4>
##     const optimizer& maxopt = maxopt_default, // <5>
##     const optimizer& minopt = minopt_default, // <6>
## 	const midpoint& mid = midpoint_default    // <7>
## );
## 
## RealConstRegion(
## 	double a,                                 // <1>
## 	const dfdb& w,                            // <3>
## 	const UnivariateHelper& helper            // <4>
##     const optimizer& maxopt = maxopt_default, // <5>
##     const optimizer& minopt = minopt_default, // <6>
## 	const midpoint& mid = midpoint_default    // <7>
## );

## double RealConstRegion::midpoint() const

## std::pair<RealConstRegion,RealConstRegion> bifurcate() const;
## std::pair<RealConstRegion,RealConstRegion> bifurcate(const double& x) const;

## bool is_bifurcatable() const;

## RealConstRegion singleton(const double& x) const;

## bool operator<(const RealConstRegion& x) const;
## bool operator==(const RealConstRegion& x) const;

## const RealConstRegion& operator=(const RealConstRegion& x);

## class IntConstRegion : public RealConstRegion { ... };

## IntConstRegion(
## 	double a,                                 // <1>
## 	double b,                                 // <2>
## 	const dfdb& w,                            // <3>
## 	const UnivariateHelper& helper,           // <4>
##     const optimizer& maxopt = maxopt_default, // <5>
##     const optimizer& minopt = minopt_default, // <6>
## 	const midpoint& mid = midpoint_default    // <7>
## );
## 
## IntConstRegion(
## 	double a,                                 // <1>
## 	const dfdb& w,                            // <3>
## 	const UnivariateHelper& helper            // <4>
##     const optimizer& maxopt = maxopt_default, // <5>
##     const optimizer& minopt = minopt_default, // <6>
## 	const midpoint& mid = midpoint_default    // <7>
## );

## std::pair<IntConstRegion,IntConstRegion> bifurcate() const;
## std::pair<IntConstRegion,IntConstRegion> bifurcate(const double& x) const;

## bool is_bifurcatable() const;

## IntConstRegion singleton(const double& x) const;

## 
## bool operator<(const IntConstRegion& x) const;
## bool operator==(const IntConstRegion& x) const;

## const IntConstRegion& operator=(const IntConstRegion& x);

## UnivariateHelper(
## 	const fntl::density& d,  // <1>
## 	const fntl::cdf& p,      // <2>
## 	const fntl::quantile& q  // <3>
## );

## double d(double x, bool log = false) const;                     // <1>
## double p(double q, bool lower = true, bool log = false) const;  // <2>
## double q(double p, bool lower = true, bool log = false) const;  // <3>

## template <typename T, typename R>
## rejection_result<T> rejection(
## 	const FMMProposal<T,R>& h,         // <1>
## 	unsigned int n,                    // <2>
## 	const rejection_args& args         // <3>
## );
## 
## template <typename T, typename R>
## rejection_result<T> rejection(
## 	const FMMProposal<T,R>& h,         // <1>
## 	unsigned int n                     // <2>
## );

## struct rejection_args {
## 	unsigned int max_rejects = std::numeric_limits<unsigned int>::max(); // <1>
## 	unsigned int report = std::numeric_limits<unsigned int>::max();      // <2>
## 	double ratio_ub = std::exp(1e-5);                                    // <3>
## 	fntl::error_action action = fntl::error_action::STOP;                // <4>
## 
## 	rejection_args() { };                                                // <5>
## 	rejection_args(SEXP obj);                                            // <6>
## 	operator SEXP() const;                                               // <7>
## };

## template <typename T>
## struct rejection_result
## {
## 	std::vector<T> draws;               // <1>
## 	std::vector<unsigned int> rejects;  // <2>
## 
## 	operator SEXP() const;              // <3>
## };

## optimize_hybrid_result optimize_hybrid(
## 	const fntl::dfd& f,       // <1>
##     double init,              // <2>
## 	double lower,             // <3>
## 	double upper,             // <4>
## 	bool maximize,            // <5>
## 	unsigned maxiter = 100000 // <6>
## );

## struct optimize_hybrid_result {
## 	double par;            // <1>
## 	double value;          // <2>
## 	std::string method;    // <3>
## 	int status;            // <4>
## 
## 	operator SEXP() const; // <5>
## };

## double log_sum_exp(const Rcpp::NumericVector& x);

## double log_add2_exp(double x, double y);
## 
## std::vector<double> log_add2_exp(
## 	const std::vector<double>& x,
## 	const std::vector<double>& y
## );
## Rcpp::NumericVector log_add2_exp(
## 	const Rcpp::NumericVector& x,
## 	const Rcpp::NumericVector& y
## );

## double log_sub2_exp(double x, double y);
## 
## std::vector<double> log_sub2_exp(
## 	const std::vector<double>& x,
## 	const std::vector<double>& y
## );
## Rcpp::NumericVector log_sub2_exp(
## 	const Rcpp::NumericVector& x,
## 	const Rcpp::NumericVector& y
## );

## unsigned int r_categ(
## 	const Rcpp::NumericVector& p,  // <2>
## 	bool log = false,              // <3>
## 	bool one_based = false         // <4>
## );
## Rcpp::IntegerVector r_categ(
## 	unsigned int n,                // <1>
## 	const Rcpp::NumericVector& p,  // <2>
## 	bool log = false,              // <3>
## 	bool one_based = false         // <4>
## );

## double d_gumbel(
## 	double x,                      // <1>
## 	double mu = 0,                 // <5>
## 	double sigma = 1,              // <6>
## 	bool log = false               // <7>
## );
## double p_gumbel(
## 	double q,                      // <2>
## 	double mu = 0,                 // <5>
## 	double sigma = 1,              // <6>
## 	bool lower = true,             // <8>
## 	bool log = false               // <7>
## );
## double q_gumbel(
## 	double p,                      // <3>
## 	double mu = 0,                 // <5>
## 	double sigma = 1,              // <6>
## 	bool lower = true,             // <8>
## 	bool log = false               // <7>
## );
## double r_gumbel(
## 	double mu = 0,                 // <5>
## 	double sigma = 1               // <7>
## );
## 
## Rcpp::NumericVector d_gumbel(
## 	const Rcpp::NumericVector& x,  // <1>
## 	double mu = 0,                 // <5>
## 	double sigma = 1,              // <6>
## 	bool log = false               // <7>
## );
## Rcpp::NumericVector p_gumbel(
## 	const Rcpp::NumericVector& q,  // <2>
## 	double mu = 0,                 // <5>
## 	double sigma = 1,              // <6>
## 	bool lower = true,             // <8>
## 	bool log = false               // <7>
## );
## Rcpp::NumericVector q_gumbel(
## 	const Rcpp::NumericVector& p,  // <3>
## 	double mu = 0,                 // <5>
## 	double sigma = 1,              // <6>
## 	bool lower = true,             // <8>
## 	bool log = false               // <7>
## );
## Rcpp::NumericVector r_gumbel(
## 	unsigned int n,                // <4>
## 	double mu = 0,                 // <5>
## 	double sigma = 1               // <6>
## );

## -----------------------------------------------------------------------------
#| echo: false
#| prompt: true
Rcpp::sourceCpp("examples/vmf/functions.cpp")


## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v1.cpp")
out1 = r_vmf_pre_v1(n = 1000, kappa = 5, d = 4, N = 50, tol = 0.10)
head(out1$draws)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Refinement for VMF example with constant majorizer.
#| label: fig-vmf-const-refine
#| echo: false

df = data.frame(bdd = exp(out1$lbdd)) %>%
	mutate(step = row_number() - 1)
ggplot(df) +
	geom_line(aes(x = step, y = bdd)) +
	scale_x_continuous(breaks = df$step[df$step %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("Step") +
	ylab("Bound") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid) versus target (dashed) for VMF
#|   example with constant majorizer.
#| label: fig-vmf-const-draws
#| echo: false

data.frame(x = out1$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, args = list(kappa = 5, d = 4), lty = 2) +
	ylab("Empirical Density") +
	theme_minimal()


## vws::optimizer opt1 = [&](const vws::dfdb& w, double lo, double hi, bool log)
## {
## 	double x = 0;
## 	if (hi < 0) {
## 		x = hi;
## 	} else if (lo > 0){
## 		x = lo;
## 	}
## 	double out = w(x, true);
## 	return log ? out : std::exp(out);
## };

## vws::optimizer opt2 = [&](const vws::dfdb& w, double lo, double hi, bool log)
## {
## 	double w_lo = w(lo, true);
## 	double w_hi = w(hi, true);
## 	double out = std::min(w_lo, w_hi);
## 	return log ? out : std::exp(out);
## };

## vws::optimizer* maxopt;
## vws::optimizer* minopt;
## if (d >= 3) {
## 	maxopt = &opt1;
## 	minopt = &opt2;
## } else {
## 	minopt = &opt1;
## 	maxopt = &opt2;
## }

## vws::UnivariateHelper helper(df, pf, qf);
## vws::RealConstRegion supp(-1, 1, w, helper, *maxopt, *minopt);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v2.cpp")
out2 = r_vmf_pre_v2(n = 1000, kappa = 5, d = 4, N = 50, tol = 0.10)
head(out2$draws)


## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v3.cpp")
out3 = r_vmf_pre_v3(n = 1000, kappa = 5, d = 4, N = 50, tol = 0.01)
head(out3$draws)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid) versus target (dashed) for VMF
#|   example with linear majorizer.
#| label: fig-vmf-linear-draws
#| echo: false

data.frame(x = out3$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, args = list(kappa = 5, d = 4), lty = 2) +
	ylab("Empirical Density") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Refinement for VMF example with constant majorizer.
#| label: fig-vmf-refine
#| echo: false

df1 = data.frame(sampler = "v1", bdd = exp(out1$lbdd)) %>% 
	mutate(N = row_number())
df2 = data.frame(sampler = "v2", bdd = exp(out2$lbdd)) %>% 
	mutate(N = row_number())
df3 = data.frame(sampler = "v3", bdd = exp(out3$lbdd)) %>% 
	mutate(N = row_number() + 1)
df = bind_rows(df1, df2, df3) %>% mutate(N = as.integer(N))

ggplot(df) +
	geom_line(aes(x = N, y = bdd, color = sampler), lwd = 1.02, alpha = 1) +
	geom_point(aes(x = N, y = bdd, pch = sampler, color = sampler),
		cex = 3, alpha = 1) +
	scale_x_continuous(breaks = df$N[df$N %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("N") +
	ylab("Bound") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| prompt: true
mu = 5
sigma = sqrt(0.5)
lambda = 10


## -----------------------------------------------------------------------------
#| prompt: true
source("examples/ln-norm/functions.R")
y_true = rlnorm(1, mu, sigma)
z = rnorm(1, y_true, lambda)
vws::printf("y_true = %g and z = %g\n", y_true, z)


## const vws::dfdb& w = [&](double x, bool log = true) {
## 	double out = R_NegInf;
## 	if (x > 0) {
## 		double sigma2 = std::pow(sigma, 2)
## 		out = -std::log(x) - std::pow(std::log(x) - mu, 2) / (2 * sigma2);
## 	}
## 	return log ? out : std::exp(out);
## };

## fntl::density df = [&](double x, bool log = false) {
## 	return R::dnorm(x, z, lambda, log);
## };
## fntl::cdf pf = [&](double q, bool lower = true, bool log = false) {
## 	return R::pnorm(q, z, lambda, lower, log);
## };
## fntl::quantile qf = [&](double p, bool lower = true, bool log = false) {
## 	return R::qnorm(p, z, lambda, lower, log);
## };
## 
## vws::UnivariateHelper helper(df, pf, qf);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/ln-norm/ln-norm-v1.cpp")
out1 = r_ln_norm_v1(n = 1000, z, mu, sigma, lambda, N = 50, tol = 0.10)
head(out1$draws)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Refinement for lognormal-normal example with constant majorizer.
#| label: fig-ln-norm-const-refine
#| echo: false

df = data.frame(bdd = exp(out1$lbdd)) %>%
	mutate(step = row_number() - 1)
ggplot(df) +
	geom_line(aes(x = step, y = bdd)) +
	scale_x_continuous(breaks = df$step[df$step %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("Step") +
	ylab("Bound") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid black) versus target (dashed black)
#|   for lognormal-normal example with constant majorizer. Observed value of
#|   $z$ (blue line) and latent value of $y$ (red line) are displayed along
#|   with a 95\% interval (blue ribbon) based on draws from $[Y \mid Z = z]$.
#| label: fig-ln-norm-const-draws
#| echo: false

interval_lo = quantile(out1$draws, probs = 0.025)
interval_hi = quantile(out1$draws, probs = 0.975)

data.frame(x = out1$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, lty = 2,
		args = list(mu = mu, sigma = sigma, z = z, lambda = lambda)) +
	annotate("rect", xmin = interval_lo, xmax = interval_hi, ymin = 0,
		ymax = Inf, alpha = 0.1, fill = "blue") +
	geom_vline(xintercept = z, col = "blue", lwd = 1.05) +
	geom_vline(xintercept = y_true, col = "red", lwd = 1.05) +
	ylab("Empirical Density") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| prompt: true
y_star = exp(mu - sigma^2)
w_star = -mu + sigma^2 / 2
vws::printf("Maximizer y = %g obtains value log w(%g) = %g.\n",
	y_star, y_star, w_star)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Weight function for Lognormal-Normal example on the log-scale, with the
#|   maximizer $y^* = \exp(\mu - \sigma^2)$ highlighted.
#| label: fig-ln-norm-weight
#| echo: false

xlim = y_star + c(-10,10)
w_one = function(y, log = T) {
	out = -Inf
	if (y > 0) { out = -log(y) - (log(y) - mu)^2 / (2*sigma^2) }
	if (log) { return(out) } else { return(exp(out)) }
}
w = Vectorize(w_one, vectorize.args = "y")
ggplot() +
	geom_function(fun = w, xlim = xlim) +
	geom_vline(xintercept = y_star, lty = 2) +
	geom_hline(yintercept = w_star, lty = 2) +
	xlab("y") +
	ylab(expression("log w(y)")) +
	theme_minimal()


## const vws::optimizer& maxopt = [&](const vws::dfdb& w, double lo,
## 	double hi, bool log)
## {
## 	double y_star = exp(mu - sigma2);
## 
## 	double y = y_star;
## 	if (y_star > hi) {
## 		y = hi;
## 	} else if (y_star < lo) {
## 		y = lo;
## 	}
## 
## 	double out = w(y, true);
## 	return log ? out : exp(out);
## };

## const vws::optimizer& minopt = [&](const vws::dfdb& w, double lo,
## 	double hi, bool log)
## {
## 	double lwa = w(lo, true);
## 	double lwb = w(hi, true);
## 	double out = std::min(lwa, lwb);
## 	return log ? out : exp(out);
## };

## vws::RealConstRegion supp(0, R_PosInf, w, helper, maxopt, minopt);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/ln-norm/ln-norm-v2.cpp")
out2 = r_ln_norm_v2(n = 1000, z, mu, sigma, lambda, N = 50, tol = 0.10)
head(out2$draws)


## -----------------------------------------------------------------------------
#| prompt: true
y0 = exp(mu - sigma^2 + 1)
printf("Convexity changes at y = %g.\n", y0)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Weight function for Lognormal-Normal example on the log-scale,
#|   highlighting $y_0 = \exp(\mu - \sigma^2 + 1)$ where there is a change in
#|   convexity.
#| label: fig-ln-norm-weight-convexity
#| echo: false

xlim = c(0, 2000)
ggplot() +
	geom_function(fun = w, xlim = xlim) +
	geom_vline(xintercept = y0, lty = 2) +
	xlab("y") +
	ylab(expression("log w(y)")) +
	theme_minimal()


## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/ln-norm/ln-norm-v3.cpp")
out3 = r_ln_norm_v3(n = 1000, z, mu, sigma, lambda, lo = 1e-8, 1e8,
	N = 50, tol = 0.10)
head(out3$draws)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid) versus target (dashed) for
#|   lognormal-normal example with constant majorizer.
#| label: fig-ln-norm-linear-draws
#| echo: false

data.frame(x = out1$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, lty = 2,
		args = list(mu = mu, sigma = sigma, z = z, lambda = lambda)) +
	ylab("Empirical Density") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Refinement for lognormal-normal example with constant majorizer.
#| label: fig-ln-norm-refine
#| echo: false

df1 = data.frame(sampler = "v1", bdd = exp(out1$lbdd)) %>% 
	mutate(N = row_number())
df2 = data.frame(sampler = "v2", bdd = exp(out2$lbdd)) %>% 
	mutate(N = row_number())
df3 = data.frame(sampler = "v3", bdd = exp(out3$lbdd)) %>% 
	mutate(N = row_number() + 1)
df = bind_rows(df1, df2, df3) %>% mutate(N = as.integer(N))

ggplot(df) +
	geom_line(aes(x = N, y = bdd, color = sampler), lwd = 1.02, alpha = 1) +
	geom_point(aes(x = N, y = bdd, pch = sampler, color = sampler),
		cex = 3, alpha = 1) +
	scale_x_continuous(breaks = df$N[df$N %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("N") +
	ylab("Bound") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| prompt: true
source("examples/bessel/bessel.R")
lambda = 10
nu = 2


## -----------------------------------------------------------------------------
#| echo: false
source("examples/common/plots.R")
xseq = seq(0, 15)
fseq = d_bessel(xseq, lambda, nu)
lwseq = w(xseq, log = T)


## const vws::dfdb& w = [&](double x, bool log = true) {
## 	double out = -std::lgamma(x + nu + 1);
## 	return log ? out : std::exp(out);
## };

## fntl::density df = [&](double x, bool log = false) {
## 	return R::dpois(x, mean, log);
## };
## fntl::cdf pf = [&](double q, bool lower = true, bool log = false) {
## 	return R::ppois(q, mean, lower, log);
## };
## fntl::quantile qf = [&](double p, bool lower = true, bool log = false) {
## 	return R::qpois(p, mean, lower, log);
## };
## 
## vws::UnivariateHelper helper(df, pf, qf);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/bessel/bessel-v1.cpp")
out1 = r_bessel_v1(n = 1000, lambda, nu, N = 50, tol = 0.10)
head(out1$draws)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Refinement for lognormal-normal example with constant majorizer.
#| label: fig-bessel-const-refine
#| echo: false

df = data.frame(bdd = exp(out1$lbdd)) %>%
	mutate(step = row_number() - 1)
ggplot(df) +
	geom_line(aes(x = step, y = bdd)) +
	scale_x_continuous(breaks = df$step[df$step %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("Step") +
	ylab("Bound") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (bars) versus target (points) for Bessel
#|   example with constant majorizer.
#| label: fig-bessel-const-draws
#| echo: false

plot_pmf(out1$draws) +
	geom_point(data = data.frame(x = xseq, y = fseq), aes(x,y)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0])


## const vws::optimizer& maxopt = [](const vws::dfdb& w, double lo,
## 	double hi, bool log)
## {
## 	if (lo < 0 && hi < 0) { Rcpp::stop("Did not code this case"); }
## 	double x = (lo < 0) ? 0 : std::floor(lo) + 1;
## 	double out = w(x, true);
## 	return log ? out : std::exp(out);
## };
## 
## const vws::optimizer& minopt = [](const vws::dfdb& w, double lo,
## 	double hi, bool log)
## {
## 	if (lo < 0 && hi < 0) { Rcpp::stop("Did not code this case"); }
## 	double x = std::isinf(hi) ? hi : std::floor(hi);
## 	double out = w(x, true);
## 	return log ? out : std::exp(out);
## };

## vws::IntConstRegion supp(-0.1, R_PosInf, w, helper, maxopt, minopt);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/bessel/bessel-v2.cpp")
out2 = r_bessel_v2(n = 1000, lambda, nu, N = 50, tol = 0.10)
head(out2$draws)


## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/bessel/bessel-v3.cpp")
out3 = r_bessel_v3(n = 1000, lambda, nu, lo = -0.1, hi = 1e5, N = 50,
	tol = 0.10)
head(out3$draws)


## -----------------------------------------------------------------------------
#| out-width: 60%
#| fig-cap: |
#|   Refinement for Bessel example with three majorizers.
#| label: fig-bessel-refine
#| echo: false

df1 = data.frame(sampler = "v1", bdd = exp(out1$lbdd)) %>% 
	mutate(N = row_number())
df2 = data.frame(sampler = "v2", bdd = exp(out2$lbdd)) %>% 
	mutate(N = row_number())
df3 = data.frame(sampler = "v3", bdd = exp(out3$lbdd)) %>% 
	mutate(N = row_number())
df = bind_rows(df1, df2, df3) %>% mutate(N = as.integer(N)) %>% filter(N <= 18)

ggplot(df) +
	geom_line(aes(x = N, y = bdd, color = sampler), lwd = 1.02, alpha = 1) +
	geom_point(aes(x = N, y = bdd, pch = sampler, color = sampler),
		cex = 3, alpha = 1) +
	scale_x_continuous(breaks = df$N[df$N %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("N") +
	ylab("Bound") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| fig-cap: |
#|   Majorizers (on the log-scale) for Bessel example for three sampler
#|   versions. The blue curve with is the majorizer with a blue dot marking the
#|   right endpoint of each region. The black triangle displays the value of
#|   $\log w(x)$ at integers $x \in \mathbb{N}$.
#| fig-subcap:
#|   - v1.
#|   - v2.
#|   - v3.
#| label: fig-bessel-majorizers
#| layout-ncol: 3
#| fig-width: 3
#| fig-height: 2.5
#| echo: false

ggplot(out1$df_weight) +
	geom_segment(aes(x = lo, xend = hi, y = lwmax), col = "blue") +
	geom_point(aes(x = hi, y = lwmax), col = "blue") +
	geom_point(data = data.frame(x = xseq, w = lwseq), aes(x, w), pch = 2) +
	coord_cartesian(xlim = c(NA, 15), ylim = c(min(lwseq), NA)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Logarithm of Weight") +
	theme_minimal()

ggplot(out2$df_weight) +
	geom_segment(aes(x = lo, xend = hi, y = lwmax), col = "blue") +
	geom_point(aes(x = hi, y = lwmax), col = "blue") +
	geom_point(data = data.frame(x = xseq, w = lwseq), aes(x, w), pch = 2) +
	coord_cartesian(xlim = c(NA, 15), ylim = c(min(lwseq), NA)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Logarithm of Weight") +
	theme_minimal()

out3$df_weight %>%
	mutate(w_lo = beta0_max + beta1_max * lo) %>%
	mutate(w_hi = beta0_max + beta1_max * hi) %>%
	ggplot() +
	geom_segment(aes(x = lo, xend = hi, y = w_lo, yend = w_hi), col = "blue") +
	geom_point(aes(x = hi, y = w_hi), col = "blue") +
	geom_point(data = data.frame(x = xseq, y = lwseq), aes(x, y), pch = 2) +
	coord_cartesian(xlim = c(NA, 15), ylim = c(min(lwseq), NA)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Logarithm of Weight") +
	theme_minimal()


## -----------------------------------------------------------------------------
#| fig-cap: |
#|   Proposal density (bars) for Bessel example for three sampler versions.
#|   Triangles display the values of $f(x)$, $x \in \mathbb{N}$.
#| fig-subcap:
#|   - v1.
#|   - v2.
#|   - v3.
#| label: fig-bessel-proposals
#| layout-ncol: 3
#| fig-width: 3
#| fig-height: 2
#| echo: false

out1 = r_bessel_v1(0, lambda, nu, N = 50, tol = 0.1, x = xseq)
out2 = r_bessel_v2(0, lambda, nu, N = 50, tol = 0.1, x = xseq)
out3 = r_bessel_v3(0, lambda, nu, lo = -0.1, hi = 1e5, N = 50, tol = 0.1, x = xseq)

data.frame(x = xseq, h = out1$hx) %>%
	ggplot() +
	geom_bar(aes(x, h), stat = "identity", fill = NA, col = "black") +
	geom_point(aes(x, fseq), pch = 2) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Density") +
	theme_minimal()

data.frame(x = xseq, h = out2$hx) %>%
	ggplot() +
	geom_bar(aes(x, h), stat = "identity", fill = NA, col = "black") +
	geom_point(aes(x, fseq), pch = 2) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Density") +
	theme_minimal()

data.frame(x = xseq, h = out3$hx) %>%
	ggplot() +
	geom_bar(aes(x, h), stat = "identity", fill = NA, col = "black") +
	geom_point(aes(x, fseq), pch = 2) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Density") +
	theme_minimal()