aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CMakeLists.txt6
-rw-r--r--src/centroid.sc66
-rw-r--r--src/extent-sampling.sc22
-rw-r--r--src/macros/macros.sc15
-rw-r--r--src/sample-on-ellipsoid.sc154
-rw-r--r--src/whitening.sc145
6 files changed, 408 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b4fae28..6701515 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -75,6 +75,12 @@ add_executable(integral integral.c)
target_link_libraries(integral nsmc -lm)
add_executable(spheroid spheroid.c)
target_link_libraries(spheroid nsmc)
+add_executable(centroid centroid.c)
+target_link_libraries(centroid nsmc)
+add_executable(whitening whitening.c)
+target_link_libraries(whitening nsmc -lm)
+add_executable(sample-on-ellipsoid sample-on-ellipsoid.c)
+target_link_libraries(sample-on-ellipsoid nsmc -lm)
# Build and install scheme wrapper.
if(${GUILE_FOUND})
diff --git a/src/centroid.sc b/src/centroid.sc
new file mode 100644
index 0000000..64f2df3
--- /dev/null
+++ b/src/centroid.sc
@@ -0,0 +1,66 @@
+(sc-include "macros/macros")
+
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_vector.h")
+(pre-include "extent-sampling.h")
+(pre-include "nd-random.h")
+(pre-include "oracles.h")
+
+(pre-define DIMENSION 100)
+(pre-define CUBE-EDGE 1.0)
+(pre-define ELLIPSOID-AXES-START 0.25)
+(pre-define ELLIPSOID-AXES-END 0.5)
+
+(pre-define SAMPLES-PER-ITERATION 5000)
+(pre-define ITERATIONS 10)
+
+(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
+ ((struct-get extent-oracle oracle) r x (struct-get extent-oracle params)))
+
+(define (find-centroid oracle center filename) (void extent-oracle-t gsl-vector* (const char*))
+ (with-rng rng
+ (with-vector new-center DIMENSION
+ (with-vector x DIMENSION
+ (with-file fp filename "w"
+ (fprintf fp "%g\n" (gsl-blas-dnrm2 center))
+ (for-i iteration ITERATIONS
+ (gsl-vector-set-zero new-center)
+ (for-i trial SAMPLES-PER-ITERATION
+ (random-direction-vector rng x)
+ (gsl-vector-scale x (invoke-extent-oracle oracle rng x))
+ (gsl-vector-add new-center x)
+
+ (gsl-vector-scale x (/ -1.0 (gsl-blas-dnrm2 x)))
+ (gsl-vector-scale x (invoke-extent-oracle oracle rng x))
+ (gsl-vector-add new-center x))
+ (gsl-vector-scale new-center (/ 1.0 2 SAMPLES-PER-ITERATION))
+ (gsl-vector-memcpy center new-center)
+ (fprintf fp "%g\n" (gsl-blas-dnrm2 center))))))))
+
+(define (main) (int)
+ (with-vector center DIMENSION
+ (gsl-vector-set-basis center 0)
+ (gsl-vector-scale center (* 0.25 CUBE-EDGE))
+ (declare cube-params (struct-variable cube-params
+ (edge CUBE-EDGE)
+ (center center)))
+ (declare cube-oracle (struct-variable extent-oracle-t
+ (oracle cube-extent-oracle-with-center)
+ (params (address-of cube-params))))
+ (find-centroid cube-oracle center "cube.dat")
+
+ (gsl-vector-set-basis center 0)
+ (gsl-vector-scale center (* 0.5 ELLIPSOID-AXES-START))
+ (with-vector axes DIMENSION
+ (for-i i DIMENSION
+ (gsl-vector-set axes i (+ ELLIPSOID-AXES-START (* i (/ (- ELLIPSOID-AXES-END
+ ELLIPSOID-AXES-START)
+ (- DIMENSION 1))))))
+ (declare ellipsoid-params (struct-variable ellipsoid-params
+ (axes axes)
+ (center center)))
+ (declare ellipsoid-oracle (struct-variable extent-oracle-t
+ (oracle ellipsoid-extent-oracle-with-center)
+ (params (address-of ellipsoid-params))))
+ (find-centroid ellipsoid-oracle center "ellipsoid.dat")))
+ (return 0))
diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc
index d46f31d..b74051f 100644
--- a/src/extent-sampling.sc
+++ b/src/extent-sampling.sc
@@ -52,6 +52,28 @@
(else
(set accurate-estimates 0))))))))
+;; (pre-let* (WINDOW-LENGTH 1000 COVARIANCE-SAMPLES 1000 HARDCODED-ECCENTRICITY 2)
+;; (define (volume-whitening extent-oracle true-volume r dimension rtol stats)
+;; (void (const extent-oracle-t*) double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
+;; (define accurate-estimates int 0)
+;; (let* ((vn double (ln-volume-of-ball dimension)))
+;; (with-square-matrix transform dimension
+;; (gsl-matrix-set-identity transform)
+;; ;; TODO: Replace hardcoded eccentricity.
+;; (gsl-matrix-set transform 0 0 HARDCODED-ECCENTRICITY)
+;; (with-vector x dimension
+;; (with-vector y dimension
+;; (do-while (< accurate-estimates WINDOW-LENGTH)
+;; (random-direction-vector r x)
+;; (cblas-dgemv CblasRowMajor CblasNoTrans dimension dimension 1 transform dimension x 1 0 y 1)
+;; (gsl-rstat-add (exp (+ vn (* dimension (log (gsl-blas-dnrm2 y)))))
+;; stats)
+;; (cond
+;; ((rtol? (gsl-rstat-mean stats) true-volume rtol)
+;; (set+ accurate-estimates 1))
+;; (else
+;; (set accurate-estimates 0))))))))))
+
(sc-define-syntax (invoke-integrand integrand r x)
((: integrand integrand) r x (: integrand params)))
diff --git a/src/macros/macros.sc b/src/macros/macros.sc
index 1ccfeab..4832887 100644
--- a/src/macros/macros.sc
+++ b/src/macros/macros.sc
@@ -26,6 +26,9 @@
(when (not condition)
body ...))
+(sc-define-syntax (incr var)
+ (set+ var 1))
+
(sc-define-syntax (for-i index limit body ...)
(for ((define index int 0)
(< index limit)
@@ -49,6 +52,18 @@
gsl-vector-free
body ...))
+(sc-define-syntax (with-matrix var m n body ...)
+ (with-alloc var gsl-matrix*
+ (gsl-matrix-alloc m n)
+ gsl-matrix-free
+ body ...))
+
+(sc-define-syntax (with-square-matrix var n body ...)
+ (with-alloc var gsl-matrix*
+ (gsl-matrix-alloc n n)
+ gsl-matrix-free
+ body ...))
+
(sc-define-syntax (with-rng var body ...)
(with-alloc var gsl-rng*
(gsl-rng-alloc gsl-rng-default) gsl-rng-free
diff --git a/src/sample-on-ellipsoid.sc b/src/sample-on-ellipsoid.sc
new file mode 100644
index 0000000..a3012ed
--- /dev/null
+++ b/src/sample-on-ellipsoid.sc
@@ -0,0 +1,154 @@
+(sc-include "macros/macros")
+
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_cblas.h")
+(pre-include "gsl/gsl_eigen.h")
+(pre-include "extent-sampling.h")
+(pre-include "nd-random.h")
+(pre-include "oracles.h")
+(pre-include "utils.h")
+
+(pre-define ECCENTRICITY 2)
+(pre-define DIMENSION 4)
+
+(pre-define COVARIANCE_SAMPLES 1000)
+(pre-define SAMPLES 10000)
+
+;; We are generalizing ellipsoid sampling with estimation of
+;; covariance matrix.
+;;
+;; We need a fw(alpha) for a non-diagonal but symmetric A. With that,
+;; we will know how to normalize the volume estimate. This new
+;; fw(alpha) must involve some simple rotation to the eigen basis of
+;; A.
+
+(declare covariance-accumulator
+ (type (struct (square gsl-matrix*)
+ (sum gsl-vector*)
+ (n int))))
+
+(define (covariance-accumulator-alloc dimension) (covariance-accumulator int)
+ (declare accumulator covariance-accumulator)
+ (set (struct-get accumulator square) (gsl-matrix-alloc dimension dimension))
+ (set (struct-get accumulator sum) (gsl-vector-alloc dimension))
+ (set (struct-get accumulator n) 0)
+ (return accumulator))
+
+(define (covariance-accumulator-add accumulator sample) (void covariance-accumulator* gsl-vector*)
+ (define square gsl-matrix* (: accumulator square))
+ (cblas-dsyr CblasRowMajor CblasLower (: sample size) 1
+ (: sample data) (: sample stride)
+ (: square data) (: square tda))
+ (gsl-vector-add (: accumulator sum) sample)
+ (incr (: accumulator n)))
+
+(define (covariance-accumulator-covariance accumulator result) (void covariance-accumulator gsl-matrix*)
+ (let* ((dim int (: (struct-get accumulator sum) size))
+ (n int (struct-get accumulator n))
+ (sum gsl-vector* (struct-get accumulator sum))
+ (square gsl-matrix* (struct-get accumulator square)))
+ (gsl-matrix-memcpy result square)
+ (gsl-matrix-scale result (/ 1.0 (- n 1)))
+ (cblas-dsyr CblasRowMajor CblasLower dim (/ -1.0 n n)
+ (: sum data) (: sum stride)
+ (: result data) (: result tda))))
+
+(define (covariance-accumulator-free accumulator) (void covariance-accumulator)
+ (gsl-matrix-free (struct-get accumulator square))
+ (gsl-vector-free (struct-get accumulator sum)))
+
+(sc-define-syntax (with-covariance-accumulator var dimension body ...)
+ (with-alloc var covariance-accumulator
+ (covariance-accumulator-alloc dimension)
+ covariance-accumulator-free
+ body ...))
+
+(define (matrix-vector-multiply A x y)
+ (void (const gsl-matrix*) (const gsl-vector*) gsl-vector*)
+ (cblas-dgemv CblasRowMajor CblasNoTrans (: A size1) (: A size2)
+ 1 (: A data) (: A tda)
+ (: x data) (: x stride) 0 (: y data) (: y stride)))
+
+;; (define (product-of-diagonal-entries mat)
+;; (double (const gsl-matrix*))
+;; (define result double 1)
+;; ;; TODO: Do not assume matrix is square.
+;; (for-i i (: mat size1)
+;; (set* result (gsl-matrix-get mat i i)))
+;; (return result))
+
+(define (product-of-entries vec) (double (const gsl-vector*))
+ (define result double 1)
+ (for-i i (: vec size)
+ (set* result (gsl-vector-get vec i)))
+ (return result))
+
+(sc-define-syntax (with-eigen-symmv workspace dimension body ...)
+ (with-alloc workspace gsl-eigen-symmv-workspace*
+ (gsl-eigen-symmv-alloc dimension)
+ gsl-eigen-symmv-free
+ body ...))
+
+(define (eigen-symmetric A eval evec) (void gsl-matrix* gsl-vector* gsl-matrix*)
+ (with-eigen-symmv workspace (: A size1)
+ (gsl-eigen-symmv A eval evec workspace)))
+
+(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
+ ((: extent-oracle oracle) r x (: extent-oracle params)))
+
+(define (print-vector x) (void gsl-vector*)
+ (for-i i (: x size)
+ (printf "%g\t" (gsl-vector-get x i)))
+ (printf "\n"))
+
+(define (print-matrix A) (void gsl-matrix*)
+ (for-i i (: A size1)
+ (for-i j (: A size2)
+ (printf "%g\t" (gsl-matrix-get A i j)))
+ (printf "\n")))
+
+(define (volume-by-ellipsoid extent-oracle rng dimension stats)
+ (void (const extent-oracle-t*) (const gsl-rng*) (unsigned int) gsl-rstat-workspace*)
+ (let* ((vn double (ln-volume-of-ball dimension)))
+ (with-square-matrix transform dimension
+ (with-covariance-accumulator accumulator dimension
+ (with-vector x dimension
+ (for-i i COVARIANCE-SAMPLES
+ (random-direction-vector rng x)
+ (gsl-vector-scale x (invoke-extent-oracle extent-oracle rng x))
+ (covariance-accumulator-add (address-of accumulator) x)))
+ (covariance-accumulator-covariance accumulator transform))
+ (with-vector eval dimension
+ (with-square-matrix evec dimension
+ (eigen-symmetric transform eval evec)
+ (print-vector eval)
+ (printf "\n")
+ (print-matrix evec))
+ ;; (gsl-matrix-set-identity transform)
+ ;; (gsl-matrix-set transform 0 0 (* 0.5 ECCENTRICITY))
+ (let* ((ln-determinant double (log (product-of-entries eval))))
+ (with-vector x dimension
+ (with-vector y dimension
+ (for-i i SAMPLES
+ (random-direction-vector rng x)
+ (matrix-vector-multiply transform x y)
+ (let* ((norm double (gsl-blas-dnrm2 y)))
+ (gsl-vector-scale y (/ norm))
+ (gsl-rstat-add (exp (+ vn
+ (* dimension (- (log (invoke-extent-oracle extent-oracle rng y))
+ (log norm)))
+ ln-determinant))
+ stats))))))))))
+
+(define (main) (int)
+ (declare params (struct-variable spheroid-params
+ (eccentricity ECCENTRICITY)))
+ (declare oracle (struct-variable extent-oracle-t
+ (oracle spheroid-extent-oracle)
+ (params (address-of params))))
+ (printf "True volume: %g\n" (spheroid-true-volume DIMENSION (address-of params)))
+ (with-rng rng
+ (with-rstats stats
+ (volume-by-ellipsoid (address-of oracle) rng DIMENSION stats)
+ (printf "Estimated volume: %g\n" (gsl-rstat-mean stats))))
+ (return 0))
diff --git a/src/whitening.sc b/src/whitening.sc
new file mode 100644
index 0000000..c330ed4
--- /dev/null
+++ b/src/whitening.sc
@@ -0,0 +1,145 @@
+(sc-include "macros/macros")
+
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_cblas.h")
+(pre-include "gsl/gsl_eigen.h")
+(pre-include "gsl/gsl_matrix.h")
+(pre-include "extent-sampling.h")
+(pre-include "nd-random.h")
+(pre-include "oracles.h")
+(pre-include "utils.h")
+
+(pre-define ECCENTRICITY 2)
+(pre-define DIMENSION 2)
+
+(pre-define SAMPLES 100)
+
+(declare covariance-accumulator
+ (type (struct (square gsl-matrix*)
+ (sum gsl-vector*)
+ (n int))))
+
+(define (covariance-accumulator-alloc dimension) (covariance-accumulator int)
+ (declare accumulator covariance-accumulator)
+ (set (struct-get accumulator square) (gsl-matrix-alloc dimension dimension))
+ (set (struct-get accumulator sum) (gsl-vector-alloc dimension))
+ (set (struct-get accumulator n) 0)
+ (return accumulator))
+
+(define (covariance-accumulator-add accumulator sample) (void covariance-accumulator* gsl-vector*)
+ (define square gsl-matrix* (: accumulator square))
+ (cblas-dsyr CblasRowMajor CblasLower (: sample size) 1
+ (: sample data) (: sample stride)
+ (: square data) (: square tda))
+ (gsl-vector-add (: accumulator sum) sample)
+ (set+ (: accumulator n) 1))
+
+(define (covariance-accumulator-covariance accumulator result) (void covariance-accumulator gsl-matrix*)
+ (let* ((dim int (: (struct-get accumulator sum) size))
+ (n int (struct-get accumulator n))
+ (sum gsl-vector* (struct-get accumulator sum))
+ (square gsl-matrix* (struct-get accumulator square)))
+ (gsl-matrix-memcpy result square)
+ (gsl-matrix-scale result (/ 1.0 (- n 1)))
+ (cblas-dsyr CblasRowMajor CblasLower dim (/ -1.0 n n)
+ (: sum data) (: sum stride)
+ (: result data) (: result tda))))
+
+(define (covariance-accumulator-free accumulator) (void covariance-accumulator)
+ (gsl-matrix-free (struct-get accumulator square))
+ (gsl-vector-free (struct-get accumulator sum)))
+
+(sc-define-syntax (with-covariance-accumulator var dimension body ...)
+ (with-alloc var covariance-accumulator
+ (covariance-accumulator-alloc dimension)
+ covariance-accumulator-free
+ body ...))
+
+(define (matrix-vector-multiply A x y)
+ (void (const gsl-matrix*) (const gsl-vector*) gsl-vector*)
+ (cblas-dgemv CblasRowMajor CblasNoTrans (: A size1) (: A size2)
+ 1 (: A data) (: A tda)
+ (: x data) (: x stride) 0 (: y data) (: y stride)))
+
+(sc-define-syntax (with-eigen-symmv workspace dimension body ...)
+ (with-alloc workspace gsl-eigen-symmv-workspace*
+ (gsl-eigen-symmv-alloc dimension)
+ gsl-eigen-symmv-free
+ body ...))
+
+(define (eigen-symmetric A eval evec) (void gsl-matrix* gsl-vector* gsl-matrix*)
+ (with-eigen-symmv workspace (: A size1)
+ (gsl-eigen-symmv A eval evec workspace)))
+
+(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
+ ((: extent-oracle oracle) r x (: extent-oracle params)))
+
+(define (print-vector x) (void gsl-vector*)
+ (for-i i (: x size)
+ (printf "%g\t" (gsl-vector-get x i)))
+ (printf "\n"))
+
+(define (print-matrix A) (void gsl-matrix*)
+ (for-i i (: A size1)
+ (for-i j (: A size2)
+ (printf "%g\t" (gsl-matrix-get A i j)))
+ (printf "\n")))
+
+(pre-let* (COVARIANCE-SAMPLES 10000)
+ (define (volume-whitening extent-oracle rng dimension stats)
+ (void (const extent-oracle-t*) (const gsl-rng*) (unsigned int) gsl-rstat-workspace*)
+ (let* ((vn double (ln-volume-of-ball dimension)))
+ (with-square-matrix transform dimension
+ ;; (with-covariance-accumulator accumulator dimension
+ ;; (with-vector x dimension
+ ;; (for-i i COVARIANCE-SAMPLES
+ ;; (random-direction-vector rng x)
+ ;; (gsl-vector-scale x (invoke-extent-oracle extent-oracle rng x))
+ ;; (covariance-accumulator-add (address-of accumulator) x)))
+ ;; (covariance-accumulator-covariance accumulator transform))
+ ;; (with-vector eval dimension
+ ;; (with-square-matrix evec dimension
+ ;; (eigen-symmetric transform eval evec)
+ ;; (print-vector eval)
+ ;; (printf "\n")
+ ;; (print-matrix evec)))
+ ;; TODO: Replace hardcoded eccentricity.
+ (gsl-matrix-set-identity transform)
+ (gsl-matrix-set transform 0 0 (/ 1.0 ECCENTRICITY))
+ ;; find surface vector as usual using the extent function
+ ;; transform this surface vector to a virtual body
+ ;; find the volume of this virtual body as usual
+ ;; use this volume to find the volume of the real body
+ ;;
+ ;; The idea is to transform the real body to a virtual body
+ ;; whose volume can be estimated very accurately. Then, scale
+ ;; up the virtual volume to get the real volume. Hopefully,
+ ;; the enlargement of error will be small.
+
+ ;; So, the simple determinant relation between the volumes of
+ ;; the real and virtual bodies is incorrect. What then is
+ ;; correct?
+ (with-vector x dimension
+ (with-vector y dimension
+ (for-i i SAMPLES
+ (random-direction-vector rng x)
+ (gsl-vector-scale x (invoke-extent-oracle extent-oracle rng x))
+ (matrix-vector-multiply transform x y)
+ (gsl-rstat-add (exp (+ vn (* dimension (log (gsl-blas-dnrm2 y)))))
+ stats)
+ (printf "%g (%g/%g)\t" (gsl-rstat-mean stats)
+ (gsl-blas-dnrm2 x)
+ (gsl-blas-dnrm2 y)))))))
+ (printf "\n")))
+
+(define (main) (int)
+ (declare params (struct-variable spheroid-params
+ (eccentricity ECCENTRICITY)))
+ (declare oracle (struct-variable extent-oracle-t
+ (oracle spheroid-extent-oracle)
+ (params (address-of params))))
+ (printf "True volume: %g\n" (spheroid-true-volume DIMENSION (address-of params)))
+ (with-rng rng
+ (with-rstats stats
+ (volume-whitening (address-of oracle) rng DIMENSION stats)))
+ (return 0))