diff options
| -rw-r--r-- | include/oracles.h | 2 | ||||
| -rw-r--r-- | src/oracles.sc | 24 |
2 files changed, 26 insertions, 0 deletions
diff --git a/include/oracles.h b/include/oracles.h index 6bd4d39..2ae8cc9 100644 --- a/include/oracles.h +++ b/include/oracles.h @@ -57,9 +57,11 @@ double cube_true_volume (unsigned int dimension, void *params); typedef struct { const gsl_vector* axes; + const gsl_vector* center; } ellipsoid_params; double ellipsoid_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params); +double ellipsoid_extent_oracle_with_center (const gsl_rng *r, const gsl_vector *x, void *params); double ellipsoid_true_volume (unsigned int dimension, void *params); typedef struct { diff --git a/src/oracles.sc b/src/oracles.sc index 6f0ecbf..689ac9d 100644 --- a/src/oracles.sc +++ b/src/oracles.sc @@ -21,6 +21,7 @@ (pre-include "math.h") (pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_poly.h") (pre-include "gsl/gsl_randist.h") (pre-include "oracles.h") (pre-include "utils.h") @@ -101,6 +102,29 @@ (gsl-vector-get (: params axes) i))))) (return (/ (sqrt k))))) +(define (ellipsoid-extent-oracle-with-center r x -params) (double (const gsl-rng*) (const gsl-vector*) void*) + (let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*)) + ;; a, b and c are the coefficients of a quadratic equation + ;; ax^2 + bx + c = 0. + (a double 0) + (b double 0) + (c double 0)) + (for-i i (: (: params axes) size) + (set+ a (gsl-pow-2 (/ (gsl-vector-get x i) + (gsl-vector-get (: params axes) i)))) + (set+ b (/ (* (gsl-vector-get x i) + (gsl-vector-get (: params center) i)) + (gsl-pow-2 (gsl-vector-get (: params axes) i)))) + (set+ c (gsl-pow-2 (/ (gsl-vector-get (: params center) i) + (gsl-vector-get (: params axes) i))))) + (set* b 2) + (set- c 1) + (declare k1 double + k2 double) + (gsl-poly-solve-quadratic a b c (address-of k1) (address-of k2)) + ;; k2 > k1, k1 < 0, k2 > 0. k2 is our desired root. + (return k2))) + (define (ellipsoid-true-volume dimension -params) (double (unsigned int) void*) (let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*)) (vol double (volume-of-ball (: (: params axes) size)))) |