From d0f5bc85b53e7450ba8cb90e964779688460a22b Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Mon, 17 Jul 2006 02:45:58 +0000
Subject: [PATCH] Quasi-random number generators and Gaussian distributions.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3124 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
---
 gsll.asd               |   7 +-
 random/gaussian.lisp   | 165 +++++++++++++++++++++++++++++++++++++++++
 random/generators.lisp |  34 +++++----
 random/quasi.lisp      | 138 ++++++++++++++++++++++++++++++++++
 random/rng-types.lisp  |   7 +-
 5 files changed, 333 insertions(+), 18 deletions(-)
 create mode 100644 random/gaussian.lisp
 create mode 100644 random/quasi.lisp

diff --git a/gsll.asd b/gsll.asd
index d1777250..f496756b 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Sat Jul 15 2006 - 22:02
+; modified:    Sun Jul 16 2006 - 22:44
 ;********************************************************
 ;;; $Id: $
 
@@ -92,4 +92,7 @@
 	    :depends-on (init)
 	    :components
 	    ((:file "rng-types")
-	     (:file "generators" :depends-on (rng-types))))))
+	     (:file "generators" :depends-on (rng-types))
+	     (:file "quasi" :depends-on (rng-types))
+	     (:file "gaussian" :depends-on (rng-types))))))
+
diff --git a/random/gaussian.lisp b/random/gaussian.lisp
new file mode 100644
index 00000000..c43bcbf2
--- /dev/null
+++ b/random/gaussian.lisp
@@ -0,0 +1,165 @@
+;********************************************************
+; file:        gaussian.lisp                             
+; description: Gaussian distribution                     
+; date:        Sun Jul 16 2006 - 22:09                   
+; author:      Liam M. Healy                             
+; modified:    Sun Jul 16 2006 - 22:42
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl gaussian (generator sigma)
+  "gsl_ran_gaussian"
+  (((generator generator) :pointer) (sigma :double))
+  :c-return :double
+  :documentation
+  "A Gaussian random variate, with mean zero and
+  standard deviation @var{sigma}.  The probability distribution for
+  Gaussian random variates is
+  p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx
+  for @math{x} in the range @math{-\infty} to @math{+\infty}.  Use the
+  transformation @math{z = \mu + x} on the numbers returned by
+  @code{gsl_ran_gaussian} to obtain a Gaussian distribution with mean
+  @math{\mu}.  This function uses the Box-Mueller algorithm which requires two
+  calls to the random number generator @var{r}.")
+
+(defun-gsl gaussian-pdf (x sigma)
+  "gsl_ran_gaussian_pdf" ((x :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "Compute the probability density @math{p(x)} at @var{x}
+   for a Gaussian distribution with standard deviation @var{sigma}.")
+
+(defun-gsl gaussian-ziggurat (generator sigma)
+  "gsl_ran_gaussian_ziggurat"
+  (((generator generator) :pointer) (sigma :double))
+  :c-return :double
+  :documentation
+  "Compute a Gaussian random variate using the alternative
+   Marsaglia-Tsang ziggurat method. The Ziggurat algorithm
+   is the fastest available algorithm in most cases.")
+
+(defun-gsl gaussian-ratio-method (generator sigma)
+  "gsl_ran_gaussian_ratio_method"
+  (((generator generator) :pointer) (sigma :double))
+  :c-return :double
+  :documentation
+  "Compute a Gaussian random variate using the Kinderman-Monahan-Leva
+   ratio method.")
+
+(defun-gsl ugaussian (generator)
+  "gsl_ran_ugaussian" (((generator generator) :pointer))
+  :c-return :double
+  :documentation
+  "Compute results for the unit Gaussian distribution,
+   equivalent to the #'gaussian with a standard deviation of one,
+   @var{sigma} = 1.")
+
+(defun-gsl ugaussian-pdf (x)
+  "gsl_ran_ugaussian_pdf" ((x :double))
+  :c-return :double
+  :documentation
+  "Compute results for the unit Gaussian distribution,
+   equivalent to the #'gaussian-pdf with a standard deviation of one,
+   @var{sigma} = 1.")
+
+(defun-gsl ugaussian-ratio-method (generator)
+  "gsl_ran_ugaussian_ratio_method"
+  (((generator generator) :pointer))
+  :c-return :double
+  :documentation
+  "Compute results for the unit Gaussian distribution,
+   equivalent to the #'gaussian-ration-method with a
+   standard deviation of one, @var{sigma} = 1.")
+
+(defun-gsl gaussian-P (x sigma)
+  "gsl_cdf_gaussian_P" ((x :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "The cumulative distribution function @math{P(x)} for the Gaussian
+   distribution with standard deviation @var{sigma}.")
+
+(defun-gsl gaussian-Q (x sigma)
+  "gsl_cdf_gaussian_Q" ((x :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "The cumulative distribution function @math{Q(x)} for the Gaussian
+   distribution with standard deviation @var{sigma}.")
+
+(defun-gsl gaussian-Pinv (P sigma)
+  "gsl_cdf_gaussian_Pinv" ((P :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "The inverse cumulative distribution function @math{P(x)} for the Gaussian
+   distribution with standard deviation @var{sigma}.")
+
+(defun-gsl gaussian-Qinv (Q sigma)
+  "gsl_cdf_gaussian_Qinv" ((Q :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "The inverse cumulative distribution function @math{Q(x)} for the Gaussian
+   distribution with standard deviation @var{sigma}.")
+
+(defun-gsl ugaussian-P (x)
+  "gsl_cdf_ugaussian_P" ((x :double))
+  :c-return :double
+  :documentation
+  "The cumulative distribution function @math{P(x)} for the Gaussian
+   distribution with unit standard deviation.")
+
+(defun-gsl ugaussian-Q (x)
+  "gsl_cdf_ugaussian_Q" ((x :double))
+  :c-return :double
+  :documentation
+  "The cumulative distribution function @math{Q(x)} for the Gaussian
+   distribution with unit standard deviation.")
+
+(defun-gsl ugaussian-Pinv (P)
+  "gsl_cdf_ugaussian_Pinv" ((P :double))
+  :c-return :double
+  :documentation
+  "The inverse cumulative distribution function @math{P(x)} for the Gaussian
+   distribution with unit standard deviation.")
+
+(defun-gsl ugaussian-Qinv (Q)
+  "gsl_cdf_ugaussian_Qinv" ((Q :double))
+  :c-return :double
+  :documentation
+  "The inverse cumulative distribution function @math{Q(x)} for the Gaussian
+   distribution with unit standard deviation.")
+
+;;; Examples and unit test
+(lisp-unit:define-test gaussian
+  (lisp-unit:assert-equal
+   '("0.133918608119d+01" "-0.881009918314d+00" "0.167440840625d+02"
+     "0.733641107293d+01" "0.997524631602d+01" "-0.127750208100d+02"
+     "-0.239671528273d+02" "-0.679280164729d+01" "-0.390913184336d+00"
+     "0.893555545521d+01" "-0.176477945898d+00")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	 collect
+	 (gaussian *rng-mt19937* 10.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.398942280401d-01"
+   (gaussian-pdf 0.0d0 10.0d0))
+  (lisp-unit:assert-equal
+   '("0.610792663776d+01" "-0.157649178949d+02" "-0.720097139253d+01"
+     "0.195115024682d+02" "0.524820292158d+01" "0.104085102623d+02"
+     "-0.143332423262d+02" "0.945508871923d+00" "-0.111203439859d+02"
+     "-0.307919574280d+01" "-0.209068077307d+01")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	 collect
+	 (gaussian-ziggurat *rng-mt19937* 10.0d0)))))
+  ;; Given in examples in GSL documentation
+  (lisp-unit:assert-first-fp-equal "0.977249868052d+00" (ugaussian-p 2.0d0))
+  (lisp-unit:assert-first-fp-equal "0.227501319482d-01" (ugaussian-q 2.0d0))
+  (lisp-unit:assert-first-fp-equal "0.200000000000d+01"
+				   (ugaussian-pinv 0.9772498680518208d0))
+  (lisp-unit:assert-first-fp-equal "0.200000000000d+01"
+				   (ugaussian-qinv 0.02275013194817921d0)))
diff --git a/random/generators.lisp b/random/generators.lisp
index cb20d584..d48f7d0a 100644
--- a/random/generators.lisp
+++ b/random/generators.lisp
@@ -3,21 +3,25 @@
 ; description: Generators of random numbers.             
 ; date:        Sat Jul 15 2006 - 14:43                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Jul 15 2006 - 21:54
+; modified:    Sun Jul 16 2006 - 17:41
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
+;;; Would it be useful to have a make-data like macro?
 
 (export '(random-number-generator make-random-number-generator))
 
-(defclass random-number-generator ()
+(defclass gsl-random ()
   ((type :initarg :type :reader rng-type)
-   (generator :initarg :generator :accessor generator))
+   (generator :initarg :generator :accessor generator)))
+
+(defclass random-number-generator (gsl-random)
+  ()
   (:documentation "A generator of random numbers."))
 
-(defmethod print-object ((object random-number-generator) stream)
+(defmethod print-object ((object gsl-random) stream)
   (print-unreadable-object (object stream :type t :identity t) 
     (format stream "~a" (when (generator object) (rng-name object)))))
 
@@ -132,8 +136,9 @@
 ;;;; Information functions about instances
 ;;;;****************************************************************************
 
-(defun-gsl rng-name (rng-instance)
+(defun-gsl rng-name ((rng-instance random-number-generator))
   "gsl_rng_name" (((generator rng-instance) :pointer))
+  :type :method
   :c-return :string
   :documentation "The name of the random number generator.")
 
@@ -144,23 +149,25 @@
    can return.")
 
 (defun-gsl rng-min (rng-instance)
-  "gsl_rng_min" ((generator rng-instance :pointer))
+  "gsl_rng_min" (((generator rng-instance) :pointer))
   :c-return :unsigned-long
   :documentation "The smallest value that @code{gsl_rng_get}
    can return.  Usually this value is zero.  There are some generators with
    algorithms that cannot return zero, and for these generators the minimum
    value is 1.")
 
-(defun-gsl rng-state (rng-instance)
-  "gsl_rng_state" ((rng-instance :pointer))
+(defun-gsl rng-state ((rng-instance random-number-generator))
+  "gsl_rng_state" (((generator rng-instance) :pointer))
   :c-return :pointer
+  :type :method
   :export nil
   :index gsl-random-state
   :documentation "A pointer to the state of generator.")
 
-(defun-gsl rng-size (rng-instance)
-  "gsl_rng_size" ((rng-instance :pointer))
+(defun-gsl rng-size ((rng-instance random-number-generator))
+  "gsl_rng_size" (((generator rng-instance) :pointer))
   :c-return :size
+  :type :method
   :export nil
   :index gsl-random-state
   :documentation "The size of the generator.")
@@ -169,7 +176,7 @@
 (defun gsl-random-state (rng-instance)
   "The complete state of a given random number generator, specified
    as a vector of bytes."
-  (let* ((gen (generator rng-instance))
+  (let* ((gen rng-instance)
 	 (ans
 	  (make-array (rng-size gen)
 		      :element-type '(unsigned-byte 8))))
@@ -194,9 +201,10 @@
    making @var{destination} into an exact copy
    of @var{source}.  The two generators must be of the same type.")
 
-(defun-gsl clone-generator (rng-instance)
-  "gsl_rng_clone" ((rng-instance :pointer))
+(defun-gsl clone-generator ((instance random-number-generator))
+  "gsl_rng_clone" (((generator instance) :pointer))
   :c-return :pointer
+  :type :method
   :documentation
   "Create a new generator which is an exact copy of the original.
    Dont' use; use #'make-random-number-generator, #'copy instead.")
diff --git a/random/quasi.lisp b/random/quasi.lisp
new file mode 100644
index 00000000..906251b9
--- /dev/null
+++ b/random/quasi.lisp
@@ -0,0 +1,138 @@
+;********************************************************
+; file:        quasi.lisp                                
+; description: Quasi-Random Number Generators            
+; date:        Sun Jul 16 2006 - 15:54                   
+; author:      Liam M. Healy                             
+; modified:    Sun Jul 16 2006 - 17:40
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defclass quasi-random-number-generator (gsl-random)
+  ((dimension :initarg :dimension :accessor qr-dimension))
+  (:documentation "A generator of quasi-random numbers."))
+
+(defun make-quasi-random-number-generator 
+    (dimension &optional (type *niederreiter2*) (generator t))
+  "Make a random number generator; by default it is allocated on creation."
+  (let ((instance
+	 (make-instance
+	  'quasi-random-number-generator
+	  :type type
+	  :dimension dimension
+	  :generator generator)))
+    (if (eq generator t) (alloc instance) instance)))
+
+(defun-gsl alloc ((generator quasi-random-number-generator))
+  "gsl_qrng_alloc"
+  (((rng-type generator) :pointer) ((qr-dimension generator) :uint))
+  :type :method
+  :c-return (ptr :pointer)
+  :return ((progn (setf (generator generator) ptr) generator))
+  :documentation
+  "Instatiate a random number generator of specified type.
+   For example, create an instance of the Tausworthe
+   generator: (rng-alloc *taus*).
+   The generator is automatically initialized with the default seed,
+   @code{gsl_rng_default_seed}.  This is zero by default but can be changed
+   either directly or by using the environment variable @code{GSL_RNG_SEED}")
+
+(defun-gsl free ((generator quasi-random-number-generator))
+  "gsl_qrng_free" (((generator generator) :pointer))
+  :type :method
+  :c-return :void
+  :after ((setf (generator generator) nil))
+  :documentation "Free all the memory associated with the generator.")
+
+(defun-gsl init (generator)
+  "gsl_qrng_init" (((generator generator) :pointer))
+  :c-return :void
+  :documentation
+  "Reinitializes the generator @var{q} to its starting point.
+   Note that quasi-random sequences do not use a seed and always produce
+   the same set of values.")
+
+(defun-gsl qrng-get (generator return-vector)
+  "gsl_qrng_get"
+  (((generator generator) :pointer) ((gsl-array return-vector) :pointer))
+  :documentation
+  "Store the next point from the sequence generator @var{q}
+   in the array @var{x}.  The space available for @var{x} must match the
+   dimension of the generator.  The point @var{x} will lie in the range
+   @math{0 < x_i < 1} for each @math{x_i}."
+  :invalidate (return-vector))
+
+(defun-gsl rng-name ((instance quasi-random-number-generator))
+  "gsl_qrng_name" (((generator instance) :pointer))
+  :type :method
+  :c-return :string
+  :documentation "The name of the random number generator.")
+
+(defun-gsl rng-state ((instance quasi-random-number-generator))
+  "gsl_qrng_state" (((generator instance) :pointer))
+  :c-return :pointer
+  :type :method
+  :export nil
+  :index gsl-random-state
+  :documentation "A pointer to the state of generator.")
+
+(defun-gsl rng-size ((instance quasi-random-number-generator))
+  "gsl_qrng_size" (((generator instance) :pointer))
+  :c-return :size
+  :type :method
+  :export nil
+  :index gsl-random-state
+  :documentation "The size of the generator.")
+
+(defun-gsl copy
+    ((destination quasi-random-number-generator)
+     (source quasi-random-number-generator))
+  "gsl_qrng_memcpy"
+  (((generator destination) :pointer) ((generator source) :pointer))
+  :type :method
+  :documentation
+  "Copy the quasi-random sequence generator @var{src} into the
+   pre-existing generator @var{dest}, making @var{dest} into an exact copy
+   of @var{src}.  The two generators must be of the same type.")
+
+(defun-gsl clone-generator ((instance quasi-random-number-generator))
+  "gsl_qrng_clone" (((generator instance) :pointer))
+  :c-return :pointer
+  :type :method
+  :documentation
+  "Create a new generator which is an exact copy of the original.
+   Dont' use; use #'make-random-number-generator, #'copy instead.")
+
+(def-rng-type *niederreiter2*
+    "Described in Bratley, Fox, Niederreiter,
+     @cite{ACM Trans. Model. Comp. Sim.} 2, 195 (1992). It is
+     valid up to 12 dimensions."
+  "gsl_qrng_niederreiter_2")
+
+(def-rng-type *sobol*
+    "This generator uses the Sobol sequence described in Antonov, Saleev,
+    @cite{USSR Comput. Maths. Math. Phys.} 19, 252 (1980). It is valid up to
+    40 dimensions."
+  "gsl_qrng_sobol")
+
+(with-data (vec vector-double 2)
+  (let ((gen (make-quasi-random-number-generator 2 *sobol*)))
+    (init gen)
+    (prog1
+	(loop repeat 5 times
+	   append (coerce (data (qrng-get gen vec)) 'list))
+      (free gen))))
+
+;;; Examples and unit test
+(lisp-unit:define-test quasi-random-number-generators
+  ;; This example is given in the GSL documentation
+  (lisp-unit:assert-equal
+   '(0.5d0 0.5d0 0.75d0 0.25d0 0.25d0 0.75d0 0.375d0 0.375d0 0.875d0 0.875d0)
+   (with-data (vec vector-double 2)
+     (let ((gen (make-quasi-random-number-generator 2 *sobol*)))
+       (init gen)
+       (prog1
+	   (loop repeat 5 
+	      append (coerce (data (qrng-get gen vec)) 'list))
+	 (free gen))))))
diff --git a/random/rng-types.lisp b/random/rng-types.lisp
index 6c94a581..65cb016d 100644
--- a/random/rng-types.lisp
+++ b/random/rng-types.lisp
@@ -3,7 +3,7 @@
 ; description: Random number generation                  
 ; date:        Tue Jul 11 2006 - 23:39                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Jul 15 2006 - 16:29
+; modified:    Sun Jul 16 2006 - 16:28
 ;********************************************************
 ;;; $Id: $
 
@@ -61,8 +61,9 @@
   (let ((cname
 	 (or gsl-name
 	     (remove #\*
-		     (substitute #\_ #\- 
-				 (format nil "gsl_rng_~(~a~)" lisp-name))))))
+		     (substitute
+		      #\_ #\- 
+		      (format nil "gsl_rng_~(~a~)" lisp-name))))))
     `(progn
        (cffi:defcvar (,cname ,lisp-name) :pointer :read-only t)
        (setf (documentation ',lisp-name 'variable) ,documentation)
-- 
GitLab