aivika (empty) → 0.1

raw patch · 15 files changed

+2625/−0 lines, 15 filesdep +arraydep +basedep +mtlsetup-changedbinary-added

Dependencies added: array, base, mtl

Files

• LICENSE +30/−0
• Setup.lhs +3/−0
• Simulation/Aivika/Dynamics.hs +1628/−0
• Simulation/Aivika/PriorityQueue.hs +173/−0
• Simulation/Aivika/Queue.hs +112/−0
• aivika.cabal +47/−0
• doc/aivika.pdf binary
• examples/BassDiffusion.hs +106/−0
• examples/ChemicalReaction.hs +26/−0
• examples/FishBank.hs +57/−0
• examples/MachRep1.hs +69/−0
• examples/MachRep1EventDriven.hs +75/−0
• examples/MachRep1TimeDriven.hs +116/−0
• examples/MachRep2.hs +96/−0
• examples/MachRep3.hs +87/−0

@@ -0,0 +1,30 @@+Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.

@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain

@@ -0,0 +1,1628 @@++-- Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>+-- +-- All rights reserved.+-- +-- Redistribution and use in source and binary forms, with or without+-- modification, are permitted provided that the following conditions+-- are met:+-- +-- 1. Redistributions of source code must retain the above copyright+--    notice, this list of conditions and the following disclaimer.+-- +-- 2. Redistributions in binary form must reproduce the above copyright+--    notice, this list of conditions and the following disclaimer in the+--    documentation and/or other materials provided with the distribution.+-- +-- 3. Neither the name of the author nor the names of his contributors+--    may be used to endorse or promote products derived from this software+--    without specific prior written permission.+-- +-- THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+-- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+-- ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+-- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+-- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+-- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+-- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+-- SUCH DAMAGE.++{-# LANGUAGE FlexibleContexts, FlexibleInstances, UndecidableInstances #-}++-- |+-- Module     : Simulation.Aivika.Dynamics+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 6.12.1+--+-- Aivika is a multi-paradigm simulation library. It allows us to integrate +-- a system of ordinary differential equations. Also it can be applied to+-- the Discrete Event Simulation. It supports the event-oriented, +-- process-oriented and activity-oriented paradigms. Aivika also supports +-- the Agent-based Modeling. Finally, it can be applied to System Dynamics.+--+module Simulation.Aivika.Dynamics +       (-- * Dynamics+        Dynamics,+        DynamicsTrans(..),+        Specs(..),+        Method(..),+        runDynamics1,+        runDynamics,+        runDynamicsIO,+        -- ** Time parameters+        starttime,+        stoptime,+        dt,+        time,+        -- ** Maximum and Minimum+        maxD,+        minD,+        -- ** Integrals+        Integ,+        newInteg,+        integInit,+        integValue,+        integDiff,+        -- ** Table Functions+        lookupD,+        lookupStepwiseD,+        -- ** Interpolation+        initD,+        discrete,+        interpolate,+        -- ** Memoization and Sequential Calculations+        Memo,+        UMemo,+        memo,+        umemo,+        memo0,+        umemo0,+        -- ** Utility+        once,+        -- * Event Queue+        DynamicsQueue,+        newQueue,+        enqueueDC,+        enqueueD,+        runQueue,+        -- * References+        DynamicsRef,+        newRef,+        refQueue,+        readRef,+        writeRef,+        writeRef',+        modifyRef,+        modifyRef',+        -- * Discontinuous Processes+        DynamicsPID,+        DynamicsProc,+        newPID,+        pidQueue,+        holdProcD,+        holdProc,+        passivateProc,+        procPassive,+        reactivateProc,+        procPID,+        runProc,+        -- * Resources+        DynamicsResource,+        newResource,+        resourceQueue,+        resourceInitCount,+        resourceCount,+        requestResource,+        releaseResource,+        -- * Agent-based Modeling+        Agent,+        AgentState,+        newAgent,+        newState,+        newSubstate,+        agentQueue,+        agentState,+        activateState,+        initState,+        stateAgent,+        stateParent,+        addTimeoutD,+        addTimeout,+        addTimerD,+        addTimer,+        stateActivation,+        stateDeactivation) where++import Data.Array+import Data.Array.IO+import Data.IORef+import Control.Monad+import Control.Monad.Trans++import qualified Simulation.Aivika.Queue as Q+import qualified Simulation.Aivika.PriorityQueue as PQ++--+-- The Dynamics Monad+--+-- A value of the Dynamics monad represents an abstract dynamic +-- process, i.e. a time varying polymorphic function. This is +-- a key point of the Aivika simulation library.+--++-- | A value in the 'Dynamics' monad represents a dynamic process, i.e.+-- a polymorphic time varying function.+newtype Dynamics a = Dynamics (Parameters -> IO a)++-- | It defines the simulation time appended with additional information.+data Parameters = Parameters { parSpecs :: Specs,    -- ^ the simulation specs+                               parTime :: Double,    -- ^ the current time+                               parIteration :: Int,  -- ^ the current iteration+                               parPhase :: Int }     -- ^ the current phase++-- | It defines the simulation specs.+data Specs = Specs { spcStartTime :: Double,    -- ^ the start time+                     spcStopTime :: Double,     -- ^ the stop time+                     spcDT :: Double,           -- ^ the integration time step+                     spcMethod :: Method        -- ^ the integration method+                   } deriving (Eq, Ord, Show)++-- | It defines the integration method.+data Method = Euler          -- ^ Euler's method+            | RungeKutta2    -- ^ the 2nd order Runge-Kutta method+            | RungeKutta4    -- ^ the 4th order Runge-Kutta method+            deriving (Eq, Ord, Show)++iterations :: Specs -> [Int]+iterations sc = [i1 .. i2] where+  i1 = 0+  i2 = round ((spcStopTime sc - +               spcStartTime sc) / spcDT sc)++iterationBnds :: Specs -> (Int, Int)+iterationBnds sc = (0, round ((spcStopTime sc - +                               spcStartTime sc) / spcDT sc))++iterationLoBnd :: Specs -> Int+iterationLoBnd sc = 0++iterationHiBnd :: Specs -> Int+iterationHiBnd sc = round ((spcStopTime sc - +                            spcStartTime sc) / spcDT sc)++phases :: Specs -> [Int]+phases sc = +  case spcMethod sc of+    Euler -> [0]+    RungeKutta2 -> [0, 1]+    RungeKutta4 -> [0, 1, 2, 3]++phaseBnds :: Specs -> (Int, Int)+phaseBnds sc = +  case spcMethod sc of+    Euler -> (0, 0)+    RungeKutta2 -> (0, 1)+    RungeKutta4 -> (0, 3)++phaseLoBnd :: Specs -> Int+phaseLoBnd sc = 0+                  +phaseHiBnd :: Specs -> Int+phaseHiBnd sc = +  case spcMethod sc of+    Euler -> 0+    RungeKutta2 -> 1+    RungeKutta4 -> 3++basicTime :: Specs -> Int -> Int -> Double+basicTime sc n ph =+  if ph < 0 then +    error "Incorrect phase: basicTime"+  else+    spcStartTime sc + n' * spcDT sc + delta (spcMethod sc) ph +      where n' = fromInteger (toInteger n)+            delta Euler       0 = 0+            delta RungeKutta2 0 = 0+            delta RungeKutta2 1 = spcDT sc+            delta RungeKutta4 0 = 0+            delta RungeKutta4 1 = spcDT sc / 2+            delta RungeKutta4 2 = spcDT sc / 2+            delta RungeKutta4 3 = spcDT sc++neighborhood :: Specs -> Double -> Double -> Bool+neighborhood sc t t' = +  abs (t - t') <= spcDT sc / 1.0e6++instance Monad Dynamics where+  return  = returnD+  m >>= k = bindD m k++returnD :: a -> Dynamics a+returnD a = Dynamics (\ps -> return a)++bindD :: Dynamics a -> (a -> Dynamics b) -> Dynamics b+bindD (Dynamics m) k = +  Dynamics $ \ps -> +  do a <- m ps+     let Dynamics m' = k a+     m' ps++subrunDynamics1 :: Dynamics a -> Specs -> IO a+subrunDynamics1 (Dynamics m) sc =+  do let n = iterationHiBnd sc+         t = basicTime sc n 0+     m Parameters { parSpecs = sc,+                    parTime = t,+                    parIteration = n,+                    parPhase = 0 }++subrunDynamics :: Dynamics a -> Specs -> [IO a]+subrunDynamics (Dynamics m) sc =+  do let (nl, nu) = iterationBnds sc+         parameterise n = Parameters { parSpecs = sc,+                                       parTime = basicTime sc n 0,+                                       parIteration = n,+                                       parPhase = 0 }+     map (m . parameterise) [nl .. nu]++-- | Run the simulation and return the result in the last +-- time point using the specified simulation specs.+runDynamics1 :: Dynamics (Dynamics a) -> Specs -> IO a+runDynamics1 (Dynamics m) sc = +  do d <- m Parameters { parSpecs = sc,+                         parTime = spcStartTime sc,+                         parIteration = 0,+                         parPhase = 0 }+     subrunDynamics1 d sc++-- | Run the simulation and return the results in all +-- integration time points using the specified simulation specs.+runDynamics :: Dynamics (Dynamics a) -> Specs -> IO [a]+runDynamics (Dynamics m) sc = +  do d <- m Parameters { parSpecs = sc,+                         parTime = spcStartTime sc,+                         parIteration = 0,+                         parPhase = 0 }+     sequence $ subrunDynamics d sc++-- | Run the simulation and return the results in all +-- integration time points using the specified simulation specs.+runDynamicsIO :: Dynamics (Dynamics a) -> Specs -> IO [IO a]+runDynamicsIO (Dynamics m) sc =+  do d <- m Parameters { parSpecs = sc,+                         parTime = spcStartTime sc,+                         parIteration = 0,+                         parPhase = 0 }+     return $ subrunDynamics d sc++instance Functor Dynamics where+  fmap f (Dynamics m) = +    Dynamics $ \ps -> do { a <- m ps; return $ f a }++instance Eq (Dynamics a) where+  x == y = error "Can't compare dynamics." ++instance Show (Dynamics a) where+  showsPrec _ x = showString "<< Dynamics >>"++liftMD :: (a -> b) -> Dynamics a -> Dynamics b+liftMD f (Dynamics x) =+  Dynamics $ \ps -> do { a <- x ps; return $ f a }++liftM2D :: (a -> b -> c) -> Dynamics a -> Dynamics b -> Dynamics c+liftM2D f (Dynamics x) (Dynamics y) =+  Dynamics $ \ps -> do { a <- x ps; b <- y ps; return $ f a b }++instance (Num a) => Num (Dynamics a) where+  x + y = liftM2D (+) x y+  x - y = liftM2D (-) x y+  x * y = liftM2D (*) x y+  negate = liftMD negate+  abs = liftMD abs+  signum = liftMD signum+  fromInteger i = return $ fromInteger i++instance (Fractional a) => Fractional (Dynamics a) where+  x / y = liftM2D (/) x y+  recip = liftMD recip+  fromRational t = return $ fromRational t++instance (Floating a) => Floating (Dynamics a) where+  pi = return pi+  exp = liftMD exp+  log = liftMD log+  sqrt = liftMD sqrt+  x ** y = liftM2D (**) x y+  sin = liftMD sin+  cos = liftMD cos+  tan = liftMD tan+  asin = liftMD asin+  acos = liftMD acos+  atan = liftMD atan+  sinh = liftMD sinh+  cosh = liftMD cosh+  tanh = liftMD tanh+  asinh = liftMD asinh+  acosh = liftMD acosh+  atanh = liftMD atanh++instance MonadIO Dynamics where+  liftIO m = Dynamics $ const m++-- | The 'DynamicsTrans' class defines a type which the 'Dynamics' +-- computation can be lifted to.+class DynamicsTrans m where+  -- | Lift the computation.+  liftD :: Dynamics a -> m a++--+-- Integration Parameters and Time+--++-- | Return the start simulation time.+starttime :: Dynamics Double+starttime = Dynamics $ return . spcStartTime . parSpecs++-- | Return the stop simulation time.+stoptime :: Dynamics Double+stoptime = Dynamics $ return . spcStopTime . parSpecs++-- | Return the integration time step.+dt :: Dynamics Double+dt = Dynamics $ return . spcDT . parSpecs++-- | Return the current simulation time.+time :: Dynamics Double+time = Dynamics $ return . parTime ++--+-- Maximum and Minimum+--++-- | Return the maximum.+maxD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a+maxD = liftM2D max++-- | Return the minimum.+minD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a+minD = liftM2D min+           +--+-- System Dynamics+--++-- | The 'Integ' type represents an integral.+data Integ = Integ { integInit     :: Dynamics Double,   -- ^ The initial value.+                     integExternal :: IORef (Dynamics Double),+                     integInternal :: IORef (Dynamics Double) }++-- | Create a new integral with the specified initial value.+newInteg :: Dynamics Double -> Dynamics Integ+newInteg i = +  do r1 <- liftIO $ newIORef $ initD i +     r2 <- liftIO $ newIORef $ initD i +     let integ = Integ { integInit     = i, +                         integExternal = r1,+                         integInternal = r2 }+         z = Dynamics $ \ps -> +           do (Dynamics m) <- readIORef (integInternal integ)+              m ps+     y <- umemo interpolate z+     liftIO $ writeIORef (integExternal integ) y+     return integ++-- | Return the integral's value.+integValue :: Integ -> Dynamics Double+integValue integ = +  Dynamics $ \ps ->+  do (Dynamics m) <- readIORef (integExternal integ)+     m ps++-- | Set the derivative for the integral.+integDiff :: Integ -> Dynamics Double -> Dynamics ()+integDiff integ diff =+  do let z = Dynamics $ \ps ->+           do y <- readIORef (integExternal integ)+              let i = integInit integ+              case spcMethod (parSpecs ps) of+                Euler -> integEuler diff i y ps+                RungeKutta2 -> integRK2 diff i y ps+                RungeKutta4 -> integRK4 diff i y ps+     liftIO $ writeIORef (integInternal integ) z++integEuler :: Dynamics Double+             -> Dynamics Double +             -> Dynamics Double +             -> Parameters -> IO Double+integEuler (Dynamics f) (Dynamics i) (Dynamics y) ps = +  case parIteration ps of+    0 -> +      i ps+    n -> do +      let sc  = parSpecs ps+          ty  = basicTime sc (n - 1) 0+          psy = ps { parTime = ty, parIteration = n - 1, parPhase = 0 }+      a <- y psy+      b <- f psy+      let !v = a + spcDT (parSpecs ps) * b+      return v++integRK2 :: Dynamics Double+           -> Dynamics Double+           -> Dynamics Double+           -> Parameters -> IO Double+integRK2 (Dynamics f) (Dynamics i) (Dynamics y) ps =+  case parPhase ps of+    0 -> case parIteration ps of+      0 ->+        i ps+      n -> do+        let sc = parSpecs ps+            ty = basicTime sc (n - 1) 0+            t1 = ty+            t2 = basicTime sc (n - 1) 1+            psy = ps { parTime = ty, parIteration = n - 1, parPhase = 0 }+            ps1 = psy+            ps2 = ps { parTime = t2, parIteration = n - 1, parPhase = 1 }+        vy <- y psy+        k1 <- f ps1+        k2 <- f ps2+        let !v = vy + spcDT sc / 2.0 * (k1 + k2)+        return v+    1 -> do+      let sc = parSpecs ps+          n  = parIteration ps+          ty = basicTime sc n 0+          t1 = ty+          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }+          ps1 = psy+      vy <- y psy+      k1 <- f ps1+      let !v = vy + spcDT sc * k1+      return v+    _ -> +      error "Incorrect phase: integ"++integRK4 :: Dynamics Double+           -> Dynamics Double+           -> Dynamics Double+           -> Parameters -> IO Double+integRK4 (Dynamics f) (Dynamics i) (Dynamics y) ps =+  case parPhase ps of+    0 -> case parIteration ps of+      0 -> +        i ps+      n -> do+        let sc = parSpecs ps+            ty = basicTime sc (n - 1) 0+            t1 = ty+            t2 = basicTime sc (n - 1) 1+            t3 = basicTime sc (n - 1) 2+            t4 = basicTime sc (n - 1) 3+            psy = ps { parTime = ty, parIteration = n - 1, parPhase = 0 }+            ps1 = psy+            ps2 = ps { parTime = t2, parIteration = n - 1, parPhase = 1 }+            ps3 = ps { parTime = t3, parIteration = n - 1, parPhase = 2 }+            ps4 = ps { parTime = t4, parIteration = n - 1, parPhase = 3 }+        vy <- y psy+        k1 <- f ps1+        k2 <- f ps2+        k3 <- f ps3+        k4 <- f ps4+        let !v = vy + spcDT sc / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)+        return v+    1 -> do+      let sc = parSpecs ps+          n  = parIteration ps+          ty = basicTime sc n 0+          t1 = ty+          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }+          ps1 = psy+      vy <- y psy+      k1 <- f ps1+      let !v = vy + spcDT sc / 2.0 * k1+      return v+    2 -> do+      let sc = parSpecs ps+          n  = parIteration ps+          ty = basicTime sc n 0+          t2 = basicTime sc n 1+          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }+          ps2 = ps { parTime = t2, parIteration = n, parPhase = 1 }+      vy <- y psy+      k2 <- f ps2+      let !v = vy + spcDT sc / 2.0 * k2+      return v+    3 -> do+      let sc = parSpecs ps+          n  = parIteration ps+          ty = basicTime sc n 0+          t3 = basicTime sc n 2+          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }+          ps3 = ps { parTime = t3, parIteration = n, parPhase = 2 }+      vy <- y psy+      k3 <- f ps3+      let !v = vy + spcDT sc * k3+      return v+    _ -> +      error "Incorrect phase: integ"++-- smoothI :: Dynamics Double -> Dynamics Double -> Dynamics Double +--           -> Dynamics Double+-- smoothI x t i = y where+--   y = integ ((x - y) / t) i++-- smooth :: Dynamics Double -> Dynamics Double -> Dynamics Double+-- smooth x t = smoothI x t x++-- smooth3I :: Dynamics Double -> Dynamics Double -> Dynamics Double +--            -> Dynamics Double+-- smooth3I x t i = y where+--   y  = integ ((s1 - y) / t') i+--   s1 = integ ((s0 - s1) / t') i+--   s0 = integ ((x - s0) / t') i+--   t' = t / 3.0++-- smooth3 :: Dynamics Double -> Dynamics Double -> Dynamics Double+-- smooth3 x t = smooth3I x t x++-- smoothNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double +--            -> Dynamics Double+-- smoothNI x t n i = s ! n where+--   s   = array (1, n) [(k, f k) | k <- [1 .. n]]+--   f 0 = integ ((x - s ! 0) / t') i+--   f k = integ ((s ! (k - 1) - s ! k) / t') i+--   t'  = t / fromIntegral n++-- smoothN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double+-- smoothN x t n = smoothNI x t n x++-- delay1I :: Dynamics Double -> Dynamics Double -> Dynamics Double +--           -> Dynamics Double+-- delay1I x t i = y where+--   y = integ (x - y) (i * t) / t++-- delay1 :: Dynamics Double -> Dynamics Double -> Dynamics Double+-- delay1 x t = delay1I x t x++-- delay3I :: Dynamics Double -> Dynamics Double -> Dynamics Double +--           -> Dynamics Double+-- delay3I x t i = y where+--   y  = integ (s1 - y) (i * t') / t'+--   s1 = integ (s0 - s1) (i * t') / t'+--   s0 = integ (x - s0) (i * t') / t'+--   t' = t / 3.0++-- delay3 :: Dynamics Double -> Dynamics Double -> Dynamics Double+-- delay3 x t = delay3I x t x++-- delayNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double +--           -> Dynamics Double+-- delayNI x t n i = s ! n where+--   s   = array (1, n) [(k, f k) | k <- [1 .. n]]+--   f 0 = integ (x - s ! 0) (i * t') / t'+--   f k = integ (s ! (k - 1) - s ! k) (i * t') / t'+--   t'  = t / fromIntegral n++-- delayN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double+-- delayN x t n = delayNI x t n x++-- forecast :: Dynamics Double -> Dynamics Double -> Dynamics Double +--            -> Dynamics Double+-- forecast x at hz =+--   x * (1.0 + (x / smooth x at - 1.0) / at * hz)++-- trend :: Dynamics Double -> Dynamics Double -> Dynamics Double +--         -> Dynamics Double+-- trend x at i =+--   (x / smoothI x at (x / (1.0 + i * at)) - 1.0) / at++--+-- Table Functions+--++-- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.+lookupD :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double+lookupD (Dynamics m) tbl =+  Dynamics (\ps -> do a <- m ps; return $ find first last a) where+    (first, last) = bounds tbl+    find left right x =+      if left > right then+        error "Incorrect index: table"+      else+        let index = (left + 1 + right) `div` 2+            x1    = fst $ tbl ! index+        in if x1 <= x then +             let y | index < right = find index right x+                   | right == last  = snd $ tbl ! right+                   | otherwise     = +                     let x2 = fst $ tbl ! (index + 1)+                         y1 = snd $ tbl ! index+                         y2 = snd $ tbl ! (index + 1)+                     in y1 + (y2 - y1) * (x - x1) / (x2 - x1) +             in y+           else+             let y | left < index = find left (index - 1) x+                   | left == first = snd $ tbl ! left+                   | otherwise    = error "Incorrect index: table"+             in y++-- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function.+lookupStepwiseD :: Dynamics Double -> Array Int (Double, Double)+                  -> Dynamics Double+lookupStepwiseD (Dynamics m) tbl =+  Dynamics (\ps -> do a <- m ps; return $ find first last a) where+    (first, last) = bounds tbl+    find left right x =+      if left > right then+        error "Incorrect index: table"+      else+        let index = (left + 1 + right) `div` 2+            x1    = fst $ tbl ! index+        in if x1 <= x then +             let y | index < right = find index right x+                   | right == last  = snd $ tbl ! right+                   | otherwise     = snd $ tbl ! right+             in y+           else+             let y | left < index = find left (index - 1) x+                   | left == first = snd $ tbl ! left+                   | otherwise    = error "Incorrect index: table"+             in y++-- --+-- -- Discrete Functions+-- --+    +-- delayTrans :: Dynamics a -> Dynamics Double -> Dynamics a +--               -> (Dynamics a -> Dynamics a) -> Dynamics a+-- delayTrans (Dynamics x) (Dynamics d) (Dynamics i) tr = tr $ Dynamics r +--   where+--     r ps = do +--       let t  = parTime ps+--           sc = parSpecs ps+--           n  = parIteration ps+--       a <- d ps+--       let t' = (t - a) - spcStartTime sc+--           n' = fromInteger $ toInteger $ floor $ t' / spcDT sc+--           y | n' < 0    = i $ ps { parTime = spcStartTime sc,+--                                    parIteration = 0, +--                                    parPhase = 0 }+--             | n' < n    = x $ ps { parTime = t',+--                                    parIteration = n',+--                                    parPhase = -1 }+--             | n' > n    = error "Cannot return the future data: delay"+--             | otherwise = error "Cannot return the current data: delay"+--       y    ++-- delay :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a+-- delay x d = delayTrans x d x $ memo0 discrete++-- delay' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a+-- delay' x d = delayTrans x d x $ memo0' discrete++-- delayI :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a+-- delayI x d i = delayTrans x d i $ memo0 discrete         ++-- delayI' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a+-- delayI' x d i = delayTrans x d i $ memo0' discrete         ++--+-- Interpolation and Initial Value+--+-- These functions complement the memoization, possibly except for +-- the initial function which can be also useful to get an initial +-- value of any dynamic process. See comments to the Memoization +-- section.+--      ++-- | Return the initial value.+initD :: Dynamics a -> Dynamics a+initD (Dynamics m) =+  Dynamics $ \ps ->+  if parIteration ps == 0 && parPhase ps == 0 then+    m ps+  else+    let sc = parSpecs ps+    in m $ ps { parTime = basicTime sc 0 0,+                parIteration = 0,+                parPhase = 0 } ++-- | Discretize the computation in the integration time points.+discrete :: Dynamics a -> Dynamics a+discrete (Dynamics m) =+  Dynamics $ \ps ->+  let ph = parPhase ps+      r | ph == 0    = m ps+        | ph > 0    = let sc = parSpecs ps+                          n  = parIteration ps+                      in m $ ps { parTime = basicTime sc n 0,+                                  parPhase = 0 }+        | otherwise = let sc = parSpecs ps+                          t  = parTime ps+                          n  = parIteration ps+                          t' = spcStartTime sc + fromIntegral (n + 1) * spcDT sc+                          n' = if neighborhood sc t t' then n + 1 else n+                      in m $ ps { parTime = basicTime sc n' 0,+                                  parIteration = n',+                                  parPhase = 0 }+  in r++-- | Interpolate the computation based on the integration time points only.+interpolate :: Dynamics Double -> Dynamics Double+interpolate (Dynamics m) = +  Dynamics $ \ps -> +  if parPhase ps >= 0 then +    m ps+  else +    let sc = parSpecs ps+        t  = parTime ps+        x  = (t - spcStartTime sc) / spcDT sc+        n1 = max (floor x) (iterationLoBnd sc)+        n2 = min (ceiling x) (iterationHiBnd sc)+        t1 = basicTime sc n1 0+        t2 = basicTime sc n2 0+        z1 = m $ ps { parTime = t1, +                      parIteration = n1, +                      parPhase = 0 }+        z2 = m $ ps { parTime = t2,+                      parIteration = n2,+                      parPhase = 0 }+        r | t == t1   = z1+          | t == t2   = z2+          | otherwise = +            do y1 <- z1+               y2 <- z2+               return $ y1 + (y2 - y1) * (t - t1) / (t2 - t1)+    in r++-- --+-- -- Memoization+-- --+-- -- The memoization creates such processes, which values are +-- -- defined and then stored in the cache for the points of+-- -- integration. You should use some kind of interpolation +-- -- like the interpolate function to process all other time +-- -- points that don't coincide with the integration points:+-- --+-- --   x = memo interpolate y    -- a linear interpolation+-- --   x = memo discrete y       -- a discrete process+-- --++-- | The 'Memo' class specifies a type for which an array can be created.+class (MArray IOArray e IO) => Memo e where+  newMemoArray_ :: Ix i => (i, i) -> IO (IOArray i e)++-- | The 'UMemo' class specifies a type for which an unboxed array exists.+class (MArray IOUArray e IO) => UMemo e where+  newMemoUArray_ :: Ix i => (i, i) -> IO (IOUArray i e)++instance Memo e where+  newMemoArray_ = newArray_+    +instance (MArray IOUArray e IO) => UMemo e where+  newMemoUArray_ = newArray_++-- | Memoize and order the computation in the integration time points using +-- the specified interpolation and being aware of the Runge-Kutta method.+memo :: Memo e => (Dynamics e -> Dynamics e) -> Dynamics e +       -> Dynamics (Dynamics e)+memo tr (Dynamics m) = +  Dynamics $ \ps ->+  do let sc = parSpecs ps+         (phl, phu) = phaseBnds sc+         (nl, nu)   = iterationBnds sc+     arr   <- newMemoArray_ ((phl, nl), (phu, nu))+     nref  <- newIORef 0+     phref <- newIORef 0+     let r ps = +           do let sc  = parSpecs ps+                  n   = parIteration ps+                  ph  = parPhase ps+                  phu = phaseHiBnd sc +                  loop n' ph' = +                    if (n' > n) || ((n' == n) && (ph' > ph)) +                    then +                      readArray arr (ph, n)+                    else +                      let ps' = ps { parIteration = n', parPhase = ph',+                                     parTime = basicTime sc n' ph' }+                      in do a <- m ps'+                            a `seq` writeArray arr (ph', n') a+                            if ph' >= phu +                              then do writeIORef phref 0+                                      writeIORef nref (n' + 1)+                                      loop (n' + 1) 0+                              else do writeIORef phref (ph' + 1)+                                      loop n' (ph' + 1)+              n'  <- readIORef nref+              ph' <- readIORef phref+              loop n' ph'+     return $ tr $ Dynamics r++-- | Memoize and order the computation in the integration time points using +-- the specified interpolation and being aware of the Runge-Kutta method.+umemo :: UMemo e => (Dynamics e -> Dynamics e) -> Dynamics e +        -> Dynamics (Dynamics e)+umemo tr (Dynamics m) = +  Dynamics $ \ps ->+  do let sc = parSpecs ps+         (phl, phu) = phaseBnds sc+         (nl, nu)   = iterationBnds sc+     arr   <- newMemoUArray_ ((phl, nl), (phu, nu))+     nref  <- newIORef 0+     phref <- newIORef 0+     let r ps =+           do let sc  = parSpecs ps+                  n   = parIteration ps+                  ph  = parPhase ps+                  phu = phaseHiBnd sc +                  loop n' ph' = +                    if (n' > n) || ((n' == n) && (ph' > ph)) +                    then +                      readArray arr (ph, n)+                    else +                      let ps' = ps { parIteration = n', +                                     parPhase = ph',+                                     parTime = basicTime sc n' ph' }+                      in do a <- m ps'+                            a `seq` writeArray arr (ph', n') a+                            if ph' >= phu +                              then do writeIORef phref 0+                                      writeIORef nref (n' + 1)+                                      loop (n' + 1) 0+                              else do writeIORef phref (ph' + 1)+                                      loop n' (ph' + 1)+              n'  <- readIORef nref+              ph' <- readIORef phref+              loop n' ph'+     return $ tr $ Dynamics r++-- | Memoize and order the computation in the integration time points using +-- the specified interpolation and without knowledge of the Runge-Kutta method.+memo0 :: Memo e => (Dynamics e -> Dynamics e) -> Dynamics e +        -> Dynamics (Dynamics e)+memo0 tr (Dynamics m) = +  Dynamics $ \ps ->+  do let sc   = parSpecs ps+         bnds = iterationBnds sc+     arr   <- newMemoArray_ bnds+     nref  <- newIORef 0+     let r ps =+           do let sc = parSpecs ps+                  n  = parIteration ps+                  loop n' = +                    if n' > n+                    then +                      readArray arr n+                    else +                      let ps' = ps { parIteration = n', parPhase = 0,+                                     parTime = basicTime sc n' 0 }+                      in do a <- m ps'+                            a `seq` writeArray arr n' a+                            writeIORef nref (n' + 1)+                            loop (n' + 1)+              n' <- readIORef nref+              loop n'+     return $ tr $ Dynamics r++-- | Memoize and order the computation in the integration time points using +-- the specified interpolation and without knowledge of the Runge-Kutta method.+umemo0 :: UMemo e => (Dynamics e -> Dynamics e) -> Dynamics e +         -> Dynamics (Dynamics e)+umemo0 tr (Dynamics m) = +  Dynamics $ \ps ->+  do let sc   = parSpecs ps+         bnds = iterationBnds sc+     arr   <- newMemoUArray_ bnds+     nref  <- newIORef 0+     let r ps =+           do let sc = parSpecs ps+                  n  = parIteration ps+                  loop n' = +                    if n' > n+                    then +                      readArray arr n+                    else +                      let ps' = ps { parIteration = n', parPhase = 0,+                                     parTime = basicTime sc n' 0 }+                      in do a <- m ps'+                            a `seq` writeArray arr n' a+                            writeIORef nref (n' + 1)+                            loop (n' + 1)+              n' <- readIORef nref+              loop n'+     return $ tr $ Dynamics r++--+-- Once+--++-- | Call the computation only once.+once :: Dynamics a -> Dynamics (Dynamics a)+once (Dynamics m) =+  Dynamics $ \ps ->+  do x <- newIORef Nothing+     let r ps =+           do a <- readIORef x+              case a of+                Just b -> +                  return b+                Nothing ->+                  do b <- m ps+                     writeIORef x $ Just b+                     return $! b+     return $ Dynamics r++--+-- The DynamicsCont Monad+--+-- It looks somewhere like the ContT monad transformer parameterized by +-- the Dynamics monad, although this analogy is not strong. The main +-- idea is to represent the continuation as a dynamic process varying +-- in time.+--++newtype DynamicsCont a = DynamicsCont (Dynamics (a -> IO ()) -> Dynamics ())++instance Monad DynamicsCont where+  return  = returnDC+  m >>= k = bindDC m k++returnDC :: a -> DynamicsCont a+returnDC a = +  DynamicsCont $ \(Dynamics c) -> +  Dynamics $ \ps -> +  do cont' <- c ps+     cont' a+                          +bindDC :: DynamicsCont a -> (a -> DynamicsCont b) -> DynamicsCont b+bindDC (DynamicsCont m) k =+  DynamicsCont $ \c ->+  m $ Dynamics $ \ps -> +  let cont' a = let (DynamicsCont m') = k a+                    (Dynamics u) = m' c+                in u ps+  in return cont'++runCont :: DynamicsCont a -> IO (a -> IO ()) -> Dynamics ()+runCont (DynamicsCont m) f = m $ Dynamics $ const f++instance DynamicsTrans DynamicsCont where+  liftD (Dynamics m) =+    DynamicsCont $ \(Dynamics c) ->+    Dynamics $ \ps ->+    do cont' <- c ps+       a <- m ps+       cont' a++instance Functor DynamicsCont where+  fmap = liftM++instance MonadIO DynamicsCont where+  liftIO m =+    DynamicsCont $ \(Dynamics c) ->+    Dynamics $ \ps ->+    do cont' <- c ps+       a <- m+       cont' a++--+-- The Event Queue (The Dynamics Queue)+--+-- The most exciting thing is that any event is just some value in +-- the Dynamics monad, i.e. a computation, or saying differently, +-- a dynamic process that has a single purpose to perform some +-- side effect. To pass the message, we actually use a closure.+--++-- | The 'DynamicsQueue' type represents the event queue.+data DynamicsQueue = DynamicsQueue { +  queuePQ   :: PQ.PriorityQueue (Dynamics (() -> IO ())),+  queueBusy :: IORef Bool,+  queueTime :: IORef Double, +  queueRun  :: Dynamics () }++-- | Create a new event queue.+newQueue :: Dynamics DynamicsQueue+newQueue = +  Dynamics $ \ps ->+  do let sc = parSpecs ps+     f <- newIORef False+     t <- newIORef $ spcStartTime sc+     let cont () = return ()+     pq <- PQ.newQueue $ return cont+     let q = DynamicsQueue { queuePQ   = pq,+                             queueBusy = f,+                             queueTime = t, +                             queueRun  = subrunQueue q }+     return q+             +-- | Enqueue the event which must be actuated at the specified time.+enqueueDC :: DynamicsQueue -> Dynamics Double -> Dynamics (() -> IO ()) +            -> Dynamics ()+enqueueDC q (Dynamics t) c = Dynamics r where+  r ps =+    do t' <- t ps+       let pq = queuePQ q+       PQ.enqueue pq t' c+    +-- | Enqueue the event which must be actuated at the specified time.+enqueueD :: DynamicsQueue -> Dynamics Double -> Dynamics () -> Dynamics ()+enqueueD q t (Dynamics m) = enqueueDC q t (Dynamics c) where+  c ps = let f () = m ps in return f+    +subrunQueue :: DynamicsQueue -> Dynamics ()+subrunQueue q = Dynamics r where+  r ps =+    do let f = queueBusy q+       f' <- readIORef f+       unless f' $+         do writeIORef f True+            call q ps+            writeIORef f False+  call q ps =+    do let pq = queuePQ q+       f <- PQ.queueNull pq+       unless f $+         do (t2, Dynamics c2) <- PQ.queueFront pq+            let t = queueTime q+            t' <- readIORef t+            when (t2 < t') $ +              error "The time value is too small: subrunQueue"+            when (t2 <= parTime ps) $+              do writeIORef t t2+                 PQ.dequeue pq+                 let sc  = parSpecs ps+                     t0  = spcStartTime sc+                     dt  = spcDT sc+                     n2  = fromInteger $ toInteger $ floor ((t2 - t0) / dt)+                 k <- c2 $ ps { parTime = t2,+                               parIteration = n2,+                               parPhase = -1 }+                 k ()    -- raise the event+                 call q ps++-- | Run the event queue processing its events.+runQueue :: DynamicsQueue -> Dynamics ()+runQueue = queueRun++--+-- DynamicsPID and DynamicsProc+--+-- A value in the DynamicsProc monad represents a control process that can be +-- suspended and resumed at any time. It behaves like a dynamic process too. +-- Any value in the Dynamics monad can be lifted to the DynamicsProc monad. +-- Moreover, a value in the DynamicsProc monad can be run in the Dynamics monad.+--+-- A value of the DynamicsPID type is just an identifier of such a process.+--++-- Public functions:+--+--   pidQueue+--   holdProcD+--   holdProc+--   passivateProc+--   procPassive+--   reactivateProc+--   runProc+--   procPID+--   newPID++-- | Represents a process handler, its PID.+data DynamicsPID = +  DynamicsPID { pidQueue   :: DynamicsQueue,  -- ^ Return the bound event queue.+                pidStarted :: IORef Bool,+                pidCont    :: IORef (Maybe (Dynamics (() -> IO ()))) }++-- | Specifies a discontinuous process that can be suspended at any time+-- and then resumed later.+newtype DynamicsProc a = DynamicsProc (DynamicsPID -> DynamicsCont a)++-- | Hold the process for the specified time period.+holdProcD :: Dynamics Double -> DynamicsProc ()+holdProcD t =+  DynamicsProc $ \pid ->+  DynamicsCont $ \c ->+  enqueueDC (pidQueue pid) (t + time) c++-- | Hold the process for the specified time period.+holdProc :: Double -> DynamicsProc ()+holdProc t = holdProcD $ return t++-- | Passivate the process.+passivateProc :: DynamicsProc ()+passivateProc =+  DynamicsProc $ \pid ->+  DynamicsCont $ \c ->+  Dynamics $ \ps ->+  do let x = pidCont pid+     a <- readIORef x+     case a of+       Nothing -> writeIORef x $ Just c+       Just _  -> error "Cannot passivate the process twice: passivate"++-- | Test whether the process with the specified PID is passivated.+procPassive :: DynamicsPID -> DynamicsProc Bool+procPassive pid =+  DynamicsProc $ \_ ->+  DynamicsCont $ \(Dynamics c) ->+  Dynamics $ \ps ->+  do cont' <- c ps +     let x = pidCont pid+     a <- readIORef x+     case a of+       Nothing -> cont' False+       Just _  -> cont' True++-- | Reactivate a process with the specified PID.+reactivateProc :: DynamicsPID -> DynamicsProc ()+reactivateProc pid =+  DynamicsProc $ \pid' ->+  DynamicsCont $ \c@(Dynamics cont) ->+  Dynamics $ \ps ->+  do let x = pidCont pid+     a <- readIORef x+     case a of+       Nothing ->+         do cont' <- cont ps+            cont' ()+       Just (Dynamics cont2) ->+         do writeIORef x Nothing+            let Dynamics m = enqueueDC (pidQueue pid') time c+            m ps+            cont2' <- cont2 ps+            cont2' ()++-- | Start the process with the specified PID at the desired time.+runProc :: DynamicsProc () -> DynamicsPID -> Dynamics Double -> Dynamics ()+runProc (DynamicsProc p) pid t =+  runCont m r+    where m = do y <- liftIO $ readIORef (pidStarted pid)+                 if y +                   then error $+                        "A process with such PID " +++                        "has been started already: runProc"+                   else liftIO $ writeIORef (pidStarted pid) True+                 DynamicsCont $ \c -> enqueueDC (pidQueue pid) t c+                 p pid+          r = let f () = return () in return f++-- | Return the current process PID.+procPID :: DynamicsProc DynamicsPID+procPID = DynamicsProc $ \pid -> return pid++-- | Create a new process PID.+newPID :: DynamicsQueue -> Dynamics DynamicsPID+newPID q =+  do x <- liftIO $ newIORef Nothing+     y <- liftIO $ newIORef False+     return DynamicsPID { pidQueue   = q,+                          pidStarted = y,+                          pidCont    = x }++instance Eq DynamicsPID where+  x == y = pidCont x == pidCont y    -- for the references are unique++instance Monad DynamicsProc where+  return  = returnDP+  m >>= k = bindDP m k++returnDP :: a -> DynamicsProc a+returnDP a = DynamicsProc (\pid -> return a)++bindDP :: DynamicsProc a -> (a -> DynamicsProc b) -> DynamicsProc b+bindDP (DynamicsProc m) k = +  DynamicsProc $ \pid -> +  do a <- m pid+     let DynamicsProc m' = k a+     m' pid++instance Functor DynamicsProc where+  fmap = liftM++instance DynamicsTrans DynamicsProc where+  liftD m = DynamicsProc $ \pid -> liftD m+  +instance MonadIO DynamicsProc where+  liftIO m = DynamicsProc $ \pid -> liftIO m++--+-- DynamicsResource+--+  +-- Public functions:  +--+--   resourceQueue+--   resourceInitCount+--   resourceCount+--   requestResource+--   releaseResource+--   newResource+  +-- | Represents a limited resource.+data DynamicsResource = +  DynamicsResource { resourceQueue     :: DynamicsQueue,  +                     -- ^ Return the bound event queue.+                     resourceInitCount :: Int,+                     -- ^ Return the initial count of the resource.+                     resourceCountRef  :: IORef Int, +                     resourceWaitQueue :: Q.Queue (Dynamics (() -> IO ()))}++instance Eq DynamicsResource where+  x == y = resourceCountRef x == resourceCountRef y  -- unique references++-- | Create a new resource with the specified initial count.+newResource :: DynamicsQueue -> Int -> Dynamics DynamicsResource+newResource q initCount =+  Dynamics $ \ps ->+  do countRef  <- newIORef initCount+     waitQueue <- Q.newQueue+     return DynamicsResource { resourceQueue     = q,+                               resourceInitCount = initCount,+                               resourceCountRef  = countRef,+                               resourceWaitQueue = waitQueue }++-- | Return the current count of the resource.+resourceCount :: DynamicsResource -> DynamicsProc Int+resourceCount r =+  DynamicsProc $ \_ ->+  DynamicsCont $ \(Dynamics c) ->+  Dynamics $ \ps ->+  do cont' <- c ps +     a <- readIORef (resourceCountRef r)+     cont' a++-- | Request for the resource decreasing its count in case of success,+-- otherwise suspending the discontinuous process until some other +-- process releases the resource.+requestResource :: DynamicsResource -> DynamicsProc ()+requestResource r =+  DynamicsProc $ \_ ->+  DynamicsCont $ \c@(Dynamics cont) ->+  Dynamics $ \ps ->+  do a <- readIORef (resourceCountRef r)+     if a == 0 +       then Q.enqueue (resourceWaitQueue r) c+       else do let a' = a - 1+               a' `seq` writeIORef (resourceCountRef r) a'+               cont' <- cont ps+               cont' ()++-- | Release the resource increasing its count and resuming one of the+-- previously suspended processes as possible.+releaseResource :: DynamicsResource -> DynamicsProc ()+releaseResource r =+  DynamicsProc $ \_ ->+  DynamicsCont $ \(Dynamics c) ->+  Dynamics $ \ps ->+  do a <- readIORef (resourceCountRef r)+     let a' = a + 1+     when (a' > resourceInitCount r) $+       error $+       "The resource count cannot be greater than " +++       "its initial value: releaseResource."+     f <- Q.queueNull (resourceWaitQueue r)+     if f +       then a' `seq` writeIORef (resourceCountRef r) a'+       else do c2 <- Q.queueFront (resourceWaitQueue r)+               Q.dequeue (resourceWaitQueue r)+               let Dynamics m = enqueueDC (resourceQueue r) time c2+               m ps+     cont' <- c ps+     cont' ()++--+-- DynamicsRef+--++-- | The 'DynamicsRef' type represents a mutable variable similar to +-- the 'IORef' variable but only bound to some event queue, which makes +-- the variable coordinated with that queue.+data DynamicsRef a = +  DynamicsRef { refQueue  :: DynamicsQueue,    -- ^ Return the bound event queue.+                refRunner :: Dynamics (),+                refValue  :: IORef a }++-- | Create a new reference bound to the specified event queue.+newRef :: DynamicsQueue -> a -> Dynamics (DynamicsRef a)+newRef q a =+  do x <- liftIO $ newIORef a+     return DynamicsRef { refQueue  = q,+                          refRunner = runQueue q,+                          refValue  = x }+     +-- | Read the value of a reference, forcing the bound event queue to raise +-- the events in case of need.+readRef :: DynamicsRef a -> Dynamics a+readRef r = Dynamics $ \ps -> +  do let Dynamics m = refRunner r+     m ps+     readIORef (refValue r)++-- | Write a new value into the reference.+writeRef :: DynamicsRef a -> a -> Dynamics ()+writeRef r a = Dynamics $ \ps -> +  do writeIORef (refValue r) a+     let Dynamics m = refRunner r +     m ps++-- | Mutate the contents of the reference, forcing the bound event queue to+-- raise all pending events in case of need.+modifyRef :: DynamicsRef a -> (a -> a) -> Dynamics ()+modifyRef r f = Dynamics $ \ps -> +  do let Dynamics m = refRunner r +     m ps+     modifyIORef (refValue r) f++-- | A strict version of the 'writeRef' function.+writeRef' :: DynamicsRef a -> a -> Dynamics ()+writeRef' r a = a `seq` writeRef r a++-- | A strict version of the 'modifyRef' function.+modifyRef' :: DynamicsRef a -> (a -> a) -> Dynamics ()+modifyRef' r f = Dynamics $ \ps ->+  do let Dynamics m = refRunner r+     m ps+     a <- readIORef (refValue r)+     let b = f a+     b `seq` writeIORef (refValue r) b++--+-- Agent-based Modeling+--++-- Public functions:+--+--   agentQueue+--   agentState+--   activateState+--   initState+--   stateAgent     +--   stateParent+--   addTimeoutD+--   addTimeout+--   addTimerD+--   addTimer+--   stateActivation+--   stateDeactivation+--   newState     +--   newSubstate+--   newAgent++-- | Represents an agent.+data Agent = Agent { agentQueue :: DynamicsQueue,+                     -- ^ Return the bound event queue.+                     agentModeRef :: IORef AgentMode,+                     agentStateRef :: IORef (Maybe AgentState) }++-- | Represents the agent state.+data AgentState = AgentState { stateAgent :: Agent,+                               -- ^ Return the corresponded agent.+                               stateParent :: Maybe AgentState,+                               -- ^ Return the parent state or 'Nothing'.+                               stateActivateRef :: IORef (Dynamics ()),+                               stateDeactivateRef :: IORef (Dynamics ()), +                               stateVersionRef :: IORef Int }+                  +data AgentMode = CreationMode+               | InitialMode+               | TransientMode+               | ProcessingMode+                      +instance Eq Agent where+  x == y = agentStateRef x == agentStateRef y      -- unique references+  +instance Eq AgentState where+  x == y = stateVersionRef x == stateVersionRef y  -- unique references++findPath :: AgentState -> AgentState -> ([AgentState], [AgentState])+findPath source target = +  if stateAgent source == stateAgent target +  then+    partitionPath path1 path2+  else+    error "Different agents: findPath."+      where+        path1 = fullPath source []+        path2 = fullPath target []+        fullPath st acc =+          case stateParent st of+            Nothing  -> st : acc+            Just st' -> fullPath st' (st : acc)+        partitionPath path1 path2 =+          case (path1, path2) of+            (h1 : t1, [h2]) | h1 == h2 -> +              (reverse path1, path2)+            (h1 : t1, h2 : t2) | h1 == h2 -> +              partitionPath t1 t2+            _ -> +              (reverse path1, path2)+            +traversePath :: AgentState -> AgentState -> Dynamics ()+traversePath source target =+  let (path1, path2) = findPath source target+      agent = stateAgent source+      activate st ps   =+        do Dynamics m <- readIORef (stateActivateRef st)+           m ps+      deactivate st ps =+        do Dynamics m <- readIORef (stateDeactivateRef st)+           m ps+  in Dynamics $ \ps ->+       do writeIORef (agentModeRef agent) TransientMode+          forM_ path1 $ \st ->+            do writeIORef (agentStateRef agent) (Just st)+               deactivate st ps+               -- it makes all timeout and timer handlers obsolete+               modifyIORef (stateVersionRef st) (1 +)+          forM_ path2 $ \st ->+            do when (st == target) $+                 writeIORef (agentModeRef agent) InitialMode+               writeIORef (agentStateRef agent) (Just st)+               activate st ps+               when (st == target) $+                 writeIORef (agentModeRef agent) ProcessingMode++-- | Add to the state a timeout handler that will be actuated +-- in the specified time period, while the state remains active.+addTimeoutD :: AgentState -> Dynamics Double -> Dynamics () -> Dynamics ()+addTimeoutD st t (Dynamics action) =+  Dynamics $ \ps ->+  do v <- readIORef (stateVersionRef st)+     let m1 = Dynamics $ \ps ->+           do v' <- readIORef (stateVersionRef st)+              when (v == v') $ action ps+         q = agentQueue (stateAgent st)+         Dynamics m2 = enqueueD q (t + time) m1+     m2 ps++-- | Add to the state a timer handler that will be actuated+-- in the specified time period and then repeated again many times,+-- while the state remains active.+addTimerD :: AgentState -> Dynamics Double -> Dynamics () -> Dynamics ()+addTimerD st t (Dynamics action) =+  Dynamics $ \ps ->+  do v <- readIORef (stateVersionRef st)+     let m1 = Dynamics $ \ps ->+           do v' <- readIORef (stateVersionRef st)+              when (v == v') $ do { m2 ps; action ps }+         q = agentQueue (stateAgent st)+         Dynamics m2 = enqueueD q (t + time) m1+     m2 ps++-- | Add to the state a timeout handler that will be actuated +-- in the specified time period, while the state remains active.+addTimeout :: AgentState -> Double -> Dynamics () -> Dynamics ()+addTimeout st t = addTimeoutD st (return t)++-- | Add to the state a timer handler that will be actuated+-- in the specified time period and then repeated again many times,+-- while the state remains active.+addTimer :: AgentState -> Double -> Dynamics () -> Dynamics ()+addTimer st t = addTimerD st (return t)++-- | Create a new state.+newState :: Agent -> Dynamics AgentState+newState agent =+  Dynamics $ \ps ->+  do aref <- newIORef $ return ()+     dref <- newIORef $ return ()+     vref <- newIORef 0+     return AgentState { stateAgent = agent,+                         stateParent = Nothing,+                         stateActivateRef = aref,+                         stateDeactivateRef = dref,+                         stateVersionRef = vref }++-- | Create a child state.+newSubstate :: AgentState -> Dynamics AgentState+newSubstate parent =+  Dynamics $ \ps ->+  do let agent = stateAgent parent +     aref <- newIORef $ return ()+     dref <- newIORef $ return ()+     vref <- newIORef 0+     return AgentState { stateAgent = agent,+                         stateParent = Just parent,+                         stateActivateRef= aref,+                         stateDeactivateRef = dref,+                         stateVersionRef = vref }++-- | Create an agent bound with the specified event queue.+newAgent :: DynamicsQueue -> Dynamics Agent+newAgent queue =+  Dynamics $ \ps ->+  do modeRef    <- newIORef CreationMode+     stateRef   <- newIORef Nothing+     return Agent { agentQueue = queue,+                    agentModeRef = modeRef,+                    agentStateRef = stateRef }++-- | Return the selected downmost active state.+agentState :: Agent -> Dynamics (Maybe AgentState)+agentState agent =+  Dynamics $ \ps -> +  do let Dynamics m = queueRun $ agentQueue agent +     m ps    -- ensure that the agent state is actual+     readIORef (agentStateRef agent)+                   +-- | Select the next downmost active state.       +activateState :: AgentState -> Dynamics ()+activateState st =+  Dynamics $ \ps ->+  do let agent = stateAgent st+         Dynamics m = queueRun $ agentQueue agent +     m ps    -- ensure that the agent state is actual+     mode <- readIORef (agentModeRef agent)+     case mode of+       CreationMode ->+         case stateParent st of+           Just _ ->+             error $ +             "To run the agent for the first time, an initial state " +++             "must be top-level: activateState."+           Nothing ->+             do writeIORef (agentModeRef agent) InitialMode+                writeIORef (agentStateRef agent) (Just st)+                Dynamics m <- readIORef (stateActivateRef st)+                m ps+                writeIORef (agentModeRef agent) ProcessingMode+       InitialMode ->+         error $ +         "Use the initState function during " +++         "the state activation: activateState."+       TransientMode ->+         error $+         "Use the initState function during " +++         "the state activation: activateState."+       ProcessingMode ->+         do Just st0 <- readIORef (agentStateRef agent)+            let Dynamics m = traversePath st0 st+            m ps+              +-- | Activate the child state during the direct activation of +-- the parent state. This call is ignored in other cases.+initState :: AgentState -> Dynamics ()+initState st =+  Dynamics $ \ps ->+  do let agent = stateAgent st+         Dynamics m = queueRun $ agentQueue agent +     m ps    -- ensure that the agent state is actual+     mode <- readIORef (agentModeRef agent)+     case mode of+       CreationMode ->+         error $+         "To run the agent for the fist time, use " +++         "the activateState function: initState."+       InitialMode ->+         do Just st0 <- readIORef (agentStateRef agent)+            let Dynamics m = traversePath st0 st+            m ps+       TransientMode -> +         return ()+       ProcessingMode ->+         error $+         "Use the activateState function everywhere outside " +++         "the state activation: initState."++-- | Set the activation computation for the specified state.+stateActivation :: AgentState -> Dynamics () -> Dynamics ()+stateActivation st action =+  Dynamics $ \ps ->+  writeIORef (stateActivateRef st) action+  +-- | Set the deactivation computation for the specified state.+stateDeactivation :: AgentState -> Dynamics () -> Dynamics ()+stateDeactivation st action =+  Dynamics $ \ps ->+  writeIORef (stateDeactivateRef st) action+  +          

@@ -0,0 +1,173 @@++-- Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>+-- +-- All rights reserved.+-- +-- Redistribution and use in source and binary forms, with or without+-- modification, are permitted provided that the following conditions+-- are met:+-- +-- 1. Redistributions of source code must retain the above copyright+--    notice, this list of conditions and the following disclaimer.+-- +-- 2. Redistributions in binary form must reproduce the above copyright+--    notice, this list of conditions and the following disclaimer in the+--    documentation and/or other materials provided with the distribution.+-- +-- 3. Neither the name of the author nor the names of his contributors+--    may be used to endorse or promote products derived from this software+--    without specific prior written permission.+-- +-- THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+-- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+-- ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+-- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+-- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+-- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+-- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+-- SUCH DAMAGE.++-- This is an imperative version of the heap-based priority queue in monad IO.++module Simulation.Aivika.PriorityQueue +       (PriorityQueue, +        queueNull, +        newQueue, +        enqueue, +        dequeue, +        queueFront) where ++import Data.Array+import Data.Array.MArray+import Data.Array.IO+import Data.IORef+import Control.Monad++data PriorityQueue a = +  PriorityQueue { pqKeys  :: IORef (IOUArray Int Double),+                  pqVals  :: IORef (IOArray Int a),+                  pqSize  :: IORef Int,+                  pqNoVal :: a    -- to release references                 +                }++increase :: PriorityQueue a -> Int -> IO ()+increase pq capacity = +  do let keyRef = pqKeys pq+         valRef = pqVals pq+     keys <- readIORef keyRef+     vals <- readIORef valRef+     (il, iu)  <- getBounds keys+     let len = (iu - il) + 1+         capacity' | len < 64  = max capacity ((len + 1) * 2)+                   | otherwise = max capacity ((len `div` 2) * 3)+         il' = il+         iu' = il + capacity' - 1+     keys' <- newArray_ (il', iu')+     vals' <- newArray_ (il', iu')+     mapM_ (\i -> do { k <- readArray keys i; writeArray keys' i k }) [il..iu]+     mapM_ (\i -> do { v <- readArray vals i; writeArray vals' i v }) [il..iu]+     writeIORef keyRef keys'+     writeIORef valRef vals'++siftUp :: IOUArray Int Double +         -> IOArray Int a+         -> Int -> Double -> a +         -> IO ()+siftUp keys vals i k v =+  if i == 0 +  then do writeArray keys i k+          writeArray vals i v+  else do let n = (i - 1) `div` 2+          kn <- readArray keys n+          if k >= kn +            then do writeArray keys i k+                    writeArray vals i v+            else do vn <- readArray vals n+                    writeArray keys i kn+                    writeArray vals i vn+                    siftUp keys vals n k v++siftDown :: IOUArray Int Double +           -> IOArray Int a -> Int+           -> Int -> Double -> a +           -> IO ()+siftDown keys vals size i k v =+  if i >= (size `div` 2)+  then do writeArray keys i k+          writeArray vals i v+  else do let n  = 2 * i + 1+              n' = n + 1+          kn  <- readArray keys n+          if n' >= size +            then if k <= kn+                 then do writeArray keys i k+                         writeArray vals i v+                 else do vn <- readArray vals n+                         writeArray keys i kn+                         writeArray vals i vn+                         siftDown keys vals size n k v+            else do kn' <- readArray keys n'+                    let n''  = if kn > kn' then n' else n+                        kn'' = min kn' kn+                    if k <= kn''+                      then do writeArray keys i k+                              writeArray vals i v+                      else do vn'' <- readArray vals n''+                              writeArray keys i kn''+                              writeArray vals i vn''+                              siftDown keys vals size n'' k v++queueNull :: PriorityQueue a -> IO Bool+queueNull pq =+  do size <- readIORef (pqSize pq)+     return $ size == 0++newQueue :: a -> IO (PriorityQueue a)+newQueue defaultValue =+  do keys <- newArray_ (0, 10)+     vals <- newArray_ (0, 10)+     keyRef  <- newIORef keys+     valRef  <- newIORef vals+     sizeRef <- newIORef 0+     return PriorityQueue { pqKeys = keyRef, +                            pqVals = valRef, +                            pqSize = sizeRef,+                            pqNoVal = defaultValue }++enqueue :: PriorityQueue a -> Double -> a -> IO ()+enqueue pq k v =+  do i <- readIORef (pqSize pq)+     keys <- readIORef (pqKeys pq)+     (il, iu) <- getBounds keys+     when (i >= iu - il + 1) $ increase pq (i + 1)+     writeIORef (pqSize pq) (i + 1)+     keys <- readIORef (pqKeys pq)  -- it can be another! (side-effect)+     vals <- readIORef (pqVals pq)+     siftUp keys vals i k v++dequeue :: PriorityQueue a -> IO ()+dequeue pq =+  do size <- readIORef (pqSize pq)+     when (size == 0) $ error "Empty priority queue: dequeue"+     let i = size - 1+     writeIORef (pqSize pq) i+     keys <- readIORef (pqKeys pq)+     vals <- readIORef (pqVals pq)+     k <- readArray keys i+     v <- readArray vals i+     writeArray keys i 0.0+     writeArray vals i (pqNoVal pq)    -- to release the reference!+     siftDown keys vals i 0 k v++queueFront :: PriorityQueue a -> IO (Double, a)+queueFront pq =+  do size <- readIORef (pqSize pq)+     when (size == 0) $ error "Empty priority queue: front"+     keys <- readIORef (pqKeys pq)+     vals <- readIORef (pqVals pq)+     k <- readArray keys 0+     v <- readArray vals 0+     return (k, v)

@@ -0,0 +1,112 @@++-- Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>+-- +-- All rights reserved.+-- +-- Redistribution and use in source and binary forms, with or without+-- modification, are permitted provided that the following conditions+-- are met:+-- +-- 1. Redistributions of source code must retain the above copyright+--    notice, this list of conditions and the following disclaimer.+-- +-- 2. Redistributions in binary form must reproduce the above copyright+--    notice, this list of conditions and the following disclaimer in the+--    documentation and/or other materials provided with the distribution.+-- +-- 3. Neither the name of the author nor the names of his contributors+--    may be used to endorse or promote products derived from this software+--    without specific prior written permission.+-- +-- THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+-- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+-- ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+-- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+-- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+-- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+-- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+-- SUCH DAMAGE.++-- This is an imperative version of the queue in monad IO.++module Simulation.Aivika.Queue +       (Queue, +        queueNull, +        newQueue, +        enqueue, +        dequeue, +        queueFront) where ++import Data.IORef+import Control.Monad++data QueueItem a = +  QueueItem { qiVal  :: a,+              qiPrev :: IORef (Maybe (QueueItem a)),+              qiNext :: IORef (Maybe (QueueItem a)) }+  +data Queue a =  +  Queue { qHead :: IORef (Maybe (QueueItem a)),+          qTail :: IORef (Maybe (QueueItem a)) }++queueNull :: Queue a -> IO Bool+queueNull q =+  do head <- readIORef (qHead q) +     case head of+       Nothing -> return True+       Just _  -> return False+    +newQueue :: IO (Queue a)+newQueue =+  do head <- newIORef Nothing +     tail <- newIORef Nothing+     return Queue { qHead = head, qTail = tail }++enqueue :: Queue a -> a -> IO ()+enqueue q v =+  do head <- readIORef (qHead q)+     case head of+       Nothing ->+         do prev <- newIORef Nothing+            next <- newIORef Nothing+            let item = Just QueueItem { qiVal = v, +                                        qiPrev = prev, +                                        qiNext = next }+            writeIORef (qHead q) item+            writeIORef (qTail q) item+       Just h ->+         do prev <- newIORef Nothing+            next <- newIORef head+            let item = Just QueueItem { qiVal = v,+                                        qiPrev = prev,+                                        qiNext = next }+            writeIORef (qiPrev h) item+            writeIORef (qHead q) item++dequeue :: Queue a -> IO ()+dequeue q =+  do tail <- readIORef (qTail q) +     case tail of+       Nothing ->+         error "Empty queue: dequeue"+       Just t ->+         do tail' <- readIORef (qiPrev t)+            case tail' of+              Nothing ->+                do writeIORef (qHead q) Nothing+                   writeIORef (qTail q) Nothing+              Just t' ->+                do writeIORef (qiNext t') Nothing+                   writeIORef (qTail q) tail'++queueFront :: Queue a -> IO a+queueFront q =+  do tail <- readIORef (qTail q)+     case tail of+       Nothing ->+         error "Empty queue: front"+       Just t ->+         return $ qiVal t

