#include <Rcpp.h>
#include "../include/brownian_bridge_min.hpp"
#include "../include/inverse_gauss.hpp"

using namespace Rcpp;

// [[Rcpp::plugins(cpp17)]]

// [[Rcpp::export]]
double M_function(const double &a, const double &x, const double &y, const double &s, const double &t) {
  // M function that is used to simulate a minimum of a Brownian bridge
  return exp(-2.0 * (a-x) * (a-y) / (t-s));
}

// [[Rcpp::export]]
std::array<double, 2> min_sampler(const double &x, const double &y,
                                  const double &s, const double &t,
                                  const double &low_bound, const double &up_bound,
                                  std::mt19937 &RNG)
{
  // function simulates a minimum of a Brownian bridge between (low_bound) and (up_bound)
  // first element returned is the simulated minimum
  // second element returned is the simulated time which the minimum occurs

  // calculate bounds for u1
  double low_M {M_function(low_bound, x, y, s, t)};
  double up_M {M_function(up_bound, x, y, s, t)};

  // simulate uniform random variables
  std::uniform_real_distribution<double> u1(low_M, up_M);
  std::uniform_real_distribution<double> u2(0.0, 1.0);

  // set simulated minimum value
  double min {x - (0.5*(sqrt((y-x)*(y-x) - 2.0*(t-s)*log(u1(RNG))) - y + x))};

  // condition for setting V
  double condition {(x-min)/(x+y-(2.0*min))};
  // simulating from Inverse Gaussian
  double mu, lambda, V;
  if (u2(RNG) < condition) {
    mu = (y-min)/(x-min);
    lambda = (y-min)*(y-min)/(t-s);
    V = inv_gauss_sampler(mu, lambda, RNG);
  } else {
    mu = (x-min)/(y-min);
    lambda = (x-min)*(x-min)/(t-s);
    V = 1.0 / inv_gauss_sampler(mu, lambda, RNG);
  }

  // set tau (time of simualted minimum)
  double tau{((s*V)+t)/(1.0+V)};
  // setting simulated minimum and tau in array
  std::array<double, 2> simulated_min {min, tau};

  return simulated_min;
}