diff --git a/Simulation/Aivika.hs b/Simulation/Aivika.hs --- a/Simulation/Aivika.hs +++ b/Simulation/Aivika.hs @@ -14,6 +14,7 @@ (-- * Modules module Simulation.Aivika.Agent, module Simulation.Aivika.Arrival, + module Simulation.Aivika.Circuit, module Simulation.Aivika.Cont, module Simulation.Aivika.Dynamics, module Simulation.Aivika.Dynamics.Interpolate, @@ -38,10 +39,12 @@ module Simulation.Aivika.Stream, module Simulation.Aivika.Stream.Random, module Simulation.Aivika.Task, + module Simulation.Aivika.Transform, module Simulation.Aivika.Var.Unboxed) where import Simulation.Aivika.Agent import Simulation.Aivika.Arrival +import Simulation.Aivika.Circuit import Simulation.Aivika.Cont import Simulation.Aivika.Dynamics import Simulation.Aivika.Dynamics.Interpolate @@ -66,4 +69,5 @@ import Simulation.Aivika.Stream import Simulation.Aivika.Stream.Random import Simulation.Aivika.Task +import Simulation.Aivika.Transform import Simulation.Aivika.Var.Unboxed diff --git a/Simulation/Aivika/Arrival.hs b/Simulation/Aivika/Arrival.hs --- a/Simulation/Aivika/Arrival.hs +++ b/Simulation/Aivika/Arrival.hs @@ -7,8 +7,11 @@ -- Stability : experimental -- Tested with: GHC 7.6.3 -- --- This module defines the types and functions for working with --- the external events that usually arrive from outside the model. +-- This module defines the types and functions for working with the events +-- that can represent something that arrive from outside the model, or +-- represent other things which computation is delayed and hence is not synchronized. +-- +-- Therefore, the additional information is provided about the time and delay of arrival. module Simulation.Aivika.Arrival (Arrival(..), @@ -27,22 +30,7 @@ import Simulation.Aivika.Stream import Simulation.Aivika.Statistics import Simulation.Aivika.Ref - --- | It defines when an external event has arrived, usually generated by --- some random 'Stream'. --- --- Such events should arrive one by one without time lag in the following sense --- that the model should start awaiting the next event exactly in that time --- when the previous event has arrived. It usually happens automatically when --- using a 'queueProcessor' with an ability to lost the arrived event if the queue --- is full. -data Arrival = - Arrival { arrivalTime :: Double, - -- ^ the simulation time at which the event has arrived - arrivalDelay :: Double - -- ^ the delay time which has passed from the time of - -- arriving the previous event - } deriving (Eq, Ord, Show) +import Simulation.Aivika.Internal.Arrival -- | Accumulates the statistics about that how long the arrived events are processed. data ArrivalTimer = @@ -60,7 +48,7 @@ -- | Return a processor that actually measures how much time has passed from -- the time of arriving the events. -arrivalTimerProcessor :: ArrivalTimer -> Processor Arrival Arrival +arrivalTimerProcessor :: ArrivalTimer -> Processor (Arrival a) (Arrival a) arrivalTimerProcessor timer = Processor $ \xs -> Cons $ loop xs where loop xs = diff --git a/Simulation/Aivika/Circuit.hs b/Simulation/Aivika/Circuit.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Circuit.hs @@ -0,0 +1,378 @@ + +{-# LANGUAGE RecursiveDo, Arrows #-} + +-- | +-- Module : Simulation.Aivika.Circuit +-- Copyright : Copyright (c) 2009-2014, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.6.3 +-- +-- It represents a circuit synchronized with the event queue. +-- Also it allows creating the recursive links with help of +-- the proc-notation. +-- +-- The implementation is based on the . +-- +module Simulation.Aivika.Circuit + (-- * Circuit Arrow + Circuit(..), + -- * Circuit Primitives + arrCircuit, + accumCircuit, + -- * Arrival Circuit + arrivalCircuit, + -- * Delaying Circuit + delayCircuit, + -- * Time Circuit + timeCircuit, + -- * Conditional Computation + (?>), + filterCircuit, + filterCircuitM, + neverCircuit, + -- * Converting to Signals and Processors + circuitSignaling, + circuitProcessor, + -- * Integrals and Difference Equations + integCircuit, + sumCircuit, + -- * Circuit Transform + circuitTransform) where + +import qualified Control.Category as C +import Control.Arrow +import Control.Monad.Fix + +import Data.IORef + +import Simulation.Aivika.Internal.Arrival +import Simulation.Aivika.Internal.Specs +import Simulation.Aivika.Internal.Simulation +import Simulation.Aivika.Internal.Dynamics +import Simulation.Aivika.Internal.Event +import Simulation.Aivika.Dynamics.Memo +import Simulation.Aivika.Transform +import Simulation.Aivika.SystemDynamics +import Simulation.Aivika.Signal +import Simulation.Aivika.Stream +import Simulation.Aivika.Processor + +-- | Represents a circuit synchronized with the event queue. +-- Besides, it allows creating the recursive links with help of +-- the proc-notation. +-- +newtype Circuit a b = + Circuit { runCircuit :: a -> Event (Circuit a b, b) + -- ^ Run the circuit. + } + +instance C.Category Circuit where + + id = Circuit $ \a -> return (C.id, a) + + (.) = dot + where + (Circuit g) `dot` (Circuit f) = + Circuit $ \a -> + Event $ \p -> + do (cir1, b) <- invokeEvent p (f a) + (cir2, c) <- invokeEvent p (g b) + return (cir2 `dot` cir1, c) + +instance Arrow Circuit where + + arr f = Circuit $ \a -> return (arr f, f a) + + first (Circuit f) = + Circuit $ \(b, d) -> + Event $ \p -> + do (cir, c) <- invokeEvent p (f b) + return (first cir, (c, d)) + + second (Circuit f) = + Circuit $ \(d, b) -> + Event $ \p -> + do (cir, c) <- invokeEvent p (f b) + return (second cir, (d, c)) + + (Circuit f) *** (Circuit g) = + Circuit $ \(b, b') -> + Event $ \p -> + do (cir1, c) <- invokeEvent p (f b) + (cir2, c') <- invokeEvent p (g b') + return (cir1 *** cir2, (c, c')) + + (Circuit f) &&& (Circuit g) = + Circuit $ \b -> + Event $ \p -> + do (cir1, c) <- invokeEvent p (f b) + (cir2, c') <- invokeEvent p (g b) + return (cir1 &&& cir2, (c, c')) + +instance ArrowLoop Circuit where + + loop (Circuit f) = + Circuit $ \b -> + Event $ \p -> + do rec (cir, (c, d)) <- invokeEvent p (f (b, d)) + return (loop cir, c) + +instance ArrowChoice Circuit where + + left x@(Circuit f) = + Circuit $ \ebd -> + Event $ \p -> + case ebd of + Left b -> + do (cir, c) <- invokeEvent p (f b) + return (left cir, Left c) + Right d -> + return (left x, Right d) + + right x@(Circuit f) = + Circuit $ \edb -> + Event $ \p -> + case edb of + Right b -> + do (cir, c) <- invokeEvent p (f b) + return (right cir, Right c) + Left d -> + return (right x, Left d) + + x@(Circuit f) +++ y@(Circuit g) = + Circuit $ \ebb' -> + Event $ \p -> + case ebb' of + Left b -> + do (cir1, c) <- invokeEvent p (f b) + return (cir1 +++ y, Left c) + Right b' -> + do (cir2, c') <- invokeEvent p (g b') + return (x +++ cir2, Right c') + + x@(Circuit f) ||| y@(Circuit g) = + Circuit $ \ebc -> + Event $ \p -> + case ebc of + Left b -> + do (cir1, d) <- invokeEvent p (f b) + return (cir1 ||| y, d) + Right b' -> + do (cir2, d) <- invokeEvent p (g b') + return (x ||| cir2, d) + +-- | Get a signal transform by the specified circuit. +circuitSignaling :: Circuit a b -> Signal a -> Signal b +circuitSignaling (Circuit cir) sa = + Signal { handleSignal = \f -> + Event $ \p -> + do r <- newIORef cir + invokeEvent p $ + handleSignal sa $ \a -> + Event $ \p -> + do cir <- readIORef r + (Circuit cir', b) <- invokeEvent p (cir a) + writeIORef r cir' + invokeEvent p (f b) } + +-- | Transform the circuit to a processor. +circuitProcessor :: Circuit a b -> Processor a b +circuitProcessor (Circuit cir) = Processor $ \sa -> + Cons $ + do (a, xs) <- runStream sa + (cir', b) <- liftEvent (cir a) + let f = runProcessor (circuitProcessor cir') + return (b, f xs) + +-- | Lift the 'Event' function to a curcuit. +arrCircuit :: (a -> Event b) -> Circuit a b +arrCircuit f = + let x = + Circuit $ \a -> + Event $ \p -> + do b <- invokeEvent p (f a) + return (x, b) + in x + +-- | Accumulator that outputs a value determined by the supplied function. +accumCircuit :: (acc -> a -> Event (acc, b)) -> acc -> Circuit a b +accumCircuit f acc = + Circuit $ \a -> + Event $ \p -> + do (acc', b) <- invokeEvent p (f acc a) + return (accumCircuit f acc', b) + +-- | A circuit that adds the information about the time points at which +-- the values were received. +arrivalCircuit :: Circuit a (Arrival a) +arrivalCircuit = + let loop t0 = + Circuit $ \a -> + Event $ \p -> + let t = pointTime p + b = Arrival { arrivalValue = a, + arrivalTime = t, + arrivalDelay = + case t0 of + Nothing -> 0 + Just t0 -> t - t0 } + in return (loop $ Just t, b) + in loop Nothing + +-- | Delay the input by one step using the specified initial value. +delayCircuit :: a -> Circuit a a +delayCircuit a0 = + Circuit $ \a -> + return (delayCircuit a, a0) + +-- | A circuit that returns the current modeling time. +timeCircuit :: Circuit a Double +timeCircuit = + Circuit $ \a -> + Event $ \p -> + return (timeCircuit, pointTime p) + +-- | Like '>>>' but processes only the represented events. +(>?>) :: Circuit a (Maybe b) + -- ^ whether there is an event + -> Circuit b c + -- ^ process the event if it presents + -> Circuit a (Maybe c) + -- ^ the resulting circuit that processes only the represented events +whether >?> process = + Circuit $ \a -> + Event $ \p -> + do (whether', b) <- invokeEvent p (runCircuit whether a) + case b of + Nothing -> + return (whether' >?> process, Nothing) + Just b -> + do (process', c) <- invokeEvent p (runCircuit process b) + return (whether' >?> process', Just c) + +-- | Like '<<<' but processes only the represented events. +( Circuit a (Maybe b) + -- ^ whether there is an event + -> Circuit a (Maybe c) + -- ^ the resulting circuit that processes only the represented events +(?>) + +-- | Filter the circuit, calculating only those parts of the circuit that satisfy +-- the specified predicate. +filterCircuit :: (a -> Bool) -> Circuit a b -> Circuit a (Maybe b) +filterCircuit pred = filterCircuitM (return . pred) + +-- | Filter the circuit within the 'Event' computation, calculating only those parts +-- of the circuit that satisfy the specified predicate. +filterCircuitM :: (a -> Event Bool) -> Circuit a b -> Circuit a (Maybe b) +filterCircuitM pred cir = + Circuit $ \a -> + Event $ \p -> + do x <- invokeEvent p (pred a) + if x + then do (cir', b) <- invokeEvent p (runCircuit cir a) + return (filterCircuitM pred cir', Just b) + else return (filterCircuitM pred cir, Nothing) + +-- | The source of events that never occur. +neverCircuit :: Circuit a (Maybe b) +neverCircuit = + Circuit $ \a -> return (neverCircuit, Nothing) + +-- | An approximation of the integral using Euler's method. +-- +-- This function can be rather inaccurate as it depends on +-- the time points at wich the 'Circuit' computation is actuated. +-- Also Euler's method per se is not most accurate, although simple +-- enough for implementation. +-- +-- Consider using the 'integ' function whenever possible. +-- That function can integrate with help of the Runge-Kutta method by +-- the specified integration time points that are passed in the simulation +-- specs to every 'Simulation', when running the model. +-- +-- At the same time, the 'integCircuit' function has no mutable state +-- unlike the former. The latter consumes less memory but at the cost +-- of inaccuracy and relatively more slow simulation, had we requested +-- the integral in the same time points. +-- +-- Regarding the recursive equations, the both functions allow defining them +-- but whithin different computations (either with help of the recursive +-- do-notation or the proc-notation). +integCircuit :: Double + -- ^ the initial value + -> Circuit Double Double + -- ^ map the derivative to an integral +integCircuit init = start + where + start = + Circuit $ \a -> + Event $ \p -> + do let t = pointTime p + return (next t init a, init) + next t0 v0 a0 = + Circuit $ \a -> + Event $ \p -> + do let t = pointTime p + dt = t - t0 + v = v0 + a0 * dt + v `seq` return (next t v a, v) + +-- | A sum of differences starting from the specified initial value. +-- +-- Consider using the more accurate 'diffsum' function whener possible as +-- it is calculated in every integration time point specified by specs +-- passed in to every 'Simulation', when running the model. +-- +-- At the same time, the 'sumCircuit' function has no mutable state and +-- it consumes less memory than the former. +-- +-- Regarding the recursive equations, the both functions allow defining them +-- but whithin different computations (either with help of the recursive +-- do-notation or the proc-notation). +sumCircuit :: Num a => + a + -- ^ the initial value + -> Circuit a a + -- ^ map the difference to a sum +sumCircuit init = start + where + start = + Circuit $ \a -> + Event $ \p -> + return (next init a, init) + next v0 a0 = + Circuit $ \a -> + Event $ \p -> + do let v = v0 + a0 + v `seq` return (next v a, v) + +-- | Approximate the circuit as a transform of time varying function, +-- calculating the values in the integration time points and then +-- interpolating in all other time points. The resulting transform +-- computation is synchronized with the event queue. +-- +-- This procedure consumes memory as the underlying memoization allocates +-- an array to store the calculated values. +circuitTransform :: Circuit a b -> Transform a b +circuitTransform cir = Transform start + where + start m = + Simulation $ \r -> + do ref <- newIORef cir + invokeSimulation r $ + memo0Dynamics (next ref m) + next ref m = + Dynamics $ \p -> + do a <- invokeDynamics p m + cir <- readIORef ref + (cir', b) <- + invokeDynamics p $ + runEvent (runCircuit cir a) + writeIORef ref cir' + return b diff --git a/Simulation/Aivika/Event.hs b/Simulation/Aivika/Event.hs --- a/Simulation/Aivika/Event.hs +++ b/Simulation/Aivika/Event.hs @@ -24,6 +24,7 @@ enqueueEventWithCancellation, enqueueEventWithTimes, enqueueEventWithIntegTimes, + yieldEvent, eventQueueCount, -- * Cancelling Event EventCancellation, @@ -35,6 +36,7 @@ finallyEvent, throwEvent, -- * Memoization - memoEvent) where + memoEvent, + memoEventInTime) where import Simulation.Aivika.Internal.Event diff --git a/Simulation/Aivika/Internal/Arrival.hs b/Simulation/Aivika/Internal/Arrival.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Internal/Arrival.hs @@ -0,0 +1,39 @@ + +-- | +-- Module : Simulation.Aivika.Internal.Arrival +-- Copyright : Copyright (c) 2009-2014, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.6.3 +-- +-- This module defines the types and functions for working with the events +-- that can represent something that arrive from outside the model, or +-- represent other things which computation is delayed and hence is not synchronized. +-- +-- Therefore, the additional information is provided about the time and delay of arrival. + +module Simulation.Aivika.Internal.Arrival + (Arrival(..)) where + +import Simulation.Aivika.Event + +-- | It defines when an event has arrived, usually generated by some random stream. +-- +-- Such events should arrive one by one without time lag in the following sense +-- that the model should start awaiting the next event exactly in that time +-- when the previous event has arrived. +-- +-- Another use case is a situation when the actual event is not synchronized with +-- the 'Event' computation, being synchronized with the event queue, nevertheless. +-- Then the arrival is used for providing the additional information about the time +-- at which the event had been actually arrived. +data Arrival a = + Arrival { arrivalValue :: a, + -- ^ the data we received with the event + arrivalTime :: Double, + -- ^ the simulation time at which the event has arrived + arrivalDelay :: Double + -- ^ the delay time which has passed from the time of + -- arriving the previous event + } deriving (Eq, Ord, Show) diff --git a/Simulation/Aivika/Internal/Event.hs b/Simulation/Aivika/Internal/Event.hs --- a/Simulation/Aivika/Internal/Event.hs +++ b/Simulation/Aivika/Internal/Event.hs @@ -28,6 +28,7 @@ enqueueEventWithTimes, enqueueEventWithPoints, enqueueEventWithIntegTimes, + yieldEvent, eventQueueCount, -- * Cancelling Event EventCancellation, @@ -39,7 +40,8 @@ finallyEvent, throwEvent, -- * Memoization - memoEvent) where + memoEvent, + memoEventInTime) where import Data.IORef @@ -360,3 +362,31 @@ do v <- invokeEvent p m writeIORef ref (Just v) return v + +-- | Memoize the 'Event' computation, always returning the same value +-- in the same modeling time. After the time changes, the value is +-- recalculated by demand. +-- +-- It is possible to implement this function efficiently, for the 'Event' +-- computation is always synchronized with the event queue which time +-- flows in one direction only. This synchronization is a key difference +-- between the 'Event' and 'Dynamics' computations. +memoEventInTime :: Event a -> Simulation (Event a) +memoEventInTime m = + do ref <- liftIO $ newIORef Nothing + return $ Event $ \p -> + do x <- readIORef ref + case x of + Just (t, v) | t == pointTime p -> + return v + _ -> + do v <- invokeEvent p m + writeIORef ref (Just (pointTime p, v)) + return v + +-- | Enqueue the event which must be actuated with the current modeling time but later. +yieldEvent :: Event () -> Event () +yieldEvent m = + Event $ \p -> + invokeEvent p $ + enqueueEvent (pointTime p) m diff --git a/Simulation/Aivika/Internal/Process.hs b/Simulation/Aivika/Internal/Process.hs --- a/Simulation/Aivika/Internal/Process.hs +++ b/Simulation/Aivika/Internal/Process.hs @@ -50,8 +50,12 @@ cancelProcessWithId, cancelProcess, processCancelled, + processCancelling, + whenCancellingProcess, -- * Awaiting Signal processAwait, + -- * Yield of Process + processYield, -- * Process Timeout timeoutProcess, timeoutProcessUsingId, @@ -69,7 +73,9 @@ zip3ProcessParallel, unzipProcess, -- * Memoizing Process - memoProcess) where + memoProcess, + -- * Never Ending Process + neverProcess) where import Data.Maybe import Data.IORef @@ -194,7 +200,7 @@ "Another process with the specified identifier " ++ "has been started already: processIdPrepare" else writeIORef (processStarted pid) True - let signal = (contCancellationInitiating $ processCancelSource pid) + let signal = processCancelling pid invokeEvent p $ handleSignal_ signal $ \_ -> do interruptProcess pid @@ -295,6 +301,18 @@ processCancelled :: ProcessId -> Event Bool processCancelled pid = contCancellationInitiated (processCancelSource pid) +-- | Return a signal that notifies about cancelling the process with +-- the specified identifier. +processCancelling :: ProcessId -> Signal () +processCancelling pid = contCancellationInitiating (processCancelSource pid) + +-- | Register a handler that will be invoked in case of cancelling the current process. +whenCancellingProcess :: Event () -> Process () +whenCancellingProcess h = + Process $ \pid -> + liftEvent $ + handleSignal_ (processCancelling pid) $ \() -> h + instance Eq ProcessId where x == y = processReactCont x == processReactCont y -- for the references are unique @@ -443,7 +461,7 @@ -- should be cancelled in case of need. spawnProcess :: ContCancellation -> Process () -> Process () spawnProcess cancellation x = - do pid <- liftSimulation $ newProcessId + do pid <- liftSimulation newProcessId spawnProcessUsingId cancellation pid x -- | Spawn the child process specifying how the child and parent processes @@ -578,3 +596,26 @@ Nothing -> return Nothing Just (Right a) -> return (Just a) Just (Left e) -> throwProcess e + +-- | Yield to allow other 'Process' and 'Event' computations to run +-- at the current simulation time point. +processYield :: Process () +processYield = + Process $ \pid -> + Cont $ \c -> + Event $ \p -> + invokeEvent p $ + enqueueEvent (pointTime p) $ + resumeCont c () + +-- | A computation that never computes the result. It behaves like a black hole for +-- the discontinuous process, although such a process can still be canceled outside +-- (see 'cancelProcessWithId'), but then only its finalization parts (see 'finallyProcess') +-- will be called, usually, to release the resources acquired before. +neverProcess :: Process a +neverProcess = + Process $ \pid -> + Cont $ \c -> + let signal = processCancelling pid + in handleSignal_ signal $ \_ -> + resumeCont c $ error "It must never be computed: neverProcess" diff --git a/Simulation/Aivika/Internal/Signal.hs b/Simulation/Aivika/Internal/Signal.hs --- a/Simulation/Aivika/Internal/Signal.hs +++ b/Simulation/Aivika/Internal/Signal.hs @@ -32,6 +32,8 @@ merge3Signals, merge4Signals, merge5Signals, + -- * Signal Arriving + arrivalSignal, -- * Creating Signal in Time Points newSignalInTimes, newSignalInIntegTimes, @@ -61,6 +63,7 @@ import Simulation.Aivika.Internal.Parameter import Simulation.Aivika.Internal.Simulation import Simulation.Aivika.Internal.Event +import Simulation.Aivika.Internal.Arrival import qualified Simulation.Aivika.Vector as V import qualified Simulation.Aivika.Vector.Unboxed as UV @@ -352,3 +355,22 @@ appendSignalable m1 m2 = Signalable { readSignalable = liftM2 (<>) (readSignalable m1) (readSignalable m2), signalableChanged_ = (signalableChanged_ m1) <> (signalableChanged_ m2) } + +-- | Transform a signal so that the resulting signal returns a sequence of arrivals +-- saving the information about the time points at which the original signal was received. +arrivalSignal :: Signal a -> Signal (Arrival a) +arrivalSignal m = + Signal { handleSignal = \h -> + Event $ \p -> + do r <- newIORef (pointTime p) + invokeEvent p $ + handleSignal m $ \a -> + Event $ \p -> + do t0 <- readIORef r + let t = pointTime p + writeIORef r t + invokeEvent p $ + h Arrival { arrivalValue = a, + arrivalTime = t, + arrivalDelay = t - t0 } + } diff --git a/Simulation/Aivika/Process.hs b/Simulation/Aivika/Process.hs --- a/Simulation/Aivika/Process.hs +++ b/Simulation/Aivika/Process.hs @@ -52,8 +52,12 @@ cancelProcessWithId, cancelProcess, processCancelled, + processCancelling, + whenCancellingProcess, -- * Awaiting Signal processAwait, + -- * Yield of Process + processYield, -- * Process Timeout timeoutProcess, timeoutProcessUsingId, @@ -71,6 +75,8 @@ zip3ProcessParallel, unzipProcess, -- * Memoizing Process - memoProcess) where + memoProcess, + -- * Never Ending Process + neverProcess) where import Simulation.Aivika.Internal.Process diff --git a/Simulation/Aivika/Processor.hs b/Simulation/Aivika/Processor.hs --- a/Simulation/Aivika/Processor.hs +++ b/Simulation/Aivika/Processor.hs @@ -12,7 +12,10 @@ module Simulation.Aivika.Processor (-- * Processor Type Processor(..), - -- * Creating Simple Processor + -- * Processor Primitives + emptyProcessor, + arrProcessor, + accumProcessor, simpleProcessor, statefulProcessor, -- * Specifying Identifier @@ -27,12 +30,19 @@ queueProcessorLoopMerging, queueProcessorLoopSeq, queueProcessorLoopParallel, + -- * Sequencing Processors + processorSeq, -- * Parallelizing Processors processorParallel, processorQueuedParallel, processorPrioritisingOutputParallel, processorPrioritisingInputParallel, - processorPrioritisingInputOutputParallel) where + processorPrioritisingInputOutputParallel, + -- * Arrival Processor + arrivalProcessor, + -- * Integrating with Signals + signalProcessor, + processorSignaling) where import qualified Control.Category as C import Control.Arrow @@ -44,6 +54,8 @@ import Simulation.Aivika.Process import Simulation.Aivika.Stream import Simulation.Aivika.QueueStrategy +import Simulation.Aivika.Signal +import Simulation.Aivika.Internal.Arrival -- | Represents a processor of simulation data. newtype Processor a b = @@ -63,6 +75,7 @@ -- already depend on the Process monad, -- while the pure streams were considered in the -- mentioned article. + instance Arrow Processor where arr = Processor . mapStream @@ -85,39 +98,6 @@ do (xs, ys) <- liftSimulation $ unzipStream xys runStream $ zipStreamSeq (f xs) (g ys) --- N.B. --- Very probably, Processor is not ArrowLoop, --- which would be natural as Process is not MonadFix, --- for the discontinuous process is not irreversible --- and the time flows in one direction only. --- --- -- The implementation is based on article --- -- A New Notation for Arrows by Ross Paterson, --- -- although my streams are different and they --- -- already depend on the Process monad, --- -- while the pure streams were considered in the --- -- mentioned article. --- instance ArrowLoop Processor where --- --- loop (Processor f) = --- Processor $ \xs -> --- Cons $ --- do Cons zs <- liftSimulation $ --- simulationLoop (\(xs, ys) -> --- unzipStream $ f $ zipStreamSeq xs ys) xs --- zs --- --- simulationLoop :: ((b, d) -> Simulation (c, d)) -> b -> Simulation c --- simulationLoop f b = --- mdo (c, d) <- f (b, d) --- return c - --- The implementation is based on article --- A New Notation for Arrows by Ross Paterson, --- although my streams are different and they --- already depend on the Process monad, --- while the pure streams were considered in the --- mentioned article. instance ArrowChoice Processor where left (Processor f) = @@ -144,34 +124,35 @@ do [xs1, xs2] <- liftSimulation $ splitStream 2 xs runStream $ mergeStreams (f xs1) (g xs2) --- These instances are meaningless: --- --- instance SimulationLift (Processor a) where --- liftSimulation = Processor . mapStreamM . const . liftSimulation --- --- instance DynamicsLift (Processor a) where --- liftDynamics = Processor . mapStreamM . const . liftDynamics --- --- instance EventLift (Processor a) where --- liftEvent = Processor . mapStreamM . const . liftEvent --- --- instance ProcessLift (Processor a) where --- liftProcess = Processor . mapStreamM . const -- data first! +-- | A processor that never finishes its work producing an 'emptyStream'. +emptyProcessor :: Processor a b +emptyProcessor = Processor $ const emptyStream -- | Create a simple processor by the specified handling function -- that runs the discontinuous process for each input value to get the output. +arrProcessor :: (a -> Process b) -> Processor a b +arrProcessor = Processor . mapStreamM + +-- | Accumulator that outputs a value determined by the supplied function. +accumProcessor :: (acc -> a -> Process (acc, b)) -> acc -> Processor a b +accumProcessor f acc = + Processor $ \xs -> Cons $ loop xs acc where + loop xs acc = + do (a, xs') <- runStream xs + (acc', b) <- f acc a + return (b, Cons $ loop xs' acc') + +-- | Create a simple processor by the specified handling function +-- that runs the discontinuous process for each input value to get the output. simpleProcessor :: (a -> Process b) -> Processor a b +{-# DEPRECATED simpleProcessor "Use arrProcessor instead" #-} simpleProcessor = Processor . mapStreamM -- | Like 'simpleProcessor' but allows creating a processor that has a state -- which is passed in to every new iteration. statefulProcessor :: s -> ((s, a) -> Process (s, b)) -> Processor a b -statefulProcessor s f = - Processor $ \xs -> Cons $ loop s xs where - loop s xs = - do (a, xs') <- runStream xs - (s', b) <- f (s, a) - return (b, Cons $ loop s' xs') +{-# DEPRECATED statefulProcessor "Use accumProcessor instead" #-} +statefulProcessor s f = accumProcessor (\acc a -> f (s, a)) s -- | Create a processor that will use the specified process identifier. -- It can be useful to refer to the underlying 'Process' computation which @@ -201,7 +182,7 @@ Processor $ \xs -> Cons $ do let n = length ps - input <- liftSimulation $ splitStreamQueuing si n xs + input <- liftSimulation $ splitStreamQueueing si n xs let results = flip map (zip input ps) $ \(input, p) -> runProcessor p input output = concatQueuedStreams so results @@ -222,7 +203,7 @@ Processor $ \xs -> Cons $ do let n = length ps - input <- liftSimulation $ splitStreamQueuing si n xs + input <- liftSimulation $ splitStreamQueueing si n xs let results = flip map (zip input ps) $ \(input, p) -> runProcessor p input output = concatPriorityStreams so results @@ -277,6 +258,13 @@ processorParallel :: [Processor a b] -> Processor a b processorParallel = processorQueuedParallel FCFS FCFS +-- | Launches the processors sequentially using the 'prefetchProcessor' between them +-- to model an autonomous work of each of the processors specified. +processorSeq :: [Processor a a] -> Processor a a +processorSeq [] = emptyProcessor +processorSeq [p] = p +processorSeq (p : ps) = p >>> prefetchProcessor >>> processorSeq ps + -- | Create a buffer processor, where the process from the first argument -- consumes the input stream but the stream passed in as the second argument -- and produced usually by some other process is returned as an output. @@ -436,3 +424,44 @@ -- for modeling a sequence of separate and independent work places. prefetchProcessor :: Processor a a prefetchProcessor = Processor prefetchStream + +-- | Convert the specified signal transform to a processor. +-- +-- The processor may return data with delay as the values are requested by demand. +-- Consider using the 'arrivalSignal' function to provide with the information +-- about the time points at which the signal was actually triggered. +-- +-- The point is that the 'Stream' used in the 'Processor' is requested outside, +-- while the 'Signal' is triggered inside. They are different by nature. +-- The former is passive, while the latter is active. +-- +-- Cancel the processor's process to unsubscribe from the signals provided. +signalProcessor :: (Signal a -> Signal b) -> Processor a b +signalProcessor f = + Processor $ \xs -> + Cons $ + do sa <- streamSignal xs + sb <- signalStream (f sa) + runStream sb + +-- | Convert the specified processor to a signal transform. +-- +-- The processor may return data with delay as the values are requested by demand. +-- Consider using the 'arrivalSignal' function to provide with the information +-- about the time points at which the signal was actually triggered. +-- +-- The point is that the 'Stream' used in the 'Processor' is requested outside, +-- while the 'Signal' is triggered inside. They are different by nature. +-- The former is passive, while the latter is active. +-- +-- Cancel the returned process to unsubscribe from the signal specified. +processorSignaling :: Processor a b -> Signal a -> Process (Signal b) +processorSignaling (Processor f) sa = + do xs <- signalStream sa + let ys = f xs + streamSignal ys + +-- | A processor that adds the information about the time points at which +-- the original stream items were received by demand. +arrivalProcessor :: Processor a (Arrival a) +arrivalProcessor = Processor arrivalStream diff --git a/Simulation/Aivika/Queue.hs b/Simulation/Aivika/Queue.hs --- a/Simulation/Aivika/Queue.hs +++ b/Simulation/Aivika/Queue.hs @@ -117,8 +117,6 @@ import Simulation.Aivika.Resource import Simulation.Aivika.QueueStrategy import Simulation.Aivika.Statistics -import Simulation.Aivika.Stream -import Simulation.Aivika.Processor import qualified Simulation.Aivika.DoubleLinkedList as DLL import qualified Simulation.Aivika.Vector as V diff --git a/Simulation/Aivika/Queue/Infinite.hs b/Simulation/Aivika/Queue/Infinite.hs --- a/Simulation/Aivika/Queue/Infinite.hs +++ b/Simulation/Aivika/Queue/Infinite.hs @@ -81,8 +81,6 @@ import Simulation.Aivika.Resource import Simulation.Aivika.QueueStrategy import Simulation.Aivika.Statistics -import Simulation.Aivika.Stream -import Simulation.Aivika.Processor import qualified Simulation.Aivika.DoubleLinkedList as DLL import qualified Simulation.Aivika.Vector as V diff --git a/Simulation/Aivika/Ref/Light.hs b/Simulation/Aivika/Ref/Light.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Ref/Light.hs @@ -0,0 +1,53 @@ + +-- | +-- Module : Simulation.Aivika.Ref.Light +-- Copyright : Copyright (c) 2009-2014, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.6.3 +-- +-- This module defines a light-weight and more fast version of an updatable reference +-- that depends on the event queue but that doesn't supply with the signal notification. +-- +module Simulation.Aivika.Ref.Light + (Ref, + newRef, + readRef, + writeRef, + modifyRef) where + +import Data.IORef +import Control.Monad +import Control.Monad.Trans + +import Simulation.Aivika.Internal.Simulation +import Simulation.Aivika.Internal.Event + +-- | The 'Ref' type represents a mutable variable similar to the 'IORef' variable +-- but only dependent on the event queue, which allows synchronizing the reference +-- with the model explicitly through the 'Event' monad. +newtype Ref a = + Ref { refValue :: IORef a } + +-- | Create a new reference. +newRef :: a -> Simulation (Ref a) +newRef a = + do x <- liftIO $ newIORef a + return Ref { refValue = x } + +-- | Read the value of a reference. +readRef :: Ref a -> Event a +readRef r = Event $ \p -> readIORef (refValue r) + +-- | Write a new value into the reference. +writeRef :: Ref a -> a -> Event () +writeRef r a = Event $ \p -> + a `seq` writeIORef (refValue r) a + +-- | Mutate the contents of the reference. +modifyRef :: Ref a -> (a -> a) -> Event () +modifyRef r f = Event $ \p -> + do a <- readIORef (refValue r) + let b = f a + b `seq` writeIORef (refValue r) b diff --git a/Simulation/Aivika/Server.hs b/Simulation/Aivika/Server.hs --- a/Simulation/Aivika/Server.hs +++ b/Simulation/Aivika/Server.hs @@ -12,6 +12,7 @@ (-- * Server Server, newServer, + newStateServer, newServerWithState, -- * Processing serverProcessor, @@ -80,7 +81,7 @@ -- ^ The initial state of the server. serverStateRef :: IORef s, -- ^ The current state of the server. - serverProcess :: (s, a) -> Process (s, b), + serverProcess :: s -> a -> Process (s, b), -- ^ Provide @b@ by specified @a@. serverTotalInputWaitTimeRef :: IORef Double, -- ^ The counted total time spent in awating the input. @@ -117,13 +118,13 @@ -- | Create a new server that can provide output @b@ by input @a@ -- starting from state @s@. Also it returns the corresponded processor -- that being applied updates the server state. -newServerWithState :: s - -- ^ the initial state - -> ((s, a) -> Process (s, b)) - -- ^ provide an output by the specified input - -- and update the state - -> Simulation (Server s a b) -newServerWithState state provide = +newStateServer :: (s -> a -> Process (s, b)) + -- ^ provide a new state and output by the specified + -- old state and input + -> s + -- ^ the initial state + -> Simulation (Server s a b) +newStateServer provide state = do r0 <- liftIO $ newIORef state r1 <- liftIO $ newIORef 0 r2 <- liftIO $ newIORef 0 @@ -148,6 +149,18 @@ serverOutputProvidedSource = s3 } return server +-- | Create a new server that can provide output @b@ by input @a@ +-- starting from state @s@. Also it returns the corresponded processor +-- that being applied updates the server state. +newServerWithState :: s + -- ^ the initial state + -> ((s, a) -> Process (s, b)) + -- ^ provide an output by the specified input + -- and update the state + -> Simulation (Server s a b) +{-# DEPRECATED newServerWithState "Use newStateServer instead" #-} +newServerWithState state provide = newStateServer (curry provide) state + -- | Return a processor for the specified server. -- -- The processor updates the internal state of the server. The usual case is when @@ -195,7 +208,7 @@ addSamplingStats (t1 - t0) triggerSignal (serverInputReceivedSource server) a -- provide the service - (s', b) <- serverProcess server (s, a) + (s', b) <- serverProcess server s a t2 <- liftDynamics time liftEvent $ do liftIO $ @@ -204,7 +217,7 @@ modifyIORef' (serverProcessingTimeRef server) $ addSamplingStats (t2 - t1) triggerSignal (serverTaskProcessedSource server) (a, b) - return (b, loop s' (Just $ (t2, a, b)) xs') + return (b, loop s' (Just (t2, a, b)) xs') -- | Return the current state of the server. -- diff --git a/Simulation/Aivika/Signal.hs b/Simulation/Aivika/Signal.hs --- a/Simulation/Aivika/Signal.hs +++ b/Simulation/Aivika/Signal.hs @@ -31,6 +31,8 @@ merge3Signals, merge4Signals, merge5Signals, + -- * Signal Arriving + arrivalSignal, -- * Creating Signal in Time Points newSignalInTimes, newSignalInIntegTimes, diff --git a/Simulation/Aivika/Stream.hs b/Simulation/Aivika/Stream.hs --- a/Simulation/Aivika/Stream.hs +++ b/Simulation/Aivika/Stream.hs @@ -22,11 +22,14 @@ concatPriorityStreams, splitStream, splitStreamQueuing, + splitStreamQueueing, splitStreamPrioritising, -- * Specifying Identifier streamUsingId, -- * Prefetching Stream prefetchStream, + -- * Stream Arriving + arrivalStream, -- * Memoizing, Zipping and Uzipping Stream memoStream, zipStreamSeq, @@ -48,6 +51,9 @@ apStreamParallel, filterStream, filterStreamM, + -- * Integrating with Signals + signalStream, + streamSignal, -- * Utilities leftStream, rightStream, @@ -63,10 +69,15 @@ import Control.Monad.Trans import Simulation.Aivika.Simulation +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Event import Simulation.Aivika.Cont import Simulation.Aivika.Process +import Simulation.Aivika.Signal import Simulation.Aivika.Resource import Simulation.Aivika.QueueStrategy +import Simulation.Aivika.Queue.Infinite +import Simulation.Aivika.Internal.Arrival -- | Represents an infinite stream of data in time, -- some kind of the cons cell. @@ -149,7 +160,7 @@ -- This is a generalization of 'zipStreamSeq'. streamSeq :: [Stream a] -> Stream [a] streamSeq xs = Cons y where - y = do ps <- forM xs $ runStream + y = do ps <- forM xs runStream return (map fst ps, streamSeq $ map snd ps) -- | To form each new portion of data for the output stream, @@ -263,23 +274,23 @@ -- | Split the input stream into the specified number of output streams -- after applying the 'FCFS' strategy for enqueuing the output requests. splitStream :: Int -> Stream a -> Simulation [Stream a] -splitStream = splitStreamQueuing FCFS +splitStream = splitStreamQueueing FCFS -- | Split the input stream into the specified number of output streams. -- -- If you don't know what the strategy to apply, then you probably -- need the 'FCFS' strategy, or function 'splitStream' that -- does namely this. -splitStreamQueuing :: EnqueueStrategy s q - => s - -- ^ the strategy applied for enqueuing the output requests - -> Int - -- ^ the number of output streams - -> Stream a - -- ^ the input stream - -> Simulation [Stream a] - -- ^ the splitted output streams -splitStreamQueuing s n x = +splitStreamQueueing :: EnqueueStrategy s q + => s + -- ^ the strategy applied for enqueuing the output requests + -> Int + -- ^ the number of output streams + -> Stream a + -- ^ the input stream + -> Simulation [Stream a] + -- ^ the splitted output streams +splitStreamQueueing s n x = do ref <- liftIO $ newIORef x res <- newResource s 1 let reader = @@ -290,6 +301,19 @@ return a return $ map (\i -> repeatProcess reader) [1..n] +-- | It was renamed to 'splitStreamQueueing'. +{-# DEPRECATED splitStreamQueuing "Use splitStreamQueueing instead" #-} +splitStreamQueuing :: EnqueueStrategy s q + => s + -- ^ the strategy applied for enqueuing the output requests + -> Int + -- ^ the number of output streams + -> Stream a + -- ^ the input stream + -> Simulation [Stream a] + -- ^ the splitted output streams +splitStreamQueuing = splitStreamQueueing + -- | Split the input stream into a list of output streams -- using the specified priorities. splitStreamPrioritising :: PriorityQueueStrategy s q p @@ -419,20 +443,14 @@ -- | An empty stream that never returns data. emptyStream :: Stream a -emptyStream = Cons z where - z = do pid <- liftSimulation newProcessId - -- use the generated identifier so that - -- nobody could reactivate the process, - -- although it can be still canceled - processUsingId pid passivateProcess - error "It should never happen: emptyStream." +emptyStream = Cons neverProcess -- | Consume the stream. It returns a process that infinitely reads data -- from the stream and then redirects them to the provided function. -- It is useful for modeling the process of enqueueing data in the queue -- from the input stream. consumeStream :: (a -> Process ()) -> Stream a -> Process () -consumeStream f s = p s where +consumeStream f = p where p (Cons s) = do (a, xs) <- s f a p xs @@ -442,7 +460,7 @@ -- to simulate the whole system of the interconnected streams and -- processors. sinkStream :: Stream a -> Process () -sinkStream s = p s where +sinkStream = p where p (Cons s) = do (a, xs) <- s p xs @@ -472,3 +490,53 @@ return a spawnProcess CancelTogether $ writer s runStream $ repeatProcess reader + +-- | Return a stream of values triggered by the specified signal. +-- +-- Since the time at which the values of the stream are requested for may differ from +-- the time at which the signal is triggered, it can be useful to apply the 'arrivalSignal' +-- function to add the information about the time points at which the signal was +-- actually received. +-- +-- The point is that the 'Stream' is requested outside, while the 'Signal' is triggered +-- inside. They are different by nature. The former is passive, while the latter is active. +-- +-- The resulting stream may be a root of space leak as it uses an internal queue to store +-- the values received from the signal. The oldest value is dequeued each time we request +-- the stream and it is returned within the computation. +-- +-- Cancel the stream's process to unsubscribe from the specified signal. +signalStream :: Signal a -> Process (Stream a) +signalStream s = + do q <- liftSimulation newFCFSQueue + h <- liftEvent $ + handleSignal s $ + enqueue q + whenCancellingProcess h + return $ repeatProcess $ dequeue q + +-- | Return a computation of the signal that triggers values from the specified stream, +-- each time the next value of the stream is received within the underlying 'Process' +-- computation. +-- +-- Cancel the returned process to stop reading from the specified stream. +streamSignal :: Stream a -> Process (Signal a) +streamSignal z = + do s <- liftSimulation newSignalSource + spawnProcess CancelTogether $ + consumeStream (liftEvent . triggerSignal s) z + return $ publishSignal s + +-- | Transform a stream so that the resulting stream returns a sequence of arrivals +-- saving the information about the time points at which the original stream items +-- were received by demand. +arrivalStream :: Stream a -> Stream (Arrival a) +arrivalStream s = Cons z where + z = do t <- liftDynamics time + loop s t + loop s t0 = do (a, xs) <- runStream s + t <- liftDynamics time + let b = Arrival { arrivalValue = a, + arrivalTime = t, + arrivalDelay = t - t0 } + return (b, Cons $ loop xs t) diff --git a/Simulation/Aivika/Stream/Random.hs b/Simulation/Aivika/Stream/Random.hs --- a/Simulation/Aivika/Stream/Random.hs +++ b/Simulation/Aivika/Stream/Random.hs @@ -39,7 +39,10 @@ import Simulation.Aivika.Arrival -- | Return a sream of random events that arrive with the specified delay. -randomStream :: Parameter Double -> Stream Arrival +randomStream :: Parameter (Double, a) + -- ^ compute a pair of the delay and event of type @a@ + -> Stream (Arrival a) + -- ^ a stream of delayed events randomStream delay = Cons z0 where z0 = do t0 <- liftDynamics time @@ -53,10 +56,11 @@ "contains a logical error. The random events should be requested permanently. " ++ "At least, they can be lost, for example, when trying to enqueue them, but " ++ "the random stream itself must always work: randomStream." - delay <- liftParameter delay + (delay, a) <- liftParameter delay holdProcess delay t2 <- liftDynamics time - let arrival = Arrival { arrivalTime = t2, + let arrival = Arrival { arrivalValue = a, + arrivalTime = t2, arrivalDelay = delay } return (arrival, Cons $ loop t2) @@ -65,29 +69,35 @@ -- ^ the minimum delay -> Double -- ^ the maximum delay - -> Stream Arrival - -- ^ the stream of random events + -> Stream (Arrival Double) + -- ^ the stream of random events with the delays generated randomUniformStream min max = - randomStream $ randomUniform min max + randomStream $ + randomUniform min max >>= \x -> + return (x, x) -- | Create a new stream with delays distributed normally. randomNormalStream :: Double -- ^ the mean delay -> Double -- ^ the delay deviation - -> Stream Arrival - -- ^ the stream of random events + -> Stream (Arrival Double) + -- ^ the stream of random events with the delays generated randomNormalStream mu nu = - randomStream $ randomNormal mu nu + randomStream $ + randomNormal mu nu >>= \x -> + return (x, x) -- | Return a new stream with delays distibuted exponentially with the specified mean -- (the reciprocal of the rate). randomExponentialStream :: Double -- ^ the mean delay (the reciprocal of the rate) - -> Stream Arrival - -- ^ the stream of random events + -> Stream (Arrival Double) + -- ^ the stream of random events with the delays generated randomExponentialStream mu = - randomStream $ randomExponential mu + randomStream $ + randomExponential mu >>= \x -> + return (x, x) -- | Return a new stream with delays having the Erlang distribution with the specified -- scale (the reciprocal of the rate) and shape parameters. @@ -95,19 +105,23 @@ -- ^ the scale (the reciprocal of the rate) -> Int -- ^ the shape - -> Stream Arrival - -- ^ the stream of random events + -> Stream (Arrival Double) + -- ^ the stream of random events with the delays generated randomErlangStream beta m = - randomStream $ randomErlang beta m + randomStream $ + randomErlang beta m >>= \x -> + return (x, x) -- | Return a new stream with delays having the Poisson distribution with -- the specified mean. randomPoissonStream :: Double -- ^ the mean delay - -> Stream Arrival - -- ^ the stream of random events + -> Stream (Arrival Int) + -- ^ the stream of random events with the delays generated randomPoissonStream mu = - randomStream $ fmap fromIntegral $ randomPoisson mu + randomStream $ + randomPoisson mu >>= \x -> + return (fromIntegral x, x) -- | Return a new stream with delays having the binomial distribution with the specified -- probability and trials. @@ -115,7 +129,9 @@ -- ^ the probability -> Int -- ^ the number of trials - -> Stream Arrival - -- ^ the stream of random events + -> Stream (Arrival Int) + -- ^ the stream of random events with the delays generated randomBinomialStream prob trials = - randomStream $ fmap fromIntegral $ randomBinomial prob trials + randomStream $ + randomBinomial prob trials >>= \x -> + return (fromIntegral x, x) diff --git a/Simulation/Aivika/SystemDynamics.hs b/Simulation/Aivika/SystemDynamics.hs --- a/Simulation/Aivika/SystemDynamics.hs +++ b/Simulation/Aivika/SystemDynamics.hs @@ -119,13 +119,13 @@ -- integEuler :: Dynamics Double - -> Dynamics Double + -> Double -> Dynamics Double -> Point -> IO Double -integEuler (Dynamics f) (Dynamics i) (Dynamics y) p = +integEuler (Dynamics f) i (Dynamics y) p = case pointIteration p of - 0 -> - i p + 0 -> + return i n -> do let sc = pointSpecs p ty = basicTime sc (n - 1) 0 @@ -136,14 +136,14 @@ return v integRK2 :: Dynamics Double - -> Dynamics Double + -> Double -> Dynamics Double -> Point -> IO Double -integRK2 (Dynamics f) (Dynamics i) (Dynamics y) p = +integRK2 (Dynamics f) i (Dynamics y) p = case pointPhase p of 0 -> case pointIteration p of 0 -> - i p + return i n -> do let sc = pointSpecs p ty = basicTime sc (n - 1) 0 @@ -172,14 +172,14 @@ error "Incorrect phase: integRK2" integRK4 :: Dynamics Double - -> Dynamics Double + -> Double -> Dynamics Double -> Point -> IO Double -integRK4 (Dynamics f) (Dynamics i) (Dynamics y) p = +integRK4 (Dynamics f) i (Dynamics y) p = case pointPhase p of 0 -> case pointIteration p of 0 -> - i p + return i n -> do let sc = pointSpecs p ty = basicTime sc (n - 1) 0 @@ -251,8 +251,11 @@ -- kb = 1 -- runDynamicsInStopTime $ sequence [a, b, c] -- @ +-- +-- To receive the initial value for an abitrary 'Dynamics' computation, +-- you can always use the 'runDynamicsInStartTime' function. integ :: Dynamics Double -- ^ the derivative - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the integral integ diff i = mdo y <- MU.memoDynamics z @@ -275,7 +278,7 @@ -- @ smoothI :: Dynamics Double -- ^ the value to smooth over time -> Dynamics Double -- ^ time - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the first order exponential smooth smoothI x t i = mdo y <- integ ((x - y) / t) i @@ -288,7 +291,9 @@ smooth :: Dynamics Double -- ^ the value to smooth over time -> Dynamics Double -- ^ time -> Simulation (Dynamics Double) -- ^ the first order exponential smooth -smooth x t = smoothI x t x +smooth x t = + do i <- runDynamicsInStartTime x + smoothI x t i -- | Return the third order exponential smooth. -- @@ -305,7 +310,7 @@ -- @ smooth3I :: Dynamics Double -- ^ the value to smooth over time -> Dynamics Double -- ^ time - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the third order exponential smooth smooth3I x t i = mdo y <- integ ((s2 - y) / t') i @@ -321,7 +326,9 @@ smooth3 :: Dynamics Double -- ^ the value to smooth over time -> Dynamics Double -- ^ time -> Simulation (Dynamics Double) -- ^ the third order exponential smooth -smooth3 x t = smooth3I x t x +smooth3 x t = + do i <- runDynamicsInStartTime x + smooth3I x t i -- | Return the n'th order exponential smooth. -- @@ -332,7 +339,7 @@ smoothNI :: Dynamics Double -- ^ the value to smooth over time -> Dynamics Double -- ^ time -> Int -- ^ the order - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the n'th order exponential smooth smoothNI x t n i = mdo s <- forM [1 .. n] $ \k -> @@ -351,7 +358,9 @@ -> Dynamics Double -- ^ time -> Int -- ^ the order -> Simulation (Dynamics Double) -- ^ the n'th order exponential smooth -smoothN x t n = smoothNI x t n x +smoothN x t n = + do i <- runDynamicsInStartTime x + smoothNI x t n i -- | Return the first order exponential delay. -- @@ -360,15 +369,17 @@ -- -- @ -- delay1I x t i = --- mdo y <- integ (x - y \/ t) (i * t) +-- mdo t0 <- runDynamicsInStartTime t +-- y <- integ (x - y \/ t) (i * t0) -- return $ y \/ t -- @ delay1I :: Dynamics Double -- ^ the value to conserve -> Dynamics Double -- ^ time - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the first order exponential delay delay1I x t i = - mdo y <- integ (x - y / t) (i * t) + mdo t0 <- runDynamicsInStartTime t + y <- integ (x - y / t) (i * t0) return $ y / t -- | Return the first order exponential delay. @@ -378,17 +389,20 @@ delay1 :: Dynamics Double -- ^ the value to conserve -> Dynamics Double -- ^ time -> Simulation (Dynamics Double) -- ^ the first order exponential delay -delay1 x t = delay1I x t x +delay1 x t = + do i <- runDynamicsInStartTime x + delay1I x t i -- | Return the third order exponential delay. delay3I :: Dynamics Double -- ^ the value to conserve -> Dynamics Double -- ^ time - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the third order exponential delay delay3I x t i = - mdo y <- integ (s2 / t' - y / t') (i * t') - s2 <- integ (s1 / t' - s2 / t') (i * t') - s1 <- integ (x - s1 / t') (i * t') + mdo t0' <- runDynamicsInStartTime t' + y <- integ (s2 / t' - y / t') (i * t0') + s2 <- integ (s1 / t' - s2 / t') (i * t0') + s1 <- integ (x - s1 / t') (i * t0') let t' = t / 3.0 return $ y / t' @@ -399,19 +413,22 @@ delay3 :: Dynamics Double -- ^ the value to conserve -> Dynamics Double -- ^ time -> Simulation (Dynamics Double) -- ^ the third order exponential delay -delay3 x t = delay3I x t x +delay3 x t = + do i <- runDynamicsInStartTime x + delay3I x t i -- | Return the n'th order exponential delay. delayNI :: Dynamics Double -- ^ the value to conserve -> Dynamics Double -- ^ time -> Int -- ^ the order - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the n'th order exponential delay delayNI x t n i = - mdo s <- forM [1 .. n] $ \k -> + mdo t0' <- runDynamicsInStartTime t' + s <- forM [1 .. n] $ \k -> if k == 1 - then integ (x - (a ! 1) / t') (i * t') - else integ ((a ! (k - 1)) / t' - (a ! k) / t') (i * t') + then integ (x - (a ! 1) / t') (i * t0') + else integ ((a ! (k - 1)) / t' - (a ! k) / t') (i * t0') let a = listArray (1, n) s t' = t / fromIntegral n return $ (a ! n) / t' @@ -424,7 +441,9 @@ -> Dynamics Double -- ^ time -> Int -- ^ the order -> Simulation (Dynamics Double) -- ^ the n'th order exponential delay -delayN x t n = delayNI x t n x +delayN x t n = + do i <- runDynamicsInStartTime x + delayNI x t n i -- | Return the forecast. -- @@ -449,16 +468,20 @@ -- -- @ -- trend x at i = --- do y <- smoothI x at (x \/ (1.0 + i * at)) --- return $ (x \/ y - 1.0) \/ at +-- mdo x0 <- runDynamicsInStartTime x +-- at0 <- runDynamicsInStartTime at +-- y <- smoothI x at (x0 \/ (1.0 + i * at0)) +-- return $ (x \/ y - 1.0) \/ at -- @ trend :: Dynamics Double -- ^ the value for which the trend is calculated -> Dynamics Double -- ^ the average time - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Simulation (Dynamics Double) -- ^ the fractional change rate trend x at i = - do y <- smoothI x at (x / (1.0 + i * at)) - return $ (x / y - 1.0) / at + mdo x0 <- runDynamicsInStartTime x + at0 <- runDynamicsInStartTime at + y <- smoothI x at (x0 / (1.0 + i * at0)) + return $ (x / y - 1.0) / at -- -- Difference Equations @@ -471,14 +494,15 @@ -- As usual, to create a loopback, you should use the recursive do-notation. diffsum :: (Num a, Unboxed a) => Dynamics a -- ^ the difference - -> Dynamics a -- ^ the initial value + -> a -- ^ the initial value -> Simulation (Dynamics a) -- ^ the sum -diffsum (Dynamics diff) (Dynamics i) = +diffsum (Dynamics diff) i = mdo y <- MU.memo0Dynamics $ Dynamics $ \p -> case pointIteration p of - 0 -> i p + 0 -> + return i n -> do let Dynamics m = y sc = pointSpecs p @@ -545,9 +569,9 @@ -- Because of the latter, it allows creating a loop back. delayI :: Dynamics a -- ^ the value to delay -> Dynamics Double -- ^ the lag time - -> Dynamics a -- ^ the initial value + -> a -- ^ the initial value -> Simulation (Dynamics a) -- ^ the delayed value -delayI (Dynamics x) (Dynamics d) (Dynamics i) = M.memo0Dynamics $ Dynamics r +delayI (Dynamics x) (Dynamics d) i = M.memo0Dynamics $ Dynamics r where r p = do let t = pointTime p @@ -556,9 +580,7 @@ a <- d p let t' = t - a n' = fromIntegral $ floor $ (t' - spcStartTime sc) / spcDT sc - y | n' < 0 = i $ p { pointTime = spcStartTime sc, - pointIteration = 0, - pointPhase = 0 } + y | n' < 0 = return i | n' < n = x $ p { pointTime = t', pointIteration = n', pointPhase = -1 } @@ -582,18 +604,18 @@ -- @ -- npv stream rate init factor = -- mdo let dt' = liftParameter dt --- df <- integ (- df * rate) 1 +-- df <- integ (- df * rate) 1 -- accum <- integ (stream * df) init -- return $ (accum + dt' * stream * df) * factor -- @ npv :: Dynamics Double -- ^ the stream -> Dynamics Double -- ^ the discount rate - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Dynamics Double -- ^ factor -> Simulation (Dynamics Double) -- ^ the Net Present Value (NPV) npv stream rate init factor = mdo let dt' = liftParameter dt - df <- integ (- df * rate) 1 + df <- integ (- df * rate) 1 accum <- integ (stream * df) init return $ (accum + dt' * stream * df) * factor @@ -605,18 +627,22 @@ -- @ -- npve stream rate init factor = -- mdo let dt' = liftParameter dt --- df <- integ (- df * rate \/ (1 + rate * dt')) (1 \/ (1 + rate * dt')) +-- rate0 <- runDynamicsInStartTime rate +-- dt0 <- liftParameter dt +-- df <- integ (- df * rate \/ (1 + rate * dt')) (1 \/ (1 + rate0 * dt0)) -- accum <- integ (stream * df) init -- return $ (accum + dt' * stream * df) * factor -- @ npve :: Dynamics Double -- ^ the stream -> Dynamics Double -- ^ the discount rate - -> Dynamics Double -- ^ the initial value + -> Double -- ^ the initial value -> Dynamics Double -- ^ factor -> Simulation (Dynamics Double) -- ^ the Net Present Value End (NPVE) npve stream rate init factor = mdo let dt' = liftParameter dt - df <- integ (- df * rate / (1 + rate * dt')) (1 / (1 + rate * dt')) + rate0 <- runDynamicsInStartTime rate + dt0 <- liftParameter dt + df <- integ (- df * rate / (1 + rate * dt')) (1 / (1 + rate0 * dt0)) accum <- integ (stream * df) init return $ (accum + dt' * stream * df) * factor diff --git a/Simulation/Aivika/Task.hs b/Simulation/Aivika/Task.hs --- a/Simulation/Aivika/Task.hs +++ b/Simulation/Aivika/Task.hs @@ -80,9 +80,7 @@ do x <- liftIO $ readIORef (taskResultRef t) case x of Just x -> return x - Nothing -> - do x <- processAwait (taskResultReceived t) - return x + Nothing -> processAwait (taskResultReceived t) -- | Cancel the task. cancelTask :: Task a -> Event () diff --git a/Simulation/Aivika/Transform.hs b/Simulation/Aivika/Transform.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Transform.hs @@ -0,0 +1,30 @@ + +-- | +-- Module : Simulation.Aivika.Transform +-- Copyright : Copyright (c) 2009-2014, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.6.3 +-- +-- The module defines a transform of one time varying function to another +-- usually specified in the integration time points and then interpolated in +-- other time points with help of one of the memoization functions +-- like 'memo0Dynamics'. +-- +module Simulation.Aivika.Transform + (Transform(..)) where + +import Simulation.Aivika.Simulation +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Memo + +-- | The transform of one time varying function to another usually +-- specified in the integration time points and then interpolated in +-- other time points with help of one of the memoization functions +-- like 'memo0Dynamics'. +-- +newtype Transform a b = + Transform { runTransform :: Dynamics a -> Simulation (Dynamics b) + -- ^ Run the transform. + } diff --git a/aivika.cabal b/aivika.cabal --- a/aivika.cabal +++ b/aivika.cabal @@ -1,5 +1,5 @@ name: aivika -version: 1.1 +version: 1.2 synopsis: A multi-paradigm simulation library description: Aivika is a multi-paradigm simulation library with a strong emphasis @@ -87,6 +87,9 @@ . \[3] . + P.S. Aivika is actually a genuine female Mari name which is pronounced + with stress on the last syllable as in French, but the Russians usually + pronounce it wrong :) category: Simulation license: BSD3 license-file: LICENSE @@ -94,12 +97,13 @@ author: David Sorokin maintainer: David Sorokin homepage: http://github.com/dsorokin/aivika -cabal-version: >= 1.2.0 +cabal-version: >= 1.6 build-type: Simple tested-with: GHC == 7.6.3 extra-source-files: examples/BassDiffusion.hs examples/ChemicalReaction.hs + examples/ChemicalReactionCircuit.hs examples/FishBank.hs examples/MachRep1.hs examples/MachRep1EventDriven.hs @@ -107,11 +111,11 @@ examples/MachRep2.hs examples/MachRep3.hs examples/Furnace.hs - examples/SimpleWorkflow.hs + examples/InspectionAdjustmentStations.hs + examples/WorkStationsInSeries.hs examples/TimeOut.hs examples/TimeOutInt.hs examples/TimeOutWait.hs - examples/WorkflowLoop.hs examples/README data-files: doc/aivika.pdf @@ -121,6 +125,7 @@ exposed-modules: Simulation.Aivika Simulation.Aivika.Agent Simulation.Aivika.Arrival + Simulation.Aivika.Circuit Simulation.Aivika.Cont Simulation.Aivika.DoubleLinkedList Simulation.Aivika.Dynamics @@ -141,6 +146,7 @@ Simulation.Aivika.Queue.Infinite Simulation.Aivika.QueueStrategy Simulation.Aivika.Ref + Simulation.Aivika.Ref.Light Simulation.Aivika.Resource Simulation.Aivika.Server Simulation.Aivika.Signal @@ -153,6 +159,7 @@ Simulation.Aivika.SystemDynamics Simulation.Aivika.Table Simulation.Aivika.Task + Simulation.Aivika.Transform Simulation.Aivika.Unboxed Simulation.Aivika.Var Simulation.Aivika.Var.Unboxed @@ -167,6 +174,7 @@ Simulation.Aivika.Internal.Signal Simulation.Aivika.Internal.Simulation Simulation.Aivika.Internal.Specs + Simulation.Aivika.Internal.Arrival build-depends: base >= 4.5.0.0 && < 6, mtl >= 2.1.1, @@ -177,7 +185,13 @@ extensions: FlexibleContexts, BangPatterns, RecursiveDo, + Arrows, MultiParamTypeClasses, FunctionalDependencies ghc-options: -O2 + +source-repository head + + type: git + location: https://github.com/dsorokin/aivika diff --git a/examples/ChemicalReactionCircuit.hs b/examples/ChemicalReactionCircuit.hs new file mode 100644 --- /dev/null +++ b/examples/ChemicalReactionCircuit.hs @@ -0,0 +1,43 @@ + +-- Note that the integCircut function uses Euler's method regardless of +-- the simulation specs specified. Therefore, to receieve almost the same +-- results in the old example based on using the integ function, you should +-- specify Euler's method in their specs in that file, although the Runge-Kutta +-- method gives similar results too, which is expected. +-- +-- Finally, the integ function can be significantly faster than integCircuit, +-- although they have different purposes. + +{-# LANGUAGE Arrows #-} + +import Control.Arrow + +import Simulation.Aivika + +specs = Specs { spcStartTime = 0, + spcStopTime = 13, + spcDT = 0.01, + spcMethod = RungeKutta4, + spcGeneratorType = SimpleGenerator } + +circuit :: Circuit () [Double] +circuit = + let ka = 1 + kb = 1 + in proc () -> do + rec let da = - ka * a + db = ka * a - kb * b + dc = kb * b + a <- integCircuit 100 -< da + b <- integCircuit 0 -< db + c <- integCircuit 0 -< dc + returnA -< [a, b, c] + +model :: Simulation [Double] +model = + do results <- + runTransform (circuitTransform circuit) $ + return () + runDynamicsInStopTime results + +main = runSimulation model specs >>= print diff --git a/examples/InspectionAdjustmentStations.hs b/examples/InspectionAdjustmentStations.hs new file mode 100644 --- /dev/null +++ b/examples/InspectionAdjustmentStations.hs @@ -0,0 +1,215 @@ + +{-# LANGUAGE RecursiveDo, Arrows #-} + +-- Example: Inspection and Adjustment Stations on a Production Line +-- +-- This is a model of the workflow with a loop. Also there are two infinite queues. +-- +-- It is described in different sources [1, 2]. So, this is chapter 8 of [2] and section 5.15 of [1]. +-- +-- [1] A. Alan B. Pritsker, Simulation with Visual SLAM and AweSim, 2nd ed. +-- +-- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 + +-- CAUTION: +-- +-- This model is not yet fully tested and it may contain logical errors but it seems to be working, +-- although some results may differ slightly but it can be related to a great value of the deviation +-- for some variables as well as to a small number of samples in [1]. +-- +-- The results for the queue sizes in [2] seem doubtful for me, while my results for these queue sizes +-- are similar to [1] but I also made 1000 runs (see the aivika-experiment-chart package) versus 1 run +-- in [1]. In comparison with [1] I see a difference in the queue size for the adjustment station and +-- it can be realized as there was a too small number of samples (= 13) in [1], for the TV settings must +-- fail when inspecting to be directed to the adjustor. +-- +-- Also I have received more small values for the wait time in comparison with [1] but they have +-- a relatively great deviation, which may be acceptable (??), taking into account a small number of +-- samples used in [1]. +-- +-- At the same time, all my other results except for these queue sizes correspond to [2], where the author +-- launched 1000 simulation runs too. +-- +-- Some new things that I have added the past summer (2013), i.e. Streams / Processors / Queues / Servers, +-- should be yet verified for other models but, as I wrote, they seem to be working. + +import Prelude hiding (id, (.)) + +import Control.Monad +import Control.Monad.Trans +import Control.Arrow +import Control.Category (id, (.)) + +import Simulation.Aivika +import Simulation.Aivika.Queue.Infinite + +-- | The simulation specs. +specs = Specs { spcStartTime = 0.0, + spcStopTime = 480.0, + spcDT = 0.1, + spcMethod = RungeKutta4, + spcGeneratorType = SimpleGenerator } + +-- the minimum delay of arriving the next TV set +minArrivalDelay = 3.5 + +-- the maximum delay of arriving the next TV set +maxArrivalDelay = 7.5 + +-- the minimum time to inspect the TV set +minInspectionTime = 6 + +-- the maximum time to inspect the TV set +maxInspectionTime = 12 + +-- the probability of passing the inspection phase +inspectionPassingProb = 0.85 + +-- how many are inspection stations? +inspectionStationCount = 2 + +-- the minimum time to adjust an improper TV set +minAdjustmentTime = 20 + +-- the maximum time to adjust an improper TV set +maxAdjustmentTime = 40 + +-- how many are adjustment stations? +adjustmentStationCount = 1 + +-- create an accumulator to gather the queue size statistics +newQueueSizeAccumulator queue = + newTimingStatsAccumulator $ + Signalable (queueCount queue) (queueCountChanged_ queue) + +-- create an inspection station (server) +newInspectionStation = + newServer $ \a -> + do holdProcess =<< + (liftParameter $ + randomUniform minInspectionTime maxInspectionTime) + passed <- + liftParameter $ + randomTrue inspectionPassingProb + if passed + then return $ Right a + else return $ Left a + +-- create an adjustment station (server) +newAdjustmentStation = + newServer $ \a -> + do holdProcess =<< + (liftParameter $ + randomUniform minAdjustmentTime maxAdjustmentTime) + return a + +model :: Simulation () +model = mdo + -- to count the arrived TV sets for inspecting and adjusting + inputArrivalTimer <- newArrivalTimer + -- it will gather the statistics of the processing time + outputArrivalTimer <- newArrivalTimer + -- define a stream of input events + let inputStream = + randomUniformStream minArrivalDelay maxArrivalDelay + -- create a queue before the inspection stations + inspectionQueue <- newFCFSQueue + -- create a queue before the adjustment stations + adjustmentQueue <- newFCFSQueue + -- the inspection stations' queue size statistics + inspectionQueueSizeAcc <- + runEventInStartTime $ + newQueueSizeAccumulator inspectionQueue + -- the adjustment stations' queue size statistics + adjustmentQueueSizeAcc <- + runEventInStartTime $ + newQueueSizeAccumulator adjustmentQueue + -- create the inspection stations (servers) + inspectionStations <- + forM [1 .. inspectionStationCount] $ \_ -> + newInspectionStation + -- create the adjustment stations (servers) + adjustmentStations <- + forM [1 .. adjustmentStationCount] $ \_ -> + newAdjustmentStation + -- a processor loop for the inspection stations' queue + let inspectionQueueProcessorLoop = + queueProcessorLoopSeq + (liftEvent . enqueue inspectionQueue) + (dequeue inspectionQueue) + inspectionProcessor + (adjustmentQueueProcessor >>> adjustmentProcessor) + -- a processor for the adjustment stations' queue + let adjustmentQueueProcessor = + queueProcessor + (liftEvent . enqueue adjustmentQueue) + (dequeue adjustmentQueue) + -- a parallel work of the inspection stations + let inspectionProcessor = + processorParallel (map serverProcessor inspectionStations) + -- a parallel work of the adjustment stations + let adjustmentProcessor = + processorParallel (map serverProcessor adjustmentStations) + -- the entire processor from input to output + let entireProcessor = + arrivalTimerProcessor inputArrivalTimer >>> + inspectionQueueProcessorLoop >>> + arrivalTimerProcessor outputArrivalTimer + -- start simulating the model + runProcessInStartTime $ + sinkStream $ runProcessor entireProcessor inputStream + -- show the results in the final time + runEventInStopTime $ + do let indent = 2 + inspectionQueueSum <- queueSummary inspectionQueue indent + adjustmentQueueSum <- queueSummary adjustmentQueue indent + inspectionStationSums <- + forM inspectionStations $ \x -> serverSummary x indent + adjustmentStationSums <- + forM adjustmentStations $ \x -> serverSummary x indent + inputProcessingTime <- arrivalProcessingTime inputArrivalTimer + outputProcessingTime <- arrivalProcessingTime outputArrivalTimer + inspectionQueueSize <- timingStatsAccumulated inspectionQueueSizeAcc + adjustmentQueueSize <- timingStatsAccumulated adjustmentQueueSizeAcc + liftIO $ + do putStrLn "" + putStrLn "--- the inspection stations' queue summary (in the final time) ---" + putStrLn "" + putStrLn $ inspectionQueueSum [] + putStrLn "" + forM_ (zip [1..] inspectionStationSums) $ \(i, x) -> + do putStrLn $ "--- the inspection station no. " + ++ show i ++ " (in the final time) ---" + putStrLn "" + putStrLn $ x [] + putStrLn "" + putStrLn "--- the adjustment stations' queue summary (in the final time) ---" + putStrLn "" + putStrLn $ adjustmentQueueSum [] + putStrLn "" + forM_ (zip [1..] adjustmentStationSums) $ \(i, x) -> + do putStrLn $ "--- the adjustment station no. " + ++ show i ++ " (in the final time) ---" + putStrLn "" + putStrLn $ x [] + putStrLn "" + putStrLn "--- the input arrival time summary (we are interested in their count) ---" + putStrLn "" + putStrLn $ samplingStatsSummary inputProcessingTime indent [] + putStrLn "" + putStrLn "--- the arrival processing time summary ---" + putStrLn "" + putStrLn $ samplingStatsSummary outputProcessingTime indent [] + putStrLn "" + putStrLn $ "--- the inspection stations' queue size summary " + ++ "(updated when enqueueing and dequeueing) ---" + putStrLn "" + putStrLn $ timingStatsSummary inspectionQueueSize indent [] + putStrLn "" + putStrLn $ "--- the adjustment stations' queue size summary " + ++ "(updated when enqueueing and dequeueing) ---" + putStrLn "" + putStrLn $ timingStatsSummary adjustmentQueueSize indent [] + putStrLn "" + +main = runSimulation model specs diff --git a/examples/SimpleWorkflow.hs b/examples/SimpleWorkflow.hs deleted file mode 100644 --- a/examples/SimpleWorkflow.hs +++ /dev/null @@ -1,157 +0,0 @@ - --- This is a model of the workflow with originally two work places. Also there are two finite queues. --- --- It is described in different sources [1, 2]. So, this is chapter 7 of [2]. --- --- [1] { add a foreign source in English } --- --- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 - -import Prelude hiding (id, (.)) - -import Control.Monad -import Control.Monad.Trans -import Control.Arrow -import Control.Category (id, (.)) - -import Simulation.Aivika -import Simulation.Aivika.Queue - --- | The simulation specs. -specs = Specs { spcStartTime = 0.0, - spcStopTime = 300.0, - spcDT = 0.1, - spcMethod = RungeKutta4, - spcGeneratorType = SimpleGenerator } - --- the mean delay of the input arrivals distributed exponentially -meanOrderDelay = 0.4 - --- the capacity of the queue before the first work places -queueMaxCount1 = 4 - --- the capacity of the queue before the second work places -queueMaxCount2 = 2 - --- the mean processing time distributed exponentially in --- the first work places -meanProcessingTime1 = 0.25 - --- the mean processing time distributed exponentially in --- the second work places -meanProcessingTime2 = 0.5 - --- the number of the first work places --- (in parallel but the commented code allocates them sequentially) -workplaceCount1 = 1 - --- the number of the second work places --- (in parallel but the commented code allocates them sequentially) -workplaceCount2 = 1 - --- create an accumulator to gather the queue size statistics -newQueueSizeAccumulator queue = - newTimingStatsAccumulator $ - Signalable (queueCount queue) (queueCountChanged_ queue) - --- create a workflow with the exponential processing time -newWorkplaceExponential meanTime = - newServer $ \a -> - do holdProcess =<< - (liftParameter $ - randomExponential meanTime) - return a - --- interpose the prefetch processor between two processors -interposePrefetchProcessor x y = - x >>> prefetchProcessor >>> y - -model :: Simulation () -model = do - -- it will gather the statistics of the processing time - arrivalTimer <- newArrivalTimer - -- define a stream of input events - let inputStream = randomExponentialStream meanOrderDelay - -- create a queue before the first work place - queue1 <- newFCFSQueue queueMaxCount1 - -- create a queue before the second work place - queue2 <- newFCFSQueue queueMaxCount2 - -- the first queue size statistics - queueSizeAcc1 <- - runEventInStartTime $ - newQueueSizeAccumulator queue1 - -- the second queue size statistics - queueSizeAcc2 <- - runEventInStartTime $ - newQueueSizeAccumulator queue2 - -- create the first work places, i.e. the "servers" - workplace1s <- forM [1 .. workplaceCount1] $ \_ -> - newWorkplaceExponential meanProcessingTime1 - -- create the second work places, i.e. the "servers" - workplace2s <- forM [1 .. workplaceCount2] $ \_ -> - newWorkplaceExponential meanProcessingTime2 - -- processor for the queue before the first work place - let queueProcessor1 = - queueProcessor - (\a -> liftEvent $ enqueueOrLost_ queue1 a) - (dequeue queue1) - -- processor for the queue before the second work place - let queueProcessor2 = - queueProcessor - (enqueue queue2) - (dequeue queue2) - -- the entire processor from input to output - let entireProcessor = - queueProcessor1 >>> - processorParallel (map serverProcessor workplace1s) >>> - -- foldr1 interposePrefetchProcessor (map serverProcessor workplace1s) >>> - queueProcessor2 >>> - processorParallel (map serverProcessor workplace2s) >>> - -- foldr1 interposePrefetchProcessor (map serverProcessor workplace2s) >>> - arrivalTimerProcessor arrivalTimer - -- start simulating the model - runProcessInStartTime $ - sinkStream $ runProcessor entireProcessor inputStream - -- show the results in the final time - runEventInStopTime $ - do queueSum1 <- queueSummary queue1 2 - queueSum2 <- queueSummary queue2 2 - workplaceSum1s <- forM workplace1s $ \x -> serverSummary x 2 - workplaceSum2s <- forM workplace2s $ \x -> serverSummary x 2 - processingTime <- arrivalProcessingTime arrivalTimer - queueSize1 <- timingStatsAccumulated queueSizeAcc1 - queueSize2 <- timingStatsAccumulated queueSizeAcc2 - liftIO $ - do putStrLn "" - putStrLn "--- the first queue summary (in the final time) ---" - putStrLn "" - putStrLn $ queueSum1 [] - putStrLn "" - forM_ (zip [1..] workplaceSum1s) $ \(i, x) -> - do putStrLn $ "--- the first work place no." ++ show i ++ " (in the final time) ---" - putStrLn "" - putStrLn $ x [] - putStrLn "" - putStrLn "--- the second queue summary (in the final time) ---" - putStrLn "" - putStrLn $ queueSum2 [] - putStrLn "" - forM_ (zip [1..] workplaceSum2s) $ \(i, x) -> - do putStrLn $ "--- the second work place no. " ++ show i ++ " (in the final time) ---" - putStrLn "" - putStrLn $ x [] - putStrLn "" - putStrLn "--- the processing time summary ---" - putStrLn "" - putStrLn $ samplingStatsSummary processingTime 2 [] - putStrLn "" - putStrLn "--- the first queue size summary ---" - putStrLn "" - putStrLn $ timingStatsSummary queueSize1 2 [] - putStrLn "" - putStrLn "--- the second queue size summary ---" - putStrLn "" - putStrLn $ timingStatsSummary queueSize2 2 [] - putStrLn "" - -main = runSimulation model specs diff --git a/examples/WorkStationsInSeries.hs b/examples/WorkStationsInSeries.hs new file mode 100644 --- /dev/null +++ b/examples/WorkStationsInSeries.hs @@ -0,0 +1,159 @@ + +-- Example: Work Stations in Series +-- +-- This is a model of two work stations connected in a series and separated by finite queues. +-- +-- It is described in different sources [1, 2]. So, this is chapter 7 of [2] and section 5.14 of [1]. +-- +-- [1] A. Alan B. Pritsker, Simulation with Visual SLAM and AweSim, 2nd ed. +-- +-- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 + +import Prelude hiding (id, (.)) + +import Control.Monad +import Control.Monad.Trans +import Control.Arrow +import Control.Category (id, (.)) + +import Simulation.Aivika +import Simulation.Aivika.Queue + +-- | The simulation specs. +specs = Specs { spcStartTime = 0.0, + spcStopTime = 300.0, + spcDT = 0.1, + spcMethod = RungeKutta4, + spcGeneratorType = SimpleGenerator } + +-- the mean delay of the input arrivals distributed exponentially +meanOrderDelay = 0.4 + +-- the capacity of the queue before the first work places +queueMaxCount1 = 4 + +-- the capacity of the queue before the second work places +queueMaxCount2 = 2 + +-- the mean processing time distributed exponentially in +-- the first work stations +meanProcessingTime1 = 0.25 + +-- the mean processing time distributed exponentially in +-- the second work stations +meanProcessingTime2 = 0.5 + +-- the number of the first work stations +-- (in parallel but the commented code allocates them sequentially) +workStationCount1 = 1 + +-- the number of the second work stations +-- (in parallel but the commented code allocates them sequentially) +workStationCount2 = 1 + +-- create an accumulator to gather the queue size statistics +newQueueSizeAccumulator queue = + newTimingStatsAccumulator $ + Signalable (queueCount queue) (queueCountChanged_ queue) + +-- create a work station (server) with the exponential processing time +newWorkStationExponential meanTime = + newServer $ \a -> + do holdProcess =<< + (liftParameter $ + randomExponential meanTime) + return a + +-- interpose the prefetch processor between two processors +interposePrefetchProcessor x y = + x >>> prefetchProcessor >>> y + +model :: Simulation () +model = do + -- it will gather the statistics of the processing time + arrivalTimer <- newArrivalTimer + -- define a stream of input events + let inputStream = randomExponentialStream meanOrderDelay + -- create a queue before the first work stations + queue1 <- newFCFSQueue queueMaxCount1 + -- create a queue before the second work stations + queue2 <- newFCFSQueue queueMaxCount2 + -- the first queue size statistics + queueSizeAcc1 <- + runEventInStartTime $ + newQueueSizeAccumulator queue1 + -- the second queue size statistics + queueSizeAcc2 <- + runEventInStartTime $ + newQueueSizeAccumulator queue2 + -- create the first work stations (servers) + workStation1s <- forM [1 .. workStationCount1] $ \_ -> + newWorkStationExponential meanProcessingTime1 + -- create the second work stations (servers) + workStation2s <- forM [1 .. workStationCount2] $ \_ -> + newWorkStationExponential meanProcessingTime2 + -- processor for the queue before the first work station + let queueProcessor1 = + queueProcessor + (\a -> liftEvent $ enqueueOrLost_ queue1 a) + (dequeue queue1) + -- processor for the queue before the second work station + let queueProcessor2 = + queueProcessor + (enqueue queue2) + (dequeue queue2) + -- the entire processor from input to output + let entireProcessor = + queueProcessor1 >>> + processorParallel (map serverProcessor workStation1s) >>> + -- foldr1 interposePrefetchProcessor (map serverProcessor workStation1s) >>> + queueProcessor2 >>> + processorParallel (map serverProcessor workStation2s) >>> + -- foldr1 interposePrefetchProcessor (map serverProcessor workStation2s) >>> + arrivalTimerProcessor arrivalTimer + -- start simulating the model + runProcessInStartTime $ + sinkStream $ runProcessor entireProcessor inputStream + -- show the results in the final time + runEventInStopTime $ + do queueSum1 <- queueSummary queue1 2 + queueSum2 <- queueSummary queue2 2 + workStationSum1s <- forM workStation1s $ \x -> serverSummary x 2 + workStationSum2s <- forM workStation2s $ \x -> serverSummary x 2 + processingTime <- arrivalProcessingTime arrivalTimer + queueSize1 <- timingStatsAccumulated queueSizeAcc1 + queueSize2 <- timingStatsAccumulated queueSizeAcc2 + liftIO $ + do putStrLn "" + putStrLn "--- the first queue summary (in the final time) ---" + putStrLn "" + putStrLn $ queueSum1 [] + putStrLn "" + forM_ (zip [1..] workStationSum1s) $ \(i, x) -> + do putStrLn $ "--- the first work station no. " ++ show i ++ " (in the final time) ---" + putStrLn "" + putStrLn $ x [] + putStrLn "" + putStrLn "--- the second queue summary (in the final time) ---" + putStrLn "" + putStrLn $ queueSum2 [] + putStrLn "" + forM_ (zip [1..] workStationSum2s) $ \(i, x) -> + do putStrLn $ "--- the second work station no. " ++ show i ++ " (in the final time) ---" + putStrLn "" + putStrLn $ x [] + putStrLn "" + putStrLn "--- the processing time summary ---" + putStrLn "" + putStrLn $ samplingStatsSummary processingTime 2 [] + putStrLn "" + putStrLn "--- the first queue size summary ---" + putStrLn "" + putStrLn $ timingStatsSummary queueSize1 2 [] + putStrLn "" + putStrLn "--- the second queue size summary ---" + putStrLn "" + putStrLn $ timingStatsSummary queueSize2 2 [] + putStrLn "" + +main = runSimulation model specs diff --git a/examples/WorkflowLoop.hs b/examples/WorkflowLoop.hs deleted file mode 100644 --- a/examples/WorkflowLoop.hs +++ /dev/null @@ -1,191 +0,0 @@ - -{-# LANGUAGE RecursiveDo, Arrows #-} - --- This is a model of the workflow with a loop. Also there are two infinite queues. --- --- It is described in different sources [1, 2]. So, this is chapter 8 of [2]. --- --- [1] { add a foreign source in English } --- --- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 - --- CAUTION: this is model is not yet fully tested and it may contain logical errors. - -import Prelude hiding (id, (.)) - -import Control.Monad -import Control.Monad.Trans -import Control.Arrow -import Control.Category (id, (.)) - -import Simulation.Aivika -import Simulation.Aivika.Queue.Infinite - --- | The simulation specs. -specs = Specs { spcStartTime = 0.0, - spcStopTime = 480.0, - spcDT = 0.1, - spcMethod = RungeKutta4, - spcGeneratorType = SimpleGenerator } - --- the minimum delay for arriving the next TV set -minArrivalDelay = 3.5 - --- the maximum delay for arriving the next TV set -maxArrivalDelay = 7.5 - --- the minimum test time -minTestTime = 6 - --- the maximum test time -maxTestTime = 12 - --- the probability of passing the test -testPassingProb = 0.85 - --- how many testers are there? -testerWorkplaceCount = 2 - --- the minimum time of tuning the TV set --- that has not passed the test -minTuningTime = 20 - --- the maximum time of tuning the TV set --- that has not passed the test -maxTuningTime = 40 - --- how many persons perform a tuning of TV sets? -tunerWorkplaceCount = 1 - --- create an accumulator to gather the queue size statistics -newQueueSizeAccumulator queue = - newTimingStatsAccumulator $ - Signalable (queueCount queue) (queueCountChanged_ queue) - --- create a tester's workplace -newTesterWorkplace = - newServer $ \a -> - do holdProcess =<< - (liftParameter $ - randomUniform minTestTime maxTestTime) - passed <- - liftParameter $ - randomTrue testPassingProb - if passed - then return $ Right a - else return $ Left a - --- create a tuner's workplace -newTunerWorkplace = - newServer $ \a -> - do holdProcess =<< - (liftParameter $ - randomUniform minTuningTime maxTuningTime) - return a - -model :: Simulation () -model = mdo - -- it will gather the statistics of the processing time - inputArrivalTimer <- newArrivalTimer - outputArrivalTimer <- newArrivalTimer - -- define a stream of input events - let inputStream = - randomUniformStream minArrivalDelay maxArrivalDelay - -- create a queue before the tester's work place - testerQueue <- newFCFSQueue - -- create a queue before the tuner's work place - tunerQueue <- newFCFSQueue - -- the tester's queue size statistics - testerQueueSizeAcc <- - runEventInStartTime $ - newQueueSizeAccumulator testerQueue - -- the tuner's queue size statistics - tunerQueueSizeAcc <- - runEventInStartTime $ - newQueueSizeAccumulator tunerQueue - -- create the tester's work places, i.e. the "servers" - testerWorkplaces <- - forM [1 .. testerWorkplaceCount] $ \_ -> - newTesterWorkplace - -- create the tuner's work places, i.e. the "servers" - tunerWorkplaces <- - forM [1 .. tunerWorkplaceCount] $ \_ -> - newTunerWorkplace - -- a processor loop for the tester's queue - let testerQueueProcessorLoop = - queueProcessorLoopSeq - (liftEvent . enqueue testerQueue) - (dequeue testerQueue) - testerProcessor - (tunerQueueProcessor >>> tunerProcessor) - -- a processor for the tuner's queue - let tunerQueueProcessor = - queueProcessor - (liftEvent . enqueue tunerQueue) - (dequeue tunerQueue) - -- the parallel work of all the testers - let testerProcessor = - processorParallel (map serverProcessor testerWorkplaces) - -- the parallel work of all the tuners - let tunerProcessor = - processorParallel (map serverProcessor tunerWorkplaces) - -- the entire processor from input to output - let entireProcessor = - arrivalTimerProcessor inputArrivalTimer >>> - testerQueueProcessorLoop >>> - arrivalTimerProcessor outputArrivalTimer - -- start simulating the model - runProcessInStartTime $ - sinkStream $ runProcessor entireProcessor inputStream - -- show the results in the final time - runEventInStopTime $ - do testerQueueSum <- queueSummary testerQueue 2 - tunerQueueSum <- queueSummary tunerQueue 2 - testerWorkplaceSums <- - forM testerWorkplaces $ \x -> serverSummary x 2 - tunerWorkplaceSums <- - forM tunerWorkplaces $ \x -> serverSummary x 2 - inputProcessingTime <- arrivalProcessingTime inputArrivalTimer - outputProcessingTime <- arrivalProcessingTime outputArrivalTimer - testerQueueSize <- timingStatsAccumulated testerQueueSizeAcc - tunerQueueSize <- timingStatsAccumulated tunerQueueSizeAcc - liftIO $ - do putStrLn "" - putStrLn "--- the tester's queue summary (in the final time) ---" - putStrLn "" - putStrLn $ testerQueueSum [] - putStrLn "" - forM_ (zip [1..] testerWorkplaceSums) $ \(i, x) -> - do putStrLn $ "--- the tester's work place no." - ++ show i ++ " (in the final time) ---" - putStrLn "" - putStrLn $ x [] - putStrLn "" - putStrLn "--- the tuner's queue summary (in the final time) ---" - putStrLn "" - putStrLn $ tunerQueueSum [] - putStrLn "" - forM_ (zip [1..] tunerWorkplaceSums) $ \(i, x) -> - do putStrLn $ "--- the tuner's work place no. " - ++ show i ++ " (in the final time) ---" - putStrLn "" - putStrLn $ x [] - putStrLn "" - putStrLn "--- the arrival receiving time summary (we are interested in their count) ---" - putStrLn "" - putStrLn $ samplingStatsSummary inputProcessingTime 2 [] - putStrLn "" - putStrLn "--- the arrival processing time summary ---" - putStrLn "" - putStrLn $ samplingStatsSummary outputProcessingTime 2 [] - putStrLn "" - putStrLn "--- the tester's queue size summary (updated when enqueueing and dequeueing) ---" - putStrLn "" - putStrLn $ timingStatsSummary testerQueueSize 2 [] - putStrLn "" - putStrLn "--- the tuner's queue size summary (updated when enqueueing and dequeueing) ---" - putStrLn "" - putStrLn $ timingStatsSummary tunerQueueSize 2 [] - putStrLn "" - -main = runSimulation model specs