@@ -0,0 +1,47 @@+name:            aivika+version:         0.1+synopsis:        A multi-paradigm simulation library+description:+    Aivika is a multi-paradigm simulation library. It allows us to integrate +    a system of ordinary differential equations. Also it can be applied to+    the Discrete Event Simulation. It supports the event-oriented, +    process-oriented and activity-oriented paradigms. Aivika also supports +    the Agent-based Modeling. Finally, it can be applied to System Dynamics.+    .+    The library widely uses monads. The dynamic system is represented as +    a computation in the Dynamics monad. There is also the DynamicsProc+    monad to represent the discontinuous processes which can be suspended+    at any time and then resumed later. Everything else is expressed through +    these two monads, including the events, agent handlers and even integrals.+    .+category:        Simulation+license:         BSD3+license-file:    LICENSE+copyright:       (c) 2009-2011. David Sorokin <david.sorokin@gmail.com>+author:          David Sorokin+maintainer:      David Sorokin <david.sorokin@gmail.com>+cabal-version:   >= 1.2.0+build-type:      Simple+tested-with:     GHC == 6.12.1++extra-source-files:  examples/BassDiffusion.hs+                     examples/ChemicalReaction.hs+                     examples/FishBank.hs+                     examples/MachRep1.hs+                     examples/MachRep1EventDriven.hs+                     examples/MachRep1TimeDriven.hs+                     examples/MachRep2.hs+                     examples/MachRep3.hs++data-files:          doc/aivika.pdf++library+    exposed-modules: Simulation.Aivika.Dynamics+    other-modules:   Simulation.Aivika.Queue+                     Simulation.Aivika.PriorityQueue+                     +    build-depends:   base >= 3 && < 5,+                     mtl >= 1.1.0.2,+                     array >= 0.3.0.0+                     +    ghc-options:     -O2

