diff options
| -rw-r--r-- | CMakeLists.txt | 5 | ||||
| -rw-r--r-- | src/integral.sc | 140 |
2 files changed, 144 insertions, 1 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 4e363c9..e56757d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,13 +61,16 @@ endforeach() include_directories("include") # Build and install shared library with headers. -add_library(nsmc SHARED ${C_SOURCES}) +add_library(nsmc SHARED extent-sampling.c gaussian-nd-random.c integrands.c nd-random.c oracles.c utils.c) target_link_libraries(nsmc -lgsl -lgslcblas) set_target_properties(nsmc PROPERTIES PUBLIC_HEADER "include/extent-sampling.h;include/gaussian-nd-random.h;include/nd-random.h;include/oracles.h") install(TARGETS nsmc LIBRARY PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/nsmc) +# Build executables +add_executable(integral integral.c) +target_link_libraries(integral nsmc -lm) # Build and install scheme wrapper. if(${GUILE_FOUND}) configure_file(scm/nsmc/load-libs.scm.in scm/nsmc/load-libs.scm) diff --git a/src/integral.sc b/src/integral.sc new file mode 100644 index 0000000..fa3d7a6 --- /dev/null +++ b/src/integral.sc @@ -0,0 +1,140 @@ +(sc-include "macros/macros") + +(pre-include "gsl/gsl_poly.h") +(pre-include "gsl/gsl_sf.h") +(pre-include "extent-sampling.h") +(pre-include "oracles.h") +(pre-include "integrands.h") +(pre-include "utils.h") + +(pre-define TRIALS 1000) +(pre-define MAX-EXTENT 1.0) + +(sc-define-syntax* (define-polynomial var) + (let ((coefficients (list -0.09375 0.6875 -1.5 1)) + (coeff (sc-gensym))) + `(begin + (declare ,coeff (array double ,(length coefficients))) + (array-set* ,coeff ,@coefficients) + (declare ,var (struct-variable polynomial-integrand-params + (coefficients ,coeff) + (degree ,(1- (length coefficients)))))))) + +(sc-define-syntax (polynomial-map! polynomial i coeff form) + (for-i i (+ 1 (struct-get polynomial degree)) + (let* ((coeff double (array-get (struct-get polynomial coefficients) + i))) + (array-set (convert-type (struct-get polynomial coefficients) + double*) + i form)))) + +(define (polyval polynomial value) (double polynomial-integrand-params double) + (return (gsl-poly-eval (struct-get polynomial coefficients) + (+ 1 (struct-get polynomial degree)) + value))) + +(sc-define-syntax* (define-polynomial-integrand integrand) + (let ((params (sc-gensym))) + `(begin + (define-polynomial ,params) + (declare ,integrand (struct-variable (const integrand-t) + (integrand polynomial-integrand) + (params (address-of ,params))))))) + +(sc-define-syntax* (define-uniform-extent-oracle oracle a b) + (let ((oracle-params (sc-gensym))) + `(begin + (declare ,oracle-params (struct-variable uniform-params + (a ,a) + (b ,b))) + (declare ,oracle (struct-variable extent-oracle-t + (oracle uniform-extent-oracle) + (params (address-of ,oracle-params))))))) + +(sc-define-syntax* (experiment integrand true-integral filename) + (let ((rtols (list 0.2 0.1 0.05))) + `(begin + ,@(map (lambda (i rtol) + `(begin + (with-data-file fp ,(format #f "integral-~a-~a.dat" filename rtol) "w" + (fprintf fp "dimension\tsamples\n")))) + (iota (length rtols)) + rtols) + (for-i-step dimension 10 100 10 + (printf "Dimension %d: " dimension) + (with-rstats* samples ,(length rtols) + (for-i trial TRIALS + (with-rstats stats + ,@(map (lambda (i rtol) + `(begin + (integral (address-of ,integrand) + (address-of oracle) + (,true-integral dimension) + rng + dimension + ,rtol + stats) + (gsl-rstat-add (gsl-rstat-n stats) + (array-get samples ,i)))) + (iota (length rtols)) + rtols)) + (when (= (modulo trial (/ TRIALS 10)) 0) + (printf "%d " trial)) + (fflush NULL)) + (printf "\n") + ,@(map (lambda (i rtol) + `(with-data-file fp ,(format #f "integral-~a-~a.dat" filename rtol) "a" + (fprintf fp "%d\t%g\n" dimension (gsl-rstat-mean (array-get samples ,i))))) + (iota (length rtols)) + rtols)))))) + +(define (true-polynomial-integral dimension) + (double (unsigned int)) + (define-polynomial polynomial) + (polynomial-map! polynomial k coeff + (/ coeff (+ dimension k) (+ dimension k 1))) + (return (* (surface-area-of-ball dimension) + (gsl-pow-uint MAX-EXTENT dimension) + (polyval polynomial MAX-EXTENT)))) + +(define (lower-incomplete-gamma s x) + (double double double) + (return (* (gsl-sf-gamma s) + (gsl-sf-gamma-inc-P s x)))) + +(define (true-gaussian-integral dimension) + (double (unsigned int)) + (return (/ (* (pow 2.0 (- (* dimension 0.5) 1)) + (surface-area-of-ball dimension) + (- (* MAX-EXTENT + (lower-incomplete-gamma (* 0.5 dimension) + (* 0.5 (gsl-pow-2 MAX-EXTENT)))) + (* (sqrt 2.0) + (lower-incomplete-gamma (* 0.5 (+ dimension 1)) + (* 0.5 (gsl-pow-2 MAX-EXTENT)))))) + MAX-EXTENT))) + +(define (true-x-coordinate-integral dimension) + (double (unsigned int)) + (return (/ (* (gsl-pow-uint MAX-EXTENT (1+ dimension)) + (surface-area-of-ball (+ dimension 3))) + 2.0 + (pow M-PI 2) + (+ dimension 2)))) + +(define (main) (int) + (define-polynomial-integrand -polynomial-integrand) + (declare -gaussian-integrand (struct-variable (const integrand-t) + gaussian-integrand + NULL)) + (declare -x-coordinate-integrand (struct-variable (const integrand-t) + x-coordinate-integrand + NULL)) + (define-uniform-extent-oracle oracle 0 MAX-EXTENT) + (with-rng rng + (experiment -polynomial-integrand true-polynomial-integral "polynomial") + (experiment -gaussian-integrand true-gaussian-integral "gaussian") + (experiment -x-coordinate-integrand true-x-coordinate-integral "x-coordinate")) + (with-data-file fp "integral-parameters.tex" "w" + (fprintf fp "\\\\newcommand{\\\\integraltrials}{%d}\n" TRIALS)) + (return 0)) |