#include <math.h>
#include <stdio.h>

#include <gsl/gsl_randist.h>

#include "extent-sampling.h"
#include "oracles.h"
#include "utils.h"

#define RTOLS 3
#define RTOL1 0.2
#define RTOL2 0.1
#define RTOL3 0.05

typedef double (*true_volume_t) (unsigned int);
gsl_rng* r;

#define MAX_EXTENT 1.0
double uniform_oracle (const gsl_vector* x) {
  return gsl_ran_flat(r, 0, MAX_EXTENT);
}

double uniform_true_volume_dim (unsigned int dim)
{
  return uniform_true_volume(0, MAX_EXTENT, dim);
}
#undef MAX_EXTENT

#define ALPHA 2
#define BETA 2
double beta_oracle (const gsl_vector* x) {
  return beta_extent_generator(r, ALPHA, BETA);
}

double beta_true_volume_dim (unsigned int dim)
{
  return beta_true_volume(ALPHA, BETA, dim);
}
#undef ALPHA
#undef BETA

double arcsine_oracle (const gsl_vector* x) {
  return beta_extent_generator(r, 0.5, 0.5);
}

double arcsine_true_volume_dim (unsigned int dim)
{
  return beta_true_volume(0.5, 0.5, dim);
}

#define EDGE 1.0
double cube_oracle (const gsl_vector* x) {
  return cube_extent_oracle(x, EDGE);
}

double cube_true_volume_dim (unsigned int dim)
{
  return cube_true_volume(EDGE, dim);
}
#undef EDGE

void run_experiment (extent_oracle_t oracle, true_volume_t true_volume, const char* filename_prefix,
                     int dim_start, int dim_step, int dim_stop)
{
  gsl_rstat_workspace* stats = gsl_rstat_alloc();
  printf("dimension\trtol=%g\trtol=%g\trtol=%g\n", RTOL1, RTOL2, RTOL3);
  double rtol[RTOLS] = {RTOL1, RTOL2, RTOL3};
  FILE *fp[RTOLS];
  for (int i=0; i<RTOLS; i++) {
    char *filename;
    asprintf(&filename, "%s-%g.dat", filename_prefix, rtol[i]);
    fp[i] = fopen(filename, "w");
    fprintf(fp[i], "dimension\tsamples\n");
    free(filename);
  }
  for (int dim=dim_start; dim<=dim_stop; dim+=dim_step) {
    printf("%d\t", dim);
    for (int i=0; i<RTOLS; i++) {
      volume(oracle, true_volume(dim), r, dim, rtol[i], stats);
      fprintf(fp[i], "%d\t%zd\n", dim, gsl_rstat_n(stats));
      printf("%zd\t", gsl_rstat_n(stats));
    }
    printf("\n");
    gsl_rstat_reset(stats);
  }
  for (int i=0; i<RTOLS; i++) fclose(fp[i]);
  gsl_rstat_free(stats);
}

int main (int argc, char** argv)
{
  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  printf("# uniform\n");
  run_experiment(uniform_oracle, uniform_true_volume_dim, "volume-uniform", 5, 5, 100);
  printf("# beta\n");
  run_experiment(beta_oracle, beta_true_volume_dim, "volume-beta", 5, 5, 100);
  printf("# arcsine\n");
  run_experiment(arcsine_oracle, arcsine_true_volume_dim, "volume-arcsine", 5, 5, 100);
  printf("# cube\n");
  run_experiment(cube_oracle, cube_true_volume_dim, "volume-cube", 10, 10, 50);

  gsl_rng_free(r);
  return 0;
}