binary file changed (absent → 275619 bytes)

@@ -0,0 +1,106 @@++import Random+import Data.Array+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics++n = 500    -- the number of agents++advertisingEffectiveness = 0.011+contactRate = 100.0+adoptionFraction = 0.015++specs = Specs { spcStartTime = 0.0, +                spcStopTime = 8.0,+                spcDT = 0.1,+                spcMethod = RungeKutta4 }++exprnd :: Double -> IO Double+exprnd lambda =+  do x <- getStdRandom random+     return (- log x / lambda)+     +boolrnd :: Double -> IO Bool+boolrnd p =+  do x <- getStdRandom random+     return (x <= p)++data Person = Person { personAgent :: Agent,+                       personPotentialAdopter :: AgentState,+                       personAdopter :: AgentState }+              +createPerson :: DynamicsQueue -> Dynamics Person              +createPerson q =    +  do agent <- newAgent q+     potentialAdopter <- newState agent+     adopter <- newState agent+     return Person { personAgent = agent,+                     personPotentialAdopter = potentialAdopter,+                     personAdopter = adopter }+       +createPersons :: DynamicsQueue -> Dynamics (Array Int Person)+createPersons q =+  do list <- forM [1 .. n] $ \i ->+       do p <- createPerson q+          return (i, p)+     return $ array (1, n) list+     +definePerson :: Person -> Array Int Person +               -> DynamicsRef Int -> DynamicsRef Int+               -> Dynamics ()+definePerson p ps potentialAdopters adopters =+  do stateActivation (personPotentialAdopter p) $+       do modifyRef' potentialAdopters $ \a -> a + 1+          -- add a timeout+          t <- liftIO $ exprnd advertisingEffectiveness +          let st  = personPotentialAdopter p+              st' = personAdopter p+          addTimeout st t $ activateState st'+     stateActivation (personAdopter p) $ +       do modifyRef' adopters  $ \a -> a + 1+          -- add a timer that works while the state is active+          let t = liftIO $ exprnd contactRate    -- many times!+          addTimerD (personAdopter p) t $+            do i <- liftIO $ getStdRandom $ randomR (1, n)+               let p' = ps ! i+               st <- agentState (personAgent p')+               when (st == Just (personPotentialAdopter p')) $+                 do b <- liftIO $ boolrnd adoptionFraction+                    when b $ activateState (personAdopter p')+     stateDeactivation (personPotentialAdopter p) $+       modifyRef' potentialAdopters $ \a -> a - 1+     stateDeactivation (personAdopter p) $+       modifyRef' adopters $ \a -> a - 1+        +definePersons :: Array Int Person +                -> DynamicsRef Int +                -> DynamicsRef Int +                -> Dynamics ()+definePersons ps potentialAdopters adopters =+  forM_ (elems ps) $ \p -> +  definePerson p ps potentialAdopters adopters+                               +activatePerson :: Person -> Dynamics ()+activatePerson p = activateState (personPotentialAdopter p)++activatePersons :: Array Int Person -> Dynamics ()+activatePersons ps =+  forM_ (elems ps) $ \p -> activatePerson p++model :: Dynamics (Dynamics [Int])+model =+  do q <- newQueue+     potentialAdopters <- newRef q 0+     adopters <- newRef q 0+     ps <- createPersons q+     definePersons ps potentialAdopters adopters+     activatePersons ps+     return $ do i1 <- readRef potentialAdopters+                 i2 <- readRef adopters+                 return [i1, i2]++main =+  do xs <- runDynamics model specs+     print xs

@@ -0,0 +1,26 @@++import Simulation.Aivika.Dynamics++specs = Specs { spcStartTime = 0, +                spcStopTime = 13, +                spcDT = 0.01,+                spcMethod = RungeKutta4 }++model :: Dynamics (Dynamics [Double])+model =+  do integA <- newInteg 100+     integB <- newInteg 0+     integC <- newInteg 0+     let a = integValue integA+         b = integValue integB+         c = integValue integC+     let ka = 1+         kb = 1+     integDiff integA (- ka * a)+     integDiff integB (ka * a - kb * b)+     integDiff integC (kb * b)+     return $ sequence [a, b, c]++main = +  do a <- runDynamics1 model specs+     print a

@@ -0,0 +1,57 @@++import Data.Array++import Simulation.Aivika.Dynamics++specs = Specs { spcStartTime = 0, +                spcStopTime = 13, +                spcDT = 0.01,+                -- spcDT = 0.000005,+                spcMethod = RungeKutta4 }++model :: Dynamics (Dynamics Double)+model =+  do fishInteg <- newInteg 1000+     shipsInteg <- newInteg 10+     totalProfitInteg <- newInteg 0+     -- integral values --+     let fish = integValue fishInteg+         ships = integValue shipsInteg+         totalProfit = integValue totalProfitInteg+     -- auxiliary values --+     let annualProfit = profit+         area = 100+         carryingCapacity = 1000+         catchPerShip = +           lookupD density $+           listArray (1, 11) [(0.0, -0.048), (1.2, 10.875), (2.4, 17.194), +                              (3.6, 20.548), (4.8, 22.086), (6.0, 23.344), +                              (7.2, 23.903), (8.4, 24.462), (9.6, 24.882), +                              (10.8, 25.301), (12.0, 25.86)]+         deathFraction = +           lookupD (fish / carryingCapacity) $+           listArray (1, 11) [(0.0, 5.161), (0.1, 5.161), (0.2, 5.161), +                              (0.3, 5.161), (0.4, 5.161), (0.5, 5.161), +                              (0.6, 5.118), (0.7, 5.247), (0.8, 5.849), +                              (0.9, 6.151), (10.0, 6.194)]+         density = fish / area+         fishDeathRate = maxD 0 (fish * deathFraction)+         fishHatchRate = maxD 0 (fish * hatchFraction)+         fishPrice = 20+         fractionInvested = 0.2+         hatchFraction = 6+         operatingCost = ships * 250+         profit = revenue - operatingCost+         revenue = totalCatchPerYear * fishPrice+         shipBuildingRate = maxD 0 (profit * fractionInvested / shipCost)+         shipCost = 300+         totalCatchPerYear = maxD 0 (ships * catchPerShip)+     -- derivatives --+     integDiff fishInteg (fishHatchRate - fishDeathRate - totalCatchPerYear)+     integDiff shipsInteg shipBuildingRate+     integDiff totalProfitInteg annualProfit+     -- results --+     return annualProfit++main = do a <- runDynamics1 model specs+          print a

@@ -0,0 +1,69 @@++-- It corresponds to model MachRep1 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/].+--   +-- The model description is as follows.+--+-- Two machines, which sometimes break down.+-- Up time is exponentially distributed with mean 1.0, and repair time is+-- exponentially distributed with mean 0.5. There are two repairpersons,+-- so the two machines can be repaired simultaneously if they are down+-- at the same time.+--+-- Output is long-run proportion of up time. Should get value of about+-- 0.66.++import Random+import Control.Monad.Trans++import Simulation.Aivika.Dynamics++upRate = 1.0 / 1.0       -- reciprocal of mean up time+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time++specs = Specs { spcStartTime = 0.0,+                spcStopTime = 1000.0,+                spcDT = 1.0,+                spcMethod = RungeKutta4 }+        +exprnd :: Double -> IO Double+exprnd lambda =+  do x <- getStdRandom random+     return (- log x / lambda)+     +model :: Dynamics (Dynamics Double)+model =+  do queue <- newQueue+     totalUpTime <- newRef queue 0.0+     +     pid1 <- newPID queue+     pid2 <- newPID queue+     +     let machine :: DynamicsProc ()+         machine =+           do startUpTime <- liftD time+              upTime <- liftIO $ exprnd upRate+              holdProc upTime+              finishUpTime <- liftD time+              liftD $ modifyRef' totalUpTime+                (+ (finishUpTime - startUpTime))+              repairTime <- liftIO $ exprnd repairRate+              holdProc repairTime+              machine+         +     runProc machine pid1 starttime+     runProc machine pid2 starttime+     +     let system :: Dynamics Double+         system =+           do x <- readRef totalUpTime+              y <- stoptime+              return $ x / (2 * y)+     +     return system+  +main =         +  do a <- runDynamics1 model specs+     print a

