## -----------------------------------------------------------------------------
#| 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::fmm_proposal<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>
## using FMMProposal = fmm_proposal<T,R>;
## 
## template<class T>
## using Region = region<T>;
## 
## using RealConstRegion = real_const_region;
## using IntConstRegion = int_const_region;
## using UnivariateHelper = univariate_helper;

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

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

## fmm_proposal(const R& region)

## fmm_proposal(const fmm_proposal& p)

## 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::IntegerVector mergeable(unsigned int i) const;        // <4>
## Rcpp::NumericVector pi(bool log = false) const;             // <5>
## Rcpp::NumericVector bound_contrib(bool log = false) const;  // <6>
## double bound(bool log = false) const;                       // <7>
## double nc(bool log = false) const;                          // <8>
## unsigned int size() const;                                  // <9>

## std::set<R>::const_iterator begin() const;
## std::set<R>::const_iterator end() const;
## std::set<R>::iterator begin() const;
## std::set<R>::iterator 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>
## )

## bool is_bifurcatable() const;                             // <1>
## std::pair<MyRegion,MyRegion> bifurcate() const;           // <2>
## std::pair<MyRegion,MyRegion> bifurcate(const T& x) const; // <3>

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

## double merge(unsigned int i, unsigned int j, bool log = true)

## 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 double xi_upper(bool log = true) const = 0;            // <6>
## virtual double xi_lower(bool log = true) const = 0;            // <7>
## virtual double bound(bool log) const                           // <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>
##     virtual bool is_bifurcatable() const;                        // <3>
##     virtual bool is_mergeable(const region<T>& x) const;         // <4>
## 
##     MyRegion singleton(const type& x) const;                     // <5>
##     bool operator<(const MyRegion& x) const;                     // <6>
##     bool operator==(const MyRegion& x) const;                    // <7>
##     const MyRegion& operator=(const MyRegion& x);                // <8>
## ...
## };

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

## double real_const_region::midpoint() const

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

## bool is_bifurcatable() const;

## real_const_region merge(const real_const_region& x) const;

## bool is_mergeable(const real_const_region& x) const;

## void set_w(const dfdb& w);                         // <1>
## void set_helper(const univariate_helper& helper);  // <2>
## void set_maxopt(const optimizer& maxopt);          // <3>
## void set_minopt(const optimizer& minopt);          // <4>
## void set_mid(const vws::midpoint& mid);            // <5>
## void init();                                       // <6>

## real_const_region singleton(const double& x) const

## bool operator<(const real_const_region& x) const
## bool operator==(const real_const_region& x) const

## const real_const_region& operator=(const real_const_region& x)

## class int_const_region : public real_const_region { ... };

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

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

## bool is_bifurcatable() const;

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

## 
## bool operator<(const int_const_region& x) const
## bool operator==(const int_const_region& x) const

## const int_const_region& operator=(const int_const_region& x)

## univariate_helper(
## 	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 fmm_proposal<T,R>& h,        // <1>
## 	unsigned int n,                    // <2>
## 	const rejection_args& args         // <3>
## )
## 
## template <typename T, typename R>
## rejection_result<T> rejection(
## 	const fmm_proposal<T,R>& h,        // <1>
## 	unsigned int n = 1                 // <2>
## )

## template <typename T, typename R>
## rejection_result<T> rejection_tune(
## 	fmm_proposal<T,R>& h,       // <1>
## 	unsigned int n,             // <2>
## 	const rejection_args& args  // <3>
## )
## 
## template <typename T, typename R>
## rejection_result<T> rejection_tune(
## 	const fmm_proposal<T,R>& h, // <1>
## 	unsigned int n = 1          // <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>
## 	double tol_suff = 0.85;                                              // <5>
## 	double tol_merge = 0.01;                                             // <6>
## 	bool metrics = true;                                                 // <7>
## 
## 	rejection_args() { };                                                // <8>
## 	rejection_args(SEXP obj);                                            // <9>
## 	operator SEXP() const;                                               // <10>
## };

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

