aivika 0.1 → 0.2
raw patch · 37 files changed
+3500/−1786 lines, 37 filesdep +containersdep +haskell98dep ~base
Dependencies added: containers, haskell98
Dependency ranges changed: base
Files
- Simulation/Aivika/Dynamics.hs +34/−1627
- Simulation/Aivika/Dynamics/Agent.hs +260/−0
- Simulation/Aivika/Dynamics/Base.hs +40/−0
- Simulation/Aivika/Dynamics/Cont.hs +20/−0
- Simulation/Aivika/Dynamics/EventQueue.hs +89/−0
- Simulation/Aivika/Dynamics/Internal/Cont.hs +84/−0
- Simulation/Aivika/Dynamics/Internal/Dynamics.hs +389/−0
- Simulation/Aivika/Dynamics/Internal/Fold.hs +91/−0
- Simulation/Aivika/Dynamics/Internal/Interpolate.hs +62/−0
- Simulation/Aivika/Dynamics/Internal/Memo.hs +195/−0
- Simulation/Aivika/Dynamics/Internal/Process.hs +162/−0
- Simulation/Aivika/Dynamics/Internal/Time.hs +35/−0
- Simulation/Aivika/Dynamics/Lift.hs +21/−0
- Simulation/Aivika/Dynamics/Parameter.hs +58/−0
- Simulation/Aivika/Dynamics/Process.hs +30/−0
- Simulation/Aivika/Dynamics/Random.hs +73/−0
- Simulation/Aivika/Dynamics/Ref.hs +63/−0
- Simulation/Aivika/Dynamics/Resource.hs +101/−0
- Simulation/Aivika/Dynamics/SystemDynamics.hs +441/−0
- Simulation/Aivika/Dynamics/UVar.hs +126/−0
- Simulation/Aivika/Dynamics/Var.hs +126/−0
- Simulation/Aivika/PriorityQueue.hs +33/−45
- Simulation/Aivika/Queue.hs +32/−37
- Simulation/Aivika/Statistics.hs +118/−0
- Simulation/Aivika/UVector.hs +130/−0
- Simulation/Aivika/Vector.hs +127/−0
- aivika.cabal +50/−14
- doc/aivika.pdf binary
- examples/BassDiffusion.hs +13/−16
- examples/ChemicalReaction.hs +1/−0
- examples/FishBank.hs +1/−0
- examples/Furnace.hs +417/−0
- examples/MachRep1.hs +15/−8
- examples/MachRep1EventDriven.hs +12/−7
- examples/MachRep1TimeDriven.hs +13/−10
- examples/MachRep2.hs +18/−10
- examples/MachRep3.hs +20/−12
Simulation/Aivika/Dynamics.hs view
@@ -1,1628 +1,35 @@ --- 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- - +-- |+-- 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 7.0.3+--+-- The module defines the 'Dynamics' monad representing an abstract dynamic +-- process, i.e. a time varying polymorphic function. +-- +-- This is a key point of the Aivika simulation library. With help of this monad +-- we can simulate the system of ordinary differential equations (ODEs) of +-- System Dynamics, define the tasks of Discrete Event Simulation (DES) supporting +-- different paradigms. Also we can use the Agent-based Modeling. Thus, +-- we can create hybrid simulation models.+--+module Simulation.Aivika.Dynamics + (Dynamics,+ Specs(..),+ Method(..),+ runDynamics1,+ runDynamics1_,+ runDynamics,+ runDynamics_,+ runDynamicsIO,+ runDynamicsSeries1,+ runDynamicsSeries1_,+ runDynamicsSeries,+ runDynamicsSeries_,+ printDynamics1,+ printDynamics) where++import Simulation.Aivika.Dynamics.Internal.Dynamics
+ Simulation/Aivika/Dynamics/Agent.hs view
@@ -0,0 +1,260 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Agent+-- 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 7.0.3+--+-- This module introduces an agent-based modeling.+--++module Simulation.Aivika.Dynamics.Agent+ (Agent,+ AgentState,+ newAgent,+ newState,+ newSubstate,+ agentQueue,+ agentState,+ activateState,+ initState,+ stateAgent,+ stateParent,+ addTimeout,+ addTimer,+ stateActivation,+ stateDeactivation) where++import Data.IORef+import Control.Monad++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.EventQueue++--+-- Agent-based Modeling+--++-- | Represents an agent.+data Agent = Agent { agentQueue :: EventQueue,+ -- ^ 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 p =+ do Dynamics m <- readIORef (stateActivateRef st)+ m p+ deactivate st p =+ do Dynamics m <- readIORef (stateDeactivateRef st)+ m p+ in Dynamics $ \p ->+ do writeIORef (agentModeRef agent) TransientMode+ forM_ path1 $ \st ->+ do writeIORef (agentStateRef agent) (Just st)+ deactivate st p+ -- 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 p+ 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.+addTimeout :: AgentState -> Double -> Dynamics () -> Dynamics ()+addTimeout st dt (Dynamics action) =+ Dynamics $ \p ->+ do v <- readIORef (stateVersionRef st)+ let m1 = Dynamics $ \p ->+ do v' <- readIORef (stateVersionRef st)+ when (v == v') $ action p+ q = agentQueue (stateAgent st)+ Dynamics m2 = enqueue q (pointTime p + dt) m1+ m2 p++-- | 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 -> Dynamics Double -> Dynamics () -> Dynamics ()+addTimer st (Dynamics dt) (Dynamics action) =+ Dynamics $ \p ->+ do v <- readIORef (stateVersionRef st)+ let m1 = Dynamics $ \p ->+ do v' <- readIORef (stateVersionRef st)+ when (v == v') $ do { m2 p; action p }+ q = agentQueue (stateAgent st)+ Dynamics m2 = + Dynamics $ \p ->+ do dt' <- dt p+ let Dynamics m3 = enqueue q (pointTime p + dt') m1+ m3 p+ m2 p++-- | Create a new state.+newState :: Agent -> Dynamics AgentState+newState agent =+ Dynamics $ \p ->+ 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 $ \p ->+ 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 :: EventQueue -> Dynamics Agent+newAgent queue =+ Dynamics $ \p ->+ 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 $ \p -> + do let Dynamics m = queueRun $ agentQueue agent + m p -- ensure that the agent state is actual+ readIORef (agentStateRef agent)+ +-- | Select the next downmost active state. +activateState :: AgentState -> Dynamics ()+activateState st =+ Dynamics $ \p ->+ do let agent = stateAgent st+ Dynamics m = queueRun $ agentQueue agent + m p -- 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 p+ 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 p+ +-- | 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 $ \p ->+ do let agent = stateAgent st+ Dynamics m = queueRun $ agentQueue agent + m p -- 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 p+ 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 $ \p ->+ writeIORef (stateActivateRef st) action+ +-- | Set the deactivation computation for the specified state.+stateDeactivation :: AgentState -> Dynamics () -> Dynamics ()+stateDeactivation st action =+ Dynamics $ \p ->+ writeIORef (stateDeactivateRef st) action+
+ Simulation/Aivika/Dynamics/Base.hs view
@@ -0,0 +1,40 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Base+-- 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 7.0.3+--+-- This module defines basic functions for the 'Dynamics' monad.+--++module Simulation.Aivika.Dynamics.Base+ (-- * Time Parameters+ starttime,+ stoptime,+ dt,+ time,+ -- * Interpolation and Initial Value+ initD,+ discrete,+ interpolate,+ -- * Memoization+ memo,+ umemo,+ memo0,+ umemo0,+ -- * Iterating+ iterateD,+ -- * Fold+ foldD1,+ foldD,+ -- * Norming+ divideD) where++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Time+import Simulation.Aivika.Dynamics.Internal.Interpolate+import Simulation.Aivika.Dynamics.Internal.Memo+import Simulation.Aivika.Dynamics.Internal.Fold
+ Simulation/Aivika/Dynamics/Cont.hs view
@@ -0,0 +1,20 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Cont+-- 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 7.0.3+--+-- The 'Cont' monad looks somewhere like the standard 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.+--+module Simulation.Aivika.Dynamics.Cont+ (Cont,+ runCont) where++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Cont
+ Simulation/Aivika/Dynamics/EventQueue.hs view
@@ -0,0 +1,89 @@++-- |+-- Module : Simulation.Aivika.Dynamics.EventQueue+-- 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 7.0.3+--+-- The module introduces the event queue. Any event is the Dynamics computation,+-- or, saying differently, a dynamic process that has a single purpose +-- to perform some side effect at the desired time. To pass the message, +-- we actually use a closure.+--+module Simulation.Aivika.Dynamics.EventQueue+ (EventQueue,+ newQueue,+ enqueueCont,+ enqueue,+ queueRun) where++import Data.IORef+import Control.Monad++import Simulation.Aivika.Dynamics.Internal.Dynamics+import qualified Simulation.Aivika.PriorityQueue as PQ++-- | The 'EventQueue' type represents the event queue.+data EventQueue = EventQueue { + queuePQ :: PQ.PriorityQueue (Dynamics (() -> IO ())),+ queueRun :: Dynamics (), -- ^ Run the event queue processing its events+ queueBusy :: IORef Bool,+ queueTime :: IORef Double }++-- | Create a new event queue.+newQueue :: Dynamics EventQueue+newQueue = + Dynamics $ \p ->+ do let sc = pointSpecs p+ f <- newIORef False+ t <- newIORef $ spcStartTime sc+ pq <- PQ.newQueue+ let q = EventQueue { queuePQ = pq,+ queueRun = runQueue q,+ queueBusy = f,+ queueTime = t }+ return q+ +-- | Enqueue the event which must be actuated at the specified time.+enqueueCont :: EventQueue -> Double -> Dynamics (() -> IO ()) -> Dynamics ()+enqueueCont q t c = Dynamics r where+ r p = let pq = queuePQ q in PQ.enqueue pq t c+ +-- | Enqueue the event which must be actuated at the specified time.+enqueue :: EventQueue -> Double -> Dynamics () -> Dynamics ()+enqueue q t (Dynamics m) = enqueueCont q t (Dynamics c) where+ c p = let f () = m p in return f+ +-- | Run the event queue processing its events.+runQueue :: EventQueue -> Dynamics ()+runQueue q = Dynamics r where+ r p =+ do let f = queueBusy q+ f' <- readIORef f+ unless f' $+ do writeIORef f True+ call q p+ writeIORef f False+ call q p =+ 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 <= pointTime p) $+ do writeIORef t t2+ PQ.dequeue pq+ let sc = pointSpecs p+ t0 = spcStartTime sc+ dt = spcDT sc+ n2 = fromInteger $ toInteger $ floor ((t2 - t0) / dt)+ k <- c2 $ p { pointTime = t2,+ pointIteration = n2,+ pointPhase = -1 }+ k () -- raise the event+ call q p
+ Simulation/Aivika/Dynamics/Internal/Cont.hs view
@@ -0,0 +1,84 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.Cont+-- 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 7.0.3+--+-- The 'Cont' monad looks somewhere like the standard 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.+--+module Simulation.Aivika.Dynamics.Internal.Cont+ (Cont(..),+ runCont) where++import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Lift++-- | The 'Cont' type is similar to the standard Cont monad but only+-- the continuation is represented as a dynamic process varying in time.+newtype Cont a = Cont (Dynamics (a -> IO ()) -> Dynamics ())++instance Monad Cont where+ return = returnC+ m >>= k = bindC m k++instance Lift Cont where+ liftD = liftC++instance Functor Cont where+ fmap = liftM++instance MonadIO Cont where+ liftIO = liftIOC ++returnC :: a -> Cont a+{-# INLINE returnC #-}+returnC a = + Cont $ \(Dynamics c) -> + Dynamics $ \p -> + do cont' <- c p+ cont' a+ +bindC :: Cont a -> (a -> Cont b) -> Cont b+{-# INLINE bindC #-}+bindC (Cont m) k =+ Cont $ \c ->+ m $ Dynamics $ \p -> + let cont' a = let (Cont m') = k a+ (Dynamics u) = m' c+ in u p+ in return cont'++-- | Run the 'Cont' computation.+runCont :: Cont a -> IO (a -> IO ()) -> Dynamics ()+{-# INLINE runCont #-}+runCont (Cont m) f = m $ Dynamics $ const f++-- | Lift the 'Dynamics' computation.+liftC :: Dynamics a -> Cont a+{-# INLINE liftC #-}+liftC (Dynamics m) =+ Cont $ \(Dynamics c) ->+ Dynamics $ \p ->+ do cont' <- c p+ a <- m p+ cont' a+ +-- | Lift the IO computation.+liftIOC :: IO a -> Cont a+{-# INLINE liftIOC #-}+liftIOC m =+ Cont $ \(Dynamics c) ->+ Dynamics $ \p ->+ do cont' <- c p+ a <- m+ cont' a+
+ Simulation/Aivika/Dynamics/Internal/Dynamics.hs view
@@ -0,0 +1,389 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.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 7.0.3+--+-- The module defines the 'Dynamics' monad representing an abstract dynamic +-- process, i.e. a time varying polymorphic function. +-- +-- This is a key point of the Aivika simulation library. With help of this monad +-- we can simulate the system of ordinary differential equations (ODEs) of +-- System Dynamics, define the tasks of Discrete Event Simulation (DES) supporting +-- different paradigms. Also we can use the Agent-based Modeling. Thus, +-- we can create hybrid simulation models.+--+module Simulation.Aivika.Dynamics.Internal.Dynamics+ (-- * Dynamics+ Dynamics(..),+ Point(..),+ Specs(..),+ Method(..),+ Run(..),+ runDynamics1,+ runDynamics1_,+ runDynamics,+ runDynamics_,+ runDynamicsIO,+ runDynamicsSeries1,+ runDynamicsSeries1_,+ runDynamicsSeries,+ runDynamicsSeries_,+ printDynamics1,+ printDynamics,+ -- * Utilities+ basicTime,+ iterationBnds,+ iterationHiBnd,+ iterationLoBnd,+ phaseBnds,+ phaseHiBnd,+ phaseLoBnd) where++import Control.Monad+import Control.Monad.Trans++--+-- 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 (Point -> IO a)++-- | It defines the simulation point appended with the additional information.+data Point = Point { pointSpecs :: Specs, -- ^ the simulation specs+ pointRun :: Run, -- ^ the simulation run+ pointTime :: Double, -- ^ the current time+ pointIteration :: Int, -- ^ the current iteration+ pointPhase :: Int -- ^ the current phase+ } deriving (Eq, Ord, Show)++-- | 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)++-- | It defined the simulation run as part of some experiment.+data Run = Run { runIndex :: Int, -- ^ the current simulation run+ runCount :: Int -- ^ the total number of runs in this experiment+ } deriving (Eq, Ord, Show)+ +-- | Returns the iterations starting from zero.+iterations :: Specs -> [Int]+iterations sc = [i1 .. i2] where+ i1 = 0+ i2 = round ((spcStopTime sc - + spcStartTime sc) / spcDT sc)++-- | Returns the first and last iterations.+iterationBnds :: Specs -> (Int, Int)+iterationBnds sc = (0, round ((spcStopTime sc - + spcStartTime sc) / spcDT sc))++-- | Returns the first iteration, i.e. zero.+iterationLoBnd :: Specs -> Int+iterationLoBnd sc = 0++-- | Returns the last iteration.+iterationHiBnd :: Specs -> Int+iterationHiBnd sc = round ((spcStopTime sc - + spcStartTime sc) / spcDT sc)++-- | Returns the phases for the specified simulation specs starting from zero.+phases :: Specs -> [Int]+phases sc = + case spcMethod sc of+ Euler -> [0]+ RungeKutta2 -> [0, 1]+ RungeKutta4 -> [0, 1, 2, 3]++-- | Returns the first and last phases.+phaseBnds :: Specs -> (Int, Int)+phaseBnds sc = + case spcMethod sc of+ Euler -> (0, 0)+ RungeKutta2 -> (0, 1)+ RungeKutta4 -> (0, 3)++-- | Returns the first phase, i.e. zero.+phaseLoBnd :: Specs -> Int+phaseLoBnd sc = 0+ +-- | Returns the last phase, 1 for Euler's method, 2 for RK2 and 4 for RK4.+phaseHiBnd :: Specs -> Int+phaseHiBnd sc = + case spcMethod sc of+ Euler -> 0+ RungeKutta2 -> 1+ RungeKutta4 -> 3++-- | Returns a simulation time for the integration point specified by +-- the specs, iteration and phase.+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++instance Monad Dynamics where+ return = returnD+ m >>= k = bindD m k++returnD :: a -> Dynamics a+returnD a = Dynamics (\p -> return a)++bindD :: Dynamics a -> (a -> Dynamics b) -> Dynamics b+bindD (Dynamics m) k = + Dynamics $ \p -> + do a <- m p+ let Dynamics m' = k a+ m' p++subrunDynamics1 :: Dynamics a -> Specs -> Run -> IO a+subrunDynamics1 (Dynamics m) sc r =+ do let n = iterationHiBnd sc+ t = basicTime sc n 0+ m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = t,+ pointIteration = n,+ pointPhase = 0 }++subrunDynamics1_ :: Dynamics a -> Specs -> Run -> IO ()+subrunDynamics1_ (Dynamics m) sc r =+ do let n = iterationHiBnd sc+ t = basicTime sc n 0+ m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = t,+ pointIteration = n,+ pointPhase = 0 }+ return ()++subrunDynamics :: Dynamics a -> Specs -> Run -> [IO a]+subrunDynamics (Dynamics m) sc r =+ do let (nl, nu) = iterationBnds sc+ point n = Point { pointSpecs = sc,+ pointRun = r,+ pointTime = basicTime sc n 0,+ pointIteration = n,+ pointPhase = 0 }+ map (m . point) [nl .. nu]++subrunDynamics_ :: Dynamics a -> Specs -> Run -> IO ()+subrunDynamics_ (Dynamics m) sc r =+ do let (nl, nu) = iterationBnds sc+ point n = Point { pointSpecs = sc,+ pointRun = r,+ pointTime = basicTime sc n 0,+ pointIteration = n,+ pointPhase = 0 }+ mapM_ (m . point) [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 let r = Run { runIndex = 1, runCount = 1 }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ subrunDynamics1 d sc r++-- | Run the simulation and return the result in the last +-- time point using the specified simulation specs.+runDynamics1_ :: Dynamics (Dynamics a) -> Specs -> IO ()+runDynamics1_ (Dynamics m) sc = + do let r = Run { runIndex = 1, runCount = 1 }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ subrunDynamics1_ d sc r++-- | 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 let r = Run { runIndex = 1, runCount = 1 }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ sequence $ subrunDynamics d sc r++-- | Run the simulation and return the results in all +-- integration time points using the specified simulation specs.+runDynamics_ :: Dynamics (Dynamics a) -> Specs -> IO ()+runDynamics_ (Dynamics m) sc = + do let r = Run { runIndex = 1, runCount = 1 }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ sequence_ $ subrunDynamics d sc r++-- | 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 let r = Run { runIndex = 1, runCount = 1 }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ return $ subrunDynamics d sc r++-- | Run an experiment consisting of the given number of simulations, where each +-- model is created and then requested in the last integration time point using +-- the specified specs.+runDynamicsSeries1_ :: Dynamics (Dynamics a) -> Specs -> Int -> [IO ()]+runDynamicsSeries1_ (Dynamics m) sc runs = map f [1 .. runs]+ where f i =+ do let r = Run { runIndex = i, runCount = runs }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ subrunDynamics1_ d sc r++-- | Run an experiment consisting of the given number of simulations, where each +-- model is created and then requested sequentially in all integration time points +-- using the specified specs.+runDynamicsSeries_ :: Dynamics (Dynamics a) -> Specs -> Int -> [IO ()]+runDynamicsSeries_ (Dynamics m) sc runs = map f [1 .. runs]+ where f i =+ do let r = Run { runIndex = i, runCount = runs }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ subrunDynamics_ d sc r++-- | Run an experiment consisting of the given number of simulations, where each +-- model is created and then requested in the last integration time point using +-- the specified specs.+runDynamicsSeries1 :: Dynamics (Dynamics a) -> Specs -> Int -> [IO a]+runDynamicsSeries1 (Dynamics m) sc runs = map f [1 .. runs]+ where f i =+ do let r = Run { runIndex = i, runCount = runs }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ subrunDynamics1 d sc r++-- | Run an experiment consisting of the given number of simulations, where each +-- model is created and then requested sequentially in all integration time points +-- using the specified specs.+runDynamicsSeries :: Dynamics (Dynamics a) -> Specs -> Int -> [IO [a]]+runDynamicsSeries (Dynamics m) sc runs = map f [1 .. runs]+ where f i =+ do let r = Run { runIndex = i, runCount = runs }+ d <- m Point { pointSpecs = sc,+ pointRun = r,+ pointTime = spcStartTime sc,+ pointIteration = 0,+ pointPhase = 0 }+ sequence $ subrunDynamics d sc r++-- | Run the simulation and print the result in the last +-- time point using the specified simulation specs.+printDynamics1 :: (Show a) => Dynamics (Dynamics a) -> Specs -> IO ()+printDynamics1 m sc = runDynamics1 m sc >>= print++-- | Run the simulation and print lazily the results in all+-- integration time points using the specified simulation specs.+printDynamics :: (Show a) => Dynamics (Dynamics a) -> Specs -> IO ()+printDynamics m sc = runDynamicsIO m sc >>= loop+ where loop [] = return ()+ loop (x : xs) = do { a <- x; print a; loop xs }++instance Functor Dynamics where+ fmap = liftMD++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+{-# INLINE liftMD #-}+liftMD f (Dynamics x) =+ Dynamics $ \p -> do { a <- x p; return $ f a }++liftM2D :: (a -> b -> c) -> Dynamics a -> Dynamics b -> Dynamics c+{-# INLINE liftM2D #-}+liftM2D f (Dynamics x) (Dynamics y) =+ Dynamics $ \p -> do { a <- x p; b <- y p; 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
+ Simulation/Aivika/Dynamics/Internal/Fold.hs view
@@ -0,0 +1,91 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.Fold+-- 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 7.0.3+--+-- This module defines the fold functions that allows traversing the values of+-- any dynamic process in the integration time points.+--+module Simulation.Aivika.Dynamics.Internal.Fold+ (foldD1,+ foldD,+ divideD) where++import Data.IORef+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Interpolate+import Simulation.Aivika.Dynamics.Internal.Memo++--+-- Fold+--++-- | Like the standard 'foldl1' function but applied to values in +-- the integration time points. The accumulator values are transformed+-- according to the first argument, which should be either function +-- 'memo0' or 'umemo0'.+foldD1 :: (Dynamics a -> Dynamics (Dynamics a))+ -> (a -> a -> a) + -> Dynamics a + -> Dynamics (Dynamics a)+foldD1 tr f (Dynamics m) =+ do r <- liftIO $ newIORef m+ let z = Dynamics $ \p ->+ case pointIteration p of+ 0 -> + m p+ n -> do + let sc = pointSpecs p+ ty = basicTime sc (n - 1) 0+ py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }+ y <- readIORef r+ s <- y py+ x <- m p+ return $! f s x+ y@(Dynamics m) <- tr z+ liftIO $ writeIORef r m+ return y++-- | Like the standard 'foldl' function but applied to values in +-- the integration time points. The accumulator values are transformed+-- according to the first argument, which should be either function+-- 'memo0' or 'umemo0'.+foldD :: (Dynamics a -> Dynamics (Dynamics a))+ -> (a -> b -> a) + -> a+ -> Dynamics b + -> Dynamics (Dynamics a)+foldD tr f acc (Dynamics m) =+ do r <- liftIO $ newIORef $ const $ return acc+ let z = Dynamics $ \p ->+ case pointIteration p of+ 0 -> do+ x <- m p+ return $! f acc x+ n -> do + let sc = pointSpecs p+ ty = basicTime sc (n - 1) 0+ py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }+ y <- readIORef r+ s <- y py+ x <- m p+ return $! f s x+ y@(Dynamics m) <- tr z+ liftIO $ writeIORef r m+ return y++-- | Divide the values in integration time points by the number of+-- the current iteration. It can be useful for statistic functions in+-- combination with the fold.+divideD :: Dynamics Double -> Dynamics Double+divideD (Dynamics m) = + discrete $ Dynamics $ \p ->+ do a <- m p+ return $ a / fromInteger (toInteger (pointIteration p + 1))
+ Simulation/Aivika/Dynamics/Internal/Interpolate.hs view
@@ -0,0 +1,62 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.Interpolate+-- 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 7.0.3+--+-- This module defines interpolation functions.+-- These functions complement the memoization, possibly except for +-- the 'initD' function which is useful to get an initial +-- value of any dynamic process.+--++module Simulation.Aivika.Dynamics.Internal.Interpolate+ (initD,+ discrete,+ interpolate) where++import Simulation.Aivika.Dynamics.Internal.Dynamics++-- | Return the initial value.+initD :: Dynamics a -> Dynamics a+{-# INLINE initD #-}+initD (Dynamics m) =+ Dynamics $ \p ->+ if pointIteration p == 0 && pointPhase p == 0 then+ m p+ else+ let sc = pointSpecs p+ in m $ p { pointTime = basicTime sc 0 0,+ pointIteration = 0,+ pointPhase = 0 } ++-- | Discretize the computation in the integration time points.+discrete :: Dynamics a -> Dynamics a+{-# INLINE discrete #-}+discrete (Dynamics m) =+ Dynamics $ \p ->+ if pointPhase p == 0 then+ m p+ else+ let sc = pointSpecs p+ n = pointIteration p+ in m $ p { pointTime = basicTime sc n 0,+ pointPhase = 0 }++-- | Interpolate the computation based on the integration time points only.+-- Unlike the 'discrete' function it knows about the intermediate time points +-- that are used in the Runge-Kutta method.+interpolate :: Dynamics a -> Dynamics a+{-# INLINE interpolate #-}+interpolate (Dynamics m) = + Dynamics $ \p -> + if pointPhase p >= 0 then + m p+ else + let sc = pointSpecs p+ n = pointIteration p+ in m $ p { pointTime = basicTime sc n 0,+ pointPhase = 0 }
+ Simulation/Aivika/Dynamics/Internal/Memo.hs view
@@ -0,0 +1,195 @@++{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.Memo+-- 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 7.0.3+--+-- This module defines memo functions. The memoization creates such dynamic processes, +-- which values are cached in the integration time points. Then these values are +-- interpolated in all other time points.+--++module Simulation.Aivika.Dynamics.Internal.Memo+ (memo,+ umemo,+ memo0,+ umemo0,+ iterateD) where++import Data.Array+import Data.Array.IO+import Data.IORef+import Control.Monad++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Interpolate++newMemoArray_ :: Ix i => (i, i) -> IO (IOArray i e)+newMemoArray_ = newArray_++newMemoUArray_ :: (MArray IOUArray e IO, Ix i) => (i, i) -> IO (IOUArray i e)+newMemoUArray_ = newArray_++-- | Memoize and order the computation in the integration time points using +-- the interpolation that knows of the Runge-Kutta method.+memo :: Dynamics e -> Dynamics (Dynamics e)+{-# INLINE memo #-}+memo (Dynamics m) = + Dynamics $ \p ->+ do let sc = pointSpecs p+ (phl, phu) = phaseBnds sc+ (nl, nu) = iterationBnds sc+ arr <- newMemoArray_ ((phl, nl), (phu, nu))+ nref <- newIORef 0+ phref <- newIORef 0+ let r p = + do let sc = pointSpecs p+ n = pointIteration p+ ph = pointPhase p+ phu = phaseHiBnd sc + loop n' ph' = + if (n' > n) || ((n' == n) && (ph' > ph)) + then + readArray arr (ph, n)+ else + let p' = p { pointIteration = n', pointPhase = ph',+ pointTime = basicTime sc n' ph' }+ in do a <- m p'+ 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 $ interpolate $ Dynamics r++-- | This is a more efficient version the 'memo' function which uses +-- an unboxed array to store the values.+umemo :: (MArray IOUArray e IO) => Dynamics e -> Dynamics (Dynamics e)+{-# INLINE umemo #-}+umemo (Dynamics m) = + Dynamics $ \p ->+ do let sc = pointSpecs p+ (phl, phu) = phaseBnds sc+ (nl, nu) = iterationBnds sc+ arr <- newMemoUArray_ ((phl, nl), (phu, nu))+ nref <- newIORef 0+ phref <- newIORef 0+ let r p =+ do let sc = pointSpecs p+ n = pointIteration p+ ph = pointPhase p+ phu = phaseHiBnd sc + loop n' ph' = + if (n' > n) || ((n' == n) && (ph' > ph)) + then + readArray arr (ph, n)+ else + let p' = p { pointIteration = n', + pointPhase = ph',+ pointTime = basicTime sc n' ph' }+ in do a <- m p'+ 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 $ interpolate $ Dynamics r++-- | Memoize and order the computation in the integration time points using +-- the 'discrete' interpolation. It consumes less memory than the 'memo'+-- function but it is not aware of the Runge-Kutta method. There is a subtle+-- difference when we request for values in the intermediate time points+-- that are used by this method to integrate. In general case you should +-- prefer the 'memo0' function above 'memo'.+memo0 :: Dynamics e -> Dynamics (Dynamics e)+{-# INLINE memo0 #-}+memo0 (Dynamics m) = + Dynamics $ \p ->+ do let sc = pointSpecs p+ bnds = iterationBnds sc+ arr <- newMemoArray_ bnds+ nref <- newIORef 0+ let r p =+ do let sc = pointSpecs p+ n = pointIteration p+ loop n' = + if n' > n+ then + readArray arr n+ else + let p' = p { pointIteration = n', pointPhase = 0,+ pointTime = basicTime sc n' 0 }+ in do a <- m p'+ a `seq` writeArray arr n' a+ writeIORef nref (n' + 1)+ loop (n' + 1)+ n' <- readIORef nref+ loop n'+ return $ discrete $ Dynamics r++-- | This is a more efficient version the 'memo0' function which uses +-- an unboxed array to store the values.+umemo0 :: (MArray IOUArray e IO) => Dynamics e -> Dynamics (Dynamics e)+{-# INLINE umemo0 #-}+umemo0 (Dynamics m) = + Dynamics $ \p ->+ do let sc = pointSpecs p+ bnds = iterationBnds sc+ arr <- newMemoUArray_ bnds+ nref <- newIORef 0+ let r p =+ do let sc = pointSpecs p+ n = pointIteration p+ loop n' = + if n' > n+ then + readArray arr n+ else + let p' = p { pointIteration = n', pointPhase = 0,+ pointTime = basicTime sc n' 0 }+ in do a <- m p'+ a `seq` writeArray arr n' a+ writeIORef nref (n' + 1)+ loop (n' + 1)+ n' <- readIORef nref+ loop n'+ return $ discrete $ Dynamics r++-- | Iterate sequentially the dynamic process with side effects in +-- the integration time points. It is equivalent to a call of the+-- 'memo0' function but significantly more efficient, for the array +-- is not created.+iterateD :: Dynamics () -> Dynamics (Dynamics ())+{-# INLINE iterateD #-}+iterateD (Dynamics m) = + Dynamics $ \p ->+ do let sc = pointSpecs p+ nref <- newIORef 0+ let r p =+ do let sc = pointSpecs p+ n = pointIteration p+ loop n' = + unless (n' > n) $+ let p' = p { pointIteration = n', pointPhase = 0,+ pointTime = basicTime sc n' 0 }+ in do a <- m p'+ a `seq` writeIORef nref (n' + 1)+ loop (n' + 1)+ n' <- readIORef nref+ loop n'+ return $ discrete $ Dynamics r
+ Simulation/Aivika/Dynamics/Internal/Process.hs view
@@ -0,0 +1,162 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.Process+-- 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 7.0.3+--+-- A value in the 'Process' monad represents a discontinuous process that +-- can suspend and resume at any time. It behaves like a dynamic process too. +-- Any value in the 'Dynamics' monad can be lifted to the Process monad. +-- Moreover, a value in the Process monad can be run in the Dynamics monad.+--+-- A value of the 'ProcessID' type is just an identifier of such a process.+--+module Simulation.Aivika.Dynamics.Internal.Process+ (ProcessID,+ Process(..),+ processQueue,+ newProcessID,+ holdProcess,+ passivateProcess,+ processPassive,+ reactivateProcess,+ processID,+ runProcess) where++import Data.IORef+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Cont+import Simulation.Aivika.Dynamics.Lift+import Simulation.Aivika.Dynamics.EventQueue++-- | Represents a process identificator.+data ProcessID = + ProcessID { processQueue :: EventQueue, -- ^ Return the event queue.+ processStarted :: IORef Bool,+ processCont :: IORef (Maybe (Dynamics (() -> IO ()))) }++-- | Specifies a discontinuous process that can suspend at any time+-- and then resume later.+newtype Process a = Process (ProcessID -> Cont a)++-- | Hold the process for the specified time period.+holdProcess :: Double -> Process ()+holdProcess dt =+ Process $ \pid ->+ Cont $ \c ->+ Dynamics $ \ps ->+ do let Dynamics m = enqueueCont (processQueue pid) (pointTime ps + dt) c+ m ps++-- | Passivate the process.+passivateProcess :: Process ()+passivateProcess =+ Process $ \pid ->+ Cont $ \c ->+ Dynamics $ \p ->+ do let x = processCont 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 ID is passivated.+processPassive :: ProcessID -> Process Bool+processPassive pid =+ Process $ \_ ->+ Cont $ \(Dynamics c) ->+ Dynamics $ \p ->+ do cont' <- c p+ let x = processCont pid+ a <- readIORef x+ case a of+ Nothing -> cont' False+ Just _ -> cont' True++-- | Reactivate a process with the specified ID.+reactivateProcess :: ProcessID -> Process ()+reactivateProcess pid =+ Process $ \pid' ->+ Cont $ \c@(Dynamics cont) ->+ Dynamics $ \p ->+ do let x = processCont pid+ a <- readIORef x+ case a of+ Nothing ->+ do cont' <- cont p+ cont' ()+ Just (Dynamics cont2) ->+ do writeIORef x Nothing+ let Dynamics m = enqueueCont (processQueue pid') (pointTime p) c+ m p+ cont2' <- cont2 p+ cont2' ()++-- | Start the process with the specified ID at the desired time.+runProcess :: Process () -> ProcessID -> Double -> Dynamics ()+runProcess (Process p) pid t =+ runCont m r+ where m = do y <- liftIO $ readIORef (processStarted pid)+ if y + then error $+ "A process with such ID " +++ "has been started already: runProc"+ else liftIO $ writeIORef (processStarted pid) True+ Cont $ \c -> enqueueCont (processQueue pid) t c+ p pid+ r = let f () = return () in return f++-- | Return the current process ID.+processID :: Process ProcessID+processID = Process $ \pid -> return pid++-- | Create a new process ID.+newProcessID :: EventQueue -> Dynamics ProcessID+newProcessID q =+ do x <- liftIO $ newIORef Nothing+ y <- liftIO $ newIORef False+ return ProcessID { processQueue = q,+ processStarted = y,+ processCont = x }++instance Eq ProcessID where+ x == y = processCont x == processCont y -- for the references are unique++instance Monad Process where+ return = returnP+ m >>= k = bindP m k++instance Functor Process where+ fmap = liftM++instance Lift Process where+ liftD = liftP+ +instance MonadIO Process where+ liftIO = liftIOP+ +returnP :: a -> Process a+{-# INLINE returnP #-}+returnP a = Process (\pid -> return a)++bindP :: Process a -> (a -> Process b) -> Process b+{-# INLINE bindP #-}+bindP (Process m) k = + Process $ \pid -> + do a <- m pid+ let Process m' = k a+ m' pid++liftP :: Dynamics a -> Process a+{-# INLINE liftP #-}+liftP m = Process $ \pid -> liftD m++liftIOP :: IO a -> Process a+{-# INLINE liftIOP #-}+liftIOP m = Process $ \pid -> liftIO m
+ Simulation/Aivika/Dynamics/Internal/Time.hs view
@@ -0,0 +1,35 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Internal.Time+-- 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 7.0.3+--+-- This module defines the time parameters.+--++module Simulation.Aivika.Dynamics.Internal.Time+ (starttime,+ stoptime,+ dt,+ time) where++import Simulation.Aivika.Dynamics.Internal.Dynamics++-- | Return the start simulation time.+starttime :: Dynamics Double+starttime = Dynamics $ return . spcStartTime . pointSpecs++-- | Return the stop simulation time.+stoptime :: Dynamics Double+stoptime = Dynamics $ return . spcStopTime . pointSpecs++-- | Return the integration time step.+dt :: Dynamics Double+dt = Dynamics $ return . spcDT . pointSpecs++-- | Return the current simulation time.+time :: Dynamics Double+time = Dynamics $ return . pointTime
+ Simulation/Aivika/Dynamics/Lift.hs view
@@ -0,0 +1,21 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Lift+-- 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 7.0.3+--+-- This module defines the 'liftD' function that allows embedding+-- the 'Dynamics' computation.+--+module Simulation.Aivika.Dynamics.Lift (Lift(..)) where++import Simulation.Aivika.Dynamics++-- | The 'Lift' class defines a type which the 'Dynamics' +-- computation can be lifted to.+class Lift m where+ -- | Lift the computation.+ liftD :: Dynamics a -> m a
+ Simulation/Aivika/Dynamics/Parameter.hs view
@@ -0,0 +1,58 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Parameter+-- 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 7.0.3+--+-- This module defines the parameters of simulation experiments.+--++module Simulation.Aivika.Dynamics.Parameter+ (newParameter,+ newTableParameter,+ newIndexedParameter) where++import Data.Array+import Data.IORef+import qualified Data.Map as M+import Control.Concurrent.MVar++import Simulation.Aivika.Dynamics.Internal.Dynamics++-- | Create a thread-safe parameter that returns always the same value during the simulation run, +-- where the value is recalculated for each new run.+newParameter :: IO a -> IO (Dynamics a)+newParameter a = newIndexedParameter $ \_ -> a++-- | Create a thread-safe parameter that returns always the same value during the simulation run,+-- where the value is taken consequently from the specified table based on the number of the +-- current run starting from zero. After all values from the table are used, it takes the first +-- value of the table, then the second one and so on.+newTableParameter :: Array Int a -> IO (Dynamics a)+newTableParameter t = newIndexedParameter (\i -> return $ t ! (((i - i1) `mod` n) + i1))+ where (i1, i2) = bounds t+ n = i2 - i1 + 1++-- | Create a thread-safe parameter that returns always the same value during the simulation run, +-- where the value depends on the number of this run starting from zero.+newIndexedParameter :: (Int -> IO a) -> IO (Dynamics a)+newIndexedParameter f = + do lock <- newMVar ()+ dict <- newIORef M.empty+ return $ Dynamics $ \p ->+ do let i = runIndex $ pointRun p+ m <- readIORef dict+ if M.member i m+ then do let Just v = M.lookup i m+ return v+ else withMVar lock $ + \() -> do { m <- readIORef dict;+ if M.member i m+ then do let Just v = M.lookup i m+ return v+ else do v <- f i+ writeIORef dict $ M.insert i v m+ return v }
+ Simulation/Aivika/Dynamics/Process.hs view
@@ -0,0 +1,30 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Process+-- 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 7.0.3+--+-- A value in the 'Process' monad represents a discontinuous process that +-- can suspend and resume at any time. It behaves like a dynamic process too. +-- Any value in the 'Dynamics' monad can be lifted to the Process monad. +-- Moreover, a value in the Process monad can be run in the Dynamics monad.+--+-- A value of the 'ProcessID' type is just an identifier of such a process.+--+module Simulation.Aivika.Dynamics.Process+ (ProcessID,+ Process,+ processQueue,+ newProcessID,+ holdProcess,+ passivateProcess,+ processPassive,+ reactivateProcess,+ processID,+ runProcess) where++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Process
+ Simulation/Aivika/Dynamics/Random.hs view
@@ -0,0 +1,73 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Random+-- 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 7.0.3+--+-- Below are defined random functions that mostly return discrete processes. +-- Literally, it means that the values are initially defined in integration +-- time points and then they are passed to the 'discrete' function.+--++module Simulation.Aivika.Dynamics.Random + (newRandom, newNormal, normalGen) where++import Random+import Data.IORef+import Control.Monad.Trans++import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base++-- | Return the uniform random numbers between 0.0 and 1.0 in+-- the integration time points.+newRandom :: Dynamics (Dynamics Double)+newRandom =+ memo0 $ liftIO $ getStdRandom random+ +-- | Return the normal random numbers with mean 0.0 and variance 1.0 in+-- the integration time points.+newNormal :: Dynamics (Dynamics Double)+newNormal =+ do g <- liftIO normalGen+ memo0 $ liftIO g++-- | Normal random number generator.+normalGen :: IO (IO Double)+normalGen =+ do nextRef <- newIORef 0.0+ flagRef <- newIORef False+ xi1Ref <- newIORef 0.0+ xi2Ref <- newIORef 0.0+ psiRef <- newIORef 0.0+ let loop =+ do psi <- readIORef psiRef+ if (psi >= 1.0) || (psi == 0.0)+ then do g1 <- getStdRandom random+ g2 <- getStdRandom random+ let xi1 = 2.0 * g1 - 1.0+ xi2 = 2.0 * g2 - 1.0+ psi = xi1 * xi1 + xi2 * xi2+ writeIORef xi1Ref xi1+ writeIORef xi2Ref xi2+ writeIORef psiRef psi+ loop+ else writeIORef psiRef $ sqrt (- 2.0 * log psi / psi)+ return $+ do flag <- readIORef flagRef+ if flag+ then do writeIORef flagRef False+ readIORef nextRef+ else do writeIORef xi1Ref 0.0+ writeIORef xi2Ref 0.0+ writeIORef psiRef 0.0+ loop+ xi1 <- readIORef xi1Ref+ xi2 <- readIORef xi2Ref+ psi <- readIORef psiRef+ writeIORef flagRef True+ writeIORef nextRef $ xi2 * psi+ return $ xi1 * psi
+ Simulation/Aivika/Dynamics/Ref.hs view
@@ -0,0 +1,63 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Ref+-- 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 7.0.3+--+-- This module defines an updatable reference that depends on the event queue.+--+module Simulation.Aivika.Dynamics.Ref+ (Ref,+ newRef,+ refQueue,+ readRef,+ writeRef,+ modifyRef) where++import Data.IORef+import Control.Monad.Trans++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.EventQueue++-- | The 'Ref' 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 Ref a = + Ref { refQueue :: EventQueue, -- ^ Return the bound event queue.+ refRun :: Dynamics (),+ refValue :: IORef a }++-- | Create a new reference bound to the specified event queue.+newRef :: EventQueue -> a -> Dynamics (Ref a)+newRef q a =+ do x <- liftIO $ newIORef a+ return Ref { refQueue = q,+ refRun = queueRun q,+ refValue = x }+ +-- | Read the value of a reference, forcing the bound event queue to raise +-- the events in case of need.+readRef :: Ref a -> Dynamics a+readRef r = Dynamics $ \p -> + do let Dynamics m = refRun r+ m p+ readIORef (refValue r)++-- | Write a new value into the reference.+writeRef :: Ref a -> a -> Dynamics ()+writeRef r a = Dynamics $ \p -> + a `seq` writeIORef (refValue r) a++-- | Mutate the contents of the reference, forcing the bound event queue to+-- raise all pending events in case of need.+modifyRef :: Ref a -> (a -> a) -> Dynamics ()+modifyRef r f = Dynamics $ \p -> + do let Dynamics m = refRun r+ m p+ a <- readIORef (refValue r)+ let b = f a+ b `seq` writeIORef (refValue r) b
+ Simulation/Aivika/Dynamics/Resource.hs view
@@ -0,0 +1,101 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Resource+-- 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 7.0.3+--+-- This module defines a limited resource which can be acquired and +-- then released by the discontinuous process 'DynamicProc'.+--+module Simulation.Aivika.Dynamics.Resource+ (Resource,+ newResource,+ resourceQueue,+ resourceInitCount,+ resourceCount,+ requestResource,+ releaseResource) where++import Data.IORef+import Control.Monad++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Internal.Cont+import Simulation.Aivika.Dynamics.Internal.Process+import Simulation.Aivika.Dynamics.EventQueue+import qualified Simulation.Aivika.Queue as Q++-- | Represents a limited resource.+data Resource = + Resource { resourceQueue :: EventQueue, + -- ^ Return the bound event queue.+ resourceInitCount :: Int,+ -- ^ Return the initial count of the resource.+ resourceCountRef :: IORef Int, + resourceWaitQueue :: Q.Queue (Dynamics (() -> IO ()))}++instance Eq Resource where+ x == y = resourceCountRef x == resourceCountRef y -- unique references++-- | Create a new resource with the specified initial count.+newResource :: EventQueue -> Int -> Dynamics Resource+newResource q initCount =+ Dynamics $ \p ->+ do countRef <- newIORef initCount+ waitQueue <- Q.newQueue+ return Resource { resourceQueue = q,+ resourceInitCount = initCount,+ resourceCountRef = countRef,+ resourceWaitQueue = waitQueue }++-- | Return the current count of the resource.+resourceCount :: Resource -> Process Int+resourceCount r =+ Process $ \_ ->+ Cont $ \(Dynamics c) ->+ Dynamics $ \p ->+ do cont' <- c p + 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 :: Resource -> Process ()+requestResource r =+ Process $ \_ ->+ Cont $ \c@(Dynamics cont) ->+ Dynamics $ \p ->+ 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 p+ cont' ()++-- | Release the resource increasing its count and resuming one of the+-- previously suspended processes as possible.+releaseResource :: Resource -> Process ()+releaseResource r =+ Process $ \_ ->+ Cont $ \(Dynamics c) ->+ Dynamics $ \p ->+ 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 = enqueueCont (resourceQueue r) (pointTime p) c2+ m p+ cont' <- c p+ cont' ()
+ Simulation/Aivika/Dynamics/SystemDynamics.hs view
@@ -0,0 +1,441 @@++{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module : Simulation.Aivika.Dynamics.SystemDynamics+-- 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 7.0.3+--+-- This module defines integrals and other functions of System Dynamics.+--++module Simulation.Aivika.Dynamics.SystemDynamics+ (-- * Maximum and Minimum+ maxD,+ minD,+ -- * Integrals+ Integ,+ newInteg,+ integInit,+ integValue,+ integDiff,+ -- * Integral Functions+ integ,+ -- * Difference Equations+ Sum,+ newSum,+ sumInit,+ sumValue,+ sumDiff,+ -- * Table Functions+ lookupD,+ lookupStepwiseD) where++import Data.Array+import Data.Array.IO+import Data.IORef+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.Base++--+-- Maximum and Minimum+--++-- | Return the maximum.+maxD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a+maxD = liftM2 max++-- | Return the minimum.+minD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a+minD = liftM2 min++--+-- Integrals+--++-- | 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 $ \p -> + do (Dynamics m) <- readIORef (integInternal integ)+ m p+ y <- umemo z+ liftIO $ writeIORef (integExternal integ) y+ return integ++-- | Return the integral's value.+integValue :: Integ -> Dynamics Double+integValue integ = + Dynamics $ \p ->+ do (Dynamics m) <- readIORef (integExternal integ)+ m p++-- | Set the derivative for the integral.+integDiff :: Integ -> Dynamics Double -> Dynamics ()+integDiff integ diff =+ do let z = Dynamics $ \p ->+ do y <- readIORef (integExternal integ)+ let i = integInit integ+ case spcMethod (pointSpecs p) of+ Euler -> integEuler diff i y p+ RungeKutta2 -> integRK2 diff i y p+ RungeKutta4 -> integRK4 diff i y p+ liftIO $ writeIORef (integInternal integ) z++integEuler :: Dynamics Double+ -> Dynamics Double + -> Dynamics Double + -> Point -> IO Double+integEuler (Dynamics f) (Dynamics i) (Dynamics y) p = + case pointIteration p of+ 0 -> + i p+ n -> do + let sc = pointSpecs p+ ty = basicTime sc (n - 1) 0+ py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }+ a <- y py+ b <- f py+ let !v = a + spcDT (pointSpecs p) * b+ return v++integRK2 :: Dynamics Double+ -> Dynamics Double+ -> Dynamics Double+ -> Point -> IO Double+integRK2 (Dynamics f) (Dynamics i) (Dynamics y) p =+ case pointPhase p of+ 0 -> case pointIteration p of+ 0 ->+ i p+ n -> do+ let sc = pointSpecs p+ ty = basicTime sc (n - 1) 0+ t1 = ty+ t2 = basicTime sc (n - 1) 1+ py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }+ p1 = py+ p2 = p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }+ vy <- y py+ k1 <- f p1+ k2 <- f p2+ let !v = vy + spcDT sc / 2.0 * (k1 + k2)+ return v+ 1 -> do+ let sc = pointSpecs p+ n = pointIteration p+ ty = basicTime sc n 0+ t1 = ty+ py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+ p1 = py+ vy <- y py+ k1 <- f p1+ let !v = vy + spcDT sc * k1+ return v+ _ -> + error "Incorrect phase: integRK2"++integRK4 :: Dynamics Double+ -> Dynamics Double+ -> Dynamics Double+ -> Point -> IO Double+integRK4 (Dynamics f) (Dynamics i) (Dynamics y) p =+ case pointPhase p of+ 0 -> case pointIteration p of+ 0 -> + i p+ n -> do+ let sc = pointSpecs p+ 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+ py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }+ p1 = py+ p2 = p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }+ p3 = p { pointTime = t3, pointIteration = n - 1, pointPhase = 2 }+ p4 = p { pointTime = t4, pointIteration = n - 1, pointPhase = 3 }+ vy <- y py+ k1 <- f p1+ k2 <- f p2+ k3 <- f p3+ k4 <- f p4+ let !v = vy + spcDT sc / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)+ return v+ 1 -> do+ let sc = pointSpecs p+ n = pointIteration p+ ty = basicTime sc n 0+ t1 = ty+ py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+ p1 = py+ vy <- y py+ k1 <- f p1+ let !v = vy + spcDT sc / 2.0 * k1+ return v+ 2 -> do+ let sc = pointSpecs p+ n = pointIteration p+ ty = basicTime sc n 0+ t2 = basicTime sc n 1+ py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+ p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }+ vy <- y py+ k2 <- f p2+ let !v = vy + spcDT sc / 2.0 * k2+ return v+ 3 -> do+ let sc = pointSpecs p+ n = pointIteration p+ ty = basicTime sc n 0+ t3 = basicTime sc n 2+ py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+ p3 = p { pointTime = t3, pointIteration = n, pointPhase = 2 }+ vy <- y py+ k3 <- f p3+ let !v = vy + spcDT sc * k3+ return v+ _ -> + error "Incorrect phase: integRK4"++-- 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++--+-- Integral Functions+--++-- | Return an integral with the specified derivative and initial value.+-- If you want to create a loopback then you should use the 'Integ' type +-- directly. The 'integ' function is just a wrapper that uses this type.+integ :: Dynamics Double -> Dynamics Double -> Dynamics (Dynamics Double)+integ diff i =+ do x <- newInteg i+ integDiff x diff+ return $ integValue x++--+-- Difference Equations+--++-- | The 'Sum' type represents a sum defined by some difference equation.+data Sum a = Sum { sumInit :: Dynamics a, -- ^ The initial value.+ sumExternal :: IORef (Dynamics a),+ sumInternal :: IORef (Dynamics a) }++-- | Create a new sum with the specified initial value.+newSum :: (MArray IOUArray a IO, Num a) => Dynamics a -> Dynamics (Sum a)+newSum i = + do r1 <- liftIO $ newIORef $ initD i + r2 <- liftIO $ newIORef $ initD i + let sum = Sum { sumInit = i, + sumExternal = r1,+ sumInternal = r2 }+ z = Dynamics $ \p -> + do (Dynamics m) <- readIORef (sumInternal sum)+ m p+ y <- umemo0 z+ liftIO $ writeIORef (sumExternal sum) y+ return sum++-- | Return the total sum defined by the difference equation.+sumValue :: Sum a -> Dynamics a+sumValue sum = + Dynamics $ \p ->+ do (Dynamics m) <- readIORef (sumExternal sum)+ m p++-- | Set the difference equation for the sum.+sumDiff :: (MArray IOUArray a IO, Num a) => Sum a -> Dynamics a -> Dynamics ()+sumDiff sum (Dynamics diff) =+ do let z = Dynamics $ \p ->+ case pointIteration p of+ 0 -> do+ let Dynamics i = sumInit sum+ i p+ n -> do + Dynamics y <- readIORef (sumExternal sum)+ let sc = pointSpecs p+ ty = basicTime sc (n - 1) 0+ py = p { pointTime = ty, + pointIteration = n - 1, + pointPhase = 0 }+ a <- y py+ b <- diff py+ let !v = a + b+ return v+ liftIO $ writeIORef (sumInternal sum) z++--+-- 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 (\p -> do a <- m p; 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 (\p -> do a <- m p; 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 p = do +-- let t = parTime p+-- sc = parSpecs p+-- n = parIteration p+-- a <- d p+-- let t' = (t - a) - spcStartTime sc+-- n' = fromInteger $ toInteger $ floor $ t' / spcDT sc+-- y | n' < 0 = i $ p { pointTime = spcStartTime sc,+-- pointIteration = 0, +-- pointPhase = 0 }+-- | n' < n = x $ p { pointTime = t',+-- pointIteration = n',+-- pointPhase = -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
+ Simulation/Aivika/Dynamics/UVar.hs view
@@ -0,0 +1,126 @@++{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module : Simulation.Aivika.Dynamics.UVar+-- 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 7.0.3+--+-- This module defines a variable that is bound to the event queue and +-- that keeps the history of changes storing the values in an unboxed array.+--+module Simulation.Aivika.Dynamics.UVar+ (UVar,+ newUVar,+ uvarQueue,+ readUVar,+ writeUVar,+ modifyUVar,+ freezeUVar) where++import Data.Array+import Data.Array.IO+import Data.IORef++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.EventQueue++import qualified Simulation.Aivika.UVector as UV++-- | A version of the 'Var' type which uses an unboxed array to store the values +-- in time points. You should prefer this type whenever possible.+data UVar a = + UVar { uvarQueue :: EventQueue, -- ^ Return the bound event queue.+ uvarRun :: Dynamics (),+ uvarXS :: UV.UVector Double, + uvarYS :: UV.UVector a}+ +-- | Create a new variable bound to the specified event queue.+newUVar :: (MArray IOUArray a IO) => EventQueue -> a -> Dynamics (UVar a)+newUVar q a =+ Dynamics $ \p ->+ do xs <- UV.newVector+ ys <- UV.newVector+ UV.appendVector xs $ spcStartTime $ pointSpecs p+ UV.appendVector ys a+ return UVar { uvarQueue = q,+ uvarRun = queueRun q,+ uvarXS = xs,+ uvarYS = ys }++-- | Read the value of a variable, forcing the bound event queue to raise +-- the events in case of need.+readUVar :: (MArray IOUArray a IO) => UVar a -> Dynamics a+readUVar v =+ Dynamics $ \p ->+ do let Dynamics m = uvarRun v+ m p+ let xs = uvarXS v+ ys = uvarYS v+ t = pointTime p+ count <- UV.vectorCount xs+ let i = count - 1+ x <- UV.readVector xs i+ if x <= t + then UV.readVector ys i+ else do i <- UV.vectorBinarySearch xs t+ if i >= 0+ then UV.readVector ys i+ else UV.readVector ys $ - (i + 1) - 1++-- | Write a new value into the variable.+writeUVar :: (MArray IOUArray a IO) => UVar a -> a -> Dynamics ()+writeUVar v a =+ Dynamics $ \p ->+ do let xs = uvarXS v+ ys = uvarYS v+ t = pointTime p+ count <- UV.vectorCount xs+ let i = count - 1+ x <- UV.readVector xs i+ if t < x + then error "Cannot update the past data: writeUVar."+ else if t == x+ then UV.writeVector ys i $! a+ else do UV.appendVector xs t+ UV.appendVector ys $! a++-- | Mutate the contents of the variable, forcing the bound event queue to+-- raise all pending events in case of need.+modifyUVar :: (MArray IOUArray a IO) => UVar a -> (a -> a) -> Dynamics ()+modifyUVar v f =+ Dynamics $ \p ->+ do let Dynamics m = uvarRun v+ m p+ let xs = uvarXS v+ ys = uvarYS v+ t = pointTime p+ count <- UV.vectorCount xs+ let i = count - 1+ x <- UV.readVector xs i+ if t < x+ then error "Cannot update the past data: modifyUVar."+ else if t == x+ then do a <- UV.readVector ys i+ UV.writeVector ys i $! f a+ else do i <- UV.vectorBinarySearch xs t+ if i >= 0+ then do a <- UV.readVector ys i+ UV.appendVector xs t+ UV.appendVector ys $! f a+ else do a <- UV.readVector ys $ - (i + 1) - 1+ UV.appendVector xs t+ UV.appendVector ys $! f a++-- | Freeze the variable and return in arrays the time points and corresponded +-- values when the variable had changed.+freezeUVar :: (MArray IOUArray a IO) => + UVar a -> Dynamics (Array Int Double, Array Int a)+freezeUVar v =+ Dynamics $ \p ->+ do xs <- UV.freezeVector (uvarXS v)+ ys <- UV.freezeVector (uvarYS v)+ return (xs, ys)
+ Simulation/Aivika/Dynamics/Var.hs view
@@ -0,0 +1,126 @@++-- |+-- Module : Simulation.Aivika.Dynamics.Var+-- 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 7.0.3+--+-- This module defines a variable that is bound to the event queue and +-- that keeps the history of changes storing the values in an array.+--+module Simulation.Aivika.Dynamics.Var+ (Var,+ newVar,+ varQueue,+ readVar,+ writeVar,+ modifyVar,+ freezeVar) where++import Data.Array+import Data.Array.IO+import Data.IORef++import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Dynamics.EventQueue++import qualified Simulation.Aivika.Vector as V+import qualified Simulation.Aivika.UVector as UV++-- | Like the 'Ref' reference but keeps the history of changes in +-- different time points. The 'Var' variable is safe in the hybrid +-- simulation and when you use different event queues, but this variable is +-- slower than references.+data Var a = + Var { varQueue :: EventQueue, -- ^ Return the bound event queue.+ varRun :: Dynamics (),+ varXS :: UV.UVector Double, + varYS :: V.Vector a}+ +-- | Create a new variable bound to the specified event queue.+newVar :: EventQueue -> a -> Dynamics (Var a)+newVar q a =+ Dynamics $ \p ->+ do xs <- UV.newVector+ ys <- V.newVector+ UV.appendVector xs $ spcStartTime $ pointSpecs p+ V.appendVector ys a+ return Var { varQueue = q,+ varRun = queueRun q,+ varXS = xs,+ varYS = ys }++-- | Read the value of a variable, forcing the bound event queue to raise +-- the events in case of need.+readVar :: Var a -> Dynamics a+readVar v =+ Dynamics $ \p ->+ do let Dynamics m = varRun v+ m p+ let xs = varXS v+ ys = varYS v+ t = pointTime p+ count <- UV.vectorCount xs+ let i = count - 1+ x <- UV.readVector xs i+ if x <= t + then V.readVector ys i+ else do i <- UV.vectorBinarySearch xs t+ if i >= 0+ then V.readVector ys i+ else V.readVector ys $ - (i + 1) - 1++-- | Write a new value into the variable.+writeVar :: Var a -> a -> Dynamics ()+writeVar v a =+ Dynamics $ \p ->+ do let xs = varXS v+ ys = varYS v+ t = pointTime p+ count <- UV.vectorCount xs+ let i = count - 1+ x <- UV.readVector xs i+ if t < x + then error "Cannot update the past data: writeVar."+ else if t == x+ then V.writeVector ys i $! a+ else do UV.appendVector xs t+ V.appendVector ys $! a++-- | Mutate the contents of the variable, forcing the bound event queue to+-- raise all pending events in case of need.+modifyVar :: Var a -> (a -> a) -> Dynamics ()+modifyVar v f =+ Dynamics $ \p ->+ do let Dynamics m = varRun v+ m p+ let xs = varXS v+ ys = varYS v+ t = pointTime p+ count <- UV.vectorCount xs+ let i = count - 1+ x <- UV.readVector xs i+ if t < x+ then error "Cannot update the past data: modifyVar."+ else if t == x+ then do a <- V.readVector ys i+ V.writeVector ys i $! f a+ else do i <- UV.vectorBinarySearch xs t+ if i >= 0+ then do a <- V.readVector ys i+ UV.appendVector xs t+ V.appendVector ys $! f a+ else do a <- V.readVector ys $ - (i + 1) - 1+ UV.appendVector xs t+ V.appendVector ys $! f a++-- | Freeze the variable and return in arrays the time points and corresponded +-- values when the variable had changed.+freezeVar :: Var a -> Dynamics (Array Int Double, Array Int a)+freezeVar v =+ Dynamics $ \p ->+ do xs <- UV.freezeVector (varXS v)+ ys <- V.freezeVector (varYS v)+ return (xs, ys)
Simulation/Aivika/PriorityQueue.hs view
@@ -1,40 +1,18 @@ --- 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+-- 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 7.0.3+--+-- An imperative heap-based priority queue.+-- module Simulation.Aivika.PriorityQueue (PriorityQueue, queueNull, + queueCount, newQueue, enqueue, dequeue, @@ -46,12 +24,12 @@ import Data.IORef import Control.Monad +-- | The 'PriorityQueue' type represents an imperative heap-based +-- priority queue. data PriorityQueue a = PriorityQueue { pqKeys :: IORef (IOUArray Int Double), pqVals :: IORef (IOArray Int a),- pqSize :: IORef Int,- pqNoVal :: a -- to release references - }+ pqSize :: IORef Int } increase :: PriorityQueue a -> Int -> IO () increase pq capacity = @@ -120,13 +98,19 @@ writeArray vals i vn'' siftDown keys vals size n'' k v +-- | Test whether the priority queue is empty. queueNull :: PriorityQueue a -> IO Bool queueNull pq = do size <- readIORef (pqSize pq) return $ size == 0 -newQueue :: a -> IO (PriorityQueue a)-newQueue defaultValue =+-- | Return the number of elements in the priority queue.+queueCount :: PriorityQueue a -> IO Int+queueCount pq = readIORef (pqSize pq)++-- | Create a new priority queue.+newQueue :: IO (PriorityQueue a)+newQueue = do keys <- newArray_ (0, 10) vals <- newArray_ (0, 10) keyRef <- newIORef keys@@ -134,20 +118,21 @@ sizeRef <- newIORef 0 return PriorityQueue { pqKeys = keyRef, pqVals = valRef, - pqSize = sizeRef,- pqNoVal = defaultValue }+ pqSize = sizeRef } +-- | Enqueue a new element with the specified priority. 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)+ when (i >= iu - il) $ increase pq (i + 2) -- plus one element on the end 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 the element with the minimal priority. dequeue :: PriorityQueue a -> IO () dequeue pq = do size <- readIORef (pqSize pq)@@ -156,12 +141,15 @@ 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!+ k <- readArray keys i+ v <- readArray vals i+ k0 <- readArray keys size+ v0 <- readArray vals size+ writeArray keys i k0+ writeArray vals i v0 siftDown keys vals i 0 k v +-- | Return the element with the minimal priority. queueFront :: PriorityQueue a -> IO (Double, a) queueFront pq = do size <- readIORef (pqSize pq)
Simulation/Aivika/Queue.hs view
@@ -1,40 +1,18 @@ --- 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+-- 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 7.0.3+--+-- An imperative double-linked queue.+-- module Simulation.Aivika.Queue (Queue, queueNull, + queueCount, newQueue, enqueue, dequeue, @@ -43,15 +21,19 @@ import Data.IORef import Control.Monad +-- | A cell of the double-linked queue. data QueueItem a = QueueItem { qiVal :: a, qiPrev :: IORef (Maybe (QueueItem a)), qiNext :: IORef (Maybe (QueueItem a)) } +-- | The 'Queue' type represents an imperative double-linked queue. data Queue a = Queue { qHead :: IORef (Maybe (QueueItem a)),- qTail :: IORef (Maybe (QueueItem a)) }+ qTail :: IORef (Maybe (QueueItem a)), + qSize :: IORef Int } +-- | Test whether the queue is empty. queueNull :: Queue a -> IO Bool queueNull q = do head <- readIORef (qHead q) @@ -59,15 +41,24 @@ Nothing -> return True Just _ -> return False +-- | Return the number of elements in the queue.+queueCount :: Queue a -> IO Int+queueCount q = readIORef (qSize q)++-- | Create a new queue. newQueue :: IO (Queue a) newQueue = do head <- newIORef Nothing tail <- newIORef Nothing- return Queue { qHead = head, qTail = tail }+ size <- newIORef 0+ return Queue { qHead = head, qTail = tail, qSize = size } +-- | Enqueue a new element. enqueue :: Queue a -> a -> IO () enqueue q v =- do head <- readIORef (qHead q)+ do size <- readIORef (qSize q)+ writeIORef (qSize q) (size + 1)+ head <- readIORef (qHead q) case head of Nothing -> do prev <- newIORef Nothing@@ -86,6 +77,7 @@ writeIORef (qiPrev h) item writeIORef (qHead q) item +-- | Dequeue the first element. dequeue :: Queue a -> IO () dequeue q = do tail <- readIORef (qTail q) @@ -93,7 +85,9 @@ Nothing -> error "Empty queue: dequeue" Just t ->- do tail' <- readIORef (qiPrev t)+ do size <- readIORef (qSize q)+ writeIORef (qSize q) (size - 1)+ tail' <- readIORef (qiPrev t) case tail' of Nothing -> do writeIORef (qHead q) Nothing@@ -102,6 +96,7 @@ do writeIORef (qiNext t') Nothing writeIORef (qTail q) tail' +-- | Return the first element. queueFront :: Queue a -> IO a queueFront q = do tail <- readIORef (qTail q)
+ Simulation/Aivika/Statistics.hs view
@@ -0,0 +1,118 @@++{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module : Simulation.Aivika.Statistics+-- 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 7.0.3+--+-- Represents statistics and results.+--+module Simulation.Aivika.Statistics+ (Statistics, + newStatistics,+ addStatistics,+ statisticsData,+ AnalysisResults(..),+ analyzeData,+ showResults) where ++import Data.Foldable+import Data.Array+import Data.Array.IO+import Control.Monad+import Control.Concurrent.MVar++import Simulation.Aivika.UVector++-- | Represents statistics. +-- +-- All functions with the statistics in this module are thread-safe. Therefore +-- you can use them in experiments when parallel simulations execute simultaneously.+data Statistics a = Statistics { statData :: UVector a, + statLock :: MVar () }++-- | Create new statistics.+newStatistics :: (MArray IOUArray a IO) => IO (Statistics a)+newStatistics = + do v <- newVector+ l <- newMVar ()+ return Statistics { statData = v, + statLock = l }++-- | Add data to the statistics. It is thread-safe.+addStatistics :: (MArray IOUArray a IO) => Statistics a -> a -> IO ()+addStatistics s x = + withMVar (statLock s) $ \() ->+ appendVector (statData s) x++-- | Return the statistics data. It is thread-safe.+statisticsData :: (MArray IOUArray a IO) => Statistics a -> IO (Array Int a)+statisticsData s =+ withMVar (statLock s) $ \() -> freezeVector (statData s)+ +-- | Represents the results of the statistic analysis.+data AnalysisResults a = + AnalysisResults { resultsData :: Array Int a,+ -- ^ Statistic data.+ resultsMean :: Double,+ -- ^ The average value.+ resultsVariance :: Double,+ -- ^ The variance.+ resultsMin :: a,+ -- ^ The minimum value.+ resultsMax :: a + -- ^ The maximum value.+ } deriving (Eq, Ord, Show)++-- | Analyze data.+analyzeData :: Real a => Array Int a -> AnalysisResults a+analyzeData xs =+ let (i1, i2) = bounds xs+ meanx = foldl' (\y i -> y * (1 - k i) + f i * k i) 0 [i1 .. i2]+ sqrx = foldl' (\y i -> y * (1 - k i) + g i * k i) 0 [i1 .. i2]+ minx = foldl' (\y i -> if i == 0 then x i else min y (x i)) 0 [i1 .. i2]+ maxx = foldl' (\y i -> if i == 0 then x i else max y (x i)) 0 [i1 .. i2]+ x i = xs ! i+ f i = fromRational (toRational (x i))+ g i = let y = f i in y * y+ k i = 1 / fromInteger (toInteger (i - i1 + 1))+ in AnalysisResults { resultsData = xs,+ resultsMean = meanx,+ resultsVariance = sqrx - meanx * meanx,+ resultsMin = minx,+ resultsMax = maxx }+ +-- | Show the results of analysis with the specified indent. +showResults :: (Show a) => AnalysisResults a -> Int -> ShowS+showResults rs indent =+ let (i1, i2) = bounds (resultsData rs)+ tab = replicate indent ' '+ in if i1 <= i2+ then+ showString tab .+ showString "mean = " . shows (resultsMean rs) . + showString "\n" . + showString tab .+ showString "deviation = " . shows (sqrt (resultsVariance rs)) . + showString "\n" .+ showString tab .+ showString "minimum = " . shows (resultsMin rs) . + showString "\n" .+ showString tab .+ showString "maximum = " . shows (resultsMax rs)+ else+ showString tab .+ showString "mean = ---" .+ showString "\n" . + showString tab .+ showString "deviation = ---" .+ showString "\n" . + showString tab .+ showString "minimum = ---" .+ showString "\n" . + showString tab .+ showString "maximum = ---"
+ Simulation/Aivika/UVector.hs view
@@ -0,0 +1,130 @@++{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module : Simulation.Aivika.UVector+-- 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 7.0.3+--+-- An imperative unboxed vector.+--+module Simulation.Aivika.UVector+ (UVector, + newVector, + copyVector, + vectorCount, + appendVector, + readVector, + writeVector, + vectorBinarySearch,+ freezeVector) where ++import Data.Array+import Data.Array.MArray+import Data.Array.IO+import Data.IORef+import Control.Monad++-- | Represents an unboxed resizable vector.+data UVector a = UVector { vectorArrayRef :: IORef (IOUArray Int a),+ vectorCountRef :: IORef Int, + vectorCapacityRef :: IORef Int }++-- | Create a new vector.+newVector :: MArray IOUArray a IO => IO (UVector a)+newVector = + do array <- newArray_ (0, 4 - 1)+ arrayRef <- newIORef array+ countRef <- newIORef 0+ capacityRef <- newIORef 4+ return UVector { vectorArrayRef = arrayRef,+ vectorCountRef = countRef,+ vectorCapacityRef = capacityRef }++-- | Copy the vector.+copyVector :: (MArray IOUArray a IO) => UVector a -> IO (UVector a)+copyVector vector =+ do array <- readIORef (vectorArrayRef vector)+ count <- readIORef (vectorCountRef vector)+ array' <- newArray_ (0, count - 1)+ arrayRef' <- newIORef array'+ countRef' <- newIORef count+ capacityRef' <- newIORef count+ forM_ [0 .. count - 1] $ \i ->+ do x <- readArray array i+ writeArray array' i x+ return UVector { vectorArrayRef = arrayRef',+ vectorCountRef = countRef',+ vectorCapacityRef = capacityRef' }++-- | Ensure that the vector has the specified capacity.+vectorEnsureCapacity :: MArray IOUArray a IO => UVector a -> Int -> IO ()+vectorEnsureCapacity vector capacity =+ do capacity' <- readIORef (vectorCapacityRef vector)+ when (capacity' < capacity) $+ do array' <- readIORef (vectorArrayRef vector)+ count' <- readIORef (vectorCountRef vector)+ let capacity'' = max (2 * capacity') capacity+ array'' <- newArray_ (0, capacity'' - 1)+ forM_ [0 .. count' - 1] $ \i ->+ do x <- readArray array' i+ writeArray array'' i x+ writeIORef (vectorArrayRef vector) array''+ writeIORef (vectorCapacityRef vector) capacity''+ +-- | Return the element count.+vectorCount :: MArray IOUArray a IO => UVector a -> IO Int+vectorCount vector = readIORef (vectorCountRef vector)+ +-- | Add the specified element to the end of the vector.+appendVector :: MArray IOUArray a IO => UVector a -> a -> IO () +appendVector vector item =+ do count <- readIORef (vectorCountRef vector)+ vectorEnsureCapacity vector (count + 1)+ array <- readIORef (vectorArrayRef vector)+ writeArray array count item+ writeIORef (vectorCountRef vector) (count + 1)+ +-- | Read a value from the vector, where indices are started from 0.+readVector :: MArray IOUArray a IO => UVector a -> Int -> IO a+readVector vector index =+ do array <- readIORef (vectorArrayRef vector)+ readArray array index+ +-- | Set an array item at the specified index which is started from 0.+writeVector :: MArray IOUArray a IO => UVector a -> Int -> a -> IO ()+writeVector vector index item =+ do array <- readIORef (vectorArrayRef vector)+ writeArray array index item+ +vectorBinarySearch' :: (MArray IOUArray a IO, Ord a) => + IOUArray Int a -> a -> Int -> Int -> IO Int+vectorBinarySearch' array item left right =+ if left > right + then return $ - (right + 1) - 1+ else+ do let index = (left + right) `div` 2+ curr <- readArray array index+ if item < curr + then vectorBinarySearch' array item left (index - 1)+ else if item == curr+ then return index+ else vectorBinarySearch' array item (index + 1) right+ +-- | Return the index of the specified element using binary search; otherwise, +-- a negated insertion index minus one: 0 -> -0 - 1, ..., i -> -i - 1, ....+vectorBinarySearch :: (MArray IOUArray a IO, Ord a) => UVector a -> a -> IO Int+vectorBinarySearch vector item =+ do array <- readIORef (vectorArrayRef vector)+ count <- readIORef (vectorCountRef vector)+ vectorBinarySearch' array item 0 (count - 1)++freezeVector :: (MArray IOUArray a IO) => UVector a -> IO (Array Int a)+freezeVector vector = + do vector' <- copyVector vector+ array <- readIORef (vectorArrayRef vector')+ freeze array+
+ Simulation/Aivika/Vector.hs view
@@ -0,0 +1,127 @@++-- |+-- Module : Simulation.Aivika.Vector+-- 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 7.0.3+--+-- An imperative vector.+--+module Simulation.Aivika.Vector+ (Vector, + newVector, + copyVector,+ vectorCount, + appendVector, + readVector, + writeVector,+ vectorBinarySearch,+ freezeVector) where ++import Data.Array+import Data.Array.MArray+import Data.Array.IO+import Data.IORef+import Control.Monad++-- | Represents a resizable vector.+data Vector a = Vector { vectorArrayRef :: IORef (IOArray Int a),+ vectorCountRef :: IORef Int, + vectorCapacityRef :: IORef Int }++-- | Create a new vector.+newVector :: IO (Vector a)+newVector = + do array <- newArray_ (0, 4 - 1)+ arrayRef <- newIORef array+ countRef <- newIORef 0+ capacityRef <- newIORef 4+ return Vector { vectorArrayRef = arrayRef,+ vectorCountRef = countRef,+ vectorCapacityRef = capacityRef }++-- | Copy the vector.+copyVector :: Vector a -> IO (Vector a)+copyVector vector =+ do array <- readIORef (vectorArrayRef vector)+ count <- readIORef (vectorCountRef vector)+ array' <- newArray_ (0, count - 1)+ arrayRef' <- newIORef array'+ countRef' <- newIORef count+ capacityRef' <- newIORef count+ forM_ [0 .. count - 1] $ \i ->+ do x <- readArray array i+ writeArray array' i x+ return Vector { vectorArrayRef = arrayRef',+ vectorCountRef = countRef',+ vectorCapacityRef = capacityRef' }++-- | Ensure that the vector has the specified capacity.+vectorEnsureCapacity :: Vector a -> Int -> IO ()+vectorEnsureCapacity vector capacity =+ do capacity' <- readIORef (vectorCapacityRef vector)+ when (capacity' < capacity) $+ do array' <- readIORef (vectorArrayRef vector)+ count' <- readIORef (vectorCountRef vector)+ let capacity'' = max (2 * capacity') capacity+ array'' <- newArray_ (0, capacity'' - 1)+ forM_ [0 .. count' - 1] $ \i ->+ do x <- readArray array' i+ writeArray array'' i x+ writeIORef (vectorArrayRef vector) array''+ writeIORef (vectorCapacityRef vector) capacity''+ +-- | Return the element count.+vectorCount :: Vector a -> IO Int+vectorCount vector = readIORef (vectorCountRef vector)+ +-- | Add the specified element to the end of the vector.+appendVector :: Vector a -> a -> IO () +appendVector vector item =+ do count <- readIORef (vectorCountRef vector)+ vectorEnsureCapacity vector (count + 1)+ array <- readIORef (vectorArrayRef vector)+ writeArray array count item+ writeIORef (vectorCountRef vector) (count + 1)+ +-- | Read a value from the vector, where indices are started from 0.+readVector :: Vector a -> Int -> IO a+readVector vector index =+ do array <- readIORef (vectorArrayRef vector)+ readArray array index+ +-- | Set an array item at the specified index which is started from 0.+writeVector :: Vector a -> Int -> a -> IO ()+writeVector vector index item =+ do array <- readIORef (vectorArrayRef vector)+ writeArray array index item++vectorBinarySearch' :: Ord a => IOArray Int a -> a -> Int -> Int -> IO Int+vectorBinarySearch' array item left right =+ if left > right + then return $ - (right + 1) - 1+ else+ do let index = (left + right) `div` 2+ curr <- readArray array index+ if item < curr + then vectorBinarySearch' array item left (index - 1)+ else if item == curr+ then return index+ else vectorBinarySearch' array item (index + 1) right+ +-- | Return the index of the specified element using binary search; otherwise, +-- a negated insertion index minus one: 0 -> -0 - 1, ..., i -> -i - 1, ....+vectorBinarySearch :: Ord a => Vector a -> a -> IO Int+vectorBinarySearch vector item =+ do array <- readIORef (vectorArrayRef vector)+ count <- readIORef (vectorCountRef vector)+ vectorBinarySearch' array item 0 (count - 1)++freezeVector :: Vector a -> IO (Array Int a)+freezeVector vector = + do vector' <- copyVector vector+ array <- readIORef (vectorArrayRef vector')+ freeze array+
aivika.cabal view
@@ -1,28 +1,34 @@ name: aivika-version: 0.1+version: 0.2 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.+ Aivika is a small simulation library that covers many paradigms. + It allows integrating 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.+ a computation in the Dynamics monad. There is also the Process+ monad to represent the discontinuous processes which can suspend+ at any time and then resume later. Everything else is expressed through + these two monads, including the event handlers, agent handlers and even + integrals. .+ The PDF documentation is available at + <https://github.com/dsorokin/aivika/blob/master/doc/aivika.pdf>+ . 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>+homepage: http://github.com/dsorokin/aivika cabal-version: >= 1.2.0 build-type: Simple-tested-with: GHC == 6.12.1+tested-with: GHC == 7.0.3 extra-source-files: examples/BassDiffusion.hs examples/ChemicalReaction.hs@@ -32,16 +38,46 @@ examples/MachRep1TimeDriven.hs examples/MachRep2.hs examples/MachRep3.hs+ examples/Furnace.hs data-files: doc/aivika.pdf library+ exposed-modules: Simulation.Aivika.Dynamics- other-modules: Simulation.Aivika.Queue+ Simulation.Aivika.Dynamics.Agent+ Simulation.Aivika.Dynamics.Base+ Simulation.Aivika.Dynamics.Cont+ Simulation.Aivika.Dynamics.EventQueue+ Simulation.Aivika.Dynamics.Lift+ Simulation.Aivika.Dynamics.Process+ Simulation.Aivika.Dynamics.Random+ Simulation.Aivika.Dynamics.Ref+ Simulation.Aivika.Dynamics.Resource+ Simulation.Aivika.Dynamics.SystemDynamics+ Simulation.Aivika.Dynamics.UVar+ Simulation.Aivika.Dynamics.Var+ Simulation.Aivika.Dynamics.Parameter Simulation.Aivika.PriorityQueue+ Simulation.Aivika.Queue+ Simulation.Aivika.Statistics++ other-modules: Simulation.Aivika.Dynamics.Internal.Dynamics+ Simulation.Aivika.Dynamics.Internal.Cont+ Simulation.Aivika.Dynamics.Internal.Process+ Simulation.Aivika.Dynamics.Internal.Time+ Simulation.Aivika.Dynamics.Internal.Memo+ Simulation.Aivika.Dynamics.Internal.Interpolate+ Simulation.Aivika.Dynamics.Internal.Fold+ Simulation.Aivika.Vector+ Simulation.Aivika.UVector - build-depends: base >= 3 && < 5,+ build-depends: base >= 3 && < 6,+ haskell98, mtl >= 1.1.0.2,- array >= 0.3.0.0+ array >= 0.3.0.0,+ containers >= 0.4.0.0++ extensions: FlexibleContexts ghc-options: -O2
doc/aivika.pdf view
binary file changed (275619 → 304809 bytes)
examples/BassDiffusion.hs view
@@ -5,6 +5,9 @@ import Control.Monad.Trans import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Agent+import Simulation.Aivika.Dynamics.Ref n = 500 -- the number of agents @@ -31,7 +34,7 @@ personPotentialAdopter :: AgentState, personAdopter :: AgentState } -createPerson :: DynamicsQueue -> Dynamics Person +createPerson :: EventQueue -> Dynamics Person createPerson q = do agent <- newAgent q potentialAdopter <- newState agent@@ -40,29 +43,27 @@ personPotentialAdopter = potentialAdopter, personAdopter = adopter } -createPersons :: DynamicsQueue -> Dynamics (Array Int Person)+createPersons :: EventQueue -> 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 :: Person -> Array Int Person -> Ref Int -> Ref Int -> Dynamics () definePerson p ps potentialAdopters adopters = do stateActivation (personPotentialAdopter p) $- do modifyRef' potentialAdopters $ \a -> a + 1+ 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+ 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 $+ addTimer (personAdopter p) t $ do i <- liftIO $ getStdRandom $ randomR (1, n) let p' = ps ! i st <- agentState (personAgent p')@@ -70,14 +71,11 @@ do b <- liftIO $ boolrnd adoptionFraction when b $ activateState (personAdopter p') stateDeactivation (personPotentialAdopter p) $- modifyRef' potentialAdopters $ \a -> a - 1+ modifyRef potentialAdopters $ \a -> a - 1 stateDeactivation (personAdopter p) $- modifyRef' adopters $ \a -> a - 1+ modifyRef adopters $ \a -> a - 1 -definePersons :: Array Int Person - -> DynamicsRef Int - -> DynamicsRef Int - -> Dynamics ()+definePersons :: Array Int Person -> Ref Int -> Ref Int -> Dynamics () definePersons ps potentialAdopters adopters = forM_ (elems ps) $ \p -> definePerson p ps potentialAdopters adopters@@ -102,5 +100,4 @@ return [i1, i2] main =- do xs <- runDynamics model specs- print xs+ printDynamics model specs
examples/ChemicalReaction.hs view
@@ -1,5 +1,6 @@ import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.SystemDynamics specs = Specs { spcStartTime = 0, spcStopTime = 13,
examples/FishBank.hs view
@@ -2,6 +2,7 @@ import Data.Array import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.SystemDynamics specs = Specs { spcStartTime = 0, spcStopTime = 13,
+ examples/Furnace.hs view
@@ -0,0 +1,417 @@++import Data.Maybe+import Random+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.Lift+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Ref+import Simulation.Aivika.Dynamics.UVar+import Simulation.Aivika.Dynamics.Process+import Simulation.Aivika.Dynamics.Random+import Simulation.Aivika.Statistics++import qualified Simulation.Aivika.Queue as Q++-- | The simulation specs.+specs = Specs { spcStartTime = 0.0,+ -- spcStopTime = 1000.0,+ spcStopTime = 300.0,+ spcDT = 0.1,+ spcMethod = RungeKutta4 }+ +-- | Return an exponentially distributed random value with mean +-- 1 / @lambda@, where @lambda@ is a parameter of the function.+exprnd :: Double -> IO Double+exprnd lambda =+ do x <- getStdRandom random+ return (- log x / lambda)+ +-- | Return a random initial temperature of the item. +temprnd :: IO Double+temprnd =+ do x <- getStdRandom random+ return (400.0 + (600.0 - 400.0) * x)++-- | Represents the furnace.+data Furnace = + Furnace { furnaceQueue :: EventQueue,+ -- ^ The event queue.+ furnaceNormalGen :: IO Double,+ -- ^ The normal random number generator.+ furnacePits :: [Pit],+ -- ^ The pits for ingots.+ furnacePitCount :: UVar Int,+ -- ^ The count of active pits with ingots.+ furnacePitCountStats :: Statistics Int,+ -- ^ The statistics about the active pits.+ furnaceAwaitingIngots :: Q.Queue Ingot,+ -- ^ The awaiting ingots in the queue.+ furnaceQueueCount :: UVar Int,+ -- ^ The queue count.+ furnaceQueueCountStats :: Statistics Int,+ -- ^ The statistics about the queue count.+ furnaceWaitCount :: Ref Int,+ -- ^ The count of awaiting ingots.+ furnaceWaitTime :: Ref Double,+ -- ^ The wait time for all loaded ingots.+ furnaceHeatingTime :: Ref Double,+ -- ^ The heating time for all unloaded ingots.+ furnaceTemp :: Ref Double,+ -- ^ The furnace temperature.+ furnaceTotalCount :: Ref Int,+ -- ^ The total count of ingots.+ furnaceLoadCount :: Ref Int,+ -- ^ The count of loaded ingots.+ furnaceUnloadCount :: Ref Int,+ -- ^ The count of unloaded ingots.+ furnaceUnloadTemps :: Ref [Double]+ -- ^ The temperatures of all unloaded ingots.+ }++-- | A pit in the furnace to place the ingots.+data Pit = + Pit { pitQueue :: EventQueue,+ -- ^ The bound dynamics queue.+ pitIngot :: Ref (Maybe Ingot),+ -- ^ The ingot in the pit.+ pitTemp :: Ref Double+ -- ^ The ingot temperature in the pit.+ }++data Ingot = + Ingot { ingotFurnace :: Furnace,+ -- ^ Return the furnace.+ ingotReceiveTime :: Double,+ -- ^ The time at which the ingot was received.+ ingotReceiveTemp :: Double,+ -- ^ The temperature with which the ingot was received.+ ingotLoadTime :: Double,+ -- ^ The time of loading in the furnace.+ ingotLoadTemp :: Double,+ -- ^ The temperature when the ingot was loaded in the furnace.+ ingotCoeff :: Double+ -- ^ The heating coefficient.+ }++-- | Create a furnace.+newFurnace :: EventQueue -> Dynamics Furnace+newFurnace queue =+ do normalGen <- liftIO normalGen+ pits <- sequence [newPit queue | i <- [1..10]]+ pitCount <- newUVar queue 0+ pitCountStats <- liftIO newStatistics+ awaitingIngots <- liftIO Q.newQueue+ queueCount <- newUVar queue 0+ queueCountStats <- liftIO newStatistics+ waitCount <- newRef queue 0+ waitTime <- newRef queue 0.0+ heatingTime <- newRef queue 0.0+ h <- newRef queue 1650.0+ totalCount <- newRef queue 0+ loadCount <- newRef queue 0+ unloadCount <- newRef queue 0+ unloadTemps <- newRef queue []+ return Furnace { furnaceQueue = queue,+ furnaceNormalGen = normalGen,+ furnacePits = pits,+ furnacePitCount = pitCount,+ furnacePitCountStats = pitCountStats,+ furnaceAwaitingIngots = awaitingIngots,+ furnaceQueueCount = queueCount,+ furnaceQueueCountStats = queueCountStats,+ furnaceWaitCount = waitCount,+ furnaceWaitTime = waitTime,+ furnaceHeatingTime = heatingTime,+ furnaceTemp = h,+ furnaceTotalCount = totalCount,+ furnaceLoadCount = loadCount, + furnaceUnloadCount = unloadCount, + furnaceUnloadTemps = unloadTemps }++-- | Create a new pit.+newPit :: EventQueue -> Dynamics Pit+newPit queue =+ do ingot <- newRef queue Nothing+ h' <- newRef queue 0.0+ return Pit { pitQueue = queue,+ pitIngot = ingot,+ pitTemp = h' }++-- | Create a new ingot.+newIngot :: Furnace -> Dynamics Ingot+newIngot furnace =+ do t <- time+ xi <- liftIO $ furnaceNormalGen furnace+ h' <- liftIO temprnd+ let c = 0.1 + (0.05 + xi * 0.01)+ return Ingot { ingotFurnace = furnace,+ ingotReceiveTime = t,+ ingotReceiveTemp = h',+ ingotLoadTime = t,+ ingotLoadTemp = h',+ ingotCoeff = c }++-- | Heat the ingot up in the pit if there is such an ingot.+heatPitUp :: Pit -> Dynamics ()+heatPitUp pit =+ do ingot <- readRef (pitIngot pit)+ case ingot of+ Nothing -> + return ()+ Just ingot -> do+ + -- update the temperature of the ingot.+ let furnace = ingotFurnace ingot+ dt' <- dt+ h' <- readRef (pitTemp pit)+ h <- readRef (furnaceTemp furnace)+ writeRef (pitTemp pit) $ + h' + dt' * (h - h') * ingotCoeff ingot++-- | Check whether there are ready ingots in the pits.+ingotsReady :: Furnace -> Dynamics Bool+ingotsReady furnace =+ fmap (not . null) $ + filterM (fmap (>= 2200.0) . readRef . pitTemp) $ + furnacePits furnace++-- | Try to unload the ready ingot from the specified pit.+tryUnloadPit :: Furnace -> Pit -> Dynamics ()+tryUnloadPit furnace pit =+ do h' <- readRef (pitTemp pit)+ when (h' >= 2000.0) $+ do Just ingot <- readRef (pitIngot pit) + unloadIngot ingot pit++-- | Try to load an awaiting ingot in the specified empty pit.+tryLoadPit :: Furnace -> Pit -> Dynamics () +tryLoadPit furnace pit =+ do let ingots = furnaceAwaitingIngots furnace+ flag <- liftIO $ Q.queueNull ingots+ unless flag $+ do ingot <- liftIO $ Q.queueFront ingots+ liftIO $ Q.dequeue ingots+ t' <- time+ modifyUVar (furnaceQueueCount furnace) (+ (-1))+ c <- readUVar (furnaceQueueCount furnace)+ liftIO $ addStatistics (furnaceQueueCountStats furnace) c+ loadIngot (ingot { ingotLoadTime = t',+ ingotLoadTemp = 400.0 }) pit+ +-- | Unload the ingot from the specified pit. +unloadIngot :: Ingot -> Pit -> Dynamics ()+unloadIngot ingot pit = + do h' <- readRef (pitTemp pit)+ writeRef (pitIngot pit) Nothing+ writeRef (pitTemp pit) 0.0+ + -- count the active pits+ let furnace = ingotFurnace ingot+ count <- readUVar (furnacePitCount furnace)+ writeUVar (furnacePitCount furnace) (count - 1)+ liftIO $ addStatistics (furnacePitCountStats furnace) (count - 1)+ + -- how long did we heat the ingot up?+ t' <- time+ modifyRef (furnaceHeatingTime furnace)+ (+ (t' - ingotLoadTime ingot))+ + -- what is the temperature of the unloaded ingot?+ modifyRef (furnaceUnloadTemps furnace) (h' :)+ + -- count the unloaded ingots+ modifyRef (furnaceUnloadCount furnace) (+ 1)+ +-- | Load the ingot in the specified pit+loadIngot :: Ingot -> Pit -> Dynamics ()+loadIngot ingot pit =+ do writeRef (pitIngot pit) $ Just ingot+ writeRef (pitTemp pit) $ ingotLoadTemp ingot+ + -- count the active pits+ let furnace = ingotFurnace ingot+ count <- readUVar (furnacePitCount furnace)+ writeUVar (furnacePitCount furnace) (count + 1)+ liftIO $ addStatistics (furnacePitCountStats furnace) (count + 1)+ + -- decrease the furnace temperature+ h <- readRef (furnaceTemp furnace)+ let h' = ingotLoadTemp ingot+ dh = - (h - h') / fromInteger (toInteger (count + 1))+ writeRef (furnaceTemp furnace) $ h + dh++ -- how long did we keep the ingot in the queue?+ t' <- time+ when (ingotReceiveTime ingot < t') $+ do modifyRef (furnaceWaitCount furnace) (+ 1) + modifyRef (furnaceWaitTime furnace)+ (+ (t' - ingotReceiveTime ingot))++ -- count the loaded ingots+ modifyRef (furnaceLoadCount furnace) (+ 1)+ +-- | Iterate the furnace processing.+iterateFurnace :: Furnace -> Dynamics (Dynamics ())+iterateFurnace furnace = + let pits = furnacePits furnace+ in iterateD $+ do ready <- ingotsReady furnace+ when ready $ + do mapM_ (tryUnloadPit furnace) pits+ pits' <- emptyPits furnace+ mapM_ (tryLoadPit furnace) pits'+ mapM_ heatPitUp pits+ + -- update the temperature of the furnace+ dt' <- dt+ h <- readRef (furnaceTemp furnace)+ writeRef (furnaceTemp furnace) $+ h + dt' * (2600.0 - h) * 0.2++-- | Return all empty pits.+emptyPits :: Furnace -> Dynamics [Pit]+emptyPits furnace =+ filterM (fmap isNothing . readRef . pitIngot) $+ furnacePits furnace++-- | Accept a new ingot.+acceptIngot :: Furnace -> Dynamics ()+acceptIngot furnace =+ do ingot <- newIngot furnace+ + -- counting+ modifyRef (furnaceTotalCount furnace) (+ 1)+ + -- check what to do with the new ingot+ count <- readUVar (furnacePitCount furnace)+ if count >= 10+ then do let ingots = furnaceAwaitingIngots furnace+ liftIO $ Q.enqueue ingots ingot+ modifyUVar (furnaceQueueCount furnace) (+ 1)+ c <- readUVar (furnaceQueueCount furnace)+ liftIO $ addStatistics (furnaceQueueCountStats furnace) c+ else do pit:_ <- emptyPits furnace+ loadIngot ingot pit+ +-- | Process the furnace.+processFurnace :: Furnace -> Process ()+processFurnace furnace =+ do delay <- liftIO $ exprnd (1.0 / 2.5)+ holdProcess delay+ -- we have got a new ingot+ liftD $ acceptIngot furnace+ -- repeat it again+ processFurnace furnace++-- | Initialize the furnace.+initializeFurnace :: Furnace -> Dynamics ()+initializeFurnace furnace =+ do x1 <- newIngot furnace+ x2 <- newIngot furnace+ x3 <- newIngot furnace+ x4 <- newIngot furnace+ x5 <- newIngot furnace+ x6 <- newIngot furnace+ let p1 : p2 : p3 : p4 : p5 : p6 : ps = + furnacePits furnace+ loadIngot (x1 { ingotLoadTemp = 550.0 }) p1+ loadIngot (x2 { ingotLoadTemp = 600.0 }) p2+ loadIngot (x3 { ingotLoadTemp = 650.0 }) p3+ loadIngot (x4 { ingotLoadTemp = 700.0 }) p4+ loadIngot (x5 { ingotLoadTemp = 750.0 }) p5+ loadIngot (x6 { ingotLoadTemp = 800.0 }) p6+ writeRef (furnaceTotalCount furnace) 6+ writeRef (furnaceTemp furnace) 1650.0+ +-- | Return a count, average and deviation.+stats :: [Double] -> (Int, Double, Double)+stats xs = (length xs, ex, sx)+ where+ n = fromInteger $ toInteger $ length xs+ ex = sum xs / n+ dx = (sum . map rho) xs / (n - 1.0)+ sx = sqrt dx+ rho x = (x - ex) ^ 2++-- | The simulation model.+model :: Dynamics (Dynamics ())+model =+ do queue <- newQueue+ furnace <- newFurnace queue+ pid <- newProcessID queue+ + initializeFurnace furnace+ + -- get the furnace iterator+ iterator <- iterateFurnace furnace+ + -- accept input ingots+ t0 <- starttime+ runProcess (processFurnace furnace) pid t0+ + let system :: Dynamics ()+ system = + do iterator -- iterate in each time point+ + -- the ingots+ c0 <- readRef (furnaceTotalCount furnace)+ c1 <- readRef (furnaceLoadCount furnace)+ c2 <- readRef (furnaceUnloadCount furnace)+ c3 <- readRef (furnaceWaitCount furnace)+ + liftIO $ do+ putStrLn "The count of ingots:"+ putStrLn $ " total = " ++ show c0+ putStrLn $ " loaded = " ++ show c1+ putStrLn $ " ready = " ++ show c2+ putStrLn $ " awaited in the queue = " ++ show c3+ putStrLn ""+ + -- the temperature of the ready ingots+ (n1, e1, d1) <- + fmap stats $ readRef (furnaceUnloadTemps furnace)+ + liftIO $ do + putStrLn "The temperature of the ready ingots:"+ putStrLn $ " average = " ++ show e1+ putStrLn $ " deviation = " ++ show d1+ putStrLn ""+ + -- the ingots in pits+ r2 <- fmap analyzeData $ liftIO $ statisticsData (furnacePitCountStats furnace)+ + liftIO $ do+ putStrLn "The ingots in pits: "+ putStrLn $ showResults r2 2 []+ putStrLn ""+ + -- the queue size+ r3 <- fmap analyzeData $ liftIO $ statisticsData (furnaceQueueCountStats furnace)+ + liftIO $ do+ putStrLn "The queue size: "+ putStrLn $ showResults r3 2 []+ putStrLn ""+ + -- the mean wait time in the queue+ t4 <- readRef (furnaceWaitTime furnace) /+ fmap (fromInteger . toInteger)+ (readRef (furnaceWaitCount furnace))+ + -- the mean heating time+ t5 <- readRef (furnaceHeatingTime furnace) /+ fmap (fromInteger . toInteger)+ (readRef (furnaceUnloadCount furnace))+ + liftIO $ do+ putStrLn $ "The mean wait time: " ++ show t4+ putStrLn $ "The mean heating time: " ++ show t5+ + return system++-- | The main program.+main = runDynamics1 model specs
examples/MachRep1.hs view
@@ -19,6 +19,11 @@ import Control.Monad.Trans import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.Lift+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Ref+import Simulation.Aivika.Dynamics.Process upRate = 1.0 / 1.0 -- reciprocal of mean up time repairRate = 1.0 / 0.5 -- reciprocal of mean repair time@@ -38,23 +43,25 @@ do queue <- newQueue totalUpTime <- newRef queue 0.0 - pid1 <- newPID queue- pid2 <- newPID queue+ pid1 <- newProcessID queue+ pid2 <- newProcessID queue - let machine :: DynamicsProc ()+ let machine :: Process () machine = do startUpTime <- liftD time upTime <- liftIO $ exprnd upRate- holdProc upTime+ holdProcess upTime finishUpTime <- liftD time- liftD $ modifyRef' totalUpTime+ liftD $ modifyRef totalUpTime (+ (finishUpTime - startUpTime)) repairTime <- liftIO $ exprnd repairRate- holdProc repairTime+ holdProcess repairTime machine - runProc machine pid1 starttime- runProc machine pid2 starttime+ t0 <- starttime+ + runProcess machine pid1 t0+ runProcess machine pid2 t0 let system :: Dynamics Double system =
examples/MachRep1EventDriven.hs view
@@ -19,6 +19,9 @@ import Control.Monad.Trans import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Ref upRate = 1.0 / 1.0 -- reciprocal of mean up time repairRate = 1.0 / 0.5 -- reciprocal of mean repair time@@ -42,12 +45,12 @@ machineBroken startUpTime = do finishUpTime <- time- modifyRef' totalUpTime (+ (finishUpTime - startUpTime))+ modifyRef totalUpTime (+ (finishUpTime - startUpTime)) repairTime <- liftIO $ exprnd repairRate -- enqueue a new event- let t = return $ finishUpTime + repairTime- enqueueD queue t machineRepaired+ let t = finishUpTime + repairTime+ enqueue queue t machineRepaired machineRepaired :: Dynamics () machineRepaired =@@ -56,11 +59,13 @@ upTime <- liftIO $ exprnd upRate -- enqueue a new event- let t = return $ startUpTime + upTime- enqueueD queue t $ machineBroken startUpTime+ let t = startUpTime + upTime+ enqueue queue t $ machineBroken startUpTime - enqueueD queue starttime machineRepaired -- start the first machine- enqueueD queue starttime machineRepaired -- start the second machine+ t0 <- starttime+ + enqueue queue t0 machineRepaired -- start the first machine+ enqueue queue t0 machineRepaired -- start the second machine let system :: Dynamics Double system =
examples/MachRep1TimeDriven.hs view
@@ -19,6 +19,9 @@ import Control.Monad.Trans import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Ref upRate = 1.0 / 1.0 -- reciprocal of mean up time repairRate = 1.0 / 0.5 -- reciprocal of mean repair time@@ -56,34 +59,34 @@ repairNum' <- readRef repairNum let untilBroken = - modifyRef' upNum $ \a -> a - 1+ modifyRef upNum $ \a -> a - 1 untilRepaired =- modifyRef' repairNum $ \a -> a - 1+ modifyRef repairNum $ \a -> a - 1 broken =- do writeRef' upNum (-1)+ do writeRef upNum (-1) -- the machine is broken startUpTime' <- readRef startUpTime finishUpTime' <- time dt' <- dt- modifyRef' totalUpTime $ + modifyRef totalUpTime $ \a -> a + (finishUpTime' - startUpTime') repairTime' <- liftIO $ exprnd repairRate- writeRef' repairNum $+ writeRef repairNum $ round (repairTime' / dt') repaired =- do writeRef' repairNum (-1)+ do writeRef repairNum (-1) -- the machine is repaired t' <- time dt' <- dt- writeRef' startUpTime t'+ writeRef startUpTime t' upTime' <- liftIO $ exprnd upRate- writeRef' upNum $+ writeRef upNum $ round (upTime' / dt') result | upNum' > 0 = untilBroken@@ -98,8 +101,8 @@ m2 <- machine -- create strictly sequential computations- c1 <- memo0 discrete m1- c2 <- memo0 discrete m2+ c1 <- iterateD m1+ c2 <- iterateD m2 let system :: Dynamics Double system =
examples/MachRep2.hs view
@@ -22,6 +22,12 @@ import Control.Monad.Trans import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.Lift+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Ref+import Simulation.Aivika.Dynamics.Resource+import Simulation.Aivika.Dynamics.Process upRate = 1.0 / 1.0 -- reciprocal of mean up time repairRate = 1.0 / 0.5 -- reciprocal of mean repair time@@ -52,33 +58,35 @@ repairPerson <- newResource queue 1 - pid1 <- newPID queue- pid2 <- newPID queue+ pid1 <- newProcessID queue+ pid2 <- newProcessID queue - let machine :: DynamicsProc ()+ let machine :: Process () machine = do startUpTime <- liftD time upTime <- liftIO $ exprnd upRate- holdProc upTime+ holdProcess upTime finishUpTime <- liftD time- liftD $ modifyRef' totalUpTime + liftD $ modifyRef totalUpTime (+ (finishUpTime - startUpTime)) -- check the resource availability- liftD $ modifyRef' nRep (+ 1)+ liftD $ modifyRef nRep (+ 1) n <- resourceCount repairPerson when (n == 1) $- liftD $ modifyRef' nImmedRep (+ 1)+ liftD $ modifyRef nImmedRep (+ 1) requestResource repairPerson repairTime <- liftIO $ exprnd repairRate- holdProc repairTime+ holdProcess repairTime releaseResource repairPerson machine - runProc machine pid1 starttime- runProc machine pid2 starttime+ t0 <- starttime+ + runProcess machine pid1 t0+ runProcess machine pid2 t0 let system :: Dynamics (Double, Double) system =
examples/MachRep3.hs view
@@ -18,6 +18,12 @@ import Control.Monad.Trans import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.Lift+import Simulation.Aivika.Dynamics.EventQueue+import Simulation.Aivika.Dynamics.Ref+import Simulation.Aivika.Dynamics.Resource+import Simulation.Aivika.Dynamics.Process upRate = 1.0 / 1.0 -- reciprocal of mean up time repairRate = 1.0 / 0.5 -- reciprocal of mean repair time@@ -44,35 +50,37 @@ repairPerson <- newResource queue 1 - pid1 <- newPID queue- pid2 <- newPID queue+ pid1 <- newProcessID queue+ pid2 <- newProcessID queue - let machine :: DynamicsPID -> DynamicsProc ()+ let machine :: ProcessID -> Process () machine pid = do startUpTime <- liftD time upTime <- liftIO $ exprnd upRate- holdProc upTime+ holdProcess upTime finishUpTime <- liftD time- liftD $ modifyRef' totalUpTime + liftD $ modifyRef totalUpTime (+ (finishUpTime - startUpTime)) - liftD $ modifyRef' nUp $ \a -> a - 1+ liftD $ modifyRef nUp $ \a -> a - 1 nUp' <- liftD $ readRef nUp if nUp' == 1- then passivateProc+ then passivateProcess else do n <- resourceCount repairPerson- when (n == 1) $ reactivateProc pid+ when (n == 1) $ reactivateProcess pid requestResource repairPerson repairTime <- liftIO $ exprnd repairRate- holdProc repairTime- liftD $ modifyRef' nUp $ \a -> a + 1+ holdProcess repairTime+ liftD $ modifyRef nUp $ \a -> a + 1 releaseResource repairPerson machine pid - runProc (machine pid2) pid1 starttime- runProc (machine pid1) pid2 starttime+ t0 <- starttime+ + runProcess (machine pid2) pid1 t0+ runProcess (machine pid1) pid2 t0 let system :: Dynamics Double system =