@@ -0,0 +1,75 @@++-- It corresponds to model MachRep1 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/].+--   +-- The model description is as follows.+--+-- Two machines, which sometimes break down.+-- Up time is exponentially distributed with mean 1.0, and repair time is+-- exponentially distributed with mean 0.5. There are two repairpersons,+-- so the two machines can be repaired simultaneously if they are down+-- at the same time.+--+-- Output is long-run proportion of up time. Should get value of about+-- 0.66.++import Random+import Control.Monad.Trans++import Simulation.Aivika.Dynamics++upRate = 1.0 / 1.0       -- reciprocal of mean up time+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time++specs = Specs { spcStartTime = 0.0,+                spcStopTime = 1000.0,+                spcDT = 1.0,+                spcMethod = RungeKutta4 }+        +exprnd :: Double -> IO Double+exprnd lambda =+  do x <- getStdRandom random+     return (- log x / lambda)+     +model :: Dynamics (Dynamics Double)+model =+  do queue <- newQueue+     totalUpTime <- newRef queue 0.0+     +     let machineBroken :: Double -> Dynamics ()+         machineBroken startUpTime =+           +           do finishUpTime <- time+              modifyRef' totalUpTime (+ (finishUpTime - startUpTime))+              repairTime <- liftIO $ exprnd repairRate+              +              -- enqueue a new event+              let t = return $ finishUpTime + repairTime+              enqueueD queue t machineRepaired+              +         machineRepaired :: Dynamics ()+         machineRepaired =+           +           do startUpTime <- time+              upTime <- liftIO $ exprnd upRate+              +              -- enqueue a new event+              let t = return $ startUpTime + upTime+              enqueueD queue t $ machineBroken startUpTime+     +     enqueueD queue starttime machineRepaired    -- start the first machine+     enqueueD queue starttime machineRepaired    -- start the second machine+     +     let system :: Dynamics Double+         system =+           do x <- readRef totalUpTime+              y <- stoptime+              return $ x / (2 * y)+              +     return system+  +main =         +  do a <- runDynamics1 model specs+     print a