## 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::univariate_helper helper(df, pf, qf)
## vws::real_const_region 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/privacy/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::univariate_helper helper(df, pf, qf);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/privacy/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-privacy-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-privacy-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-privacy-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::real_const_region supp(0, R_PosInf, w, helper, maxopt, minopt);

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/privacy/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-privacy-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/privacy/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-privacy-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-privacy-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::univariate_helper 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::int_const_region 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 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()


## typedef vws::fmm_proposal<double, vws::real_const_region> t_proposal;
## typedef Rcpp::XPtr<t_proposal> t_proposal_xptr;

## // [[Rcpp::export]]
## t_proposal_xptr vmf_pre_v1_xptr(double kappa, double d)                  // <1>
## {
## 	vws::dfdb w =
## 	[=](double x, bool log = true) {                                     // <2>
## 		double out = R_NegInf;
## 		if (std::fabs(x) < 1){
## 			out = 0.5 * (d - 3) * std::log1p(-std::pow(x, 2));
## 		}
## 		return log ? out : std::exp(out);
## 	};
## 
## 	fntl::density df = [=](double x, bool log = false) {
## 		return d_texp(x, kappa, -1, 1, log);
## 	};
## 
## 	fntl::cdf pf = [=](double q, bool lower = true, bool log = false) {
## 		return p_texp(q, kappa, -1, 1, lower, log);
## 	};
## 
## 	fntl::quantile qf = [=](double p, bool lower = true, bool log = false) {
## 		return q_texp(p, kappa, -1, 1, lower, log);
## 	};
## 
## 	vws::univariate_helper helper(df, pf, qf);
## 	vws::real_const_region supp(-1, 1, w, helper);
## 
## 	auto p = new t_proposal(supp);                                       // <3>
## 	auto out = t_proposal_xptr(p, true);
## 	set_tag(out, "t_proposal");                                          // <4>
## 	return out;
## }

## // [[Rcpp::export]]
## Rcpp::NumericVector refine(
## 	t_proposal_xptr h,                    // <1>
## 	unsigned int N, double tol = 0)
## {
## 	assert_tag(h, "t_proposal");          // <2>
##     return h->refine(N - 1, tol);         // <3>
## }
## 
## // [[Rcpp::export]]
## Rcpp::List draw(t_proposal_xptr h, unsigned int n,
## 	unsigned int max_rejects = 10000, unsigned int report = 10000)
## {
## 	assert_tag(h, "t_proposal");
## 
##     vws::rejection_args args;
##     args.max_rejects = max_rejects;
##     args.report = report;
## 
##     auto out = vws::rejection(*h, n, args);
## 
##     return Rcpp::List::create(
##         Rcpp::Named("draws") = out.draws,
##         Rcpp::Named("rejects") = out.rejects
##     );
## }

## template <typename T>
## void assert_tag(Rcpp::XPtr<T> x, const char* tag)
## {
## 	if (Rf_isNull(R_ExternalPtrTag(x))) {
## 		Rcpp::stop("XPtr does not have a tag");
## 	}
## 
## 	if (Rcpp::as<Rcpp::String>(R_ExternalPtrTag(x)) != Rcpp::String(tag)) {
## 		Rcpp::stop("XPtr is not tagged as %s", tag);
## 	}
## }
## 
## template <typename T>
## void set_tag(Rcpp::XPtr<T> x, const char* tag)
## {
## 	R_SetExternalPtrTag(x, Rcpp::wrap(Rcpp::String(tag)));
## }

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v1-xptr.cpp")
h = vmf_pre_v1_xptr(kappa = 5, d = 4)
lbdd = refine(h, N = 50, tol = 0.10)
out = draw(h, n = 1000, max_rejects = 10000)
rm(h)
print(lbdd)
head(out$draws)


