aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorArun Isaac2021-05-07 17:16:30 +0530
committerArun Isaac2021-05-07 17:16:30 +0530
commit3df13fea1c3b585910437dc2e787bfabebfed163 (patch)
tree35c53315bb1ba9b60f4c11365148445825d663a2 /src
parentf75010b875d19c8cf5e88f64a54efa033f3466b1 (diff)
downloadnsmc-3df13fea1c3b585910437dc2e787bfabebfed163.tar.gz
nsmc-3df13fea1c3b585910437dc2e787bfabebfed163.tar.lz
nsmc-3df13fea1c3b585910437dc2e787bfabebfed163.zip
Implement integral experiments.
* src/integral.sc: New file. * CMakeLists.txt: Explicitly specify sources for nsmc library. Add integral executable.
Diffstat (limited to 'src')
-rw-r--r--src/integral.sc140
1 files changed, 140 insertions, 0 deletions
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))