@@ -0,0 +1,116 @@++-- It corresponds to model MachRep1 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/].+--   +-- The model description is as follows.+--+-- Two machines, which sometimes break down.+-- Up time is exponentially distributed with mean 1.0, and repair time is+-- exponentially distributed with mean 0.5. There are two repairpersons,+-- so the two machines can be repaired simultaneously if they are down+-- at the same time.+--+-- Output is long-run proportion of up time. Should get value of about+-- 0.66.++import Random+import Control.Monad.Trans++import Simulation.Aivika.Dynamics++upRate = 1.0 / 1.0       -- reciprocal of mean up time+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time++specs = Specs { spcStartTime = 0.0,+                spcStopTime = 1000.0,+                spcDT = 0.05,+                spcMethod = RungeKutta4 }+        +exprnd :: Double -> IO Double+exprnd lambda =+  do x <- getStdRandom random+     return (- log x / lambda)+     +model :: Dynamics (Dynamics Double)+model =+  do queue <- newQueue+     totalUpTime <- newRef queue 0.0+     +     let machine :: Dynamics (Dynamics ())+         machine =+           do startUpTime <- newRef queue 0.0 +             +              -- a number of iterations when +              -- the machine works+              upNum <- newRef queue (-1)+              +              -- a number of iterations when +              -- the machine is broken+              repairNum <- newRef queue (-1)+              +              -- create a simulation model+              return $+                do upNum' <- readRef upNum+                   repairNum' <- readRef repairNum+                   +                   let untilBroken = +                         modifyRef' upNum $ \a -> a - 1+                                                  +                       untilRepaired =+                         modifyRef' repairNum $ \a -> a - 1+                                                      +                       broken =+                         do writeRef' upNum (-1)+                            -- the machine is broken+                            startUpTime' <- readRef startUpTime+                            finishUpTime' <- time+                            dt' <- dt+                            modifyRef' totalUpTime $ +                              \a -> a ++                              (finishUpTime' - startUpTime')+                            repairTime' <- +                              liftIO $ exprnd repairRate+                            writeRef' repairNum $+                              round (repairTime' / dt')+                              +                       repaired =+                         do writeRef' repairNum (-1)+                            -- the machine is repaired+                            t'  <- time+                            dt' <- dt+                            writeRef' startUpTime t'+                            upTime' <- +                              liftIO $ exprnd upRate+                            writeRef' upNum $+                              round (upTime' / dt')+                              +                       result | upNum' > 0     = untilBroken+                              | upNum' == 0     = broken+                              | repairNum' > 0 = untilRepaired+                              | repairNum' == 0 = repaired+                              | otherwise      = repaired +                   result+                            +     -- create two machines with type Dynamics ()+     m1 <- machine+     m2 <- machine+     +     -- create strictly sequential computations+     c1 <- memo0 discrete m1+     c2 <- memo0 discrete m2+       +     let system :: Dynamics Double+         system =+           do c1    -- involve in the simulation+              c2    -- involve in the simulation+              x <- readRef totalUpTime+              y <- stoptime+              return $ x / (2 * y)+     +     return system+  +main =         +  do a <- runDynamics1 model specs+     print a

