monte-carlo.lisp 9.29 KB
Newer Older
1 2
;; Monte Carlo Integration
;; Liam Healy Sat Feb  3 2007 - 17:42
3
;; Time-stamp: <2011-08-16 23:54:40EDT monte-carlo.lisp>
4
;;
5
;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
liam's avatar
liam committed
20 21 22

(in-package :gsl)

23 24 25 26 27
;;; /usr/include/gsl/gsl_monte.h
;;; /usr/include/gsl/gsl_monte_plain.h
;;; /usr/include/gsl/gsl_monte_miser.h
;;; /usr/include/gsl/gsl_monte_vegas.h

liam's avatar
liam committed
28 29 30 31
;;;;****************************************************************************
;;;; PLAIN Monte Carlo
;;;;****************************************************************************

32 33 34
(defmobject monte-carlo-plain
    "gsl_monte_plain"
  ((dim sizet))
35 36
  "plain Monte Carlo integration"
  :documentation			; FDL
37
  "Make and initialize a workspace for Monte Carlo integration in dimension dim."
38 39
  :initialize-suffix "init"
  :initialize-args nil)
liam's avatar
liam committed
40

41 42
(defparameter *monte-carlo-default-samples-per-dimension* 150000)

43
(defmfun monte-carlo-integrate-plain
44 45 46 47 48 49 50 51
    (function lower-limits upper-limits 
	      &optional
	      (number-of-samples
	       (* *monte-carlo-default-samples-per-dimension*
		  (dim0 lower-limits)))
	      (generator (make-random-number-generator +mt19937+ 0))
	      (state (make-monte-carlo-plain (dim0 lower-limits)))
	      (scalars t))
52
  "gsl_monte_plain_integrate"
53
  ((callback :pointer)
54
   ((grid:foreign-pointer lower-limits) :pointer) ((grid:foreign-pointer upper-limits) :pointer)
55
   ((dim0 lower-limits) sizet) (number-of-samples sizet)
56 57
   ((mpointer generator) :pointer)
   ((mpointer state) :pointer)
58
   (result (:pointer :double)) (abserr (:pointer :double)))
59
  :inputs (lower-limits upper-limits)
60
  :callbacks
61
  (callback fnstruct-dimension (dimension)
62 63
	    (function :double (:input :double :cvector dim0) :slug))
  :callback-dynamic (((dim0 lower-limits)) (function scalars))
64 65 66 67 68 69
  :documentation			; FDL
  "Uses the plain Monte Carlo algorithm to integrate the
   function f over the hypercubic region defined by the
   lower and upper limits in the arrays 'lower-limits and
   'upper-limits, each a gsl-vector of length dim.
   The integration uses a fixed number
70
   of function calls number-of-samples, and obtains random sampling points using
71 72 73 74
   the random number generator 'generator. A previously allocated workspace
   'state must be supplied.  The result of the integration is returned
   with an estimated absolute error.")

liam's avatar
liam committed
75 76 77 78 79 80 81 82 83
;;;;****************************************************************************
;;;; MISER
;;;;****************************************************************************

;;; The MISER algorithm of Press and Farrar is based on recursive
;;; stratified sampling.  This technique aims to reduce the overall
;;; integration error by concentrating integration points in the
;;; regions of highest variance.

84 85 86
(defmobject monte-carlo-miser
    "gsl_monte_miser"
  ((dim sizet))
87 88
  "miser Monte Carlo integration"
  :documentation			; FDL
89
  "Make and initialize a workspace for Monte Carlo integration in
90
   dimension dim.  The workspace is used to maintain
91
   the state of the integration."
92 93
  :initialize-suffix "init"
  :initialize-args nil)
94

liam's avatar
liam committed
95 96
(export 'miser-parameter)
(defmacro miser-parameter (workspace parameter)
liam's avatar
liam committed
97
  ;; FDL
liam's avatar
liam committed
98 99 100 101
  "Get or set with setf the parameter value for the MISER Monte Carlo
   integration method."
  ;; (miser-parameter ws min-calls)
  ;; (setf (miser-parameter ws min-calls) 300)
102
 `(cffi:foreign-slot-value ,workspace 'miser-state ',parameter))
liam's avatar
liam committed
103

liam's avatar
liam committed
104
(defmfun monte-carlo-integrate-miser
105 106 107 108 109 110 111 112
    (function lower-limits upper-limits 
	      &optional
	      (number-of-samples
	       (* *monte-carlo-default-samples-per-dimension*
		  (dim0 lower-limits)))
	      (generator (make-random-number-generator +mt19937+ 0))
	      (state (make-monte-carlo-miser (dim0 lower-limits)))
	      (scalars t))
liam's avatar
liam committed
113
  "gsl_monte_miser_integrate"
114
  ((callback :pointer)
115
   ((grid:foreign-pointer lower-limits) :pointer) ((grid:foreign-pointer upper-limits) :pointer)
116
   ((dim0 lower-limits) sizet) (number-of-samples sizet)
117
   ((mpointer generator) :pointer)
118
   ((mpointer state) :pointer)
119
   (result (:pointer :double)) (abserr (:pointer :double)))
120
  :inputs (lower-limits upper-limits)
121
  :callbacks
122
  (callback fnstruct-dimension (dimension)
123 124
	    (function :double (:input :double :cvector dim0) :slug))
  :callback-dynamic (((dim0 lower-limits)) (function scalars))
liam's avatar
liam committed
125
  :documentation			; FDL
liam's avatar
liam committed
126
  "Uses the miser Monte Carlo algorithm to integrate the
liam's avatar
liam committed
127
   function f over the hypercubic region defined by the
liam's avatar
liam committed
128
   lower and upper limits in the arrays 'lower-limits and
liam's avatar
liam committed
129
   'upper-limits, each a gsl-vector of the samelength
liam's avatar
liam committed
130
   The integration uses a fixed number
131
   of function calls number-of-samples, and obtains random sampling points using
liam's avatar
liam committed
132
   the random number generator 'generator. A previously allocated workspace
liam's avatar
liam committed
133 134
   'state must be supplied.  The result of the integration is returned
   with an estimated absolute error.")
liam's avatar
liam committed
135 136 137 138 139

;;;;****************************************************************************
;;;; VEGAS
;;;;****************************************************************************

liam's avatar
liam committed
140
;;; The vegas algorithm of Lepage is based on importance sampling.  It
liam's avatar
liam committed
141
;;; samples points from the probability distribution described by the
liam's avatar
liam committed
142
;;; function |f|, so that the points are concentrated in the regions
liam's avatar
liam committed
143 144
;;; that make the largest contribution to the integral.

145 146 147
(defmobject monte-carlo-vegas
    "gsl_monte_vegas"
  ((dim sizet))
148 149
  "vegas Monte Carlo integration"
  :documentation			; FDL
150
  "Make and initialize a workspace for Monte Carlo integration in
151
   dimension dim.  The workspace is used to maintain
152
   the state of the integration.  Returns a pointer to vegas-state."
153 154
  :initialize-suffix "init"
  :initialize-args nil)
155

liam's avatar
liam committed
156 157
(export 'vegas-parameter)
(defmacro vegas-parameter (workspace parameter)
liam's avatar
liam committed
158
  ;; FDL
liam's avatar
liam committed
159 160 161 162
  "Get or set with setf the parameter value for the VEGAS Monte Carlo
   integration method."
  ;; (vegas-parameter ws bins-max)
  ;; (setf (vegas-parameter ws bins-max) 300)
163
 `(cffi:foreign-slot-value ,workspace 'vegas-state ',parameter))
liam's avatar
liam committed
164

liam's avatar
liam committed
165
(defmfun monte-carlo-integrate-vegas
166 167 168 169 170 171 172 173
    (function lower-limits upper-limits 
	      &optional
	      (number-of-samples
	       (* *monte-carlo-default-samples-per-dimension*
		  (dim0 lower-limits)))
	      (generator (make-random-number-generator +mt19937+ 0))
	      (state (make-monte-carlo-vegas (dim0 lower-limits)))
	      (scalars t))
liam's avatar
liam committed
174
  "gsl_monte_vegas_integrate"
175
  ((callback :pointer)
176
   ((grid:foreign-pointer lower-limits) :pointer) ((grid:foreign-pointer upper-limits) :pointer)
177
   ((dim0 lower-limits) sizet) (number-of-samples sizet)
178
   ((mpointer generator) :pointer)
179
   ((mpointer state) :pointer)
180
   (result (:pointer :double)) (abserr (:pointer :double)))
181
  :inputs (lower-limits upper-limits)
182
  :callbacks
183
  (callback fnstruct-dimension (dimension)
184 185
	    (function :double (:input :double :cvector dim0) :slug))
  :callback-dynamic (((dim0 lower-limits)) (function scalars))
liam's avatar
liam committed
186
  :documentation			; FDL
187 188 189 190 191 192 193 194 195 196 197 198
  "Uses the vegas Monte Carlo algorithm to integrate the function f
   over the dim-dimensional hypercubic region defined by the lower and
   upper limits in the arrays x1 and xu, each of the same length.  The
   integration uses a fixed number of function calls
   number-of-samples, and obtains random sampling points using the
   random number generator r.  A previously allocated workspace s must
   be supplied.  The result of the integration is returned with an
   estimated absolute error.  The result and its error estimate are
   based on a weighted average of independent samples. The chi-squared
   per degree of freedom for the weighted average is returned via the
   state struct component, s->chisq, and must be consistent with 1 for
   the weighted average to be reliable.")
liam's avatar
liam committed
199 200 201 202 203 204 205

;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

;;; Example from Sec. 23.5

206 207
(defun mcrw (x y z)
  "Example function for Monte Carlo used in random walk studies."
208
  (* (/ (expt dpi 3))
209 210
     (/ (- 1 (* (cos x) (cos y) (cos z))))))

211 212 213
(defparameter *mc-lower*
  (grid:make-foreign-array
   'double-float :initial-contents '(0.0d0 0.0d0 0.0d0)))
214 215

(defparameter *mc-upper*
216 217
  (grid:make-foreign-array
   'double-float :initial-contents (make-list 3 :initial-element dpi)))
218

liam's avatar
liam committed
219
(defun random-walk-plain-example (&optional (nsamples 500000))
220
  (monte-carlo-integrate-plain 'mcrw *mc-lower* *mc-upper* nsamples))
liam's avatar
liam committed
221 222

(defun random-walk-miser-example (&optional (nsamples 500000))
223
  (monte-carlo-integrate-miser 'mcrw *mc-lower* *mc-upper* nsamples))
liam's avatar
liam committed
224 225

(defun random-walk-vegas-example (&optional (nsamples 500000))
226
  (monte-carlo-integrate-vegas 'mcrw *mc-lower* *mc-upper* nsamples))
liam's avatar
liam committed
227

Liam Healy's avatar
Liam Healy committed
228
(save-test monte-carlo
liam's avatar
liam committed
229 230 231
  (random-walk-plain-example)
  (random-walk-miser-example)
  (random-walk-vegas-example))