diff options
| -rw-r--r-- | CMakeLists.txt | 6 | ||||
| -rw-r--r-- | src/centroid.sc | 66 | ||||
| -rw-r--r-- | src/extent-sampling.sc | 22 | ||||
| -rw-r--r-- | src/macros/macros.sc | 15 | ||||
| -rw-r--r-- | src/sample-on-ellipsoid.sc | 154 | ||||
| -rw-r--r-- | src/whitening.sc | 145 |
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)) |