gillespie.hs

import qualified Data.Map as Map
import Random (randomRIO)

type Mol = String
type Mixture = Map.Map Mol Integer
type Rxn = (Mixture, (Mixture -> Mixture)) -- (Lhs, Action)
type State = (Double, Integer, Mixture) -- (Time, NumSteps, Mixture)

getActivities :: [Rxn] -> Mixture -> [Integer]
getActivities rxns mixture = []

selectRxn :: [Rxn] -> [Integer] -> Integer -> IO Rxn
selectRxn rxns activities totalActivity = return $ head rxns
--  do rnd <- getStdRandom $ randomR (0.0, totalActivity)

step :: [Rxn] -> State -> IO State
step rxns (t, ns, mixture) =
  let activities = getActivities rxns mixture
      totalActivity = sum activities
  in if totalActivity == 0
    then return (t, ns+1, mixture)
    else do r <- selectRxn rxns activities totalActivity
            rnd <- randomRIO (0.0, 1.0)
            let dt = log (1.0/rnd) / totalActivity
            return (t+dt, ns+1, (snd r) mixture)

simulate :: [Rxn] -> Mixture -> Integer -> IO [State]
simulate rxns mixture numSteps = take numSteps (iterate (step rxns) (0.0, 0, mixture))