diff --git a/gsll.asd b/gsll.asd index 409bcd63c19dc5881d36397cadc6c69131453a1d..d17772508bc6fbbae7390f0882f0c9132763e60a 100644 --- a/gsll.asd +++ b/gsll.asd @@ -3,7 +3,7 @@ ; description: Definition of GSLL system ; date: ; author: Liam Healy -; modified: Fri Jul 14 2006 - 21:51 +; modified: Sat Jul 15 2006 - 22:02 ;******************************************************** ;;; $Id: $ @@ -91,4 +91,5 @@ (:module random :depends-on (init) :components - ((:file "rng-types"))))) + ((:file "rng-types") + (:file "generators" :depends-on (rng-types)))))) diff --git a/random/generators.lisp b/random/generators.lisp new file mode 100644 index 0000000000000000000000000000000000000000..cb20d584488e4ec49b554cf6bb002fb75f1acca4 --- /dev/null +++ b/random/generators.lisp @@ -0,0 +1,240 @@ +;******************************************************** +; file: generators.lisp +; description: Generators of random numbers. +; date: Sat Jul 15 2006 - 14:43 +; author: Liam M. Healy +; modified: Sat Jul 15 2006 - 21:54 +;******************************************************** +;;; $Id: $ + +(in-package :gsl) + + +(export '(random-number-generator make-random-number-generator)) + +(defclass random-number-generator () + ((type :initarg :type :reader rng-type) + (generator :initarg :generator :accessor generator)) + (:documentation "A generator of random numbers.")) + +(defmethod print-object ((object random-number-generator) stream) + (print-unreadable-object (object stream :type t :identity t) + (format stream "~a" (when (generator object) (rng-name object))))) + +(defun make-random-number-generator + (&optional (type *default-type*) (generator t)) + "Make a random number generator; by default it is allocated on creation." + (let ((instance + (make-instance + 'random-number-generator + :type type + :generator generator))) + (if (eq generator t) (alloc instance) instance))) + +;;;;**************************************************************************** +;;;; Initialization +;;;;**************************************************************************** + +(defun-gsl alloc ((generator random-number-generator)) + "gsl_rng_alloc" (((rng-type generator) :pointer)) + :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 random-number-generator)) + "gsl_rng_free" (((generator generator) :pointer)) + :type :method + :c-return :void + :after ((setf (generator generator) nil)) + :documentation "Free all the memory associated with the generator.") + +;;;;**************************************************************************** +;;;; Seed +;;;;**************************************************************************** + +(cffi:defcvar ("gsl_rng_default_seed" *default-seed*) :ulong) +(setf (documentation '*default-seed* 'variable) + "The default seed for random number generators.") +(map-name '*default-seed* "gsl_rng_cmrg") +(export '*default-seed*) + +(defun-gsl rng-set (rng-instance value) + "gsl_rng_set" (((generator rng-instance) :pointer) (value :ulong)) + :c-return :void + :documentation + "Initialize (or `seeds') the random number generator. If + the generator is seeded with the same value of @var{s} on two different + runs, the same stream of random numbers will be generated by successive + calls to the routines below. If different values of @var{s} are + supplied, then the generated streams of random numbers should be + completely different. If the seed @var{s} is zero then the standard seed + from the original implementation is used instead. For example, the + original Fortran source code for the @code{ranlux} generator used a seed + of 314159265, and so choosing @var{s} equal to zero reproduces this when + using @code{gsl_rng_ranlux}.") + +;;;;**************************************************************************** +;;;; Sampling +;;;;**************************************************************************** + +(defun-gsl get-random-number (generator) + "gsl_rng_get" (((generator generator) :pointer)) + :c-return :ulong + :documentation + "Generate a random integer. The + minimum and maximum values depend on the algorithm used, but all + integers in the range [@var{min},@var{max}] are equally likely. The + values of @var{min} and @var{max} can determined using the auxiliary + functions @code{gsl_rng_max (r)} and @code{gsl_rng_min (r)}.") + +(defun-gsl uniform (generator) + "gsl_rng_uniform" (((generator generator) :pointer)) + :c-return :double + :documentation + "A double precision floating point number uniformly + distributed in the range [0,1). The range includes 0.0 but excludes 1.0. + The value is typically obtained by dividing the result of + #'rng-get by (+ (rng-max generator) 1.0) in double + precision. Some generators compute this ratio internally so that they + can provide floating point numbers with more than 32 bits of randomness + (the maximum number of bits that can be portably represented in a single + :ulong.") + +(defun-gsl uniform>0 (generator) + "gsl_rng_uniform_pos" (((generator generator) :pointer)) + :c-return :double + :documentation + "Return a positive double precision floating point number + uniformly distributed in the range (0,1), excluding both 0.0 and 1.0. + The number is obtained by sampling the generator with the algorithm of + @code{gsl_rng_uniform} until a non-zero value is obtained. You can use + this function if you need to avoid a singularity at 0.0.") + +(defun-gsl uniform-fixnum (generator upperbound) + "gsl_rng_uniform_int" (((generator generator) :pointer) (upperbound :ulong)) + :c-return :ulong + :documentation + "Generate a random integer from 0 to upperbound-1 inclusive. + All integers in the range are equally likely, regardless + of the generator used. An offset correction is applied so that zero is + always returned with the correct probability, for any minimum value of + the underlying generator. If upperbound is larger than the range + of the generator then the function signals an error.") + +;;;;**************************************************************************** +;;;; Information functions about instances +;;;;**************************************************************************** + +(defun-gsl rng-name (rng-instance) + "gsl_rng_name" (((generator rng-instance) :pointer)) + :c-return :string + :documentation "The name of the random number generator.") + +(defun-gsl rng-max (rng-instance) + "gsl_rng_max" (((generator rng-instance) :pointer)) + :c-return :unsigned-long + :documentation "The largest value that @code{gsl_rng_get} + can return.") + +(defun-gsl rng-min (rng-instance) + "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)) + :c-return :pointer + :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)) + :c-return :size + :export nil + :index gsl-random-state + :documentation "The size of the generator.") + +(export 'gsl-random-state) +(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)) + (ans + (make-array (rng-size gen) + :element-type '(unsigned-byte 8)))) + (loop for i from 0 below (length ans) + do + (setf (aref ans i) + (mem-aref (rng-state gen) :uint8 i))) + ans)) + +;;;;**************************************************************************** +;;;; Copying state +;;;;**************************************************************************** + +(defun-gsl copy + ((destination random-number-generator) (source random-number-generator)) + "gsl_rng_memcpy" + (((generator destination) :pointer) ((generator source) :pointer)) + :type :method + :documentation + "Copy the random number generator @var{source} into the + pre-existing generator @var{destination}, + 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)) + :c-return :pointer + :documentation + "Create a new generator which is an exact copy of the original. + Dont' use; use #'make-random-number-generator, #'copy instead.") + +(defun-gsl write-binary + ((object random-number-generator) stream) + "gsl_rng_fwrite" + ((stream :pointer) ((generator object) :pointer)) + :type :method) + +(defun-gsl read-binary + ((object random-number-generator) stream) + "gsl_block_fread" + ((stream :pointer) ((pointer object) :pointer)) + :type :method) + +;;;;**************************************************************************** +;;;; Examples and unit test +;;;;**************************************************************************** + +(defparameter *rng-mt19937* (make-random-number-generator *mt19937*)) +(defparameter *rng-cmrg* (make-random-number-generator *cmrg*)) + +(lisp-unit:define-test random-number-generators + (lisp-unit:assert-equal + '(999 162 282 947 231 484 957 744 540 739 759) + (progn + (rng-set *rng-mt19937* 0) + (loop for i from 0 to 10 + collect + (uniform-fixnum *rng-mt19937* 1000)))) + (lisp-unit:assert-equal + '("0.111776229978d+00" "0.959166794996d+00" "0.841526801158d+00" + "0.925403713680d+00" "0.275406984741d+00" "0.709304057392d+00" + "0.554133304187d+00" "0.880695776958d+00" "0.597139396983d+00" + "0.751874113340d+00" "0.931108462127d+00") + (lisp-unit:fp-sequence + (progn + (rng-set *rng-cmrg* 0) + (loop for i from 0 to 10 + collect (uniform *rng-cmrg*)))))) diff --git a/random/rng-types.lisp b/random/rng-types.lisp index ebab611b348bd51ca9f9d66a0410746fb13e6b3e..6c94a5810e9690113e80f4137acf9c5b0c1f3041 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: Fri Jul 14 2006 - 22:34 +; modified: Sat Jul 15 2006 - 16:29 ;******************************************************** ;;; $Id: $ @@ -11,11 +11,6 @@ (in-package :gsl) -(defun-gsl setup-random-number-generator () - "gsl_rng_env_setup" () - :c-return :pointer - :documentation "Setup random number generator.") - ;;;;**************************************************************************** ;;;; Random number generator types ;;;;**************************************************************************** @@ -34,46 +29,7 @@ ;;;; Auxiliary functions ;;;;**************************************************************************** -(export '(random-number-generator-state all-random-number-generators)) - -(defun-gsl rng-name (rng) - "gsl_rng_name" ((rng :pointer)) - :c-return :string - :documentation "The name of the random number generator.") - -(defun-gsl rng-max (rng) - "gsl_rng_max" ((rng :pointer)) - :c-return :unsigned-long - :documentation "The largest value that @code{gsl_rng_get} - can return.") - -(defun-gsl rng-min (rng) - "gsl_rng_min" ((rng :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) - "gsl_rng_state" ((rng :pointer)) - :c-return :pointer - :documentation "A pointer to the state of generator.") - -(defun-gsl rng-size (rng) - "gsl_rng_size" ((rng :pointer)) - :c-return :size - :documentation "The size of the generator.") - -(defun random-number-generator-state (rng) - "The complete state of a given random number generator, specified - as a vector of bytes." - (let ((ans (make-array (rng-size rng) :element-type '(unsigned-byte 8)))) - (loop for i from 0 below (length ans) - do - (setf (aref ans i) - (mem-aref (rng-state rng) :uint8 i))) - ans)) +(export '(all-random-number-generators)) (defun-gsl rng-types-setup () "gsl_rng_types_setup" () @@ -86,13 +42,15 @@ Users should call all-random-number-generators.") (defun all-random-number-generators () - "A list of all random number generators." + "A list of all random number generator types." (let ((start (rng-types-setup))) (loop for i from 0 - for ptr - = (cffi:inc-pointer start (* (cffi:foreign-type-size :pointer) i)) - until (cffi:null-pointer-p (cffi:mem-ref ptr :pointer)) - collect (cffi:mem-ref ptr 'random-number-generator-type)))) + for ptr + = (cffi:mem-ref + (cffi:inc-pointer start (* (cffi:foreign-type-size :pointer) i)) + :pointer) + until (cffi:null-pointer-p ptr) + collect ptr))) ;;;;**************************************************************************** ;;;; Defining RNGs and default @@ -104,18 +62,25 @@ (or gsl-name (remove #\* (substitute #\_ #\- - (format nil "gsl_rng_~(~a~)" lisp-name))))) - (sym (gensym "RNG"))) + (format nil "gsl_rng_~(~a~)" lisp-name)))))) `(progn - (cffi:defcvar (,cname ,sym) :pointer :read-only t) - (defparameter ,lisp-name (cffi:get-var-pointer ',sym) - ,documentation) - (map-name ',lisp-name ,cname)))) - -(def-rng-type *rngdefault* - "The default, set by environment variables GSL_RNG_TYPE and GSL_RNG_SEED" + (cffi:defcvar (,cname ,lisp-name) :pointer :read-only t) + (setf (documentation ',lisp-name 'variable) ,documentation) + (map-name ',lisp-name ,cname) + (export ',lisp-name)))) + +(def-rng-type *default-type* + "The default random number generator type, + set by environment variables GSL_RNG_TYPE and GSL_RNG_SEED" "gsl_rng_default") +(defun-gsl rng-environment-setup () + "gsl_rng_env_setup" () + :c-return :pointer + :documentation "Read the environment variables @code{GSL_RNG_TYPE} and + @code{GSL_RNG_SEED} and uses their values to set the corresponding + library variables *default-type* and *default-seed*") + ;;;;**************************************************************************** ;;;; Modern random number generators ;;;;****************************************************************************