## // [[Rcpp::export]]
## Rcpp::XPtr<std::vector<double>> test_vector()
## {
## 	std::vector<double>* p = new std::vector<double>();
## 	p->push_back(10);
## 	p->push_back(10);
## 	p->push_back(12);
## 	return Rcpp::XPtr<std::vector<double>>(p);
## }

## -----------------------------------------------------------------------------
#| prompt: false
tryCatch({
	h = test_vector()
	refine(h, N = 50, tol = 0.10)
}, error = function(e) {
	print(e)
})


## typedef vws::fmm_proposal<double, vws::real_const_region> t_proposal;

## class rcpp_proposal : t_proposal
## {
## public:
## 	rcpp_proposal(double kappa, double d)
## 	: t_proposal(supp(kappa, d))
## 	{
## 	}
## 
## 	Rcpp::List draw(unsigned int n, unsigned int max_rejects = 10000,
## 		unsigned int report = 10000);
## 
## 	Rcpp::NumericVector refine(unsigned int N, double tol = 0) {  // <1>
## 		return t_proposal::refine(N, tol);
## 	}
## 
## private:
## 	vws::real_const_region supp(double kappa, double d);            // <2>
## };

## RCPP_MODULE(vws_module) {
## 	Rcpp::class_<rcpp_proposal>("rcpp_proposal")
## 	.constructor<double,double>()
## 	.method("draw", &rcpp_proposal::draw)
## 	.method("refine", &rcpp_proposal::refine)
## 	;
## }

## -----------------------------------------------------------------------------
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v1-module.cpp")
h = new(rcpp_proposal, kappa = 5, d = 4)
lbdd = h$refine(N = 50 - 1, tol = 0.10)
out = h$draw(n = 1000, max_rejects = 1000, max_rejects = 10000)
print(lbdd)
head(out$draws)


## Rcpp::List draw_ln_norm(
## 	double z,
## 	double lambda,
## 	double xbeta,
## 	double sigma,
## 	double tol,
## 	unsigned int N,
## 	unsigned int max_rejects
## );

## class ln_norm_proposal : public vws::fmm_proposal<double, vws::real_const_region>
## {
## public:
## 	ln_norm_proposal(double z, double lambda, double xbeta, double sigma); // <1>
## 	void update(double xbeta, double sigma);                               // <2>
## 	vws::rejection_result<double> draw(unsigned int max_rejects) const;    // <3>
## 	vws::rejection_result<double> draw_tune(double tol_suff,
## 		double tol_merge, unsigned int max_rejects);                       // <4>
## 	unsigned int size() const;                                             // <5>
## };

## RCPP_MODULE(vws_module) {
## 	Rcpp::class_<ln_norm_proposal>("ln_norm_proposal")
## 	.constructor<double,double,double,double>()
## 	.method("draw", &ln_norm_proposal::draw)
## 	.method("draw_tune", &ln_norm_proposal::draw_tune)
## 	.method("update", &ln_norm_proposal::update)
## 	.method("size", &ln_norm_proposal::size)
## 	;
## }

## -----------------------------------------------------------------------------
#| eval: false
# inner = control_inner(tol_suff = 0.85, method = "vws-basic")
# control = control_gibbs(R = 3000, burn = 1000, report = 100, inner = inner,
# 	save_latent = 1:10)
# out = gibbs(z, lambda, X, control = control)


## -----------------------------------------------------------------------------
#| eval: false
# inner = control_inner(tol_suff = 0.85, tol_merge = 0.01, tune = 100,
# 	method = "vws-tune")
# control = control_gibbs(R = 3000, burn = 1000, report = 100, inner = inner,
# 	save_latent = 1:10)
# out = gibbs(z, lambda, X, control = control)


## class ln_norm_proposal : public vws::fmm_proposal<double, vws::real_const_region>
## {
## public:
## 	ln_norm_proposal(double z, double lambda, double xbeta, double sigma); // <1>
## 	void update(double xbeta, double sigma);                               // <2>
## };