/*
  This file contains all the core extent sampling routines. From time
  to time as research progresses, I keep refactoring this code to be
  clearer and more efficient. This mostly involves throwing out old
  code even when such old code may be required in the future. I
  figured that keeping old code around is not worth the trouble
  because it impedes clarity of thought.
*/

#include <math.h>
#include <stddef.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_rstat.h>
#include <gsl/gsl_randist.h>
#include "extent-sampling.h"
#include "nd-random.h"
#include "utils.h"

#define VOLUME_MINIMUM_SAMPLES 100
#define INTEGRAL_MINIMUM_SAMPLES 100
#define CONFIDENCE_INTERVAL_FACTOR 1.96

#define INTEGRATION_INTERVALS 1000

#define SOLID_ANGLE_START 0.0
#define SOLID_ANGLE_LAST 0.5
#define SOLID_ANGLE_STEPS 100

double volume (extent_oracle_t extent_oracle, double true_volume,
               const gsl_rng* r, unsigned int dimension, double rtol,
               gsl_rstat_workspace* stats)
{
  gsl_vector* x = gsl_vector_alloc(dimension);
  double volume, error;
  double vn = ln_volume_of_ball(dimension);
  do {
    random_direction_vector(r, x);
    double extent = extent_oracle(x);
    double volume_from_sample = exp(vn + dimension*log(extent));
    gsl_rstat_add(volume_from_sample, stats);
    volume = gsl_rstat_mean(stats);
    error = CONFIDENCE_INTERVAL_FACTOR * gsl_rstat_sd_mean(stats) / volume;
  } while ((error > rtol) || (rerror(volume, true_volume) > rtol) || (gsl_rstat_n(stats) < VOLUME_MINIMUM_SAMPLES));
  gsl_vector_free(x);
  return volume;
}

double volume_window (extent_oracle_t extent_oracle, double true_volume,
                      const gsl_rng* r, unsigned int dimension, double rtol,
                      unsigned int* number_of_samples)
{
  gsl_rstat_workspace* stats = gsl_rstat_alloc();
  gsl_vector* x = gsl_vector_alloc(dimension);
  double volume, error;
  double vn = ln_volume_of_ball(dimension);
  // This is the window length used in Volume.m of Lovasz-Vempala's
  // code.
  int window_length = 4*dimension*dimension + 500;
  int accurate_estimates = 0;
  do {
    random_direction_vector(r, x);
    double extent = extent_oracle(x);
    double volume_from_sample = exp(vn + dimension*log(extent));
    gsl_rstat_add(volume_from_sample, stats);
    volume = gsl_rstat_mean(stats);
    error = rerror(volume, true_volume);
    if (error < rtol) accurate_estimates++;
    else accurate_estimates = 0;
  } while (accurate_estimates < window_length);
  *number_of_samples = gsl_rstat_n(stats);
  gsl_vector_free(x);
  gsl_rstat_free(stats);
  return volume;
}

double integral_per_direction
(integrand_t integrand, const gsl_vector* direction, const gsl_rng* r, unsigned int n,
 double radius, double rtol, int* neval)
{
  double f (double r, void* params)
  {
    (*neval)++;
    return gsl_pow_uint(r, n-1) * integrand(r, direction);
  }

  gsl_function gsl_f = {&f};
  double result, error;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(INTEGRATION_INTERVALS);
  gsl_integration_qag(&gsl_f, 0, radius, 0, rtol, INTEGRATION_INTERVALS, GSL_INTEG_GAUSS15, w, &result, &error);
  gsl_integration_workspace_free(w);
  return surface_area_of_ball(n) * result;
}

double integral
(integrand_t integrand, extent_oracle_t extent_oracle, double true_integral,
 const gsl_rng* r, unsigned int dimension, double rtol,
 gsl_rstat_workspace* stats)
{
  gsl_vector* x = gsl_vector_alloc(dimension);
  double integral, error;
  do {
    int neval = 0;
    random_direction_vector(r, x);
    double extent = extent_oracle(x);
    double integral_from_sample = integral_per_direction(integrand, x, r, dimension, extent, rtol, &neval);
    gsl_rstat_add(integral_from_sample, stats);
    integral = gsl_rstat_mean(stats);
    error = CONFIDENCE_INTERVAL_FACTOR * gsl_rstat_sd_mean(stats) / integral;
  } while ((error > rtol) || (rerror(integral, true_integral) > rtol) || (gsl_rstat_n(stats) < INTEGRAL_MINIMUM_SAMPLES));
  if (error > rtol)
    GSL_ERROR_VAL("integral failed to converge", GSL_ETOL, integral);
  gsl_vector_free(x);
  return integral;
}

double volume_cone (extent_oracle_t extent_oracle, const gsl_rng* r,
                    const gsl_vector* mean, double omega_min, double omega_max,
                    unsigned int number_of_samples, double* variance)
{
  unsigned int n = mean->size;
  gsl_vector* x = gsl_vector_alloc(n);
  double theta_min = solid_angle_to_planar_angle(omega_min * surface_area_of_ball(n), n);
  double theta_max = solid_angle_to_planar_angle(omega_max * surface_area_of_ball(n), n);
  double vn = ln_volume_of_ball(n);
  gsl_rstat_workspace* stats = gsl_rstat_alloc();
  for (int i=0; i<number_of_samples; i++) {
    hollow_cone_random_vector(r, mean, theta_min, theta_max, x);
    double extent = extent_oracle(x);
    double volume_from_sample = exp(vn + n*log(extent));
    gsl_rstat_add(volume_from_sample, stats);
  }

  if (variance) *variance = gsl_rstat_variance(stats);

  gsl_rstat_free(stats);
  gsl_vector_free(x);
  return gsl_rstat_mean(stats) * (omega_max - omega_min);
}

double volume_experiment
(extent_oracle_t extent_oracle, const gsl_rng* r,
 const gsl_vector* mean, unsigned int samples_per_cone,
 double solid_angle_factor, double solid_angle_threshold_exponent_factor,
 unsigned int* number_of_samples)
{
  int n = mean->size;
  double volume = 0;
  double omega_max = 0.5;
  int N = 0;
  do {
    double omega_min = omega_max / solid_angle_factor;
    volume += volume_cone(extent_oracle, r, mean, omega_min, omega_max, samples_per_cone, NULL);
    N += samples_per_cone;
    omega_max = omega_min;
  } while (omega_max > pow(2, -solid_angle_threshold_exponent_factor*n));
  volume += volume_cone(extent_oracle, r, mean, 0, omega_max, samples_per_cone, NULL);
  N += samples_per_cone;
  if (number_of_samples) *number_of_samples = N;
  return volume;
}