my implementation of the gillespie algorithm

gillespie

;;;; public domain, by Bike
;;;; Gillespie, Daniel T. (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry 81 (25): 2340–2361. doi:10.1021/j100540a008.

(defun gillespie (species reactions
		  &key (stop-time 0d0 st-p) (stop-reactions 0 sr-p)
		    (rpd 0 rpd-p))
  ;; SPECIES is a vector of N positive integers, molecule counts.
  ;; REACTIONS is a vector of M lists
  ;;  The first of each list is the rate constant for that reaction.
  ;;  The second is a vector of N nonnegative integers describing the substrates.
  ;;   For example, 2A + B + C -> 2C would be #(2 1 1).
  ;;  The third is a vector of N integers describing the total stoichiometry.
  ;;   For example, 2A + B + C -> 2C would be #(-2 -1 1).
  ;; After TIME has elasped, the altered SPECIES vector is returned.
  (declare (non-negative-fixnum rpd stop-reactions))
  (setf species (coerce species '(simple-array non-negative-fixnum (*))))
  (setf reactions
	(map '(simple-array list (*))
	     (lambda (r)
	       (list (coerce (first r) 'double-float)
		     (coerce (second r) '(simple-array non-negative-fixnum (*)))
		     (coerce (third r) '(simple-array fixnum (*)))))
	     reactions))
  (setf stop-time (coerce stop-time 'double-float))
  (locally
      (declare (type (simple-array non-negative-fixnum (*)) species)
	       ;; LIST probably upgrades to T, but having a known UAET helps AREF
	       (type (simple-array list (*)) reactions)
	       #+(or)(optimize (debug 3))
	       (optimize (speed 3) (space 3)))
    (flet ((record (tau species)
	     (format t "~f" tau)
	     (map nil (lambda (s) (write-char #\Space) (write s)) species)
	     (terpri))
	   (reaction-rate (reaction)
	     (let ((constant (first reaction))
		   (substrate (second reaction)))
	       (declare (type double-float constant)
			(type (simple-array non-negative-fixnum (*)) substrate))
	       (loop for i below (length species)
		  for x of-type non-negative-fixnum
		    = (let ((order (aref substrate i)))
			(case order
			  ((0) 1) ; substrate does not participate
			  ((1) (aref species i)) ; premature optimization
			  (t (binomial-coefficient
			      (aref species i)
			      (aref substrate i)))))
		  for product of-type non-negative-fixnum = x
		  then (logand (* product x) most-positive-fixnum)
		  finally (return (* constant product)))))
	   (react (reaction)
	     (map-into species #'+ species
		       (the (simple-array fixnum (*)) (third reaction)))))
      (declare (ftype (function (list) double-float) reaction-rate))
      (with-simple-restart (finish "Stop the loop and return whatever the present species are.")
	(loop
	   for count of-type non-negative-fixnum from 0
	   for tau of-type double-float = 0d0
	   then (+ tau (/ (- (log (random 1d0))) rate-sum)) ; 21a
	   for rates of-type (simple-array double-float (*))
	     = (map '(simple-array double-float (*)) #'reaction-rate reactions)
	   then (map-into rates #'reaction-rate reactions)
	   for rate-sum of-type double-float = (reduce #'+ rates)
	   when (and rpd-p (zerop (mod count rpd)))
	   do (record tau species)
	   while (or (not sr-p) (> stop-reactions count))
	   do (let ((rn ; 21b
		     (loop with target of-type double-float = (* (random 1d0) rate-sum)
			for i from 0
			when (= (length rates) i) ; out of rxns
			do (error "No reactions can proceed.")
			summing (aref rates i) into total of-type double-float
			until (> total target)
			finally (return i))))
		(react (aref reactions rn)))
	   while (or (not st-p) (< tau stop-time))))))
  species)

;;; Example invocations
;; Exponential decay
;; (gillespie (vector 1000 0) #((5d-1 #(1 0) #(1 -1))) :stop-time 10d0)
;; -> #(2 998), or thereabouts
;; Lotka, out to file
#+(or)
(with-open-file (*standard-output* "~/gillespie_lotka.out" :direction :output
				   :if-does-not-exist :create)
  (gillespie (vector 1000 1000)
	     #((10 #(1 0) #(2 0)) (0.01 #(1 1) #(-1 1)) (10 #(0 1) #(0 -1)))
	     :stop-reactions 1000000
	     :rpd 100))
;; plus gnuplot, and bam, sick-ass cycles.
;; I'd like it to not cons in the locally, except for the initial rates.

(defun gillespie-source (n-species reactions)
  (let ((rate-names (loop for i across reactions collecting (gensym "RATE"))))
    `(lambda (species stop-time)
       (declare (type (simple-array non-negative-fixnum (*)) species)
		(type double-float stop-time)
		(optimize (speed 3)))
       (loop
	  for count of-type non-negative-fixnum from 0
	  for tau of-type double-float = 0d0
	  then (+ tau (/ (- (log (random 1d0))) rate-sum))
	    ,@(loop for name in rate-names for reaction across reactions
		 appending `(for ,name of-type double-float
				 = (* ,(first reaction)
				      ,@(loop for n below n-species
					   appending
					     (let ((order
						    (aref (second reaction) n)))
					       ;; let the compiler optimize shit.
					       `(,(/ (factorial order))
						  ,@(loop for o below order
						       collecting
							 `(- (aref species ,n) ,o))))))))
	  for rate-sum of-type double-float = (+ ,@rate-names)
	  do (let ((target (* (random 1d0) rate-sum))
		   (total 0))
	       (cond ,@(loop for name in rate-names for reaction across reactions
			  collecting `((> (incf total ,name) target)
				       ,@(loop for n below n-species
					    collecting `(incf (aref species ,n)
							      ,(aref (third reaction) n)))))
		     (t (error "No reactions can proceed."))))
	  while (< tau stop-time))
       species)))

;;; Example invocation
;; (compile nil (gillespie-source 2 #((5d-1 #(1 0) #(-1 1)))))
;; (funcall * (make-array 2 :element-type 'non-negative-fixnum :initial-contents '(1000 0)))