monte-carlo.lisp 9.01 KB
Newer Older
1 2
;; Monte Carlo Integration
;; Liam Healy Sat Feb  3 2007 - 17:42
3
;; Time-stamp: <2008-12-26 18:59:23EST monte-carlo.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
liam's avatar
liam committed
5 6 7

(in-package :gsl)

8 9 10 11 12
;;; /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
13 14 15 16
;;;;****************************************************************************
;;;; PLAIN Monte Carlo
;;;;****************************************************************************

17 18 19 20 21 22 23
(defmobject monte-carlo-plain
    "gsl_monte_plain"
  ((dim sizet))
  "plain Monte Carlo integration"				; FDL
  "Make and initialize a workspace for Monte Carlo integration in dim dimensions."
  "init"
  nil)
liam's avatar
liam committed
24

25 26 27
(cffi:defcstruct plain-state
  (dim sizet)
  (x :pointer))
liam's avatar
liam committed
28

29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
(defmfun monte-carlo-integrate-plain
    (function lower-limits upper-limits calls generator state)
  "gsl_monte_plain_integrate"
  ((function :pointer)
   ((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
   ((dim0 lower-limits) sizet) (calls sizet)
   ((mpointer generator) :pointer)
   ((mpointer state) :pointer)
   (result :double) (abserr :double))
  :inputs (lower-limits upper-limits)
  :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
   of function calls calls, and obtains random sampling points using
   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
50 51 52 53 54 55 56 57 58
;;;;****************************************************************************
;;;; 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.

59 60 61 62 63 64 65 66 67 68
(defmobject monte-carlo-miser
    "gsl_monte_miser"
  ((dim sizet))
  "miser Monte Carlo integration"				; FDL
  "Make and initialize a workspace for Monte Carlo integration in
   dim dimensions.  The workspace is used to maintain
   the state of the integration."
  "init"
  nil)

liam's avatar
liam committed
69
(cffi:defcstruct miser-state
70 71
  (min-calls sizet)
  (min-calls-per-bisection sizet)
liam's avatar
liam committed
72 73 74
  (dither :double)
  (estimate-frac :double)
  (alpha :double)
75
  (dim sizet)
liam's avatar
liam committed
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
  (estimate-style :int)
  (depth :int)
  (verbose :int)
  (x :pointer)
  (xmid :pointer)
  (sigma-l :pointer)
  (sigma-r :pointer)
  (fmax-l :pointer)
  (fmax-r :pointer)
  (fmin-l :pointer)
  (fmin-r :pointer)
  (fsum-l :pointer)
  (fsum-r :pointer)
  (fsum2-l :pointer)
  (fsum2-r :pointer)
  (hits-l :pointer)
  (hits-r :pointer))

(export 'miser-parameter)
(defmacro miser-parameter (workspace parameter)
liam's avatar
liam committed
96
  ;; FDL
liam's avatar
liam committed
97 98 99 100 101 102
  "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)
 `(foreign-slot-value ,workspace 'miser-state ',parameter))

liam's avatar
liam committed
103
(defmfun monte-carlo-integrate-miser
liam's avatar
liam committed
104 105 106
    (function lower-limits upper-limits calls generator state)
  "gsl_monte_miser_integrate"
  ((function :pointer)
107 108
   ((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
   ((dim0 lower-limits) sizet) (calls sizet)
liam's avatar
liam committed
109
   ((generator generator) :pointer)
110
   ((mpointer state) :pointer)
liam's avatar
liam committed
111
   (result :double) (abserr :double))
112
  :inputs (lower-limits upper-limits)
liam's avatar
liam committed
113
  :documentation			; FDL
liam's avatar
liam committed
114
  "Uses the miser Monte Carlo algorithm to integrate the
liam's avatar
liam committed
115
   function f over the hypercubic region defined by the
liam's avatar
liam committed
116
   lower and upper limits in the arrays 'lower-limits and
liam's avatar
liam committed
117
   'upper-limits, each a gsl-vector of the samelength
liam's avatar
liam committed
118
   The integration uses a fixed number
liam's avatar
liam committed
119
   of function calls calls, and obtains random sampling points using
liam's avatar
liam committed
120
   the random number generator 'generator. A previously allocated workspace
liam's avatar
liam committed
121 122
   'state must be supplied.  The result of the integration is returned
   with an estimated absolute error.")
liam's avatar
liam committed
123 124 125 126 127

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

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

133 134 135 136 137 138 139 140 141 142
(defmobject monte-carlo-vegas
    "gsl_monte_vegas"
  ((dim sizet))
  "vegas Monte Carlo integration"				; FDL
  "Make and initialize a workspace for Monte Carlo integration in
   dim dimensions.  The workspace is used to maintain
   the state of the integration.  Returns a pointer to vegas-state."
  "init"
  nil)

liam's avatar
liam committed
143 144
(cffi:defcstruct vegas-state
  ;; grid 
145 146
  (dim sizet)
  (bins-max sizet)
liam's avatar
liam committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
  (bins :uint)				; uint
  (boxes :uint)		       ; these are both counted along the axes
  (xi :pointer)
  (xin :pointer)
  (delx :pointer)
  (weight :pointer)
  (vol :double)
  (x :pointer)
  (bin :pointer)
  (box :pointer)
  (d :pointer)				; distribution
  ;; control variables
  (alpha :double)
  (mode :int)
  (verbose :int)
  (iterations :uint)
  (stage :int)
  ;; scratch variables preserved between calls to vegas1/2/3
  (jac :double)
  (wtd-int-sum :double)
  (sum-wgts :double)
  (chi-sum :double)
  (chisq :double)
  (result :double)
  (sigma :double)
  (it-start :uint)
  (it-num :uint)
  (samples :uint)
  (calls-per-box :uint)
  (ostream :pointer))

(export 'vegas-parameter)
(defmacro vegas-parameter (workspace parameter)
liam's avatar
liam committed
180
  ;; FDL
liam's avatar
liam committed
181 182 183 184 185 186
  "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)
 `(foreign-slot-value ,workspace 'vegas-state ',parameter))

liam's avatar
liam committed
187
(defmfun monte-carlo-integrate-vegas
liam's avatar
liam committed
188 189 190
    (function lower-limits upper-limits calls generator state)
  "gsl_monte_vegas_integrate"
  ((function :pointer)
191 192
   ((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
   ((dim0 lower-limits) sizet) (calls sizet)
liam's avatar
liam committed
193
   ((generator generator) :pointer)
194
   ((mpointer state) :pointer)
liam's avatar
liam committed
195
   (result :double) (abserr :double))
196
  :inputs (lower-limits upper-limits)
liam's avatar
liam committed
197 198 199 200
  :documentation			; FDL
  "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
201
   xu, each of the same length.  The integration uses a fixed number
liam's avatar
liam committed
202 203 204 205
   of function calls calls, 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
liam's avatar
liam committed
206 207
   and its error estimate are based on a weighted average of independent
   samples. The chi-squared per degree of freedom for the weighted average
liam's avatar
liam committed
208
   is returned via the state struct component, s->chisq, and must be
liam's avatar
liam committed
209 210
   consistent with 1 for the weighted average to be reliable.")

211 212 213 214 215 216
;;;;****************************************************************************
;;;; Callback definition
;;;;****************************************************************************

(cffi:defcstruct monte-function
  (function :pointer)
217
  (dimensions sizet)
218 219 220
  (parameters :pointer))

(export 'def-mc-function)
221
(defmacro def-mc-function (name dimensions)
222
  `(def-single-function ,name :double :pointer monte-function
223
    ((dimensions ,dimensions))))
224

liam's avatar
liam committed
225 226 227 228 229 230 231
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

;;; Example from Sec. 23.5
;;; This is a function that occurs in random walk studies.

232
(defun monte-carlo-g (arg)
233
  (with-c-double (arg x y z)
234 235 236
    (* (/ (expt pi 3))
       (/ (- 1 (* (cos x) (cos y) (cos z)))))))

237
(def-mc-function monte-carlo-g 3)
liam's avatar
liam committed
238 239

(defun random-walk-plain-example (&optional (nsamples 500000))
240
  (let ((ws (make-monte-carlo-plain 3))
241 242 243
	(lower #m(0.0d0 0.0d0 0.0d0))
	(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
	(rng (make-random-number-generator *mt19937* 0)))
244
    (monte-carlo-integrate-plain monte-carlo-g lower upper nsamples rng ws)))
liam's avatar
liam committed
245 246

(defun random-walk-miser-example (&optional (nsamples 500000))
247 248 249 250
  (let ((ws (make-monte-carlo-miser 3))
	(lower #m(0.0d0 0.0d0 0.0d0))
	(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
	(rng (make-random-number-generator *mt19937* 0)))
251
    (monte-carlo-integrate-miser monte-carlo-g lower upper nsamples rng ws)))
liam's avatar
liam committed
252 253

(defun random-walk-vegas-example (&optional (nsamples 500000))
254
  (let ((ws (make-monte-carlo-vegas 3))
255 256 257
	(lower #m(0.0d0 0.0d0 0.0d0))
	(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
	(rng (make-random-number-generator *mt19937* 0)))
258
    (monte-carlo-integrate-vegas monte-carlo-g lower upper nsamples rng ws)))
liam's avatar
liam committed
259

Liam Healy's avatar
Liam Healy committed
260
(save-test monte-carlo
liam's avatar
liam committed
261 262 263
  (random-walk-plain-example)
  (random-walk-miser-example)
  (random-walk-vegas-example))