aboutsummaryrefslogtreecommitdiff
path: root/src/oracles.c
blob: 9f6577f85ce5870ec2bfc8bac51493692dc69021 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include "oracles.h"
#include "utils.h"

double bernoulli_extent_generator
(const gsl_rng* r, double p, double r0, double r1)
{
  return gsl_ran_bernoulli(r, p) ? r1 : r0;
}

double bernoulli_true_volume (double p, double r0, double r1, unsigned int dimension)
{
  return volume_of_ball(dimension)
    * (p * gsl_pow_uint(r1, dimension)
       + (1 - p) * gsl_pow_uint(r0, dimension));
}

double uniform_extent_generator (const gsl_rng* r, double a, double b)
{
  return gsl_ran_flat(r, a, b);
}

// TODO: Verify the accuracy of this function for non-trivial a, b
double uniform_true_volume (double a, double b, unsigned int dimension)
{
  return exp(ln_volume_of_ball(dimension) + dimension*log(b) - log(dimension + 1))
    - exp(ln_volume_of_ball(dimension) + dimension*log(a) - log(dimension + 1));
}

double beta_extent_generator
(const gsl_rng* r, double alpha, double beta)
{
  return gsl_ran_beta(r, alpha, beta);
}

double beta_true_volume (double alpha, double beta, unsigned int dimension)
{
  double vol = volume_of_ball(dimension);
  for (unsigned int r=0; r<dimension; r++)
    vol *= (alpha + r) / (alpha + beta + r);
  return vol;
}

double infinity_norm (const gsl_vector* x)
{
  double max = fabs(gsl_vector_get(x, 0));
  for (int i=1; i<x->size; i++)
    max = GSL_MAX(max, fabs(gsl_vector_get(x, i)));
  return max;
}

double cube_extent_oracle (const gsl_vector* x, double edge)
{
  return edge / 2 / infinity_norm(x);
}

double cube_extent_oracle_with_center
(const gsl_vector* x, const gsl_vector* center, double edge)
{
  double min = (0.5*edge - GSL_SIGN(gsl_vector_get(x, 0))*gsl_vector_get(center, 0))
    / fabs(gsl_vector_get(x, 0));
  for (int i=1; i<center->size; i++)
    min = GSL_MIN(min, (0.5*edge - GSL_SIGN(gsl_vector_get(x, i))*gsl_vector_get(center, i))
                  / fabs(gsl_vector_get(x, i)));
  return min;
}

double cube_true_volume (double edge, unsigned int dimension)
{
  return gsl_pow_uint(edge, dimension);
}

double cube_maximum_extent (double edge, unsigned int dimension)
{
  return edge * sqrt(dimension) / 2;
}

double ellipsoid_extent_oracle (const gsl_vector* x, const gsl_vector* axes)
{
  double k = 0;
  for (int i=0; i<axes->size; i++)
    k += gsl_pow_2(gsl_vector_get(x, i) / gsl_vector_get(axes, i));
  return 1/sqrt(k);
}

double ellipsoid_true_volume (const gsl_vector* axes)
{
  unsigned int dimension = axes->size;
  double vol = volume_of_ball(dimension);
  for (int i=0; i<dimension; i++)
    vol *= gsl_vector_get(axes, i);
  return vol;
}

double spheroid_extent_oracle (const gsl_vector* x, double eccentricity)
{
  unsigned int dimension = x->size;
  gsl_vector_const_view xsub = gsl_vector_const_subvector(x, 1, dimension - 1);
  return 1/sqrt(gsl_pow_2(gsl_blas_dnrm2(&xsub.vector))
                + gsl_pow_2(gsl_vector_get(x, 0) / eccentricity));
}

double spheroid_true_volume (double eccentricity, unsigned int dimension)
{
  return volume_of_ball(dimension) * eccentricity;
}