@@ -0,0 +1,96 @@++-- It corresponds to model MachRep2 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/].+--   +-- The model description is as follows.+--   +-- Two machines, but sometimes break down. Up time is exponentially +-- distributed with mean 1.0, and repair time is exponentially distributed +-- with mean 0.5. In this example, there is only one repairperson, so +-- the two machines cannot be repaired simultaneously if they are down +-- at the same time.+--+-- In addition to finding the long-run proportion of up time as in+-- model MachRep1, let’s also find the long-run proportion of the time +-- that a given machine does not have immediate access to the repairperson +-- when the machine breaks down. Output values should be about 0.6 and 0.67. ++import Random+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics++upRate = 1.0 / 1.0       -- reciprocal of mean up time+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time++specs = Specs { spcStartTime = 0.0,+                spcStopTime = 1000.0,+                spcDT = 1.0,+                spcMethod = RungeKutta4 }+        +exprnd :: Double -> IO Double+exprnd lambda =+  do x <- getStdRandom random+     return (- log x / lambda)+     +model :: Dynamics (Dynamics (Double, Double))+model =+  do queue <- newQueue+     +     -- number of times the machines have broken down+     nRep <- newRef queue 0 +     +     -- number of breakdowns in which the machine +     -- started repair service right away+     nImmedRep <- newRef queue 0+     +     -- total up time for all machines+     totalUpTime <- newRef queue 0.0+     +     repairPerson <- newResource queue 1+     +     pid1 <- newPID queue+     pid2 <- newPID queue+     +     let machine :: DynamicsProc ()+         machine =+           do startUpTime <- liftD time+              upTime <- liftIO $ exprnd upRate+              holdProc upTime+              finishUpTime <- liftD time+              liftD $ modifyRef' totalUpTime +                (+ (finishUpTime - startUpTime))+              +              -- check the resource availability+              liftD $ modifyRef' nRep (+ 1)+              n <- resourceCount repairPerson+              when (n == 1) $+                liftD $ modifyRef' nImmedRep (+ 1)+                +              requestResource repairPerson+              repairTime <- liftIO $ exprnd repairRate+              holdProc repairTime+              releaseResource repairPerson+              +              machine+         +     runProc machine pid1 starttime+     runProc machine pid2 starttime+     +     let system :: Dynamics (Double, Double)+         system =+           do x <- readRef totalUpTime+              y <- stoptime+              n <- readRef nRep+              nImmed <- readRef nImmedRep+              return (x / (2 * y), +                      fromIntegral nImmed / fromIntegral n)+     +     return system+  +main =         +  do a <- runDynamics1 model specs+     print a

@@ -0,0 +1,87 @@++-- It corresponds to model MachRep3 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/].+--   +-- The model description is as follows.+--+-- Variation of models MachRep1, MachRep2. Two machines, but+-- sometimes break down. Up time is exponentially distributed with mean+-- 1.0, and repair time is exponentially distributed with mean 0.5. In+-- this example, there is only one repairperson, and she is not summoned+-- until both machines are down. We find the proportion of up time. It+-- should come out to about 0.45.++import Random+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics++upRate = 1.0 / 1.0       -- reciprocal of mean up time+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time++specs = Specs { spcStartTime = 0.0,+                spcStopTime = 1000.0,+                spcDT = 1.0,+                spcMethod = RungeKutta4 }+        +exprnd :: Double -> IO Double+exprnd lambda =+  do x <- getStdRandom random+     return (- log x / lambda)+     +model :: Dynamics (Dynamics Double)+model =+  do queue <- newQueue+     +     -- number of machines currently up+     nUp <- newRef queue 2+     +     -- total up time for all machines+     totalUpTime <- newRef queue 0.0+     +     repairPerson <- newResource queue 1+     +     pid1 <- newPID queue+     pid2 <- newPID queue+     +     let machine :: DynamicsPID -> DynamicsProc ()+         machine pid =+           do startUpTime <- liftD time+              upTime <- liftIO $ exprnd upRate+              holdProc upTime+              finishUpTime <- liftD time+              liftD $ modifyRef' totalUpTime +                (+ (finishUpTime - startUpTime))+                +              liftD $ modifyRef' nUp $ \a -> a - 1+              nUp' <- liftD $ readRef nUp+              if nUp' == 1+                then passivateProc+                else do n <- resourceCount repairPerson+                        when (n == 1) $ reactivateProc pid+              +              requestResource repairPerson+              repairTime <- liftIO $ exprnd repairRate+              holdProc repairTime+              liftD $ modifyRef' nUp $ \a -> a + 1+              releaseResource repairPerson+              +              machine pid++     runProc (machine pid2) pid1 starttime+     runProc (machine pid1) pid2 starttime+     +     let system :: Dynamics Double+         system =+           do x <- readRef totalUpTime+              y <- stoptime+              return $ x / (2 * y)+     +     return system+  +main =         +  do a <- runDynamics1 model specs+     print a