aivika 0.6.1 → 0.7

raw patch · 74 files changed

+5433/−5593 lines, 74 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Simulation/Aivika/Agent.hs +252/−0
• Simulation/Aivika/Cont.hs +18/−0
• Simulation/Aivika/DoubleLinkedList.hs +165/−0
• Simulation/Aivika/Dynamics.hs +13/−6
• Simulation/Aivika/Dynamics/Agent.hs +0/−343
• Simulation/Aivika/Dynamics/Base.hs +0/−46
• Simulation/Aivika/Dynamics/Buffer.hs +0/−189
• Simulation/Aivika/Dynamics/Cont.hs +0/−18
• Simulation/Aivika/Dynamics/EventQueue.hs +0/−215
• Simulation/Aivika/Dynamics/FIFO.hs +0/−207
• Simulation/Aivika/Dynamics/Fold.hs +82/−0
• Simulation/Aivika/Dynamics/Internal/Cont.hs +0/−304
• Simulation/Aivika/Dynamics/Internal/Dynamics.hs +0/−310
• Simulation/Aivika/Dynamics/Internal/Fold.hs +0/−92
• Simulation/Aivika/Dynamics/Internal/Interpolate.hs +0/−62
• Simulation/Aivika/Dynamics/Internal/Memo.hs +0/−196
• Simulation/Aivika/Dynamics/Internal/Process.hs +0/−300
• Simulation/Aivika/Dynamics/Internal/Signal.hs +0/−347
• Simulation/Aivika/Dynamics/Internal/Simulation.hs +0/−202
• Simulation/Aivika/Dynamics/Internal/Time.hs +0/−59
• Simulation/Aivika/Dynamics/Interpolate.hs +63/−0
• Simulation/Aivika/Dynamics/LIFO.hs +0/−196
• Simulation/Aivika/Dynamics/Memo.hs +125/−0
• Simulation/Aivika/Dynamics/Memo/Unboxed.hs +100/−0
• Simulation/Aivika/Dynamics/Parameter.hs +0/−83
• Simulation/Aivika/Dynamics/Process.hs +0/−42
• Simulation/Aivika/Dynamics/Random.hs +18/−54
• Simulation/Aivika/Dynamics/Ref.hs +0/−89
• Simulation/Aivika/Dynamics/Resource.hs +0/−171
• Simulation/Aivika/Dynamics/Signal.hs +0/−184
• Simulation/Aivika/Dynamics/Simulation.hs +0/−29
• Simulation/Aivika/Dynamics/SystemDynamics.hs +0/−788
• Simulation/Aivika/Dynamics/UVar.hs +0/−164
• Simulation/Aivika/Dynamics/Var.hs +0/−163
• Simulation/Aivika/Event.hs +36/−0
• Simulation/Aivika/Internal/Cont.hs +320/−0
• Simulation/Aivika/Internal/Dynamics.hs +217/−0
• Simulation/Aivika/Internal/Event.hs +349/−0
• Simulation/Aivika/Internal/Process.hs +325/−0
• Simulation/Aivika/Internal/Signal.hs +260/−0
• Simulation/Aivika/Internal/Simulation.hs +193/−0
• Simulation/Aivika/Internal/Specs.hs +199/−0
• Simulation/Aivika/Parameter.hs +59/−0
• Simulation/Aivika/Parameter/Random.hs +41/−0
• Simulation/Aivika/Process.hs +47/−0
• Simulation/Aivika/Queue.hs +294/−85
• Simulation/Aivika/QueueStrategy.hs +154/−0
• Simulation/Aivika/Random.hs +53/−0
• Simulation/Aivika/Ref.hs +69/−0
• Simulation/Aivika/Resource.hs +270/−0
• Simulation/Aivika/Signal.hs +137/−0
• Simulation/Aivika/Simulation.hs +27/−0
• Simulation/Aivika/Specs.hs +25/−0
• Simulation/Aivika/SystemDynamics.hs +620/−0
• Simulation/Aivika/UVector.hs +0/−188
• Simulation/Aivika/Unboxed.hs +43/−0
• Simulation/Aivika/Var.hs +158/−0
• Simulation/Aivika/Var/Unboxed.hs +158/−0
• Simulation/Aivika/Vector.hs +1/−0
• Simulation/Aivika/Vector/Unboxed.hs +186/−0
• aivika.cabal +44/−36
• examples/BassDiffusion.hs +18/−17
• examples/ChemicalReaction.hs +12/−15
• examples/ChemicalReactionRec.hs +0/−22
• examples/FishBank.hs +37/−43
• examples/FishBankRec.hs +0/−51
• examples/Furnace.hs +89/−117
• examples/MachRep1.hs +17/−17
• examples/MachRep1EventDriven.hs +19/−21
• examples/MachRep1TimeDriven.hs +21/−21
• examples/MachRep2.hs +24/−24
• examples/MachRep3.hs +27/−27
• examples/TimeOut.hs +27/−28
• examples/TimeOutInt.hs +21/−22

@@ -0,0 +1,252 @@++-- |+-- Module     : Simulation.Aivika.Agent+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- This module introduces basic entities for the agent-based modeling.+--+module Simulation.Aivika.Agent+       (Agent,+        AgentState,+        newAgent,+        newState,+        newSubstate,+        agentState,+        agentStateChanged,+        agentStateChanged_,+        activateState,+        stateAgent,+        stateParent,+        addTimeout,+        addTimer,+        setStateActivation,+        setStateDeactivation,+        setStateTransition) where++import Data.IORef+import Control.Monad++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Event+import Simulation.Aivika.Internal.Signal++--+-- Agent-based Modeling+--++-- | Represents an agent.+data Agent = Agent { agentModeRef            :: IORef AgentMode,+                     agentStateRef           :: IORef (Maybe AgentState), +                     agentStateChangedSource :: SignalSource (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 (Event ()),+                               stateDeactivateRef :: IORef (Event ()),+                               stateTransitRef    :: IORef (Event (Maybe AgentState)),+                               stateVersionRef    :: IORef Int }+                  +data AgentMode = CreationMode+               | 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++fullPath :: AgentState -> [AgentState] -> [AgentState]+fullPath st acc =+  case stateParent st of+    Nothing  -> st : acc+    Just st' -> fullPath st' (st : acc)++partitionPath :: [AgentState] -> [AgentState] -> ([AgentState], [AgentState])+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)++findPath :: Maybe AgentState -> AgentState -> ([AgentState], [AgentState])+findPath Nothing target = ([], fullPath target [])+findPath (Just source) target+  | stateAgent source /= stateAgent target =+    error "Different agents: findPath."+  | otherwise =+    partitionPath path1 path2+  where+    path1 = fullPath source []+    path2 = fullPath target []++traversePath :: Maybe AgentState -> AgentState -> Event ()+traversePath source target =+  let (path1, path2) = findPath source target+      agent = stateAgent target+      activate st p   = invokeEvent p =<< readIORef (stateActivateRef st)+      deactivate st p = invokeEvent p =<< readIORef (stateDeactivateRef st)+      transit st p    = invokeEvent p =<< readIORef (stateTransitRef st)+      continue st p   = invokeEvent p $ traversePath (Just target) st+  in Event $ \p ->+       unless (null path1 && null path2) $+       do writeIORef (agentModeRef agent) TransientMode+          forM_ path1 $ \st ->+            do writeIORef (agentStateRef agent) (Just st)+               deactivate st p+               -- it makes all timeout and timer handlers outdated+               modifyIORef (stateVersionRef st) (1 +)+          forM_ path2 $ \st ->+            do writeIORef (agentStateRef agent) (Just st)+               activate st p+          st' <- transit target p+          case st' of+            Nothing ->+              do writeIORef (agentModeRef agent) ProcessingMode+                 triggerAgentStateChanged p agent+            Just st' ->+              continue st' p++-- | Add to the state a timeout handler that will be actuated +-- in the specified time period, while the state remains active.+addTimeout :: AgentState -> Double -> Event () -> Event ()+addTimeout st dt action =+  Event $ \p ->+  do v <- readIORef (stateVersionRef st)+     let m1 = Event $ \p ->+           do v' <- readIORef (stateVersionRef st)+              when (v == v') $+                invokeEvent p action+         m2 = enqueueEvent (pointTime p + dt) m1+     invokeEvent p m2++-- | 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 -> Event Double -> Event () -> Event ()+addTimer st dt action =+  Event $ \p ->+  do v <- readIORef (stateVersionRef st)+     let m1 = Event $ \p ->+           do v' <- readIORef (stateVersionRef st)+              when (v == v') $+                do invokeEvent p m2+                   invokeEvent p action+         m2 = Event $ \p ->+           do dt' <- invokeEvent p dt+              invokeEvent p $ enqueueEvent (pointTime p + dt') m1+     invokeEvent p m2++-- | Create a new state.+newState :: Agent -> Simulation AgentState+newState agent =+  Simulation $ \r ->+  do aref <- newIORef $ return ()+     dref <- newIORef $ return ()+     tref <- newIORef $ return Nothing+     vref <- newIORef 0+     return AgentState { stateAgent = agent,+                         stateParent = Nothing,+                         stateActivateRef = aref,+                         stateDeactivateRef = dref,+                         stateTransitRef = tref,+                         stateVersionRef = vref }++-- | Create a child state.+newSubstate :: AgentState -> Simulation AgentState+newSubstate parent =+  Simulation $ \r ->+  do let agent = stateAgent parent +     aref <- newIORef $ return ()+     dref <- newIORef $ return ()+     tref <- newIORef $ return Nothing+     vref <- newIORef 0+     return AgentState { stateAgent = agent,+                         stateParent = Just parent,+                         stateActivateRef= aref,+                         stateDeactivateRef = dref,+                         stateTransitRef = tref,+                         stateVersionRef = vref }++-- | Create an agent.+newAgent :: Simulation Agent+newAgent =+  Simulation $ \r ->+  do modeRef  <- newIORef CreationMode+     stateRef <- newIORef Nothing+     stateChangedSource <- invokeSimulation r newSignalSource+     return Agent { agentModeRef = modeRef,+                    agentStateRef = stateRef, +                    agentStateChangedSource = stateChangedSource }++-- | Return the selected downmost active state.+agentState :: Agent -> Event (Maybe AgentState)+agentState agent =+  Event $ \p -> readIORef (agentStateRef agent)+                   +-- | Select the next downmost active state. The activation is repeated while+-- there is the transition state defined by 'setStateTransition'.+activateState :: AgentState -> Event ()+activateState st =+  Event $ \p ->+  do let agent = stateAgent st+     mode <- readIORef (agentModeRef agent)+     case mode of+       CreationMode ->+         do x0 <- readIORef (agentStateRef agent)+            invokeEvent p $ traversePath x0 st+       TransientMode ->+         error $+         "Use the setStateTransition function to define " +++         "the transition state: activateState."+       ProcessingMode ->+         do x0 @ (Just st0) <- readIORef (agentStateRef agent)+            invokeEvent p $ traversePath x0 st++-- | Set the activation computation for the specified state.+setStateActivation :: AgentState -> Event () -> Simulation ()+setStateActivation st action =+  Simulation $ \r ->+  writeIORef (stateActivateRef st) action+  +-- | Set the deactivation computation for the specified state.+setStateDeactivation :: AgentState -> Event () -> Simulation ()+setStateDeactivation st action =+  Simulation $ \r ->+  writeIORef (stateDeactivateRef st) action+  +-- | Set the transition state which will be next and which is used only+-- when activating the state directly with help of 'activateState'.+-- If the state was activated intermediately, when activating directly+-- another state, then this computation is not used.+setStateTransition :: AgentState -> Event (Maybe AgentState) -> Simulation ()+setStateTransition st action =+  Simulation $ \r ->+  writeIORef (stateTransitRef st) action+  +-- | Trigger the signal when the agent state changes.+triggerAgentStateChanged :: Point -> Agent -> IO ()+triggerAgentStateChanged p agent =+  do st <- readIORef (agentStateRef agent)+     invokeEvent p $ triggerSignal (agentStateChangedSource agent) st++-- | Return a signal that notifies about every change of the state.+agentStateChanged :: Agent -> Signal (Maybe AgentState)+agentStateChanged agent =+  publishSignal (agentStateChangedSource agent)++-- | Return a signal that notifies about every change of the state.+agentStateChanged_ :: Agent -> Signal ()+agentStateChanged_ agent =+  mapSignal (const ()) $ agentStateChanged agent

@@ -0,0 +1,18 @@++-- |+-- Module     : Simulation.Aivika.Cont+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- The 'Cont' monad is a variation of the standard Cont monad +-- and F# async workflow, where the result of applying +-- the continuations is the 'Event' computation.+--+module Simulation.Aivika.Cont+       (Cont) where++import Simulation.Aivika.Internal.Event+import Simulation.Aivika.Internal.Cont

@@ -0,0 +1,165 @@++-- |+-- Module     : Simulation.Aivika.DoubleLinkedList+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- An imperative double-linked list.+--+module Simulation.Aivika.DoubleLinkedList +       (DoubleLinkedList, +        listNull, +        listCount,+        newList, +        listInsertFirst,+        listAddLast,+        listRemoveFirst,+        listRemoveLast,+        listFirst,+        listLast) where ++import Data.IORef+import Control.Monad++-- | A cell of the double-linked list.+data DoubleLinkedItem a = +  DoubleLinkedItem { itemVal  :: a,+                     itemPrev :: IORef (Maybe (DoubleLinkedItem a)),+                     itemNext :: IORef (Maybe (DoubleLinkedItem a)) }+  +-- | The 'DoubleLinkedList' type represents an imperative double-linked list.+data DoubleLinkedList a =  +  DoubleLinkedList { listHead :: IORef (Maybe (DoubleLinkedItem a)),+                     listTail :: IORef (Maybe (DoubleLinkedItem a)), +                     listSize :: IORef Int }++-- | Test whether the list is empty.+listNull :: DoubleLinkedList a -> IO Bool+listNull x =+  do head <- readIORef (listHead x) +     case head of+       Nothing -> return True+       Just _  -> return False+    +-- | Return the number of elements in the list.+listCount :: DoubleLinkedList a -> IO Int+listCount x = readIORef (listSize x)++-- | Create a new list.+newList :: IO (DoubleLinkedList a)+newList =+  do head <- newIORef Nothing +     tail <- newIORef Nothing+     size <- newIORef 0+     return DoubleLinkedList { listHead = head,+                               listTail = tail,+                               listSize = size }++-- | Insert a new element in the beginning.+listInsertFirst :: DoubleLinkedList a -> a -> IO ()+listInsertFirst x v =+  do size <- readIORef (listSize x)+     writeIORef (listSize x) (size + 1)+     head <- readIORef (listHead x)+     case head of+       Nothing ->+         do prev <- newIORef Nothing+            next <- newIORef Nothing+            let item = Just DoubleLinkedItem { itemVal = v, +                                               itemPrev = prev, +                                               itemNext = next }+            writeIORef (listHead x) item+            writeIORef (listTail x) item+       Just h ->+         do prev <- newIORef Nothing+            next <- newIORef head+            let item = Just DoubleLinkedItem { itemVal = v,+                                               itemPrev = prev,+                                               itemNext = next }+            writeIORef (itemPrev h) item+            writeIORef (listHead x) item++-- | Add a new element to the end.+listAddLast :: DoubleLinkedList a -> a -> IO ()+listAddLast x v =+  do size <- readIORef (listSize x)+     writeIORef (listSize x) (size + 1)+     tail <- readIORef (listTail x)+     case tail of+       Nothing ->+         do prev <- newIORef Nothing+            next <- newIORef Nothing+            let item = Just DoubleLinkedItem { itemVal = v, +                                               itemPrev = prev, +                                               itemNext = next }+            writeIORef (listHead x) item+            writeIORef (listTail x) item+       Just t ->+         do prev <- newIORef tail+            next <- newIORef Nothing+            let item = Just DoubleLinkedItem { itemVal = v,+                                               itemPrev = prev,+                                               itemNext = next }+            writeIORef (itemNext t) item+            writeIORef (listTail x) item++-- | Remove the first element.+listRemoveFirst :: DoubleLinkedList a -> IO ()+listRemoveFirst x =+  do head <- readIORef (listHead x) +     case head of+       Nothing ->+         error "Empty list: listRemoveFirst"+       Just h ->+         do size  <- readIORef (listSize x)+            writeIORef (listSize x) (size - 1)+            head' <- readIORef (itemNext h)+            case head' of+              Nothing ->+                do writeIORef (listHead x) Nothing+                   writeIORef (listTail x) Nothing+              Just h' ->+                do writeIORef (itemPrev h') Nothing+                   writeIORef (listHead x) head'++-- | Remove the last element.+listRemoveLast :: DoubleLinkedList a -> IO ()+listRemoveLast x =+  do tail <- readIORef (listTail x) +     case tail of+       Nothing ->+         error "Empty list: listRemoveLast"+       Just t ->+         do size  <- readIORef (listSize x)+            writeIORef (listSize x) (size - 1)+            tail' <- readIORef (itemPrev t)+            case tail' of+              Nothing ->+                do writeIORef (listHead x) Nothing+                   writeIORef (listTail x) Nothing+              Just t' ->+                do writeIORef (itemNext t') Nothing+                   writeIORef (listTail x) tail'++-- | Return the first element.+listFirst :: DoubleLinkedList a -> IO a+listFirst x =+  do head <- readIORef (listHead x)+     case head of+       Nothing ->+         error "Empty list: listFirst"+       Just h ->+         return $ itemVal h++-- | Return the last element.+listLast :: DoubleLinkedList a -> IO a+listLast x =+  do tail <- readIORef (listTail x)+     case tail of+       Nothing ->+         error "Empty list: listLast"+       Just t ->+         return $ itemVal t

@@ -7,11 +7,10 @@ -- Stability  : experimental -- Tested with: GHC 7.6.3 ----- The module defines the 'Dynamics' monad representing an abstract dynamic --- process, i.e. a time varying polymorphic function. +-- The module defines the 'Dynamics' monad representing a time varying polymorphic function.  ---module Simulation.Aivika.Dynamics -       (-- * Dynamics+module Simulation.Aivika.Dynamics+       (-- * Dynamics Monad         Dynamics,         DynamicsLift(..),         runDynamicsInStartTime,@@ -22,6 +21,14 @@         -- * Error Handling         catchDynamics,         finallyDynamics,-        throwDynamics) where+        throwDynamics,+        -- * Time parameters+        starttime,+        stoptime,+        dt,+        time,+        isTimeInteg,+        integIteration,+        integPhase) where -import Simulation.Aivika.Dynamics.Internal.Dynamics+import Simulation.Aivika.Internal.Dynamics

@@ -1,343 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Agent--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module introduces basic entities for the agent-based modeling.------ WARNING: the module is not well tested. This caution is related mainly to--- managing the nested states.--- --- At the same time, the timer and timeout handlers seem to be well tested as--- they are just light-weight wrappers creating the event handlers that are--- already processed by the event queue.-----module Simulation.Aivika.Dynamics.Agent-       (Agent,-        AgentState,-        newAgent,-        newState,-        newSubstate,-        agentQueue,-        agentState,-        agentStateChanged,-        agentStateChanged_,-        activateState,-        initState,-        stateAgent,-        stateParent,-        addTimeout,-        addTimer,-        stateActivation,-        stateDeactivation,-        setStateActivation,-        setStateDeactivation,-        setStateTransition) where--import Data.IORef-import Control.Monad--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal------- Agent-based Modeling------- | Represents an agent.-data Agent = Agent { agentQueue :: EventQueue,-                     -- ^ Return the bound event queue.-                     agentModeRef :: IORef AgentMode,-                     agentStateRef :: IORef (Maybe AgentState), -                     agentStateChangedSource :: SignalSource (Maybe AgentState), -                     agentStateUpdatedSource :: SignalSource (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 ()),-                               stateTransitRef :: IORef (Dynamics (Maybe AgentState)),-                               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--fullPath :: AgentState -> [AgentState] -> [AgentState]-fullPath st acc =-  case stateParent st of-    Nothing  -> st : acc-    Just st' -> fullPath st' (st : acc)--partitionPath :: [AgentState] -> [AgentState] -> ([AgentState], [AgentState])-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)--findPath :: Maybe AgentState -> AgentState -> ([AgentState], [AgentState])-findPath Nothing target = ([], fullPath target [])-findPath (Just source) target-  | stateAgent source /= stateAgent target =-    error "Different agents: findPath."-  | otherwise =-    partitionPath path1 path2-  where-    path1 = fullPath source []-    path2 = fullPath target []--traversePath :: Maybe AgentState -> AgentState -> Dynamics ()-traversePath source target =-  let (path1, path2) = findPath source target-      agent = stateAgent target-      activate st p =-        do Dynamics m <- readIORef (stateActivateRef st)-           m p-      deactivate st p =-        do Dynamics m <- readIORef (stateDeactivateRef st)-           m p-      transit st p =-        do Dynamics m <- readIORef (stateTransitRef st)-           m p-      continue st p =-        do let Dynamics m = traversePath (Just target) st-           m p-  in Dynamics $ \p ->-       unless (null path1 && null path2) $-       do writeIORef (agentModeRef agent) TransientMode-          forM_ path1 $ \st ->-            do writeIORef (agentStateRef agent) (Just st)-               deactivate st p-               -- it makes all timeout and timer handlers outdated-               modifyIORef (stateVersionRef st) (1 +)-          forM_ path2 $ \st ->-            do when (st == target) $-                 writeIORef (agentModeRef agent) InitialMode-               writeIORef (agentStateRef agent) (Just st)-               activate st p-          writeIORef (agentModeRef agent) TransientMode     -          st' <- transit target p-          case st' of-            Nothing ->-              do writeIORef (agentModeRef agent) ProcessingMode-                 triggerAgentStateChanged p agent-            Just st' ->-              continue st' p---- | 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 let q = agentQueue (stateAgent st)-         Dynamics m0 = runQueueSync q-     m0 p    -- ensure that the agent state is actual-     v <- readIORef (stateVersionRef st)-     let m1 = Dynamics $ \p ->-           do -- checkTime p (stateAgent st) "addTimeout"-              v' <- readIORef (stateVersionRef st)-              when (v == v') $ action p-         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 let q = agentQueue (stateAgent st)-         Dynamics m0 = runQueueSync q-     m0 p    -- ensure that the agent state is actual-     v <- readIORef (stateVersionRef st)-     let m1 = Dynamics $ \p ->-           do -- checkTime p (stateAgent st) "addTimer"-              v' <- readIORef (stateVersionRef st)-              when (v == v') $ do { m2 p; action p }-         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 -> Simulation AgentState-newState agent =-  Simulation $ \r ->-  do aref <- newIORef $ return ()-     dref <- newIORef $ return ()-     tref <- newIORef $ return Nothing-     vref <- newIORef 0-     return AgentState { stateAgent = agent,-                         stateParent = Nothing,-                         stateActivateRef = aref,-                         stateDeactivateRef = dref,-                         stateTransitRef = tref,-                         stateVersionRef = vref }---- | Create a child state.-newSubstate :: AgentState -> Simulation AgentState-newSubstate parent =-  Simulation $ \r ->-  do let agent = stateAgent parent -     aref <- newIORef $ return ()-     dref <- newIORef $ return ()-     tref <- newIORef $ return Nothing-     vref <- newIORef 0-     return AgentState { stateAgent = agent,-                         stateParent = Just parent,-                         stateActivateRef= aref,-                         stateDeactivateRef = dref,-                         stateTransitRef = tref,-                         stateVersionRef = vref }---- | Create an agent bound with the specified event queue.-newAgent :: EventQueue -> Simulation Agent-newAgent queue =-  Simulation $ \r ->-  do modeRef    <- newIORef CreationMode-     stateRef   <- newIORef Nothing-     let Simulation m1 = newSignalSourceUnsafe-         Simulation m2 = newSignalSource queue-     stateChangedSource <- m1 r-     stateUpdatedSource <- m2 r-     return Agent { agentQueue = queue,-                    agentModeRef = modeRef,-                    agentStateRef = stateRef, -                    agentStateChangedSource = stateChangedSource, -                    agentStateUpdatedSource = stateUpdatedSource }---- | Return the selected downmost active state.-agentState :: Agent -> Dynamics (Maybe AgentState)-agentState agent =-  Dynamics $ \p -> -  do let Dynamics m = runQueueSync $ agentQueue agent -     m p    -- ensure that the agent state is actual-     readIORef (agentStateRef agent)-                   --- | Select the next downmost active state. The activation is repeated while--- there is the transition state defined by 'setStateTransition'.-activateState :: AgentState -> Dynamics ()-activateState st =-  Dynamics $ \p ->-  do let agent = stateAgent st-         Dynamics m = runQueueSync $ agentQueue agent -     m p    -- ensure that the agent state is actual-     mode <- readIORef (agentModeRef agent)-     case mode of-       CreationMode ->-         do x0 <- readIORef (agentStateRef agent)-            let Dynamics m = traversePath x0 st-            m p-       InitialMode ->-         error $ -         "Use the setStateTransition function to define " ++-         "the transition state: activateState."-       TransientMode ->-         error $-         "Use the setStateTransition function to define " ++-         "the transition state: activateState."-       ProcessingMode ->-         do x0 @ (Just st0) <- readIORef (agentStateRef agent)-            let Dynamics m = traversePath x0 st-            m p--{-# DEPRECATED initState "Rewrite using the setStateTransition function instead." #-}-              --- | 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 = runQueueSync $ 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 x0 @ (Just st0) <- readIORef (agentStateRef agent)-            let Dynamics m = traversePath x0 st-            m p-       TransientMode -> -         return ()-       ProcessingMode ->-         error $-         "Use the activateState function everywhere outside " ++-         "the state activation: initState."--{-# DEPRECATED stateActivation "Use the setStateActivation function instead" #-}-{-# DEPRECATED stateDeactivation "Use the setStateDeactivation function instead" #-}---- | Set the activation computation for the specified state.-stateActivation :: AgentState -> Dynamics () -> Simulation ()-stateActivation = setStateActivation-  --- | Set the deactivation computation for the specified state.-stateDeactivation :: AgentState -> Dynamics () -> Simulation ()-stateDeactivation = setStateDeactivation-  --- | Set the activation computation for the specified state.-setStateActivation :: AgentState -> Dynamics () -> Simulation ()-setStateActivation st action =-  Simulation $ \r ->-  writeIORef (stateActivateRef st) action-  --- | Set the deactivation computation for the specified state.-setStateDeactivation :: AgentState -> Dynamics () -> Simulation ()-setStateDeactivation st action =-  Simulation $ \r ->-  writeIORef (stateDeactivateRef st) action-  --- | Set the transition state which will be next and which is used only--- when activating the state directly with help of 'activateState'.--- If the state was activated intermediately, when activating directly--- another state, then this computation is not used.-setStateTransition :: AgentState -> Dynamics (Maybe AgentState) -> Simulation ()-setStateTransition st action =-  Simulation $ \r ->-  writeIORef (stateTransitRef st) action-  --- | Trigger the signal when the agent state changes.-triggerAgentStateChanged :: Point -> Agent -> IO ()-triggerAgentStateChanged p agent =-  do st <- readIORef (agentStateRef agent)-     let Dynamics m = triggerSignal (agentStateChangedSource agent) st-     m p---- | Return a signal that notifies about every change of the state.-agentStateChanged :: Agent -> Signal (Maybe AgentState)-agentStateChanged v = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (agentStateUpdatedSource v)-        m2 = publishSignal (agentStateChangedSource v)---- | Return a signal that notifies about every change of the state.-agentStateChanged_ :: Agent -> Signal ()-agentStateChanged_ agent =-  mapSignal (const ()) $ agentStateChanged agent

@@ -1,46 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Base--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines basic functions for the 'Dynamics' monad.-----module Simulation.Aivika.Dynamics.Base-       (-- * Time Parameters-        starttime, -        stoptime,-        dt,-        time,-        integTimes,-        isTimeInteg,-        integIteration,-        integIterationBnds,-        integIterationLoBnd,-        integIterationHiBnd,-        -- * Interpolation and Initial Value-        initDynamics,-        discrete,-        interpolate,-        -- * Memoization-        memo,-        umemo,-        memo0,-        umemo0,-        -- * Iterating-        iterateDynamics,-        -- * Fold-        foldDynamics1,-        foldDynamics,-        -- * Norming-        divideDynamics) 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

@@ -1,189 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Buffer--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the limited queue similar to 'LIFO' and 'FIFO' but where --- the items are not represented. We know only of their number in the buffer and --- how many items were lost.----module Simulation.Aivika.Dynamics.Buffer-       (Buffer,-        bufferQueue,-        bufferNull,-        bufferFull,-        bufferMaxCount,-        bufferCount,-        bufferLostCount,-        bufferEnqueue,-        bufferDequeue,-        bufferEnqueueLost,-        newBuffer,-        dequeueBuffer,-        tryDequeueBuffer,-        enqueueBuffer,-        tryEnqueueBuffer,-        enqueueBufferOrLost) where--import Data.IORef--import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics-import Simulation.Aivika.Dynamics.Simulation-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Process-import Simulation.Aivika.Dynamics.Resource-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal--import Simulation.Aivika.Dynamics.LIFO-import Simulation.Aivika.Dynamics.FIFO---- | Represents the limited queue similar to 'LIFO' and 'FIFO' but where the items are not repsented.--- So, there is no order of items but their number is strictly limited.-data Buffer =-  Buffer { bufferQueue :: EventQueue,  -- ^ Return the event queue.-           bufferMaxCount :: Int,      -- ^ The maximum available number of items.-           bufferReadRes  :: Resource,-           bufferWriteRes :: Resource,-           bufferCountRef :: IORef Int,-           bufferLostCountRef :: IORef Int, -           bufferEnqueueSource :: SignalSource (),-           bufferEnqueueLostSource :: SignalSource (),-           bufferDequeueSource :: SignalSource (),-           bufferUpdatedSource :: SignalSource () }-  --- | Create a new queue with the specified maximum available number of items.  -newBuffer :: EventQueue -> Int -> Simulation Buffer  -newBuffer q count =-  do i <- liftIO $ newIORef 0-     l <- liftIO $ newIORef 0-     r <- newResourceWithCount q count 0-     w <- newResourceWithCount q count count-     s1 <- newSignalSourceUnsafe-     s2 <- newSignalSourceUnsafe-     s3 <- newSignalSourceUnsafe-     s4 <- newSignalSource q-     return Buffer { bufferQueue = q,-                     bufferMaxCount = count,-                     bufferReadRes  = r,-                     bufferWriteRes = w,-                     bufferCountRef = i,-                     bufferLostCountRef = l, -                     bufferEnqueueSource = s1,-                     bufferEnqueueLostSource = s2,-                     bufferDequeueSource = s3,-                     bufferUpdatedSource = s4 }-  --- | Test whether the queue is empty.-bufferNull :: Buffer -> Dynamics Bool-bufferNull q =-  do a <- bufferCount q-     return (a == 0)---- | Test whether the queue is full.-bufferFull :: Buffer -> Dynamics Bool-bufferFull q =-  do a <- bufferCount q-     return (a == bufferMaxCount q)---- | Return the queue size.-bufferCount :: Buffer -> Dynamics Int-bufferCount q =-  liftIO $ readIORef (bufferCountRef q)-  --- | Return the number of lost items.-bufferLostCount :: Buffer -> Dynamics Int-bufferLostCount q =-  liftIO $ readIORef (bufferLostCountRef q)-  --- | Dequeue suspending the process if the buffer is empty.-dequeueBuffer :: Buffer -> Process ()-dequeueBuffer q =-  do requestResource (bufferReadRes q)-     liftIO $ dequeueImpl q-     releaseResource (bufferWriteRes q)-     liftDynamics $ triggerSignal (bufferDequeueSource q) ()-  --- | Try to dequeue immediately.  -tryDequeueBuffer :: Buffer -> Dynamics Bool-tryDequeueBuffer q =-  do x <- tryRequestResourceInDynamics (bufferReadRes q)-     if x -       then do liftIO $ dequeueImpl q-               releaseResourceInDynamics (bufferWriteRes q)-               triggerSignal (bufferDequeueSource q) ()-               return True-       else return False---- | Enqueue the item suspending the process --- if the buffer is full.  -enqueueBuffer :: Buffer -> Process ()-enqueueBuffer q =-  do requestResource (bufferWriteRes q)-     liftIO $ enqueueImpl q-     releaseResource (bufferReadRes q)-     liftDynamics $ triggerSignal (bufferEnqueueSource q) ()-     --- | Try to enqueue the item immediately.  -tryEnqueueBuffer :: Buffer -> Dynamics Bool-tryEnqueueBuffer q =-  do x <- tryRequestResourceInDynamics (bufferWriteRes q)-     if x -       then do liftIO $ enqueueImpl q-               releaseResourceInDynamics (bufferReadRes q)-               triggerSignal (bufferEnqueueSource q) ()-               return True-       else return False---- | Try to enqueue the item. If the buffer is full--- then the item will be lost.-enqueueBufferOrLost :: Buffer -> Dynamics ()-enqueueBufferOrLost q =-  do x <- tryRequestResourceInDynamics (bufferWriteRes q)-     if x-       then do liftIO $ enqueueImpl q-               releaseResourceInDynamics (bufferReadRes q)-               triggerSignal (bufferEnqueueSource q) ()-       else do liftIO $ modifyIORef (bufferLostCountRef q) $ (+) 1-               triggerSignal (bufferEnqueueLostSource q) ()---- | Return a signal that notifies when any item is enqueued.-bufferEnqueue :: Buffer -> Signal ()-bufferEnqueue q = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (bufferUpdatedSource q)-        m2 = publishSignal (bufferEnqueueSource q)---- | Return a signal which notifies that the item was lost when--- attempting to add it to the full queue with help of --- 'enqueueBufferOrLost'.-bufferEnqueueLost :: Buffer -> Signal ()-bufferEnqueueLost q = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (bufferUpdatedSource q)-        m2 = publishSignal (bufferEnqueueLostSource q)---- | Return a signal that notifies when any item is dequeued.-bufferDequeue :: Buffer -> Signal ()-bufferDequeue q = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (bufferUpdatedSource q)-        m2 = publishSignal (bufferDequeueSource q)---- | An implementation method.-dequeueImpl :: Buffer -> IO ()-dequeueImpl q =-  do i <- readIORef (bufferCountRef q)-     let i' = i - 1-     i' `seq` writeIORef (bufferCountRef q) i'---- | An implementation method.-enqueueImpl :: Buffer -> IO ()-enqueueImpl q =-  do i <- readIORef (bufferCountRef q)-     let i' = i + 1-     i' `seq` writeIORef (bufferCountRef q) i'

@@ -1,18 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Cont--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ The 'Cont' monad is a variation of the standard Cont monad --- and F# async workflow, where the result of applying --- the continuation is a dynamic process.----module Simulation.Aivika.Dynamics.Cont-       (Cont) where--import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Internal.Cont

@@ -1,215 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.EventQueue--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ The module introduces the event queue. An event handler is--- the Dynamics computation that has a single purpose to perform--- some side effect at the desired time. To pass in any message--- to the event, you can use a closure.----module Simulation.Aivika.Dynamics.EventQueue-       (EventQueue,-        newQueue,-        enqueue,-        enqueueWithTimes,-        enqueueWithIntegTimes,-        enqueueWithStartTime,-        enqueueWithStopTime,-        enqueueWithCurrentTime,-        runQueue,-        runQueueSync,-        runQueueBefore,-        runQueueSyncBefore,-        queueCount) where--import Data.IORef-import Control.Monad--import Simulation.Aivika.Dynamics.Internal.Simulation-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 ()),-  queueBusy :: IORef Bool,-  queueTime :: IORef Double, -  -- Optimization-  runQueue  :: Dynamics (),-  -- ^ Run the event queue processing its events.-  -- There is no restiction on the time of the queue itself. It this time-  -- is greater than the current simulation time then nothing happens.-  runQueueSync :: Dynamics (),-  -- ^ Run the event queue synchronously, i.e. the current time cannot be-  -- less than the actual time of the queue itself.-  ---  -- You will rarely need to run the event queue explicitly, but-  -- if you do want then this function is probably that one you should use.-  runQueueBefore :: Dynamics (),-  -- ^ Run the event queue processing only those events-  -- which time is less than the current simulation time.-  -- There is no restiction on the time of the queue itself. It this time-  -- is greater than the current simulation time then nothing happens.-  runQueueSyncBefore :: Dynamics ()-  -- ^ Run the event queue synchronously processing only those events-  -- which time is less than the current simulation time. But the current-  -- time cannot be less than the actual time of the queue itself.-  ---  -- This function is usually called before a handler is subscribed-  -- to the signal. Earlier 'runQueueSync' was called instead, which could-  -- lead to the lost of the signal by the handler at time of direct-  -- subscribing. Changed in version 0.6.1.-  }---- | Create a new event queue.-newQueue :: Simulation EventQueue-newQueue = -  Simulation $ \r ->-  do let sc = runSpecs r-     f <- newIORef False-     t <- newIORef $ spcStartTime sc-     pq <- PQ.newQueue-     let q = EventQueue { queuePQ   = pq,-                          queueBusy = f,-                          queueTime = t, -                          runQueue  = runQueueCore True q,-                          runQueueSync = runQueueSyncCore True q,-                          runQueueBefore = runQueueCore False q,-                          runQueueSyncBefore = runQueueSyncCore False q }-     return q-             --- | Enqueue the event which must be actuated at the specified time.-enqueue :: EventQueue -> Double -> Dynamics () -> Dynamics ()-enqueue q t c = Dynamics r where-  r p = let pq = queuePQ q in PQ.enqueue pq t c-    --- | Run the event queue processing its events.-runQueueCore :: Bool -> EventQueue -> Dynamics ()-runQueueCore includingCurrentTime 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, c2) <- PQ.queueFront pq-            let t = queueTime q-            t' <- readIORef t-            when (t2 < t') $ -              error "The time value is too small: runQueueCore"-            when ((t2 < pointTime p) ||-                  (includingCurrentTime && (t2 == pointTime p))) $-              do writeIORef t t2-                 PQ.dequeue pq-                 let sc  = pointSpecs p-                     t0  = spcStartTime sc-                     dt  = spcDT sc-                     n2  = fromIntegral $ floor ((t2 - t0) / dt)-                     Dynamics k = c2-                 k $ p { pointTime = t2,-                         pointIteration = n2,-                         pointPhase = -1 }-                 call q p---- | Run the event queue synchronously, i.e. without past.-runQueueSyncCore :: Bool -> EventQueue -> Dynamics ()-runQueueSyncCore includingCurrentTime q = Dynamics r where-  r p =-    do let t = queueTime q-       t' <- readIORef t-       if pointTime p < t'-         then error $-              "The current time is less than " ++-              "the time in the queue: runQueueSyncCore"-         else m p-  Dynamics m = if includingCurrentTime-               then runQueue q-               else runQueueBefore q-  --- | Return the number of pending events that should--- be yet actuated.-queueCount :: EventQueue -> Dynamics Int-queueCount q = Dynamics r where-  r p = -    do let Dynamics m = runQueueSync q-       m p-       PQ.queueCount $ queuePQ q-       --- | Actuate the event handler in the specified time points.-enqueueWithTimes :: EventQueue -> [Double] -> Dynamics () -> Dynamics ()-enqueueWithTimes q ts m = loop ts-  where loop []       = return ()-        loop (t : ts) = enqueue q t $ m >> loop ts-       --- | Actuate the event handler in the specified time points.-enqueueWithPoints :: EventQueue -> [Point] -> Dynamics () -> Dynamics ()-enqueueWithPoints q xs (Dynamics m) = loop xs-  where loop []       = return ()-        loop (x : xs) = enqueue q (pointTime x) $ -                        Dynamics $ \p ->-                        do m x    -- N.B. we substitute the time point!-                           let Dynamics m' = loop xs-                           m' p---- | Actuate the event handler in the integration time points.-enqueueWithIntegTimes :: EventQueue -> Dynamics () -> Dynamics ()-enqueueWithIntegTimes q m =-  Dynamics $ \p ->-  do let sc  = pointSpecs p-         (nl, nu) = integIterationBnds sc-         points  = map point [nl .. nu]-         point n = Point { pointSpecs = sc,-                           pointRun = pointRun p,-                           pointTime = basicTime sc n 0,-                           pointIteration = n,-                           pointPhase = 0 }-         Dynamics m' = enqueueWithPoints q points m-     m' p---- | Actuate the event handler in the start time.-enqueueWithStartTime :: EventQueue -> Dynamics () -> Dynamics ()-enqueueWithStartTime q m =-  Dynamics $ \p ->-  do let sc  = pointSpecs p-         (nl, nu) = integIterationBnds sc-         point n = Point { pointSpecs = sc,-                           pointRun = pointRun p,-                           pointTime = basicTime sc n 0,-                           pointIteration = n,-                           pointPhase = 0 }-         Dynamics m' = enqueueWithPoints q [point nl] m-     m' p---- | Actuate the event handler in the stop time.-enqueueWithStopTime :: EventQueue -> Dynamics () -> Dynamics ()-enqueueWithStopTime q m =-  Dynamics $ \p ->-  do let sc  = pointSpecs p-         (nl, nu) = integIterationBnds sc-         point n = Point { pointSpecs = sc,-                           pointRun = pointRun p,-                           pointTime = basicTime sc n 0,-                           pointIteration = n,-                           pointPhase = 0 }-         Dynamics m' = enqueueWithPoints q [point nu] m-     m' p---- | Actuate the event handler in the current time but --- through the event queue, which allows continuing the --- current tasks and then calling the handler after the --- tasks are finished. The simulation time will be the same.-enqueueWithCurrentTime :: EventQueue -> Dynamics () -> Dynamics ()-enqueueWithCurrentTime q m =-  Dynamics $ \p ->-  do let Dynamics m' = enqueue q (pointTime p) m-     m' p

@@ -1,207 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.FIFO--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the FIFO queue.----module Simulation.Aivika.Dynamics.FIFO-       (FIFO,-        fifoQueue,-        fifoNull,-        fifoFull,-        fifoMaxCount,-        fifoCount,-        fifoLostCount,-        fifoEnqueue,-        fifoDequeue,-        fifoEnqueueLost,-        newFIFO,-        dequeueFIFO,-        tryDequeueFIFO,-        enqueueFIFO,-        tryEnqueueFIFO,-        enqueueFIFOOrLost) where--import Data.IORef-import Data.Array-import Data.Array.IO.Safe--import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics-import Simulation.Aivika.Dynamics.Simulation-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Process-import Simulation.Aivika.Dynamics.Resource-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal---- | Represents the FIFO queue with rule: first input - first output.-data FIFO a =-  FIFO { fifoQueue :: EventQueue,  -- ^ Return the event queue.-         fifoMaxCount :: Int,      -- ^ The maximum available number of items.-         fifoReadRes  :: Resource,-         fifoWriteRes :: Resource,-         fifoCountRef :: IORef Int,-         fifoLostCountRef :: IORef Int,-         fifoStartRef :: IORef Int,-         fifoEndRef   :: IORef Int,-         fifoArray :: IOArray Int a, -         fifoEnqueueSource :: SignalSource a,-         fifoEnqueueLostSource :: SignalSource a,-         fifoDequeueSource :: SignalSource a,-         fifoUpdatedSource :: SignalSource a }-  --- | Create a new FIFO queue with the specified maximum available number of items.  -newFIFO :: EventQueue -> Int -> Simulation (FIFO a)  -newFIFO q count =-  do i <- liftIO $ newIORef 0-     l <- liftIO $ newIORef 0-     s <- liftIO $ newIORef 0-     e <- liftIO $ newIORef 0-     a <- liftIO $ newArray_ (0, count - 1)-     r <- newResourceWithCount q count 0-     w <- newResourceWithCount q count count-     s1 <- newSignalSourceUnsafe-     s2 <- newSignalSourceUnsafe-     s3 <- newSignalSourceUnsafe-     s4 <- newSignalSource q-     return FIFO { fifoQueue = q,-                   fifoMaxCount = count,-                   fifoReadRes  = r,-                   fifoWriteRes = w,-                   fifoCountRef = i,-                   fifoLostCountRef = l,-                   fifoStartRef = s,-                   fifoEndRef   = e,-                   fifoArray = a, -                   fifoEnqueueSource = s1,-                   fifoEnqueueLostSource = s2,-                   fifoDequeueSource = s3,-                   fifoUpdatedSource = s4 }-  --- | Test whether the FIFO queue is empty.-fifoNull :: FIFO a -> Dynamics Bool-fifoNull fifo =-  do a <- fifoCount fifo-     return (a == 0)---- | Test whether the FIFO queue is full.-fifoFull :: FIFO a -> Dynamics Bool-fifoFull fifo =-  do a <- fifoCount fifo-     return (a == fifoMaxCount fifo)---- | Return the queue size.-fifoCount :: FIFO a -> Dynamics Int-fifoCount fifo =-  liftIO $ readIORef (fifoCountRef fifo)-  --- | Return the number of lost items.-fifoLostCount :: FIFO a -> Dynamics Int-fifoLostCount fifo =-  liftIO $ readIORef (fifoLostCountRef fifo)-  --- | Dequeue from the FIFO queue suspending the process if--- the queue is empty.-dequeueFIFO :: FIFO a -> Process a  -dequeueFIFO fifo =-  do requestResource (fifoReadRes fifo)-     a <- liftIO $ dequeueImpl fifo-     releaseResource (fifoWriteRes fifo)-     liftDynamics $ triggerSignal (fifoDequeueSource fifo) a-     return a-  --- | Try to dequeue from the FIFO queue immediately.  -tryDequeueFIFO :: FIFO a -> Dynamics (Maybe a)-tryDequeueFIFO fifo =-  do x <- tryRequestResourceInDynamics (fifoReadRes fifo)-     if x -       then do a <- liftIO $ dequeueImpl fifo-               releaseResourceInDynamics (fifoWriteRes fifo)-               triggerSignal (fifoDequeueSource fifo) a-               return $ Just a-       else return Nothing---- | Enqueue the item in the FIFO queue suspending the process--- if the queue is full.  -enqueueFIFO :: FIFO a -> a -> Process ()-enqueueFIFO fifo a =-  do requestResource (fifoWriteRes fifo)-     liftIO $ enqueueImpl fifo a-     releaseResource (fifoReadRes fifo)-     liftDynamics $ triggerSignal (fifoEnqueueSource fifo) a-     --- | Try to enqueue the item in the FIFO queue. Return 'False' in--- the monad if the queue is full.-tryEnqueueFIFO :: FIFO a -> a -> Dynamics Bool-tryEnqueueFIFO fifo a =-  do x <- tryRequestResourceInDynamics (fifoWriteRes fifo)-     if x -       then do liftIO $ enqueueImpl fifo a-               releaseResourceInDynamics (fifoReadRes fifo)-               triggerSignal (fifoEnqueueSource fifo) a-               return True-       else return False---- | Try to enqueue the item in the FIFO queue. If the queue is full--- then the item will be lost.-enqueueFIFOOrLost :: FIFO a -> a -> Dynamics ()-enqueueFIFOOrLost fifo a =-  do x <- tryRequestResourceInDynamics (fifoWriteRes fifo)-     if x-       then do liftIO $ enqueueImpl fifo a-               releaseResourceInDynamics (fifoReadRes fifo)-               triggerSignal (fifoEnqueueSource fifo) a-       else do liftIO $ modifyIORef (fifoLostCountRef fifo) $ (+) 1-               triggerSignal (fifoEnqueueLostSource fifo) a---- | Return a signal that notifies when any item is enqueued.-fifoEnqueue :: FIFO a -> Signal a-fifoEnqueue fifo = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (fifoUpdatedSource fifo)-        m2 = publishSignal (fifoEnqueueSource fifo)---- | Return a signal which notifies that the item was lost when --- attempting to add it to the full queue with help of--- 'enqueueFIFOOrLost'.-fifoEnqueueLost :: FIFO a -> Signal a-fifoEnqueueLost fifo = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (fifoUpdatedSource fifo)-        m2 = publishSignal (fifoEnqueueLostSource fifo)---- | Return a signal that notifies when any item is dequeued.-fifoDequeue :: FIFO a -> Signal a-fifoDequeue fifo = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (fifoUpdatedSource fifo)-        m2 = publishSignal (fifoDequeueSource fifo)---- | An implementation method.-dequeueImpl :: FIFO a -> IO a-dequeueImpl fifo =-  do i <- readIORef (fifoCountRef fifo)-     s <- readIORef (fifoStartRef fifo)-     let i' = i - 1-         s' = (s + 1) `mod` fifoMaxCount fifo-     a <- readArray (fifoArray fifo) s-     writeArray (fifoArray fifo) s undefined-     i' `seq` writeIORef (fifoCountRef fifo) i'-     s' `seq` writeIORef (fifoStartRef fifo) s'-     return a---- | An implementation method.-enqueueImpl :: FIFO a -> a -> IO ()-enqueueImpl fifo a =-  do i <- readIORef (fifoCountRef fifo)-     e <- readIORef (fifoEndRef fifo)-     let i' = i + 1-         e' = (e + 1) `mod` fifoMaxCount fifo-     a `seq` writeArray (fifoArray fifo) e a-     i' `seq` writeIORef (fifoCountRef fifo) i'-     e' `seq` writeIORef (fifoEndRef fifo) e'

@@ -0,0 +1,82 @@++-- |+-- Module     : Simulation.Aivika.Dynamics.Fold+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- This module defines the fold functions that allows traversing the values of+-- any 'Dynamics' computation in the integration time points.+--+module Simulation.Aivika.Dynamics.Fold+       (foldDynamics1,+        foldDynamics) where++import Data.IORef+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Dynamics+import Simulation.Aivika.Dynamics.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 +-- 'memo0Dynamics' or its unboxed version.+foldDynamics1 :: (Dynamics a -> Simulation (Dynamics a))+                 -> (a -> a -> a) +                 -> Dynamics a +                 -> Simulation (Dynamics a)+foldDynamics1 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+-- 'memo0Dynamics' or its unboxed version.+foldDynamics :: (Dynamics a -> Simulation (Dynamics a))+                -> (a -> b -> a) +                -> a+                -> Dynamics b +                -> Simulation (Dynamics a)+foldDynamics 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

@@ -1,304 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Cont--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ The 'Cont' monad is a variation of the standard Cont monad --- and F# async workflow, where the result of applying --- the continuation is a dynamic process.----module Simulation.Aivika.Dynamics.Internal.Cont-       (Cont(..),-        ContParams,-        runCont,-        catchCont,-        finallyCont,-        throwCont,-        resumeContByParams,-        contParamsCanceled) where--import Data.IORef--import qualified Control.Exception as C-import Control.Exception (IOException, throw)--import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics---- | The 'Cont' type is similar to the standard Cont monad --- and F# async workflow but only the continuations return--- a dynamic process as a result.-newtype Cont a = Cont (ContParams a -> Dynamics ())---- | The continuation parameters.-data ContParams a = -  ContParams { contCont :: a -> Dynamics (), -               contAux  :: ContParamsAux }---- | The auxiliary continuation parameters.-data ContParamsAux =-  ContParamsAux { contECont :: IOException -> Dynamics (),-                  contCCont :: () -> Dynamics (),-                  contCancelRef :: IORef Bool, -                  contCatchFlag :: Bool }--instance Monad Cont where-  return  = returnC-  m >>= k = bindC m k--instance SimulationLift Cont where-  liftSimulation = liftSC--instance DynamicsLift Cont where-  liftDynamics = liftDC--instance Functor Cont where-  fmap = liftM--instance MonadIO Cont where-  liftIO = liftIOC --invokeC :: Cont a -> ContParams a -> Dynamics ()-{-# INLINE invokeC #-}-invokeC (Cont m) = m--invokeD :: Point -> Dynamics a -> IO a-{-# INLINE invokeD #-}-invokeD p (Dynamics m) = m p--cancelD :: Point -> ContParams a -> IO ()-{-# NOINLINE cancelD #-}-cancelD p c =-  do writeIORef (contCancelRef . contAux $ c) False-     invokeD p $ (contCCont . contAux $ c) ()--returnC :: a -> Cont a-{-# INLINE returnC #-}-returnC a = -  Cont $ \c ->-  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ contCont c a-                          --- bindC :: Cont a -> (a -> Cont b) -> Cont b--- {-# INLINE bindC #-}--- bindC m k = ---   Cont $ \c -> ---   if (contCatchFlag . contAux $ c) ---   then bindWithCatch m k c---   else bindWithoutCatch m k c-  -bindC :: Cont a -> (a -> Cont b) -> Cont b-{-# INLINE bindC #-}-bindC m k = -  Cont $ bindWithoutCatch m k  -- Another version is not tail recursive!-  -bindWithoutCatch :: Cont a -> (a -> Cont b) -> ContParams b -> Dynamics ()-{-# INLINE bindWithoutCatch #-}-bindWithoutCatch (Cont m) k c = -  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ m $ -            let cont a = invokeC (k a) c-            in c { contCont = cont }---- It is not tail recursive!-bindWithCatch :: Cont a -> (a -> Cont b) -> ContParams b -> Dynamics ()-{-# NOINLINE bindWithCatch #-}-bindWithCatch (Cont m) k c = -  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ m $ -            let cont a = catchDynamics -                         (invokeC (k a) c)-                         (contECont . contAux $ c)-            in c { contCont = cont }---- Like "bindWithoutCatch (return a) k"-callWithoutCatch :: (a -> Cont b) -> a -> ContParams b -> Dynamics ()-callWithoutCatch k a c =-  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ invokeC (k a) c---- Like "bindWithCatch (return a) k" but it is not tail recursive!-callWithCatch :: (a -> Cont b) -> a -> ContParams b -> Dynamics ()-callWithCatch k a c =-  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ catchDynamics -            (invokeC (k a) c)-            (contECont . contAux $ c)---- | Exception handling within 'Cont' computations.-catchCont :: Cont a -> (IOException -> Cont a) -> Cont a-catchCont m h = -  Cont $ \c -> -  if contCatchFlag . contAux $ c-  then catchWithCatch m h c-  else error $-       "To catch exceptions, the process must be created " ++-       "with help of newProcessIDWithCatch: catchCont."-  -catchWithCatch :: Cont a -> (IOException -> Cont a) -> ContParams a -> Dynamics ()-catchWithCatch (Cont m) h c =-  Dynamics $ \p -> -  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ m $-            -- let econt e = callWithCatch h e c   -- not tail recursive!-            let econt e = callWithoutCatch h e c-            in c { contAux = (contAux c) { contECont = econt } }-               --- | A computation with finalization part.-finallyCont :: Cont a -> Cont b -> Cont a-finallyCont m m' = -  Cont $ \c -> -  if contCatchFlag . contAux $ c-  then finallyWithCatch m m' c-  else error $-       "To finalize computation, the process must be created " ++-       "with help of newProcessIDWithCatch: finallyCont."-  -finallyWithCatch :: Cont a -> Cont b -> ContParams a -> Dynamics ()               -finallyWithCatch (Cont m) (Cont m') c =-  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z -       then cancelD p c-       else invokeD p $ m $-            let cont a   = -                  Dynamics $ \p ->-                  invokeD p $ m' $-                  let cont b = contCont c a-                  in c { contCont = cont }-                econt e  =-                  Dynamics $ \p ->-                  invokeD p $ m' $-                  let cont b = (contECont . contAux $ c) e-                  in c { contCont = cont }-                ccont () = -                  Dynamics $ \p ->-                  invokeD p $ m' $-                  let cont b  = (contCCont . contAux $ c) ()-                      econt e = (contCCont . contAux $ c) ()-                  in c { contCont = cont,-                         contAux  = (contAux c) { contECont = econt } }-            in c { contCont = cont,-                   contAux  = (contAux c) { contECont = econt,-                                            contCCont = ccont } }---- | Throw the exception with the further exception handling.--- By some reasons, the standard 'throw' function per se is not handled --- properly within 'Cont' computations, altough it will be still handled --- if it will be hidden under the 'liftIO' function. The problem arises --- namely with the @throw@ function, not 'IO' computations.-throwCont :: IOException -> Cont a-throwCont e = liftIO $ throw e---- | Run the 'Cont' computation with the specified cancelation token --- and flag indicating whether to catch exceptions.-runCont :: Cont a -> -           (a -> Dynamics ()) ->-           (IOError -> Dynamics ()) ->-           (() -> Dynamics ()) ->-           IORef Bool -> -           Bool -> -           Dynamics ()-runCont (Cont m) cont econt ccont cancelToken catchFlag = -  m ContParams { contCont = cont,-                 contAux  = -                   ContParamsAux { contECont = econt,-                                   contCCont = ccont,-                                   contCancelRef = cancelToken, -                                   contCatchFlag = catchFlag } }---- | Lift the 'Simulation' computation.-liftSC :: Simulation a -> Cont a-liftSC (Simulation m) = -  Cont $ \c ->-  Dynamics $ \p ->-  if contCatchFlag . contAux $ c-  then liftIOWithCatch (m $ pointRun p) p c-  else liftIOWithoutCatch (m $ pointRun p) p c-     --- | Lift the 'Dynamics' computation.-liftDC :: Dynamics a -> Cont a-liftDC (Dynamics m) =-  Cont $ \c ->-  Dynamics $ \p ->-  if contCatchFlag . contAux $ c-  then liftIOWithCatch (m p) p c-  else liftIOWithoutCatch (m p) p c-     --- | Lift the IO computation.-liftIOC :: IO a -> Cont a-liftIOC m =-  Cont $ \c ->-  Dynamics $ \p ->-  if contCatchFlag . contAux $ c-  then liftIOWithCatch m p c-  else liftIOWithoutCatch m p c-  -liftIOWithoutCatch :: IO a -> Point -> ContParams a -> IO ()-{-# INLINE liftIOWithoutCatch #-}-liftIOWithoutCatch m p c =-  do z <- readIORef $ (contCancelRef . contAux) c-     if z-       then cancelD p c-       else do a <- m-               invokeD p $ contCont c a--liftIOWithCatch :: IO a -> Point -> ContParams a -> IO ()-{-# NOINLINE liftIOWithCatch #-}-liftIOWithCatch m p c =-  do z <- readIORef $ (contCancelRef . contAux) c-     if z-       then cancelD p c-       else do aref <- newIORef undefined-               eref <- newIORef Nothing-               C.catch (m >>= writeIORef aref) -                 (writeIORef eref . Just)-               e <- readIORef eref-               case e of-                 Nothing -> -                   do a <- readIORef aref-                      -- tail recursive-                      invokeD p $ contCont c a-                 Just e ->-                   -- tail recursive-                   invokeD p $ (contECont . contAux) c e---- | Resume the computation by the specified parameters.-resumeContByParams :: ContParams a -> a -> Dynamics ()-{-# INLINE resumeContByParams #-}-resumeContByParams c a = -  Dynamics $ \p ->-  do z <- readIORef $ (contCancelRef . contAux) c-     if z-       then cancelD p c-       else invokeD p $ contCont c a---- | Test whether the computation is canceled-contParamsCanceled :: ContParams a -> IO Bool-{-# INLINE contParamsCanceled #-}-contParamsCanceled c = -  readIORef $ (contCancelRef . contAux) c

@@ -1,310 +0,0 @@--{-# LANGUAGE RecursiveDo #-}---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Dynamics--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ The module defines the 'Dynamics' monad representing an abstract dynamic --- process, i.e. a time varying polymorphic function. ----module Simulation.Aivika.Dynamics.Internal.Dynamics-       (-- * Dynamics-        Dynamics(..),-        DynamicsLift(..),-        Point(..),-        runDynamicsInStartTime,-        runDynamicsInStopTime,-        runDynamicsInIntegTimes,-        runDynamicsInTime,-        runDynamicsInTimes,-        -- * Error Handling-        catchDynamics,-        finallyDynamics,-        throwDynamics,-        -- * Utilities-        basicTime,-        integIterationBnds,-        integIterationHiBnd,-        integIterationLoBnd,-        integPhaseBnds,-        integPhaseHiBnd,-        integPhaseLoBnd) where--import qualified Control.Exception as C-import Control.Exception (IOException, throw, finally)--import Control.Monad-import Control.Monad.Trans-import Control.Monad.Fix--import Simulation.Aivika.Dynamics.Internal.Simulation------- 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)-           --- | Returns the integration iterations starting from zero.-integIterations :: Specs -> [Int]-integIterations sc = [i1 .. i2] where-  i1 = 0-  i2 = round ((spcStopTime sc - -               spcStartTime sc) / spcDT sc)---- | Returns the first and last integration iterations.-integIterationBnds :: Specs -> (Int, Int)-integIterationBnds sc = (0, round ((spcStopTime sc - -                                    spcStartTime sc) / spcDT sc))---- | Returns the first integration iteration, i.e. zero.-integIterationLoBnd :: Specs -> Int-integIterationLoBnd sc = 0---- | Returns the last integration iteration.-integIterationHiBnd :: Specs -> Int-integIterationHiBnd sc = round ((spcStopTime sc - -                                 spcStartTime sc) / spcDT sc)---- | Returns the phases for the specified simulation specs starting from zero.-integPhases :: Specs -> [Int]-integPhases sc = -  case spcMethod sc of-    Euler -> [0]-    RungeKutta2 -> [0, 1]-    RungeKutta4 -> [0, 1, 2, 3]---- | Returns the first and last integration phases.-integPhaseBnds :: Specs -> (Int, Int)-integPhaseBnds sc = -  case spcMethod sc of-    Euler -> (0, 0)-    RungeKutta2 -> (0, 1)-    RungeKutta4 -> (0, 3)---- | Returns the first integration phase, i.e. zero.-integPhaseLoBnd :: Specs -> Int-integPhaseLoBnd sc = 0-                  --- | Returns the last integration phase, 0 for Euler's method, 1 for RK2 and 3 for RK4.-integPhaseHiBnd :: Specs -> Int-integPhaseHiBnd 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' = fromIntegral 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---- | Run the dynamic process in the initial simulation point.-runDynamicsInStartTime :: Dynamics a -> Simulation a-runDynamicsInStartTime (Dynamics m) =-  Simulation $ \r ->-  do let sc = runSpecs r -         n  = 0-         t  = spcStartTime sc-     m Point { pointSpecs = sc,-               pointRun = r,-               pointTime = t,-               pointIteration = n,-               pointPhase = 0 }---- | Run the dynamic process in the final simulation point.-runDynamicsInStopTime :: Dynamics a -> Simulation a-runDynamicsInStopTime (Dynamics m) =-  Simulation $ \r ->-  do let sc = runSpecs r -         n  = integIterationHiBnd sc-         t  = basicTime sc n 0-     m Point { pointSpecs = sc,-               pointRun = r,-               pointTime = t,-               pointIteration = n,-               pointPhase = 0 }---- | Run the dynamic process in all integration time points-runDynamicsInIntegTimes :: Dynamics a -> Simulation [IO a]-runDynamicsInIntegTimes (Dynamics m) =-  Simulation $ \r ->-  do let sc = runSpecs r-         (nl, nu) = integIterationBnds sc-         point n = Point { pointSpecs = sc,-                           pointRun = r,-                           pointTime = basicTime sc n 0,-                           pointIteration = n,-                           pointPhase = 0 }-     return $ map (m . point) [nl .. nu]---- | Run the dynamic process in the specified time point.-runDynamicsInTime :: Double -> Dynamics a -> Simulation a-runDynamicsInTime t (Dynamics m) =-  Simulation $ \r ->-  do let sc = runSpecs r-         t0 = spcStartTime sc-         dt = spcDT sc-         n  = fromIntegral $ floor ((t - t0) / dt)-     m Point { pointSpecs = sc,-               pointRun = r,-               pointTime = t,-               pointIteration = n,-               pointPhase = -1 }---- | Run the dynamic process in the specified time points.-runDynamicsInTimes :: [Double] -> Dynamics a -> Simulation [IO a]-runDynamicsInTimes ts (Dynamics m) =-  Simulation $ \r ->-  do let sc = runSpecs r-         t0 = spcStartTime sc-         dt = spcDT sc-         point t =-           let n = fromIntegral $ floor ((t - t0) / dt)-           in Point { pointSpecs = sc,-                      pointRun = r,-                      pointTime = t,-                      pointIteration = n,-                      pointPhase = -1 }-     return $ map (m . point) ts--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--instance SimulationLift Dynamics where-  liftSimulation = liftDS-    -liftDS :: Simulation a -> Dynamics a-{-# INLINE liftDS #-}-liftDS (Simulation m) =-  Dynamics $ \p -> m $ pointRun p---- | A type class to lift the 'Dynamics' computations in other monads.-class Monad m => DynamicsLift m where-  -  -- | Lift the specified 'Dynamics' computation in another monad.-  liftDynamics :: Dynamics a -> m a-  --- | Exception handling within 'Dynamics' computations.-catchDynamics :: Dynamics a -> (IOException -> Dynamics a) -> Dynamics a-catchDynamics (Dynamics m) h =-  Dynamics $ \p -> -  C.catch (m p) $ \e ->-  let Dynamics m' = h e in m' p-                           --- | A computation with finalization part like the 'finally' function.-finallyDynamics :: Dynamics a -> Dynamics b -> Dynamics a-finallyDynamics (Dynamics m) (Dynamics m') =-  Dynamics $ \p ->-  C.finally (m p) (m' p)---- | Like the standard 'throw' function.-throwDynamics :: IOException -> Dynamics a-throwDynamics = throw---- | Invoke the 'Dynamics' computation.-invokeDynamics :: Dynamics a -> Point -> IO a-{-# INLINE invokeDynamics #-}-invokeDynamics (Dynamics m) p = m p--instance MonadFix Dynamics where-  mfix f = -    Dynamics $ \p ->-    do { rec { a <- invokeDynamics (f a) p }; return a }

@@ -1,92 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Fold--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.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-       (foldDynamics1,-        foldDynamics,-        divideDynamics) where--import Data.IORef-import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Simulation-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'.-foldDynamics1 :: (Dynamics a -> Simulation (Dynamics a))-                 -> (a -> a -> a) -                 -> Dynamics a -                 -> Simulation (Dynamics a)-foldDynamics1 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'.-foldDynamics :: (Dynamics a -> Simulation (Dynamics a))-                -> (a -> b -> a) -                -> a-                -> Dynamics b -                -> Simulation (Dynamics a)-foldDynamics 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.-divideDynamics :: Dynamics Double -> Dynamics Double-divideDynamics (Dynamics m) = -  discrete $ Dynamics $ \p ->-  do a <- m p-     return $ a / fromIntegral (pointIteration p + 1)

@@ -1,62 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Interpolate--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.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-       (initDynamics,-        discrete,-        interpolate) where--import Simulation.Aivika.Dynamics.Internal.Dynamics---- | Return the initial value.-initDynamics :: Dynamics a -> Dynamics a-{-# INLINE initDynamics #-}-initDynamics (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 }

@@ -1,196 +0,0 @@--{-# LANGUAGE FlexibleContexts #-}---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Memo--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.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,-        iterateDynamics) where--import Data.Array-import Data.Array.IO.Safe-import Data.IORef-import Control.Monad--import Simulation.Aivika.Dynamics.Internal.Simulation-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 -> Simulation (Dynamics e)-{-# INLINE memo #-}-memo (Dynamics m) = -  Simulation $ \r ->-  do let sc = runSpecs r-         (phl, phu) = integPhaseBnds sc-         (nl, nu)   = integIterationBnds 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 = integPhaseHiBnd 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 -> Simulation (Dynamics e)-{-# INLINE umemo #-}-umemo (Dynamics m) = -  Simulation $ \r ->-  do let sc = runSpecs r-         (phl, phu) = integPhaseBnds sc-         (nl, nu)   = integIterationBnds 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 = integPhaseHiBnd 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 -> Simulation (Dynamics e)-{-# INLINE memo0 #-}-memo0 (Dynamics m) = -  Simulation $ \r ->-  do let sc   = runSpecs r-         bnds = integIterationBnds 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 -> Simulation (Dynamics e)-{-# INLINE umemo0 #-}-umemo0 (Dynamics m) = -  Simulation $ \r ->-  do let sc   = runSpecs r-         bnds = integIterationBnds 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.-iterateDynamics :: Dynamics () -> Simulation (Dynamics ())-{-# INLINE iterateDynamics #-}-iterateDynamics (Dynamics m) = -  Simulation $ \r ->-  do let sc = runSpecs r-     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

@@ -1,300 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Process--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ A value in the 'Process' monad represents a discontinuous process that --- can suspend in any simulation time point and then resume later in the same --- or another time point. --- --- The process of this type behaves like a dynamic process too. So, 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,-        newProcessIDWithCatch,-        holdProcess,-        interruptProcess,-        processInterrupted,-        passivateProcess,-        processPassive,-        reactivateProcess,-        processID,-        cancelProcess,-        processCanceled,-        runProcess,-        runProcessNow,-        catchProcess,-        finallyProcess,-        throwProcess) where--import Data.Maybe-import Data.IORef-import Control.Exception (IOException, throw)-import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Internal.Cont-import Simulation.Aivika.Dynamics.EventQueue---- | Represents a process identificator.-data ProcessID = -  ProcessID { processQueue   :: EventQueue,  -- ^ Return the event queue.-              processStarted :: IORef Bool,-              processCatchFlag     :: Bool,-              processReactCont     :: IORef (Maybe (ContParams ())), -              processCancelRef     :: IORef Bool, -              processCancelToken   :: IORef Bool,-              processInterruptRef  :: IORef Bool, -              processInterruptCont :: IORef (Maybe (ContParams ())), -              processInterruptVersion :: IORef Int }---- | 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 $ \p ->-  do let x = processInterruptCont pid-     writeIORef x $ Just c-     writeIORef (processInterruptRef pid) False-     v <- readIORef (processInterruptVersion pid)-     let Dynamics m = -           enqueue (processQueue pid) (pointTime p + dt) $-           Dynamics $ \p ->-           do v' <- readIORef (processInterruptVersion pid)-              when (v == v') $ -                do writeIORef x Nothing-                   let Dynamics m = resumeContByParams c ()-                   m p-     m p---- | Interrupt a process with the specified ID if the process--- was held by computation 'holdProcess'.-interruptProcess :: ProcessID -> Dynamics ()-interruptProcess pid =-  Dynamics $ \p ->-  do let x = processInterruptCont pid-     a <- readIORef x-     case a of-       Nothing -> return ()-       Just c ->-         do writeIORef x Nothing-            writeIORef (processInterruptRef pid) True-            modifyIORef (processInterruptVersion pid) $ (+) 1-            let Dynamics m = -                  enqueue (processQueue pid) (pointTime p) $ -                  resumeContByParams c ()-            m p-            --- | Test whether the process with the specified ID was interrupted.-processInterrupted :: ProcessID -> Dynamics Bool-processInterrupted pid =-  Dynamics $ \p ->-  readIORef (processInterruptRef pid)---- | Passivate the process.-passivateProcess :: Process ()-passivateProcess =-  Process $ \pid ->-  Cont $ \c ->-  Dynamics $ \p ->-  do let x = processReactCont 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 -> Dynamics Bool-processPassive pid =-  Dynamics $ \p ->-  do let Dynamics m = runQueueSync $ processQueue pid-     m p-     let x = processReactCont pid-     a <- readIORef x-     return $ isJust a---- | Reactivate a process with the specified ID.-reactivateProcess :: ProcessID -> Dynamics ()-reactivateProcess pid =-  Dynamics $ \p ->-  do let Dynamics m = runQueueSync $ processQueue pid-     m p-     let x = processReactCont pid-     a <- readIORef x-     case a of-       Nothing -> -         return ()-       Just c ->-         do writeIORef x Nothing-            let Dynamics m  = enqueue (processQueue pid) (pointTime p) $ -                              resumeContByParams c ()-            m p---- | Start the process with the specified ID at the desired time.-runProcess :: Process () -> ProcessID -> Double -> Dynamics ()-runProcess (Process p) pid t =-  runCont m cont econt ccont (processCancelToken pid) (processCatchFlag pid)-    where cont  = return-          econt = throw-          ccont = return-          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 -> enqueue (processQueue pid) t $ -                              resumeContByParams c ()-                 p pid---- | Start the process with the specified ID at the current simulation time.-runProcessNow :: Process () -> ProcessID -> Dynamics ()-runProcessNow process pid =-  Dynamics $ \p ->-  do let Dynamics m = runProcess process pid (pointTime p)-     m p---- | Return the current process ID.-processID :: Process ProcessID-processID = Process $ \pid -> return pid---- | Create a new process ID without exception handling.-newProcessID :: EventQueue -> Simulation ProcessID-newProcessID q =-  do x <- liftIO $ newIORef Nothing-     y <- liftIO $ newIORef False-     c <- liftIO $ newIORef False-     t <- liftIO $ newIORef False-     i <- liftIO $ newIORef False-     z <- liftIO $ newIORef Nothing-     v <- liftIO $ newIORef 0-     return ProcessID { processQueue   = q,-                        processStarted = y,-                        processCatchFlag     = False,-                        processReactCont     = x, -                        processCancelRef     = c, -                        processCancelToken   = t,-                        processInterruptRef  = i,-                        processInterruptCont = z, -                        processInterruptVersion = v }---- | Create a new process ID with capabilities of catching --- the IOError exceptions and finalizing the computation. --- The corresponded process will be slower than that one--- which identifier is created with help of 'newProcessID'.-newProcessIDWithCatch :: EventQueue -> Simulation ProcessID-newProcessIDWithCatch q =-  do x <- liftIO $ newIORef Nothing-     y <- liftIO $ newIORef False-     c <- liftIO $ newIORef False-     t <- liftIO $ newIORef False-     i <- liftIO $ newIORef False-     z <- liftIO $ newIORef Nothing-     v <- liftIO $ newIORef 0-     return ProcessID { processQueue   = q,-                        processStarted = y,-                        processCatchFlag     = True,-                        processReactCont     = x, -                        processCancelRef     = c, -                        processCancelToken   = t,-                        processInterruptRef  = i,-                        processInterruptCont = z, -                        processInterruptVersion = v }---- | Cancel a process with the specified ID.-cancelProcess :: ProcessID -> Dynamics ()-cancelProcess pid =-  Dynamics $ \p ->-  do z <- readIORef (processCancelRef pid) -     unless z $-       do writeIORef (processCancelRef pid) True-          writeIORef (processCancelToken pid) True---- | Test whether the process with the specified ID is canceled.-processCanceled :: ProcessID -> Dynamics Bool-processCanceled pid =-  Dynamics $ \p ->-  readIORef (processCancelRef pid)--instance Eq ProcessID where-  x == y = processReactCont x == processReactCont y    -- for the references are unique--instance Monad Process where-  return  = returnP-  m >>= k = bindP m k--instance Functor Process where-  fmap = liftM--instance SimulationLift Process where-  liftSimulation = liftSP-  -instance DynamicsLift Process where-  liftDynamics = liftDP-  -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--liftSP :: Simulation a -> Process a-{-# INLINE liftSP #-}-liftSP m = Process $ \pid -> liftSimulation m--liftDP :: Dynamics a -> Process a-{-# INLINE liftDP #-}-liftDP m = Process $ \pid -> liftDynamics m--liftIOP :: IO a -> Process a-{-# INLINE liftIOP #-}-liftIOP m = Process $ \pid -> liftIO m---- | Exception handling within 'Process' computations.-catchProcess :: Process a -> (IOException -> Process a) -> Process a-catchProcess (Process m) h =-  Process $ \pid ->-  catchCont (m pid) $ \e ->-  let Process m' = h e in m' pid-                           --- | A computation with finalization part.-finallyProcess :: Process a -> Process b -> Process a-finallyProcess (Process m) (Process m') =-  Process $ \pid ->-  finallyCont (m pid) (m' pid)---- | Throw the exception with the further exception handling.--- By some reasons, the standard 'throw' function per se is not handled --- properly within 'Process' computations, although it will be still --- handled if it will be hidden under the 'liftIO' function. The problem --- arises namely with the @throw@ function, not 'IO' computations.-throwProcess :: IOException -> Process a-throwProcess = liftIO . throw-

@@ -1,347 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Signal--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the signal which we can subscribe handlers to. --- These handlers can be disposed. The signal is triggered in the --- current time point actuating the corresponded computations from --- the handlers. -----module Simulation.Aivika.Dynamics.Internal.Signal-       (Signal,-        SignalSource,-        newSignalSourceWithUpdate,-        newSignalSourceUnsafe,-        publishSignal,-        triggerSignal,-        handleSignal,-        handleSignal_,-        updateSignal,-        mapSignal,-        mapSignalM,-        apSignal,-        filterSignal,-        filterSignalM,-        emptySignal,-        merge2Signals,-        merge3Signals,-        merge4Signals,-        merge5Signals) where--import Data.IORef-import Data.Monoid--import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Internal.Simulation---- | The signal source that can publish its signal.-data SignalSource a =-  SignalSource { publishSignal :: Signal a,-                                  -- ^ Publish the signal.-                 triggerSignal :: a -> Dynamics ()-                                  -- ^ Trigger the signal actuating -                                  -- all its handlers at the current -                                  -- simulation time point.-               }-  --- | The signal that can have disposable handlers.  -data Signal a =-  Signal { handleSignal :: (a -> Dynamics ()) -> -                           Dynamics (Dynamics ()),-           -- ^ Subscribe the handler to the specified -           -- signal and return a nested computation -           -- that, being applied, unsubscribes the -           -- handler from this signal.-           ---           -- If the signal is bound up with the event queue-           -- then the signal in the current time is not lost-           -- by the handler any more. Changed in version 0.6.1.-           updateSignal :: Dynamics ()-           -- ^ Update the signal to its actual state.-           ---           -- You will rarely need to call this function directly-           -- as it is usually called implicitly.-           ---           -- Since version 0.6.1 it processes only those events-           -- which time is less than the current simulation time-           -- if the signal is bound up with the event queue,-           -- although you can define your own 'updateSignal'-           -- function when creating a new signal source with help of-           -- 'newSignalSourceWithUpdate'.-         }-  --- | The queue of signal handlers.-data SignalHandlerQueue a =-  SignalHandlerQueue { queueStart :: IORef (Maybe (SignalHandler a)),-                       queueEnd   :: IORef (Maybe (SignalHandler a)) }-  --- | It contains the information about the disposable queue handler.-data SignalHandler a =-  SignalHandler { handlerComp :: a -> Dynamics (),-                  handlerPrev :: IORef (Maybe (SignalHandler a)),-                  handlerNext :: IORef (Maybe (SignalHandler a)) }---- | Subscribe the handler to the specified signal.--- To subscribe the disposable handlers, use function 'handleSignal'.------ If the signal is bound up with the event queue then the signal in--- the current time is not lost by the handler any more.--- Changed in version 0.6.1.-handleSignal_ :: Signal a -> (a -> Dynamics ()) -> Dynamics ()-handleSignal_ signal h = -  do x <- handleSignal signal h-     return ()-     --- | Create a new signal source with the specified update computation.-newSignalSourceWithUpdate :: Dynamics () -> Simulation (SignalSource a)-newSignalSourceWithUpdate update =-  Simulation $ \r ->-  do start <- newIORef Nothing-     end <- newIORef Nothing-     let queue  = SignalHandlerQueue { queueStart = start,-                                       queueEnd   = end }-         signal = Signal { handleSignal = handle, -                           updateSignal = update }-         source = SignalSource { publishSignal = signal, -                                 triggerSignal = trigger }-         handle h =-           Dynamics $ \p ->-           do invokeDynamics p update-              x <- enqueueSignalHandler queue h-              return $ -                Dynamics $ \p ->-                do invokeDynamics p update-                   dequeueSignalHandler queue x-         trigger a =-           Dynamics $ \p ->-           do invokeDynamics p update -              let h = queueStart queue-              triggerSignalHandlers h a p-     return source-     --- | Create a new signal source that has no update computation.-newSignalSourceUnsafe :: Simulation (SignalSource a)-newSignalSourceUnsafe =-  Simulation $ \r ->-  do start <- newIORef Nothing-     end <- newIORef Nothing-     let queue  = SignalHandlerQueue { queueStart = start,-                                       queueEnd   = end }-         signal = Signal { handleSignal = handle, -                           updateSignal = update }-         source = SignalSource { publishSignal = signal, -                                 triggerSignal = trigger }-         handle h =-           Dynamics $ \p ->-           do x <- enqueueSignalHandler queue h-              return $ liftIO $ dequeueSignalHandler queue x-         trigger a =-           Dynamics $ \p ->-           let h = queueStart queue-           in triggerSignalHandlers h a p-         update = return ()-     return source---- | Trigger all next signal handlers.-triggerSignalHandlers :: IORef (Maybe (SignalHandler a)) -> a -> Point -> IO ()-{-# INLINE triggerSignalHandlers #-}-triggerSignalHandlers r a p =-  do x <- readIORef r-     case x of-       Nothing -> return ()-       Just h ->-         do let Dynamics m = handlerComp h a-            m p-            triggerSignalHandlers (handlerNext h) a p-            --- | Enqueue the handler and return its representative in the queue.            -enqueueSignalHandler :: SignalHandlerQueue a -> (a -> Dynamics ()) -> IO (SignalHandler a)-enqueueSignalHandler q h = -  do tail <- readIORef (queueEnd q)-     case tail of-       Nothing ->-         do prev <- newIORef Nothing-            next <- newIORef Nothing-            let handler = SignalHandler { handlerComp = h,-                                          handlerPrev = prev,-                                          handlerNext = next }-            writeIORef (queueStart q) (Just handler)-            writeIORef (queueEnd q) (Just handler)-            return handler-       Just x ->-         do prev <- newIORef tail-            next <- newIORef Nothing-            let handler = SignalHandler { handlerComp = h,-                                          handlerPrev = prev,-                                          handlerNext = next }-            writeIORef (handlerNext x) (Just handler)-            writeIORef (queueEnd q) (Just handler)-            return handler---- | Dequeue the handler representative.-dequeueSignalHandler :: SignalHandlerQueue a -> SignalHandler a -> IO ()-dequeueSignalHandler q h = -  do prev <- readIORef (handlerPrev h)-     case prev of-       Nothing ->-         do next <- readIORef (handlerNext h)-            case next of-              Nothing ->-                do writeIORef (queueStart q) Nothing-                   writeIORef (queueEnd q) Nothing-              Just y ->-                do writeIORef (handlerPrev y) Nothing-                   writeIORef (handlerNext h) Nothing-                   writeIORef (queueStart q) next-       Just x ->-         do next <- readIORef (handlerNext h)-            case next of-              Nothing ->-                do writeIORef (handlerPrev h) Nothing-                   writeIORef (handlerNext x) Nothing-                   writeIORef (queueEnd q) prev-              Just y ->-                do writeIORef (handlerPrev h) Nothing-                   writeIORef (handlerNext h) Nothing-                   writeIORef (handlerPrev y) prev-                   writeIORef (handlerNext x) next--instance Functor Signal where-  fmap = mapSignal-  -instance Monoid (Signal a) where -  -  mempty = emptySignal-  -  mappend = merge2Signals-  -  mconcat [] = emptySignal-  mconcat [x1] = x1-  mconcat [x1, x2] = merge2Signals x1 x2-  mconcat [x1, x2, x3] = merge3Signals x1 x2 x3-  mconcat [x1, x2, x3, x4] = merge4Signals x1 x2 x3 x4-  mconcat [x1, x2, x3, x4, x5] = merge5Signals x1 x2 x3 x4 x5-  mconcat (x1 : x2 : x3 : x4 : x5 : xs) = -    mconcat $ merge5Signals x1 x2 x3 x4 x5 : xs-  --- | Map the signal according the specified function.-mapSignal :: (a -> b) -> Signal a -> Signal b-mapSignal f m =-  Signal { handleSignal = \h -> -            handleSignal m $ h . f, -           updateSignal = -             updateSignal m }---- | Filter only those signal values that satisfy to --- the specified predicate.-filterSignal :: (a -> Bool) -> Signal a -> Signal a-filterSignal p m =-  Signal { handleSignal = \h ->-            handleSignal m $ \a ->-            when (p a) $ h a, -           updateSignal =-             updateSignal m }-  --- | Filter only those signal values that satisfy to --- the specified predicate.-filterSignalM :: (a -> Dynamics Bool) -> Signal a -> Signal a-filterSignalM p m =-  Signal { handleSignal = \h ->-            handleSignal m $ \a ->-            do x <- p a-               when x $ h a, -           updateSignal =-             updateSignal m }-  --- | Merge two signals.-merge2Signals :: Signal a -> Signal a -> Signal a-merge2Signals m1 m2 =-  Signal { handleSignal = \h ->-            do x1 <- handleSignal m1 h-               x2 <- handleSignal m2 h-               return $ do { x1; x2 }, -           updateSignal =-             do updateSignal m1-                updateSignal m2 }---- | Merge three signals.-merge3Signals :: Signal a -> Signal a -> Signal a -> Signal a-merge3Signals m1 m2 m3 =-  Signal { handleSignal = \h ->-            do x1 <- handleSignal m1 h-               x2 <- handleSignal m2 h-               x3 <- handleSignal m3 h-               return $ do { x1; x2; x3 },-           updateSignal =-             do updateSignal m1-                updateSignal m2 -                updateSignal m3 }---- | Merge four signals.-merge4Signals :: Signal a -> Signal a -> Signal a -> -                 Signal a -> Signal a-merge4Signals m1 m2 m3 m4 =-  Signal { handleSignal = \h ->-            do x1 <- handleSignal m1 h-               x2 <- handleSignal m2 h-               x3 <- handleSignal m3 h-               x4 <- handleSignal m4 h-               return $ do { x1; x2; x3; x4 },-           updateSignal =-             do updateSignal m1-                updateSignal m2 -                updateSignal m3 -                updateSignal m4 }-           --- | Merge five signals.-merge5Signals :: Signal a -> Signal a -> Signal a -> -                 Signal a -> Signal a -> Signal a-merge5Signals m1 m2 m3 m4 m5 =-  Signal { handleSignal = \h ->-            do x1 <- handleSignal m1 h-               x2 <- handleSignal m2 h-               x3 <- handleSignal m3 h-               x4 <- handleSignal m4 h-               x5 <- handleSignal m5 h-               return $ do { x1; x2; x3; x4; x5 },-           updateSignal =-             do updateSignal m1-                updateSignal m2 -                updateSignal m3 -                updateSignal m4-                updateSignal m5 }---- | Compose the signal.-mapSignalM :: (a -> Dynamics b) -> Signal a -> Signal b-mapSignalM f m =-  Signal { handleSignal = \h ->-            handleSignal m (f >=> h),-           updateSignal = -             updateSignal m }-  --- | Transform the signal.-apSignal :: Dynamics (a -> b) -> Signal a -> Signal b-apSignal f m =-  Signal { handleSignal = \h ->-            handleSignal m $ \a -> do { x <- f; h (x a) },-           updateSignal =-             updateSignal m }---- | An empty signal which is never triggered.-emptySignal :: Signal a-emptySignal =-  Signal { handleSignal = \h -> return $ return (),-           updateSignal = return () }-  -invokeDynamics :: Point -> Dynamics a -> IO a-{-# INLINE invokeDynamics #-}-invokeDynamics p (Dynamics m) = m p

@@ -1,202 +0,0 @@--{-# LANGUAGE RecursiveDo #-}---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Simulation--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ The module defines the 'Simulation' monad that represents a simulation run.--- -module Simulation.Aivika.Dynamics.Internal.Simulation-       (-- * Simulation-        Simulation(..),-        SimulationLift(..),-        Specs(..),-        Method(..),-        Run(..),-        runSimulation,-        runSimulations,-        -- * Error Handling-        catchSimulation,-        finallySimulation,-        throwSimulation,-        -- * Utilities-        simulationIndex,-        simulationCount,-        simulationSpecs) where--import qualified Control.Exception as C-import Control.Exception (IOException, throw, finally)--import Control.Monad-import Control.Monad.Trans-import Control.Monad.Fix------- The Simulation Monad------- | A value in the 'Simulation' monad represents something that--- doesn't change within the simulation run but may change for--- other runs.------ This monad is ideal for representing the external--- parameters for the model, when the Monte-Carlo simulation--- is used. Also this monad is useful for defining some--- actions that should occur only once within the simulation run,--- for example, setting of the integral with help of recursive--- equations.----newtype Simulation a = Simulation (Run -> IO a)---- | 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 indentifies the simulation run.-data Run = Run { runSpecs :: Specs,  -- ^ the simulation specs-                 runIndex :: Int,    -- ^ the current simulation run index-                 runCount :: Int     -- ^ the total number of runs in this experiment-               } deriving (Eq, Ord, Show)--instance Monad Simulation where-  return  = returnS-  m >>= k = bindS m k--returnS :: a -> Simulation a-returnS a = Simulation (\r -> return a)--bindS :: Simulation a -> (a -> Simulation b) -> Simulation b-bindS (Simulation m) k = -  Simulation $ \r -> -  do a <- m r-     let Simulation m' = k a-     m' r---- | Run the simulation using the specified specs.-runSimulation :: Simulation a -> Specs -> IO a-runSimulation (Simulation m) sc =-  m Run { runSpecs = sc,-          runIndex = 1,-          runCount = 1 }---- | Run the given number of simulations using the specified specs, ---   where each simulation is distinguished by its index 'simulationIndex'.-runSimulations :: Simulation a -> Specs -> Int -> [IO a]-runSimulations (Simulation m) sc runs = map f [1 .. runs]-  where f i = m Run { runSpecs = sc,-                      runIndex = i,-                      runCount = runs }---- | Return the run index for the current simulation.-simulationIndex :: Simulation Int-simulationIndex = Simulation $ return . runIndex---- | Return the number of simulations currently run.-simulationCount :: Simulation Int-simulationCount = Simulation $ return . runCount---- | Return the simulation specs-simulationSpecs :: Simulation Specs-simulationSpecs = Simulation $ return . runSpecs--instance Functor Simulation where-  fmap = liftMS--instance Eq (Simulation a) where-  x == y = error "Can't compare simulation runs." --instance Show (Simulation a) where-  showsPrec _ x = showString "<< Simulation >>"--liftMS :: (a -> b) -> Simulation a -> Simulation b-{-# INLINE liftMS #-}-liftMS f (Simulation x) =-  Simulation $ \r -> do { a <- x r; return $ f a }--liftM2S :: (a -> b -> c) -> Simulation a -> Simulation b -> Simulation c-{-# INLINE liftM2S #-}-liftM2S f (Simulation x) (Simulation y) =-  Simulation $ \r -> do { a <- x r; b <- y r; return $ f a b }--instance (Num a) => Num (Simulation a) where-  x + y = liftM2S (+) x y-  x - y = liftM2S (-) x y-  x * y = liftM2S (*) x y-  negate = liftMS negate-  abs = liftMS abs-  signum = liftMS signum-  fromInteger i = return $ fromInteger i--instance (Fractional a) => Fractional (Simulation a) where-  x / y = liftM2S (/) x y-  recip = liftMS recip-  fromRational t = return $ fromRational t--instance (Floating a) => Floating (Simulation a) where-  pi = return pi-  exp = liftMS exp-  log = liftMS log-  sqrt = liftMS sqrt-  x ** y = liftM2S (**) x y-  sin = liftMS sin-  cos = liftMS cos-  tan = liftMS tan-  asin = liftMS asin-  acos = liftMS acos-  atan = liftMS atan-  sinh = liftMS sinh-  cosh = liftMS cosh-  tanh = liftMS tanh-  asinh = liftMS asinh-  acosh = liftMS acosh-  atanh = liftMS atanh--instance MonadIO Simulation where-  liftIO m = Simulation $ const m---- | A type class to lift the simulation computations in other monads.-class Monad m => SimulationLift m where-  -  -- | Lift the specified 'Simulation' computation in another monad.-  liftSimulation :: Simulation a -> m a-    --- | Exception handling within 'Simulation' computations.-catchSimulation :: Simulation a -> (IOException -> Simulation a) -> Simulation a-catchSimulation (Simulation m) h =-  Simulation $ \r -> -  C.catch (m r) $ \e ->-  let Simulation m' = h e in m' r-                           --- | A computation with finalization part like the 'finally' function.-finallySimulation :: Simulation a -> Simulation b -> Simulation a-finallySimulation (Simulation m) (Simulation m') =-  Simulation $ \r ->-  C.finally (m r) (m' r)---- | Like the standard 'throw' function.-throwSimulation :: IOException -> Simulation a-throwSimulation = throw---- | Invoke the 'Simulation' computation.-invokeSimulation :: Simulation a -> Run -> IO a-{-# INLINE invokeSimulation #-}-invokeSimulation (Simulation m) r = m r--instance MonadFix Simulation where-  mfix f = -    Simulation $ \r ->-    do { rec { a <- invokeSimulation (f a) r }; return a }  

@@ -1,59 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Internal.Time--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the time parameters.-----module Simulation.Aivika.Dynamics.Internal.Time-       (starttime, -        stoptime, -        dt, -        time,-        integTimes, -        isTimeInteg,-        integIteration,-        integPhase) where--import Simulation.Aivika.Dynamics.Internal.Simulation-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 ---- | Return the integration time points.-integTimes :: Specs -> [Double]-integTimes sc = map t [nl .. nu]-  where (nl, nu) = integIterationBnds sc-        t n = basicTime sc n 0-     --- | Whether the current time is an integration time.-isTimeInteg :: Dynamics Bool-isTimeInteg = Dynamics $ \p -> return $ pointPhase p >= 0---- | Return the integration iteration closest to the current simulation time.-integIteration :: Dynamics Int-integIteration = Dynamics $ return . pointIteration---- | Return the integration phase for the current simulation time.--- It is @(-1)@ for non-integration time points.-integPhase :: Dynamics Int-integPhase = Dynamics $ return . pointPhase

@@ -0,0 +1,63 @@++-- |+-- Module     : Simulation.Aivika.Dynamics.Interpolate+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- This module defines interpolation functions.+-- These functions complement the memoization, possibly except for +-- the 'initDynamics' function which is useful to get an initial +-- value of any dynamic process.+--++module Simulation.Aivika.Dynamics.Interpolate+       (initDynamics,+        discreteDynamics,+        interpolateDynamics) where++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Dynamics++-- | Return the initial value.+initDynamics :: Dynamics a -> Dynamics a+{-# INLINE initDynamics #-}+initDynamics (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.+discreteDynamics :: Dynamics a -> Dynamics a+{-# INLINE discreteDynamics #-}+discreteDynamics (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 'discreteDynamics' function it knows about the intermediate +-- time points that are used in the Runge-Kutta method.+interpolateDynamics :: Dynamics a -> Dynamics a+{-# INLINE interpolateDynamics #-}+interpolateDynamics (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 }

@@ -1,196 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.LIFO--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the LIFO queue.----module Simulation.Aivika.Dynamics.LIFO-       (LIFO,-        lifoQueue,-        lifoNull,-        lifoFull,-        lifoMaxCount,-        lifoCount,-        lifoLostCount,-        lifoEnqueue,-        lifoDequeue,-        lifoEnqueueLost,-        newLIFO,-        dequeueLIFO,-        tryDequeueLIFO,-        enqueueLIFO,-        tryEnqueueLIFO,-        enqueueLIFOOrLost) where--import Data.IORef-import Data.Array-import Data.Array.IO.Safe--import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics-import Simulation.Aivika.Dynamics.Simulation-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Process-import Simulation.Aivika.Dynamics.Resource-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal---- | Represents the LIFO queue with rule: last input - first output.-data LIFO a =-  LIFO { lifoQueue :: EventQueue,  -- ^ Return the event queue.-         lifoMaxCount :: Int,      -- ^ The maximum available number of items.-         lifoReadRes  :: Resource,-         lifoWriteRes :: Resource,-         lifoCountRef :: IORef Int,-         lifoLostCountRef :: IORef Int,-         lifoArray :: IOArray Int a, -         lifoEnqueueSource :: SignalSource a,-         lifoEnqueueLostSource :: SignalSource a,-         lifoDequeueSource :: SignalSource a,-         lifoUpdatedSource :: SignalSource a }-  --- | Create a new LIFO queue with the specified maximum available number of items.  -newLIFO :: EventQueue -> Int -> Simulation (LIFO a)  -newLIFO q count =-  do i <- liftIO $ newIORef 0-     l <- liftIO $ newIORef 0-     a <- liftIO $ newArray_ (0, count - 1)-     r <- newResourceWithCount q count 0-     w <- newResourceWithCount q count count-     s1 <- newSignalSourceUnsafe-     s2 <- newSignalSourceUnsafe-     s3 <- newSignalSourceUnsafe-     s4 <- newSignalSource q-     return LIFO { lifoQueue = q,-                   lifoMaxCount = count,-                   lifoReadRes  = r,-                   lifoWriteRes = w,-                   lifoCountRef = i,-                   lifoLostCountRef = l,-                   lifoArray = a,-                   lifoEnqueueSource = s1,-                   lifoEnqueueLostSource = s2,-                   lifoDequeueSource = s3,-                   lifoUpdatedSource = s4 }-  --- | Test whether the LIFO queue is empty.-lifoNull :: LIFO a -> Dynamics Bool-lifoNull lifo =-  do a <- lifoCount lifo-     return (a == 0)---- | Test whether the LIFO queue is full.-lifoFull :: LIFO a -> Dynamics Bool-lifoFull lifo =-  do a <- lifoCount lifo-     return (a == lifoMaxCount lifo)---- | Return the queue size.-lifoCount :: LIFO a -> Dynamics Int-lifoCount lifo =-  liftIO $ readIORef (lifoCountRef lifo)-  --- | Return the number of lost items.-lifoLostCount :: LIFO a -> Dynamics Int-lifoLostCount lifo =-  liftIO $ readIORef (lifoLostCountRef lifo)-  --- | Dequeue from the LIFO queue suspending the process if--- the queue is empty.-dequeueLIFO :: LIFO a -> Process a  -dequeueLIFO lifo =-  do requestResource (lifoReadRes lifo)-     a <- liftIO $ dequeueImpl lifo-     releaseResource (lifoWriteRes lifo)-     liftDynamics $ triggerSignal (lifoDequeueSource lifo) a-     return a-  --- | Try to dequeue from the LIFO queue immediately.  -tryDequeueLIFO :: LIFO a -> Dynamics (Maybe a)-tryDequeueLIFO lifo =-  do x <- tryRequestResourceInDynamics (lifoReadRes lifo)-     if x -       then do a <- liftIO $ dequeueImpl lifo-               releaseResourceInDynamics (lifoWriteRes lifo)-               triggerSignal (lifoDequeueSource lifo) a-               return $ Just a-       else return Nothing---- | Enqueue the item in the LIFO queue suspending the process if--- the queue is full.  -enqueueLIFO :: LIFO a -> a -> Process ()-enqueueLIFO lifo a =-  do requestResource (lifoWriteRes lifo)-     liftIO $ enqueueImpl lifo a-     releaseResource (lifoReadRes lifo)-     liftDynamics $ triggerSignal (lifoEnqueueSource lifo) a-     --- | Try to enqueue the item in the LIFO queue. Return 'False' in--- the monad if the queue is full.-tryEnqueueLIFO :: LIFO a -> a -> Dynamics Bool-tryEnqueueLIFO lifo a =-  do x <- tryRequestResourceInDynamics (lifoWriteRes lifo)-     if x -       then do liftIO $ enqueueImpl lifo a-               releaseResourceInDynamics (lifoReadRes lifo)-               triggerSignal (lifoEnqueueSource lifo) a-               return True-       else return False---- | Try to enqueue the item in the LIFO queue. If the queue is full--- then the item will be lost.-enqueueLIFOOrLost :: LIFO a -> a -> Dynamics ()-enqueueLIFOOrLost lifo a =-  do x <- tryRequestResourceInDynamics (lifoWriteRes lifo)-     if x-       then do liftIO $ enqueueImpl lifo a-               releaseResourceInDynamics (lifoReadRes lifo)-               triggerSignal (lifoEnqueueSource lifo) a-       else do liftIO $ modifyIORef (lifoLostCountRef lifo) $ (+) 1-               triggerSignal (lifoEnqueueLostSource lifo) a---- | Return a signal that notifies when any item is enqueued.-lifoEnqueue :: LIFO a -> Signal a-lifoEnqueue lifo = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (lifoUpdatedSource lifo)-        m2 = publishSignal (lifoEnqueueSource lifo)---- | Return a signal which notifies that the item was lost when--- attempting to add it to the full queue with help of--- 'enqueueLIFOOrLost'.-lifoEnqueueLost :: LIFO a -> Signal a-lifoEnqueueLost lifo = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (lifoUpdatedSource lifo)-        m2 = publishSignal (lifoEnqueueLostSource lifo)---- | Return a signal that notifies when any item is dequeued.-lifoDequeue :: LIFO a -> Signal a-lifoDequeue lifo = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (lifoUpdatedSource lifo)-        m2 = publishSignal (lifoDequeueSource lifo)----- | An implementation method.-dequeueImpl :: LIFO a -> IO a-dequeueImpl lifo =-  do i <- readIORef (lifoCountRef lifo)-     let j = i - 1-     a <- j `seq` readArray (lifoArray lifo) j-     writeArray (lifoArray lifo) j undefined-     writeIORef (lifoCountRef lifo) j-     return a---- | An implementation method.-enqueueImpl :: LIFO a -> a -> IO ()-enqueueImpl lifo a =-  do i <- readIORef (lifoCountRef lifo)-     let j = i + 1-     a `seq` writeArray (lifoArray lifo) i a-     j `seq` writeIORef (lifoCountRef lifo) j

@@ -0,0 +1,125 @@++-- |+-- Module     : Simulation.Aivika.Dynamics.Memo+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- This module defines memo functions. The memoization creates such 'Dynamics'+-- computations, which values are cached in the integration time points. Then+-- these values are interpolated in all other time points.+--++module Simulation.Aivika.Dynamics.Memo+       (memoDynamics,+        memo0Dynamics,+        iterateDynamics) where++import Data.Array+import Data.Array.IO.Safe+import Data.IORef+import Control.Monad++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Dynamics+import Simulation.Aivika.Dynamics.Interpolate++-- | Create a boxed array with default values.+newBoxedArray_ :: Ix i => (i, i) -> IO (IOArray i e)+newBoxedArray_ = newArray_++-- | Memoize and order the computation in the integration time points using +-- the interpolation that knows of the Runge-Kutta method.+memoDynamics :: Dynamics e -> Simulation (Dynamics e)+{-# INLINE memoDynamics #-}+memoDynamics (Dynamics m) = +  Simulation $ \r ->+  do let sc = runSpecs r+         (phl, phu) = integPhaseBnds sc+         (nl, nu)   = integIterationBnds sc+     arr   <- newBoxedArray_ ((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 = integPhaseHiBnd 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 $ interpolateDynamics $ Dynamics r++-- | Memoize and order the computation in the integration time points using +-- the 'discreteDynamics' interpolation. It consumes less memory than the 'memoDynamics'+-- 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 'memo0Dynamics' function above 'memoDynamics'.+memo0Dynamics :: Dynamics e -> Simulation (Dynamics e)+{-# INLINE memo0Dynamics #-}+memo0Dynamics (Dynamics m) = +  Simulation $ \r ->+  do let sc   = runSpecs r+         bnds = integIterationBnds sc+     arr  <- newBoxedArray_ 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 $ discreteDynamics $ Dynamics r++-- | Iterate sequentially the dynamic process with side effects in +-- the integration time points. It is equivalent to a call of the+-- 'memo0Dynamics' function but significantly more efficient, for the array +-- is not created.+iterateDynamics :: Dynamics () -> Simulation (Dynamics ())+{-# INLINE iterateDynamics #-}+iterateDynamics (Dynamics m) = +  Simulation $ \r ->+  do let sc = runSpecs r+     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 $ discreteDynamics $ Dynamics r

@@ -0,0 +1,100 @@++{-# LANGUAGE FlexibleContexts #-}++-- |+-- Module     : Simulation.Aivika.Dynamics.Memo.Unboxed+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- This module defines the unboxed memo functions. The memoization creates such 'Dynamics'+-- computations, which values are cached in the integration time points. Then+-- these values are interpolated in all other time points.+--++module Simulation.Aivika.Dynamics.Memo.Unboxed+       (memoDynamics,+        memo0Dynamics) where++import Data.Array+import Data.Array.IO.Safe+import Data.IORef+import Control.Monad++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Dynamics+import Simulation.Aivika.Dynamics.Interpolate+import Simulation.Aivika.Unboxed++-- | Memoize and order the computation in the integration time points using +-- the interpolation that knows of the Runge-Kutta method.+memoDynamics :: Unboxed e => Dynamics e -> Simulation (Dynamics e)+{-# INLINE memoDynamics #-}+memoDynamics (Dynamics m) = +  Simulation $ \r ->+  do let sc = runSpecs r+         (phl, phu) = integPhaseBnds sc+         (nl, nu)   = integIterationBnds sc+     arr   <- newUnboxedArray_ ((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 = integPhaseHiBnd 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 $ interpolateDynamics $ Dynamics r++-- | Memoize and order the computation in the integration time points using +-- the 'discreteDynamics' interpolation. It consumes less memory than the 'memoDynamics'+-- 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 'memo0Dynamics' function above 'memoDynamics'.+memo0Dynamics :: Unboxed e => Dynamics e -> Simulation (Dynamics e)+{-# INLINE memo0Dynamics #-}+memo0Dynamics (Dynamics m) = +  Simulation $ \r ->+  do let sc   = runSpecs r+         bnds = integIterationBnds sc+     arr  <- newUnboxedArray_ 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 $ discreteDynamics $ Dynamics r

@@ -1,83 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Parameter--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the parameters of simulation experiments.-----module Simulation.Aivika.Dynamics.Parameter-       (newParameter,-        newTableParameter,-        newIndexedParameter,-        newRandomParameter,-        newNormalParameter) where--import Data.Array-import Data.IORef-import qualified Data.Map as M-import Control.Concurrent.MVar-import System.Random--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Random---- | Create a thread-safe parameter that returns always the same value within the simulation run, --- where the value is recalculated for each new run.-newParameter :: IO a -> IO (Simulation a)-newParameter a = newIndexedParameter $ \_ -> a---- | Create a thread-safe parameter that returns always the same value within 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 (Simulation 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 within the simulation run, --- where the value depends on the number of this run starting from zero.-newIndexedParameter :: (Int -> IO a) -> IO (Simulation a)-newIndexedParameter f = -  do lock <- newMVar ()-     dict <- newIORef M.empty-     return $ Simulation $ \r ->-       do let i = runIndex r-          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 }---- | Create a new random parameter distributed uniformly.--- The value doesn't change within the simulation run but--- then the value is recalculated for each new run.-newRandomParameter :: Simulation Double     -- ^ minimum-                      -> Simulation Double  -- ^ maximum-                      -> IO (Simulation Double)-newRandomParameter min max =-  do x <- newParameter $ getStdRandom random-     return $ min + x * (max - min)---- | Create a new random parameter distributed normally.--- The value doesn't change within the simulation run but--- then the value is recalculated for each new run.-newNormalParameter :: Simulation Double     -- ^ mean-                      -> Simulation Double  -- ^ variance-                      -> IO (Simulation Double)-newNormalParameter mu nu =-  do x <- normalGen >>= newParameter-     return $ mu + x * nu

@@ -1,42 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Process--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ A value in the 'Process' monad represents a discontinuous process that --- can suspend in any simulation time point and then resume later in the same --- or another time point. --- --- The process of this type behaves like a dynamic process too. So, 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,-        newProcessIDWithCatch,-        holdProcess,-        interruptProcess,-        processInterrupted,-        passivateProcess,-        processPassive,-        reactivateProcess,-        processID,-        cancelProcess,-        processCanceled,-        runProcess,-        runProcessNow,-        catchProcess,-        finallyProcess,-        throwProcess) where--import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Internal.Process

@@ -7,74 +7,38 @@ -- Stability  : experimental -- Tested with: GHC 7.6.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.+-- Below are defined random functions that return the 'Dynamics' computations. +-- The values are initially defined in the integration time points and then+-- they are passed in to the 'memo0Dynamics' function to memoize and then interpolate. --  module Simulation.Aivika.Dynamics.Random -       (newRandom, newNormal, normalGen) where+       (newRandomDynamics, newNormalDynamics) where  import System.Random import Data.IORef import Control.Monad.Trans +import Simulation.Aivika.Simulation+import Simulation.Aivika.Random import Simulation.Aivika.Dynamics-import Simulation.Aivika.Dynamics.Simulation-import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.Memo.Unboxed  -- | Return the uniform random numbers in the integration time points.-newRandom :: Dynamics Double     -- ^ minimum-             -> Dynamics Double  -- ^ maximum-             -> Simulation (Dynamics Double)-newRandom min max =-  umemo0 $ do+newRandomDynamics :: Dynamics Double     -- ^ minimum+                  -> Dynamics Double  -- ^ maximum+                  -> Simulation (Dynamics Double)+newRandomDynamics min max =+  memo0Dynamics $ do     x <- liftIO $ getStdRandom random     min + return x * (max - min)       -- | Return the normal random numbers in the integration time points.-newNormal :: Dynamics Double     -- ^ mean-             -> Dynamics Double  -- ^ variance-             -> Simulation (Dynamics Double)-newNormal mu nu =-  do g <- liftIO normalGen-     umemo0 $ do+newNormalDynamics :: Dynamics Double     -- ^ mean+                  -> Dynamics Double  -- ^ variance+                  -> Simulation (Dynamics Double)+newNormalDynamics mu nu =+  do g <- liftIO newNormalGen+     memo0Dynamics $ do        x <- liftIO g        mu + return x * nu---- | Normal random number generator with mean 0 and variance 1.-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

@@ -1,89 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Ref--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines an updatable reference that depends on the event queue.----module Simulation.Aivika.Dynamics.Ref-       (Ref,-        refQueue,-        refChanged,-        refChanged_,-        newRef,-        readRef,-        writeRef,-        modifyRef) where--import Data.IORef-import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal---- | 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, -        refChangedSource :: SignalSource a, -        refUpdatedSource :: SignalSource a }---- | Create a new reference bound to the specified event queue.-newRef :: EventQueue -> a -> Simulation (Ref a)-newRef q a =-  do x <- liftIO $ newIORef a-     s <- newSignalSourceUnsafe-     u <- newSignalSource q-     return Ref { refQueue = q,-                  refRun   = runQueueSync q,-                  refValue = x, -                  refChangedSource = s, -                  refUpdatedSource = u }-     --- | 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 invokeDynamics p $ refRun r-     readIORef (refValue r)---- | Write a new value into the reference.-writeRef :: Ref a -> a -> Dynamics ()-writeRef r a = Dynamics $ \p -> -  do a `seq` writeIORef (refValue r) a-     invokeDynamics p $ triggerSignal (refChangedSource 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 invokeDynamics p $ refRun r-     a <- readIORef (refValue r)-     let b = f a-     b `seq` writeIORef (refValue r) b-     invokeDynamics p $ triggerSignal (refChangedSource r) b---- | Return a signal that notifies about every change of the reference state.-refChanged :: Ref a -> Signal a-refChanged v = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (refUpdatedSource v)-        m2 = publishSignal (refChangedSource v)---- | Return a signal that notifies about every change of the reference state.-refChanged_ :: Ref a -> Signal ()-refChanged_ r = mapSignal (const ()) $ refChanged r--invokeDynamics :: Point -> Dynamics a -> IO a-{-# INLINE invokeDynamics #-}-invokeDynamics p (Dynamics m) = m p

@@ -1,171 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Resource--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines a limited resource which can be acquired and --- then released by the discontinuous process 'Process'.----module Simulation.Aivika.Dynamics.Resource-       (Resource,-        newResource,-        newResourceWithCount,-        resourceQueue,-        resourceInitCount,-        resourceCount,-        requestResource,-        tryRequestResourceInDynamics,-        releaseResource,-        releaseResourceInDynamics,-        usingResource) where--import Data.IORef-import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Simulation-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 (ContParams ())}--instance Eq Resource where-  x == y = resourceCountRef x == resourceCountRef y  -- unique references---- | Create a new resource with the specified initial count.-newResource :: EventQueue -> Int -> Simulation Resource-newResource q initCount =-  Simulation $ \r ->-  do countRef  <- newIORef initCount-     waitQueue <- Q.newQueue-     return Resource { resourceQueue     = q,-                       resourceInitCount = initCount,-                       resourceCountRef  = countRef,-                       resourceWaitQueue = waitQueue }---- | Create a new resource with the specified initial count.--- The third argument specifies how the resource is consumed --- at the beginning, i.e. it defines the current count, which must be --- non-negative and less or equal to the initial count.-newResourceWithCount :: EventQueue -> Int -> Int -> Simulation Resource-newResourceWithCount q initCount count = do-  when (count < 0) $-    error $-    "The resource count cannot be negative: " ++-    "newResourceWithCount."-  when (count > initCount) $-    error $-    "The resource count cannot be greater than " ++-    "its initial value: newResourceWithCount."-  Simulation $ \r ->-    do countRef  <- newIORef count-       waitQueue <- Q.newQueue-       return Resource { resourceQueue     = q,-                         resourceInitCount = initCount,-                         resourceCountRef  = countRef,-                         resourceWaitQueue = waitQueue }---- | Return the current count of the resource.-resourceCount :: Resource -> Dynamics Int-resourceCount r =-  Dynamics $ \p ->-  do invokeDynamics p $ runQueueSync (resourceQueue r)-     readIORef (resourceCountRef r)---- | 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 $ \pid ->-  Cont $ \c ->-  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'-               invokeDynamics p $ resumeContByParams c ()---- | Release the resource increasing its count and resuming one of the--- previously suspended processes as possible.-releaseResource :: Resource -> Process ()-releaseResource r = -  Process $ \_ ->-  Cont $ \c ->-  Dynamics $ \p ->-  do invokeDynamics p $ releaseResourceUnsafe r-     invokeDynamics p $ resumeContByParams c ()---- | Release the resource increasing its count and resuming one of the--- previously suspended processes as possible.-releaseResourceInDynamics :: Resource -> Dynamics ()-releaseResourceInDynamics r =-  Dynamics $ \p ->-  do invokeDynamics p $ runQueueSync (resourceQueue r)-     invokeDynamics p $ releaseResourceUnsafe r--releaseResourceUnsafe :: Resource -> Dynamics ()-{-# INLINE releaseResourceUnsafe #-}-releaseResourceUnsafe r =-  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: releaseResourceUnsafe."-     f <- Q.queueNull (resourceWaitQueue r)-     if f -       then a' `seq` writeIORef (resourceCountRef r) a'-       else do c <- Q.queueFront (resourceWaitQueue r)-               Q.dequeue (resourceWaitQueue r)-               invokeDynamics p $ enqueue (resourceQueue r) (pointTime p) $-                 Dynamics $ \p ->-                 do z <- contParamsCanceled c-                    if z-                      then do invokeDynamics p $ releaseResourceUnsafe r-                              invokeDynamics p $ resumeContByParams c ()-                      else invokeDynamics p $ resumeContByParams c ()---- | Try to request for the resource decreasing its count in case of success--- and returning 'True' in the 'Dynamics' monad; otherwise, returning 'False'.-tryRequestResourceInDynamics :: Resource -> Dynamics Bool-tryRequestResourceInDynamics r =-  Dynamics $ \p ->-  do invokeDynamics p $ runQueueSync (resourceQueue r)-     a <- readIORef (resourceCountRef r)-     if a == 0 -       then return False-       else do let a' = a - 1-               a' `seq` writeIORef (resourceCountRef r) a'-               return True-               --- | Acquire the resource, perform some action and safely release the resource               --- in the end, even if the 'IOException' was raised within the action. --- The process identifier must be created with support of exception --- handling, i.e. with help of function 'newProcessIDWithCatch'. Unfortunately,--- such processes are slower than those that are created with help of--- other function 'newProcessID'.-usingResource :: Resource -> Process a -> Process a-usingResource r m =-  do requestResource r-     finallyProcess m $ releaseResource r--invokeDynamics :: Point -> Dynamics a -> IO a-{-# INLINE invokeDynamics #-}-invokeDynamics p (Dynamics m) = m p 

@@ -1,184 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Signal--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines the signal which we can subscribe handlers to. --- These handlers can be disposed. The signal is triggered in the --- current time point actuating the corresponded computations from --- the handlers. ----module Simulation.Aivika.Dynamics.Signal-       (Signal,-        SignalSource,-        newSignalSource,-        newSignalSourceWithUpdate,-        newSignalInTimes,-        newSignalInIntegTimes,-        newSignalInStartTime,-        newSignalInStopTime,-        publishSignal,-        triggerSignal,-        handleSignal,-        handleSignal_,-        updateSignal,-        awaitSignal,-        mapSignal,-        mapSignalM,-        apSignal,-        filterSignal,-        filterSignalM,-        emptySignal,-        merge2Signals,-        merge3Signals,-        merge4Signals,-        merge5Signals,-        SignalHistory,-        signalHistorySignal,-        newSignalHistory,-        newSignalHistoryThrough,-        readSignalHistory) where--import Data.IORef-import Data.Array--import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Internal.Cont-import Simulation.Aivika.Dynamics.Internal.Process-import Simulation.Aivika.Dynamics.Base--import qualified Simulation.Aivika.Vector as V-import qualified Simulation.Aivika.UVector as UV---- | Create a new signal source when the state depends on the event queue.------ Since version 0.6.1 its 'updateSignal' function calls 'runQueueSyncBefore'--- instead of 'runQueueSync' as it was before. In case of need you can--- define your own update function with help of 'newSignalSourceWithUpdate'.------ The function has the following defintion:------ @--- newSignalSource queue = ---   newSignalSourceWithUpdate $ runQueueSyncBefore queue--- @-newSignalSource :: EventQueue -> Simulation (SignalSource a)-newSignalSource queue = -  newSignalSourceWithUpdate $ runQueueSyncBefore queue---- | Await the signal.-awaitSignal :: Signal a -> Process a-awaitSignal signal =-  Process $ \pid ->-  Cont $ \c ->-  Dynamics $ \p ->-  do r <- newIORef Nothing-     let Dynamics m = -           handleSignal signal $ -           \a -> Dynamics $ -                 \p -> do x <- readIORef r-                          case x of-                            Nothing ->-                              error "The signal was lost: awaitSignal."-                            Just x ->-                              do let Dynamics m = x-                                 m p-                                 let Dynamics m = resumeContByParams c a-                                 m p-     h <- m p-     writeIORef r $ Just h-          --- | Represents the history of the signal values.-data SignalHistory a =-  SignalHistory { signalHistorySignal :: Signal a,  -                  -- ^ The signal for which the history is created.-                  signalHistoryTimes  :: UV.UVector Double,-                  signalHistoryValues :: V.Vector a }---- | Create a history of the signal values.-newSignalHistory :: Signal a -> Dynamics (SignalHistory a)-newSignalHistory signal =-  do ts <- liftIO UV.newVector-     xs <- liftIO V.newVector-     handleSignal_ signal $ \a ->-       Dynamics $ \p ->-       do liftIO $ UV.appendVector ts (pointTime p)-          liftIO $ V.appendVector xs a-     return SignalHistory { signalHistorySignal = signal,-                            signalHistoryTimes  = ts,-                            signalHistoryValues = xs }-       --- | Create a history of the signal values with delay through the event queue.--- The history will be created at the same simulation time, just the corresponded --- handler will be subscribed to the signal after the new event will be processed --- by the queue. ------ Since version 0.6.1, this function has less meaning than before. Please use--- carefully as the behavior depends on the state of the event queue.-newSignalHistoryThrough :: EventQueue -> Signal a -> Dynamics (SignalHistory a)-newSignalHistoryThrough q signal =-  do ts <- liftIO UV.newVector-     xs <- liftIO V.newVector-     enqueueWithCurrentTime q $-       handleSignal_ signal $ \a ->-       Dynamics $ \p ->-       do liftIO $ UV.appendVector ts (pointTime p)-          liftIO $ V.appendVector xs a-     return SignalHistory { signalHistorySignal = signal,-                            signalHistoryTimes  = ts,-                            signalHistoryValues = xs }-       --- | Read the history of signal values.-readSignalHistory :: SignalHistory a -> Dynamics (Array Int Double, Array Int a)-readSignalHistory history =-  do updateSignal $ signalHistorySignal history-     xs <- liftIO $ UV.freezeVector (signalHistoryTimes history)-     ys <- liftIO $ V.freezeVector (signalHistoryValues history)-     return (xs, ys)     -     --- | Trigger the signal with the current time.-triggerSignalWithTime :: SignalSource Double -> Dynamics ()-triggerSignalWithTime s =-  Dynamics $ \p ->-  do let Dynamics m = triggerSignal s (pointTime p)-     m p---- | Return a signal that is triggered in the specified time points.-newSignalInTimes :: EventQueue -> [Double] -> Dynamics (Signal Double)-newSignalInTimes q xs =-  do s <- liftSimulation $ newSignalSource q-     enqueueWithTimes q xs $ triggerSignalWithTime s-     return $ publishSignal s-       --- | Return a signal that is triggered in the integration time points.--- It should be called with help of 'runDynamicsInStartTime'.-newSignalInIntegTimes :: EventQueue -> Dynamics (Signal Double)-newSignalInIntegTimes q =-  do s <- liftSimulation $ newSignalSource q-     enqueueWithIntegTimes q $ triggerSignalWithTime s-     return $ publishSignal s-     --- | Return a signal that is triggered in the start time.--- It should be called with help of 'runDynamicsInStartTime'.-newSignalInStartTime :: EventQueue -> Dynamics (Signal Double)-newSignalInStartTime q =-  do s <- liftSimulation $ newSignalSource q-     enqueueWithStartTime q $ triggerSignalWithTime s-     return $ publishSignal s---- | Return a signal that is triggered in the stop time.-newSignalInStopTime :: EventQueue -> Dynamics (Signal Double)-newSignalInStopTime q =-  do s <- liftSimulation $ newSignalSource q-     enqueueWithStopTime q $ triggerSignalWithTime s-     return $ publishSignal s

@@ -1,29 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Simulation--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ The module defines the 'Simulation' monad representing a simulation run.----module Simulation.Aivika.Dynamics.Simulation-       (-- * Simulation-        Simulation,-        SimulationLift(..),-        Specs(..),-        Method(..),-        runSimulation,-        runSimulations,-        -- * Error Handling-        catchSimulation,-        finallySimulation,-        throwSimulation,-        -- * Utilities-        simulationIndex,-        simulationCount,-        simulationSpecs) where--import Simulation.Aivika.Dynamics.Internal.Simulation

@@ -1,788 +0,0 @@--{-# LANGUAGE FlexibleContexts, BangPatterns, RecursiveDo #-}---- |--- Module     : Simulation.Aivika.Dynamics.SystemDynamics--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.3------ This module defines integrals and other functions of System Dynamics.-----module Simulation.Aivika.Dynamics.SystemDynamics-       (-- * Equality and Ordering-        (.==.),-        (./=.),-        (.<.),-        (.>=.),-        (.>.),-        (.<=.),-        maxDynamics,-        minDynamics,-        ifDynamics,-        -- * Integrals-        Integ,-        newInteg,-        integInit,-        integValue,-        integDiff,-        -- * Integral Functions-        integ,-        smoothI,-        smooth,-        smooth3I,-        smooth3,-        smoothNI,-        smoothN,-        delay1I,-        delay1,-        delay3I,-        delay3,-        delayNI,-        delayN,-        forecast,-        trend,-        -- * Difference Equations-        Sum,-        newSum,-        sumInit,-        sumValue,-        sumDiff,-        sumDynamics,-        -- * Table Functions-        lookupD,-        lookupStepwiseD,-        lookupDynamics,-        lookupStepwiseDynamics,-        -- * Discrete Functions-        delayTrans,-        delay,-        delayI,-        udelay,-        udelayI,-        -- * Financial Functions-        npv,-        npve) where--import Data.Array-import Data.Array.IO.Safe-import Data.IORef-import Control.Monad-import Control.Monad.Trans--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.Base------- Equality and Ordering------- | Compare for equality.-(.==.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool-(.==.) = liftM2 (==)---- | Compare for inequality.-(./=.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool-(./=.) = liftM2 (/=)---- | Compare for ordering.-(.<.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool-(.<.) = liftM2 (<)---- | Compare for ordering.-(.>=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool-(.>=.) = liftM2 (>=)---- | Compare for ordering.-(.>.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool-(.>.) = liftM2 (>)---- | Compare for ordering.-(.<=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool-(.<=.) = liftM2 (<=)---- | Return the maximum.-maxDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a-maxDynamics = liftM2 max---- | Return the minimum.-minDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a-minDynamics = liftM2 min---- | Implement the if-then-else operator.-ifDynamics :: Dynamics Bool -> Dynamics a -> Dynamics a -> Dynamics a-ifDynamics cond x y =-  do a <- cond-     if a then x else y------- Integrals-----{-# DEPRECATED Integ "Use the integ function instead" #-}-{-# DEPRECATED newInteg "Use the integ function instead" #-}-{-# DEPRECATED integInit "Use the integ function instead" #-}-{-# DEPRECATED integValue "Use the integ function instead" #-}-{-# DEPRECATED integDiff "Use the integ function instead" #-}---- | 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 -> Simulation Integ-newInteg i = -  do r1 <- liftIO $ newIORef $ initDynamics i -     r2 <- liftIO $ newIORef $ initDynamics 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 -> Simulation ()-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"---- | Return an integral with the specified derivative and initial value.------ To create a loopback, you should use the recursive do-notation.--- It allows defining the differential equations unordered as--- in mathematics:------ @--- model :: Simulation [Double]--- model = ---   mdo a <- integ (- ka * a) 100---       b <- integ (ka * a - kb * b) 0---       c <- integ (kb * b) 0---       let ka = 1---           kb = 1---       runDynamicsInStopTime $ sequence [a, b, c]--- @-integ :: Dynamics Double                  -- ^ the derivative-         -> Dynamics Double               -- ^ the initial value-         -> Simulation (Dynamics Double)  -- ^ the integral-integ diff i =-  mdo y <- umemo z-      z <- Simulation $ \r ->-        case spcMethod (runSpecs r) of-          Euler -> return $ Dynamics $ integEuler diff i y-          RungeKutta2 -> return $ Dynamics $ integRK2 diff i y-          RungeKutta4 -> return $ Dynamics $ integRK4 diff i y-      return y---- | Return the first order exponential smooth.------ To create a loopback, you should use the recursive do-notation--- with help of which the function itself is defined:------ @--- smoothI x t i =---   mdo y <- integ ((x - y) \/ t) i---       return y--- @     -smoothI :: Dynamics Double                  -- ^ the value to smooth over time-           -> Dynamics Double               -- ^ time-           -> Dynamics Double               -- ^ the initial value-           -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth-smoothI x t i =-  mdo y <- integ ((x - y) / t) i-      return y---- | Return the first order exponential smooth.------ This is a simplified version of the 'smoothI' function--- without specifing the initial value.-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---- | Return the third order exponential smooth.------ To create a loopback, you should use the recursive do-notation--- with help of which the function itself is defined:------ @--- smooth3I x t i =---   mdo y  <- integ ((s2 - y) \/ t') i---       s2 <- integ ((s1 - s2) \/ t') i---       s1 <- integ ((x - s1) \/ t') i---       let t' = t \/ 3.0---       return y--- @     -smooth3I :: Dynamics Double                  -- ^ the value to smooth over time-            -> Dynamics Double               -- ^ time-            -> Dynamics Double               -- ^ the initial value-            -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth-smooth3I x t i =-  mdo y  <- integ ((s2 - y) / t') i-      s2 <- integ ((s1 - s2) / t') i-      s1 <- integ ((x - s1) / t') i-      let t' = t / 3.0-      return y---- | Return the third order exponential smooth.--- --- This is a simplified version of the 'smooth3I' function--- without specifying the initial value.-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---- | Return the n'th order exponential smooth.------ The result is not discrete in that sense that it may change within the integration time--- interval depending on the integration method used. Probably, you should apply--- the 'discrete' function to the result if you want to achieve an effect when the value is--- not changed within the time interval, which is used sometimes.-smoothNI :: Dynamics Double                  -- ^ the value to smooth over time-            -> Dynamics Double               -- ^ time-            -> Int                           -- ^ the order-            -> Dynamics Double               -- ^ the initial value-            -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth-smoothNI x t n i =-  mdo s <- forM [1 .. n] $ \k ->-        if k == 1-        then integ ((x - a ! 1) / t') i-        else integ ((a ! (k - 1) - a ! k) / t') i-      let a  = listArray (1, n) s -          t' = t / fromIntegral n-      return $ a ! n---- | Return the n'th order exponential smooth.------ This is a simplified version of the 'smoothNI' function--- without specifying the initial value.-smoothN :: Dynamics Double                  -- ^ the value to smooth over time-           -> Dynamics Double               -- ^ time-           -> Int                           -- ^ the order-           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth-smoothN x t n = smoothNI x t n x---- | Return the first order exponential delay.------ To create a loopback, you should use the recursive do-notation--- with help of which the function itself is defined:------ @--- delay1I x t i =---   mdo y <- integ (x - y \/ t) (i * t)---       return $ y \/ t--- @     -delay1I :: Dynamics Double                  -- ^ the value to conserve-           -> Dynamics Double               -- ^ time-           -> Dynamics Double               -- ^ the initial value-           -> Simulation (Dynamics Double)  -- ^ the first order exponential delay-delay1I x t i =-  mdo y <- integ (x - y / t) (i * t)-      return $ y / t---- | Return the first order exponential delay.------ This is a simplified version of the 'delay1I' function--- without specifying the initial value.-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---- | Return the third order exponential delay.-delay3I :: Dynamics Double                  -- ^ the value to conserve-           -> Dynamics Double               -- ^ time-           -> Dynamics 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')-      let t' = t / 3.0-      return $ y / t'         ---- | Return the third order exponential delay.------ This is a simplified version of the 'delay3I' function--- without specifying the initial value.-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---- | 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-           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay-delayNI x t n i =-  mdo 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')-      let a  = listArray (1, n) s-          t' = t / fromIntegral n-      return $ (a ! n) / t'---- | Return the n'th order exponential delay.------ This is a simplified version of the 'delayNI' function--- without specifying the initial value.-delayN :: Dynamics Double                  -- ^ the value to conserve-          -> Dynamics Double               -- ^ time-          -> Int                           -- ^ the order-          -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay-delayN x t n = delayNI x t n x---- | Return the forecast.------ The function has the following definition:------ @--- forecast x at hz =---   do y <- smooth x at---      return $ x * (1.0 + (x \/ y - 1.0) \/ at * hz)--- @-forecast :: Dynamics Double                  -- ^ the value to forecast-            -> Dynamics Double               -- ^ the average time-            -> Dynamics Double               -- ^ the time horizon-            -> Simulation (Dynamics Double)  -- ^ the forecast-forecast x at hz =-  do y <- smooth x at-     return $ x * (1.0 + (x / y - 1.0) / at * hz)---- | Return the trend.------ The function has the following definition:------ @--- trend x at i =---   do y <- smoothI x at (x \/ (1.0 + i * at))---      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-         -> 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------- Difference Equations-----{-# DEPRECATED Sum "Use the sumDynamics function instead" #-}-{-# DEPRECATED newSum "Use the sumDynamics function instead" #-}-{-# DEPRECATED sumInit "Use the sumDynamics function instead" #-}-{-# DEPRECATED sumValue "Use the sumDynamics function instead" #-}-{-# DEPRECATED sumDiff "Use the sumDynamics function instead" #-}---- | 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 -> Simulation (Sum a)-newSum i =   -  do r1 <- liftIO $ newIORef $ initDynamics i -     r2 <- liftIO $ newIORef $ initDynamics 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 -> Simulation ()-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---- | Retun the sum for the difference equation.--- It is like an integral returned by the 'integ' function, only now--- the difference is used instead of derivative.------ As usual, to create a loopback, you should use the recursive do-notation.-sumDynamics :: (MArray IOUArray a IO, Num a)-               => Dynamics a               -- ^ the difference-               -> Dynamics a               -- ^ the initial value-               -> Simulation (Dynamics a)  -- ^ the sum-sumDynamics (Dynamics diff) (Dynamics i) =-  mdo y <- umemo z-      z <- Simulation $ \r ->-        return $ Dynamics $ \p ->-        case pointIteration p of-          0 -> i p-          n -> do -            let Dynamics m = y-                sc = pointSpecs p-                ty = basicTime sc (n - 1) 0-                py = p { pointTime = ty, -                         pointIteration = n - 1, -                         pointPhase = 0 }-            a <- m py-            b <- diff py-            let !v = a + b-            return v-      return y------- Table Functions-----{-# DEPRECATED lookupD "Use the lookupDynamics function instead" #-}-{-# DEPRECATED lookupStepwiseD "Use the lookupStepwiseDynamics function instead" #-}---- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.-lookupD :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double-lookupD = lookupDynamics---- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function.-lookupStepwiseD :: Dynamics Double-                   -> Array Int (Double, Double)-                   -> Dynamics Double-lookupStepwiseD = lookupStepwiseDynamics---- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.-lookupDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double-lookupDynamics (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.-lookupStepwiseDynamics :: Dynamics Double-                          -> Array Int (Double, Double)-                          -> Dynamics Double-lookupStepwiseDynamics (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------- | Return the delayed value. This is a general version using the specified transform,--- usually a memoization.-delayTrans :: Dynamics a                                  -- ^ the value to delay-              -> Dynamics Double                          -- ^ the lag time-              -> Dynamics a                               -- ^ the initial value-              -> (Dynamics a -> Simulation (Dynamics a))  -- ^ the transform (usually, a memoization)-              -> Simulation (Dynamics a)                  -- ^ the delayed value-delayTrans (Dynamics x) (Dynamics d) (Dynamics i) tr = tr $ Dynamics r -  where-    r p = do -      let t  = pointTime p-          sc = pointSpecs p-          n  = pointIteration p-      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 }-            | n' < n    = x $ p { pointTime = t',-                                  pointIteration = n',-                                  pointPhase = -1 }-            | n' > n    = error $-                          "Cannot return the future data: delayTrans. " ++-                          "The lag time cannot be negative."-            | otherwise = error $-                          "Cannot return the current data: delayTrans. " ++-                          "The lag time is too small."-      y---- | Return the delayed value.------ It is defined in the following way:------ @ delay x d = delayTrans x d x memo0 @-delay :: Dynamics a                  -- ^ the value to delay-         -> Dynamics Double          -- ^ the lag time-         -> Simulation (Dynamics a)  -- ^ the delayed value-delay x d = delayTrans x d x memo0---- | Return the delayed value.------ It is defined in the following way:------ @ delayI x d i = delayTrans x d i memo0 @-delayI :: Dynamics a                  -- ^ the value to delay-          -> Dynamics Double          -- ^ the lag time-          -> Dynamics a               -- ^ the initial value-          -> Simulation (Dynamics a)  -- ^ the delayed value-delayI x d i = delayTrans x d i memo0---- | Return the delayed value. This is a more efficient unboxed version of the 'delay' function.------ It is defined in the following way:------ @ udelay x d = delayTrans x d x umemo0 @-udelay :: (MArray IOUArray a IO, Num a)-          => Dynamics a               -- ^ the value to delay-          -> Dynamics Double          -- ^ the lag time-          -> Simulation (Dynamics a)  -- ^ the delayed value-udelay x d = delayTrans x d x umemo0---- | Return the delayed value. This is a more efficient unboxed version of the 'delayI' function.------ It is defined in the following way:------ @ udelayI x d i = delayTrans x d i umemo0 @-udelayI :: (MArray IOUArray a IO, Num a)-           => Dynamics a               -- ^ the value to delay-           -> Dynamics Double          -- ^ the lag time-           -> Dynamics a               -- ^ the initial value-           -> Simulation (Dynamics a)  -- ^ the delayed value-udelayI x d i = delayTrans x d i umemo0------- Financial Functions------- | Return the Net Present Value (NPV) of the stream computed using the specified--- discount rate, the initial value and some factor (usually 1).------ It is defined in the following way:------ @--- npv stream rate init factor =---   mdo 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-       -> Dynamics Double               -- ^ factor-       -> Simulation (Dynamics Double)  -- ^ the Net Present Value (NPV)-npv stream rate init factor =-  mdo df <- integ (- df * rate) 1-      accum <- integ (stream * df) init-      return $ (accum + dt * stream * df) * factor---- | Return the Net Present Value End of period (NPVE) of the stream computed--- using the specified discount rate, the initial value and some factor.------ It is defined in the following way:------ @--- npve stream rate init factor =---   mdo df <- integ (- df * rate \/ (1 + rate * dt)) (1 \/ (1 + rate * dt))---       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-        -> Dynamics Double               -- ^ factor-        -> Simulation (Dynamics Double)  -- ^ the Net Present Value End (NPVE)-npve stream rate init factor =-  mdo df <- integ (- df * rate / (1 + rate * dt)) (1 / (1 + rate * dt))-      accum <- integ (stream * df) init-      return $ (accum + dt * stream * df) * factor

@@ -1,164 +0,0 @@--{-# LANGUAGE FlexibleContexts #-}---- |--- Module     : Simulation.Aivika.Dynamics.UVar--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.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,-        uvarQueue,-        uvarChanged,-        uvarChanged_,-        newUVar,-        readUVar,-        writeUVar,-        modifyUVar,-        freezeUVar) where--import Control.Monad-import Data.Array-import Data.Array.IO.Safe-import Data.IORef--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal--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,-         uvarChangedSource :: SignalSource a, -         uvarUpdatedSource :: SignalSource a }-     --- | Create a new variable bound to the specified event queue.-newUVar :: (MArray IOUArray a IO) => EventQueue -> a -> Simulation (UVar a)-newUVar q a =-  Simulation $ \r ->-  do xs <- UV.newVector-     ys <- UV.newVector-     UV.appendVector xs $ spcStartTime $ runSpecs r-     UV.appendVector ys a-     s  <- invokeSimulation r newSignalSourceUnsafe-     u  <- invokeSimulation r $ newSignalSource q-     return UVar { uvarQueue = q,-                   uvarRun   = runQueue q,-                   uvarXS = xs,-                   uvarYS = ys, -                   uvarChangedSource = s, -                   uvarUpdatedSource = u }---- | 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 invokeDynamics p $ uvarRun v-     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-         s  = uvarChangedSource v-     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-     invokeDynamics p $ triggerSignal s 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 invokeDynamics p $ uvarRun v-     let xs = uvarXS v-         ys = uvarYS v-         t  = pointTime p-         s  = uvarChangedSource v-     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-                    let b = f a-                    UV.writeVector ys i $! b-                    invokeDynamics p $ triggerSignal s b-            else do i <- UV.vectorBinarySearch xs t-                    if i >= 0-                      then do a <- UV.readVector ys i-                              let b = f a-                              UV.appendVector xs t-                              UV.appendVector ys $! b-                              invokeDynamics p $ triggerSignal s b-                      else do a <- UV.readVector ys $ - (i + 1) - 1-                              let b = f a-                              UV.appendVector xs t-                              UV.appendVector ys $! b-                              invokeDynamics p $ triggerSignal s b---- | 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 invokeDynamics p $ uvarRun v-     xs <- UV.freezeVector (uvarXS v)-     ys <- UV.freezeVector (uvarYS v)-     return (xs, ys)-     --- | Return a signal that notifies about every change of the variable state.-uvarChanged :: UVar a -> Signal a-uvarChanged v = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (uvarUpdatedSource v)-        m2 = publishSignal (uvarChangedSource v)---- | Return a signal that notifies about every change of the variable state.-uvarChanged_ :: UVar a -> Signal ()-uvarChanged_ v = mapSignal (const ()) $ uvarChanged v          --invokeDynamics :: Point -> Dynamics a -> IO a-{-# INLINE invokeDynamics #-}-invokeDynamics p (Dynamics m) = m p--invokeSimulation :: Run -> Simulation a -> IO a-{-# INLINE invokeSimulation #-}-invokeSimulation r (Simulation m) = m r

@@ -1,163 +0,0 @@---- |--- Module     : Simulation.Aivika.Dynamics.Var--- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>--- License    : BSD3--- Maintainer : David Sorokin <david.sorokin@gmail.com>--- Stability  : experimental--- Tested with: GHC 7.6.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,-        varQueue,-        varChanged,-        varChanged_,-        newVar,-        readVar,-        writeVar,-        modifyVar,-        freezeVar) where--import Data.Array-import Data.Array.IO.Safe-import Data.IORef--import Simulation.Aivika.Dynamics.Internal.Simulation-import Simulation.Aivika.Dynamics.Internal.Dynamics-import Simulation.Aivika.Dynamics.EventQueue-import Simulation.Aivika.Dynamics.Internal.Signal-import Simulation.Aivika.Dynamics.Signal--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,-        varChangedSource :: SignalSource a, -        varUpdatedSource :: SignalSource a }-     --- | Create a new variable bound to the specified event queue.-newVar :: EventQueue -> a -> Simulation (Var a)-newVar q a =-  Simulation $ \r ->-  do xs <- UV.newVector-     ys <- V.newVector-     UV.appendVector xs $ spcStartTime $ runSpecs r-     V.appendVector ys a-     s  <- invokeSimulation r newSignalSourceUnsafe-     u  <- invokeSimulation r $ newSignalSource q-     return Var { varQueue = q,-                  varRun   = runQueue q,-                  varXS = xs,-                  varYS = ys, -                  varChangedSource = s, -                  varUpdatedSource = u }---- | 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 invokeDynamics p $ varRun v-     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-         s  = varChangedSource v-     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-     invokeDynamics p $ triggerSignal s 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 invokeDynamics p $ varRun v-     let xs = varXS v-         ys = varYS v-         t  = pointTime p-         s  = varChangedSource v-     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-                    let b = f a-                    V.writeVector ys i $! b-                    invokeDynamics p $ triggerSignal s b-            else do i <- UV.vectorBinarySearch xs t-                    if i >= 0-                      then do a <- V.readVector ys i-                              let b = f a-                              UV.appendVector xs t-                              V.appendVector ys $! b-                              invokeDynamics p $ triggerSignal s b-                      else do a <- V.readVector ys $ - (i + 1) - 1-                              let b = f a-                              UV.appendVector xs t-                              V.appendVector ys $! b-                              invokeDynamics p $ triggerSignal s b---- | 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 invokeDynamics p $ varRun v-     xs <- UV.freezeVector (varXS v)-     ys <- V.freezeVector (varYS v)-     return (xs, ys)-     --- | Return a signal that notifies about every change of the variable state.-varChanged :: Var a -> Signal a-varChanged v = merge2Signals m1 m2    -- N.B. The order is important!-  where m1 = publishSignal (varUpdatedSource v)-        m2 = publishSignal (varChangedSource v)---- | Return a signal that notifies about every change of the variable state.-varChanged_ :: Var a -> Signal ()-varChanged_ v = mapSignal (const ()) $ varChanged v     --invokeDynamics :: Point -> Dynamics a -> IO a-{-# INLINE invokeDynamics #-}-invokeDynamics p (Dynamics m) = m p--invokeSimulation :: Run -> Simulation a -> IO a-{-# INLINE invokeSimulation #-}-invokeSimulation r (Simulation m) = m r

@@ -0,0 +1,36 @@++-- |+-- Module     : Simulation.Aivika.Event+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- The module defines the 'Event' monad which is very similar to the 'Dynamics'+-- monad but only now the computation is strongly synchronized with the event queue.+--+module Simulation.Aivika.Event+       (-- * Event Monad+        Event,+        EventLift(..),+        EventProcessing(..),+        EventCancellation(..),+        runEvent,+        runEventInStartTime,+        runEventInStopTime,+        -- * Event Queue+        enqueueEvent,+        enqueueEventWithCancellation,+        enqueueEventWithTimes,+        enqueueEventWithIntegTimes,+        enqueueEventWithStartTime,+        enqueueEventWithStopTime,+        enqueueEventWithCurrentTime,+        eventQueueCount,+        -- * Error Handling+        catchEvent,+        finallyEvent,+        throwEvent) where++import Simulation.Aivika.Internal.Event

@@ -0,0 +1,320 @@++-- |+-- Module     : Simulation.Aivika.Internal.Cont+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- The 'Cont' monad is a variation of the standard Cont monad +-- and F# async workflow, where the result of applying +-- the continuations is the 'Event' computation.+--+module Simulation.Aivika.Internal.Cont+       (Cont(..),+        ContParams,+        invokeCont,+        runCont,+        catchCont,+        finallyCont,+        throwCont,+        resumeCont,+        contCanceled) where++import Data.IORef++import qualified Control.Exception as C+import Control.Exception (IOException, throw)++import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Dynamics+import Simulation.Aivika.Internal.Event++-- | The 'Cont' type is similar to the standard Cont monad +-- and F# async workflow but only the result of applying+-- the continuations return the 'Event' computation.+newtype Cont a = Cont (ContParams a -> Event ())++-- | The continuation parameters.+data ContParams a = +  ContParams { contCont :: a -> Event (), +               contAux  :: ContParamsAux }++-- | The auxiliary continuation parameters.+data ContParamsAux =+  ContParamsAux { contECont :: IOException -> Event (),+                  contCCont :: () -> Event (),+                  contCancelToken :: IORef Bool,+                  contCatchFlag   :: Bool }++instance Monad Cont where+  return  = returnC+  m >>= k = bindC m k++instance SimulationLift Cont where+  liftSimulation = liftSC++instance DynamicsLift Cont where+  liftDynamics = liftDC++instance EventLift Cont where+  liftEvent = liftEC++instance Functor Cont where+  fmap = liftM++instance MonadIO Cont where+  liftIO = liftIOC ++invokeCont :: ContParams a -> Cont a -> Event ()+{-# INLINE invokeCont #-}+invokeCont p (Cont m) = m p++cancelCont :: Point -> ContParams a -> IO ()+{-# NOINLINE cancelCont #-}+cancelCont p c =+  do writeIORef (contCancelToken $ contAux c) False+     invokeEvent p $ (contCCont $ contAux c) ()++returnC :: a -> Cont a+{-# INLINE returnC #-}+returnC a = +  Cont $ \c ->+  Event $ \p ->+  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ contCont c a+                          +-- bindC :: Cont a -> (a -> Cont b) -> Cont b+-- {-# INLINE bindC #-}+-- bindC m k = +--   Cont $ \c -> +--   if (contCatchFlag . contAux $ c) +--   then bindWithCatch m k c+--   else bindWithoutCatch m k c+  +bindC :: Cont a -> (a -> Cont b) -> Cont b+{-# INLINE bindC #-}+bindC m k = +  Cont $ bindWithoutCatch m k  -- Another version is not tail recursive!+  +bindWithoutCatch :: Cont a -> (a -> Cont b) -> ContParams b -> Event ()+{-# INLINE bindWithoutCatch #-}+bindWithoutCatch (Cont m) k c = +  Event $ \p ->+  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ m $ +            let cont a = invokeCont c (k a)+            in c { contCont = cont }++-- It is not tail recursive!+bindWithCatch :: Cont a -> (a -> Cont b) -> ContParams b -> Event ()+{-# NOINLINE bindWithCatch #-}+bindWithCatch (Cont m) k c = +  Event $ \p ->+  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ m $ +            let cont a = catchEvent +                         (invokeCont c (k a))+                         (contECont $ contAux c)+            in c { contCont = cont }++-- Like "bindWithoutCatch (return a) k"+callWithoutCatch :: (a -> Cont b) -> a -> ContParams b -> Event ()+callWithoutCatch k a c =+  Event $ \p ->+  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ invokeCont c (k a)++-- Like "bindWithCatch (return a) k" but it is not tail recursive!+callWithCatch :: (a -> Cont b) -> a -> ContParams b -> Event ()+callWithCatch k a c =+  Event $ \p ->+  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ catchEvent +            (invokeCont c (k a))+            (contECont $ contAux c)++-- | Exception handling within 'Cont' computations.+catchCont :: Cont a -> (IOException -> Cont a) -> Cont a+catchCont m h = +  Cont $ \c -> +  if contCatchFlag . contAux $ c+  then catchWithCatch m h c+  else error $+       "To catch exceptions, the process must be created " +++       "with help of newProcessIDWithCatch: catchCont."+  +catchWithCatch :: Cont a -> (IOException -> Cont a) -> ContParams a -> Event ()+catchWithCatch (Cont m) h c =+  Event $ \p -> +  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ m $+            -- let econt e = callWithCatch h e c   -- not tail recursive!+            let econt e = callWithoutCatch h e c+            in c { contAux = (contAux c) { contECont = econt } }+               +-- | A computation with finalization part.+finallyCont :: Cont a -> Cont b -> Cont a+finallyCont m m' = +  Cont $ \c -> +  if contCatchFlag . contAux $ c+  then finallyWithCatch m m' c+  else error $+       "To finalize computation, the process must be created " +++       "with help of newProcessIdWithCatch: finallyCont."+  +finallyWithCatch :: Cont a -> Cont b -> ContParams a -> Event ()               +finallyWithCatch (Cont m) (Cont m') c =+  Event $ \p ->+  do z <- contCanceled c+     if z +       then cancelCont p c+       else invokeEvent p $ m $+            let cont a   = +                  Event $ \p ->+                  invokeEvent p $ m' $+                  let cont b = contCont c a+                  in c { contCont = cont }+                econt e  =+                  Event $ \p ->+                  invokeEvent p $ m' $+                  let cont b = (contECont . contAux $ c) e+                  in c { contCont = cont }+                ccont () = +                  Event $ \p ->+                  invokeEvent p $ m' $+                  let cont b  = (contCCont . contAux $ c) ()+                      econt e = (contCCont . contAux $ c) ()+                  in c { contCont = cont,+                         contAux  = (contAux c) { contECont = econt } }+            in c { contCont = cont,+                   contAux  = (contAux c) { contECont = econt,+                                            contCCont = ccont } }++-- | Throw the exception with the further exception handling.+-- By some reasons, the standard 'throw' function per se is not handled +-- properly within 'Cont' computations, altough it will be still handled +-- if it will be hidden under the 'liftIO' function. The problem arises +-- namely with the @throw@ function, not 'IO' computations.+throwCont :: IOException -> Cont a+throwCont e = liftIO $ throw e++-- | Run the 'Cont' computation with the specified cancelation token +-- and flag indicating whether to catch exceptions.+runCont :: Cont a+           -- ^ the computation to run+           -> (a -> Event ())+           -- ^ the main branch +           -> (IOError -> Event ())+           -- ^ the branch for handing exceptions+           -> (() -> Event ())+           -- ^ the branch for cancellation+           -> IORef Bool+           -- ^ the cancellation token+           -> Bool+           -- ^ whether to support the exception catching+           -> Event ()+runCont (Cont m) cont econt ccont cancelToken catchFlag = +  m ContParams { contCont = cont,+                 contAux  = +                   ContParamsAux { contECont = econt,+                                   contCCont = ccont,+                                   contCancelToken = cancelToken, +                                   contCatchFlag = catchFlag } }++-- | Lift the 'Simulation' computation.+liftSC :: Simulation a -> Cont a+liftSC (Simulation m) = +  Cont $ \c ->+  Event $ \p ->+  if contCatchFlag . contAux $ c+  then liftIOWithCatch (m $ pointRun p) p c+  else liftIOWithoutCatch (m $ pointRun p) p c+     +-- | Lift the 'Dynamics' computation.+liftDC :: Dynamics a -> Cont a+liftDC (Dynamics m) =+  Cont $ \c ->+  Event $ \p ->+  if contCatchFlag . contAux $ c+  then liftIOWithCatch (m p) p c+  else liftIOWithoutCatch (m p) p c+     +-- | Lift the 'Event' computation.+liftEC :: Event a -> Cont a+liftEC (Event m) =+  Cont $ \c ->+  Event $ \p ->+  if contCatchFlag . contAux $ c+  then liftIOWithCatch (m p) p c+  else liftIOWithoutCatch (m p) p c+     +-- | Lift the IO computation.+liftIOC :: IO a -> Cont a+liftIOC m =+  Cont $ \c ->+  Event $ \p ->+  if contCatchFlag . contAux $ c+  then liftIOWithCatch m p c+  else liftIOWithoutCatch m p c+  +liftIOWithoutCatch :: IO a -> Point -> ContParams a -> IO ()+{-# INLINE liftIOWithoutCatch #-}+liftIOWithoutCatch m p c =+  do z <- contCanceled c+     if z+       then cancelCont p c+       else do a <- m+               invokeEvent p $ contCont c a++liftIOWithCatch :: IO a -> Point -> ContParams a -> IO ()+{-# NOINLINE liftIOWithCatch #-}+liftIOWithCatch m p c =+  do z <- contCanceled c+     if z+       then cancelCont p c+       else do aref <- newIORef undefined+               eref <- newIORef Nothing+               C.catch (m >>= writeIORef aref) +                 (writeIORef eref . Just)+               e <- readIORef eref+               case e of+                 Nothing -> +                   do a <- readIORef aref+                      -- tail recursive+                      invokeEvent p $ contCont c a+                 Just e ->+                   -- tail recursive+                   invokeEvent p $ (contECont . contAux) c e++-- | Resume the computation by the specified parameters.+resumeCont :: ContParams a -> a -> Event ()+{-# INLINE resumeCont #-}+resumeCont c a = +  Event $ \p ->+  do z <- contCanceled c+     if z+       then cancelCont p c+       else invokeEvent p $ contCont c a++-- | Test whether the computation is canceled.+contCanceled :: ContParams a -> IO Bool+{-# INLINE contCanceled #-}+contCanceled c = readIORef $ contCancelToken $ contAux c

@@ -0,0 +1,217 @@++{-# LANGUAGE RecursiveDo #-}++-- |+-- Module     : Simulation.Aivika.Internal.Dynamics+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- The module defines the 'Dynamics' monad representing a time varying polymorphic function. +--+module Simulation.Aivika.Internal.Dynamics+       (-- * Dynamics+        Dynamics(..),+        DynamicsLift(..),+        invokeDynamics,+        runDynamicsInStartTime,+        runDynamicsInStopTime,+        runDynamicsInIntegTimes,+        runDynamicsInTime,+        runDynamicsInTimes,+        -- * Error Handling+        catchDynamics,+        finallyDynamics,+        throwDynamics,+        -- * Time parameters+        starttime,+        stoptime,+        dt,+        time,+        isTimeInteg,+        integIteration,+        integPhase) where++import qualified Control.Exception as C+import Control.Exception (IOException, throw, finally)++import Control.Monad+import Control.Monad.Trans+import Control.Monad.Fix++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation++-- | A value in the 'Dynamics' monad represents a polymorphic time varying function.+newtype Dynamics a = Dynamics (Point -> IO a)++instance Monad Dynamics where+  return  = returnD+  m >>= k = bindD m k++returnD :: a -> Dynamics a+{-# INLINE returnD #-}+returnD a = Dynamics (\p -> return a)++bindD :: Dynamics a -> (a -> Dynamics b) -> Dynamics b+{-# INLINE bindD #-}+bindD (Dynamics m) k = +  Dynamics $ \p -> +  do a <- m p+     let Dynamics m' = k a+     m' p++-- | Run the 'Dynamics' computation in the initial time point.+runDynamicsInStartTime :: Dynamics a -> Simulation a+runDynamicsInStartTime (Dynamics m) =+  Simulation $ m . integStartPoint++-- | Run the 'Dynamics' computation in the final time point.+runDynamicsInStopTime :: Dynamics a -> Simulation a+runDynamicsInStopTime (Dynamics m) =+  Simulation $ m . integStopPoint++-- | Run the 'Dynamics' computation in all integration time points.+runDynamicsInIntegTimes :: Dynamics a -> Simulation [IO a]+runDynamicsInIntegTimes (Dynamics m) =+  Simulation $ return . map m . integPoints++-- | Run the 'Dynamics' computation in the specified time point.+runDynamicsInTime :: Double -> Dynamics a -> Simulation a+runDynamicsInTime t (Dynamics m) =+  Simulation $ \r -> m $ pointAt r t++-- | Run the 'Dynamics' computation in the specified time points.+runDynamicsInTimes :: [Double] -> Dynamics a -> Simulation [IO a]+runDynamicsInTimes ts (Dynamics m) =+  Simulation $ \r -> return $ map (m . pointAt r) ts ++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++instance SimulationLift Dynamics where+  liftSimulation = liftDS+    +liftDS :: Simulation a -> Dynamics a+{-# INLINE liftDS #-}+liftDS (Simulation m) =+  Dynamics $ \p -> m $ pointRun p++-- | A type class to lift the 'Dynamics' computations to other monads.+class Monad m => DynamicsLift m where+  +  -- | Lift the specified 'Dynamics' computation to another monad.+  liftDynamics :: Dynamics a -> m a++instance DynamicsLift Dynamics where+  liftDynamics = id+  +-- | Exception handling within 'Dynamics' computations.+catchDynamics :: Dynamics a -> (IOException -> Dynamics a) -> Dynamics a+catchDynamics (Dynamics m) h =+  Dynamics $ \p -> +  C.catch (m p) $ \e ->+  let Dynamics m' = h e in m' p+                           +-- | A computation with finalization part like the 'finally' function.+finallyDynamics :: Dynamics a -> Dynamics b -> Dynamics a+finallyDynamics (Dynamics m) (Dynamics m') =+  Dynamics $ \p ->+  C.finally (m p) (m' p)++-- | Like the standard 'throw' function.+throwDynamics :: IOException -> Dynamics a+throwDynamics = throw++-- | Invoke the 'Dynamics' computation.+invokeDynamics :: Point -> Dynamics a -> IO a+{-# INLINE invokeDynamics #-}+invokeDynamics p (Dynamics m) = m p++instance MonadFix Dynamics where+  mfix f = +    Dynamics $ \p ->+    do { rec { a <- invokeDynamics p (f a) }; return a }++-- | 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 ++-- | Whether the current time is an integration time.+isTimeInteg :: Dynamics Bool+isTimeInteg = Dynamics $ \p -> return $ pointPhase p >= 0++-- | Return the integration iteration closest to the current simulation time.+integIteration :: Dynamics Int+integIteration = Dynamics $ return . pointIteration++-- | Return the integration phase for the current simulation time.+-- It is @(-1)@ for non-integration time points.+integPhase :: Dynamics Int+integPhase = Dynamics $ return . pointPhase

@@ -0,0 +1,349 @@++{-# LANGUAGE RecursiveDo #-}++-- |+-- Module     : Simulation.Aivika.Internal.Event+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- The module defines the 'Event' monad which is very similar to the 'Dynamics'+-- monad but only now the computation is strongly synchronized with the event queue.+--+module Simulation.Aivika.Internal.Event+       (-- * Event Monad+        Event(..),+        EventLift(..),+        EventProcessing(..),+        EventCancellation(..),+        invokeEvent,+        runEvent,+        runEventInStartTime,+        runEventInStopTime,+        -- * Event Queue+        enqueueEvent,+        enqueueEventWithCancellation,+        enqueueEventWithTimes,+        enqueueEventWithPoints,+        enqueueEventWithIntegTimes,+        enqueueEventWithStartTime,+        enqueueEventWithStopTime,+        enqueueEventWithCurrentTime,+        eventQueueCount,+        -- * Error Handling+        catchEvent,+        finallyEvent,+        throwEvent) where++import Data.IORef++import qualified Control.Exception as C+import Control.Exception (IOException, throw, finally)++import Control.Monad+import Control.Monad.Trans+import Control.Monad.Fix++import qualified Simulation.Aivika.PriorityQueue as PQ++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Dynamics++-- | A value in the 'Event' monad represents a polymorphic time varying function+-- which is strongly synchronized with the event queue.+newtype Event a = Event (Point -> IO a)++instance Monad Event where+  return  = returnE+  m >>= k = bindE m k++returnE :: a -> Event a+{-# INLINE returnE #-}+returnE a = Event (\p -> return a)++bindE :: Event a -> (a -> Event b) -> Event b+{-# INLINE bindE #-}+bindE (Event m) k = +  Event $ \p -> +  do a <- m p+     let Event m' = k a+     m' p++instance Functor Event where+  fmap = liftME++liftME :: (a -> b) -> Event a -> Event b+{-# INLINE liftME #-}+liftME f (Event x) =+  Event $ \p -> do { a <- x p; return $ f a }++instance MonadIO Event where+  liftIO m = Event $ const m++instance SimulationLift Event where+  liftSimulation = liftES++instance DynamicsLift Event where+  liftDynamics = liftDS+    +liftES :: Simulation a -> Event a+{-# INLINE liftES #-}+liftES (Simulation m) =+  Event $ \p -> m $ pointRun p++liftDS :: Dynamics a -> Event a+{-# INLINE liftDS #-}+liftDS (Dynamics m) =+  Event m++-- | A type class to lift the 'Event' computation to other monads.+class Monad m => EventLift m where+  +  -- | Lift the specified 'Event' computation to another monad.+  liftEvent :: Event a -> m a++instance EventLift Event where+  liftEvent = id+  +-- | Exception handling within 'Event' computations.+catchEvent :: Event a -> (IOException -> Event a) -> Event a+catchEvent (Event m) h =+  Event $ \p -> +  C.catch (m p) $ \e ->+  let Event m' = h e in m' p+                           +-- | A computation with finalization part like the 'finally' function.+finallyEvent :: Event a -> Event b -> Event a+finallyEvent (Event m) (Event m') =+  Event $ \p ->+  C.finally (m p) (m' p)++-- | Like the standard 'throw' function.+throwEvent :: IOException -> Event a+throwEvent = throw++-- | Invoke the 'Event' computation.+invokeEvent :: Point -> Event a -> IO a+{-# INLINE invokeEvent #-}+invokeEvent p (Event m) = m p++instance MonadFix Event where+  mfix f = +    Event $ \p ->+    do { rec { a <- invokeEvent p (f a) }; return a }++-- | Defines how the events are processed.+data EventProcessing = IncludingCurrentEvents+                       -- ^ either process all earlier and then current events,+                       -- or raise an error if the current simulation time is less+                       -- than the actual time of the event queue+                     | IncludingEarlierEvents+                       -- ^ either process all earlier events not affecting+                       -- the events at the current simulation time,+                       -- or raise an error if the current simulation time is less+                       -- than the actual time of the event queue+                     | IncludingCurrentEventsOrFromPast+                       -- ^ either process all earlier and then current events,+                       -- or do nothing if the current simulation time is less+                       -- than the actual time of the event queue+                       -- (do not use unless the documentation states the opposite)+                     | IncludingEarlierEventsOrFromPast+                       -- ^ either process all earlier events,+                       -- or do nothing if the current simulation time is less+                       -- than the actual time of the event queue+                       -- (do not use unless the documentation states the opposite)+                     deriving (Eq, Ord, Show)++-- | Enqueue the event which must be actuated at the specified time.+--+-- The events are processed when calling the 'runEvent' function. So,+-- if you want to insist on their immediate execution then you can apply+-- something like+--+-- @+--   liftDynamics $ runEvent IncludingCurrentEvents $ return ()+-- @+--+-- although this is generally not good idea.  +enqueueEvent :: Double -> Event () -> Event ()+enqueueEvent t (Event m) =+  Event $ \p ->+  let pq = queuePQ $ runEventQueue $ pointRun p+  in PQ.enqueue pq t m++-- | Process the pending events.+processPendingEventsCore :: Bool -> Dynamics ()+processPendingEventsCore includingCurrentEvents = Dynamics r where+  r p =+    do let q = runEventQueue $ pointRun p+           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+           r  = pointRun p+       f <- PQ.queueNull pq+       unless f $+         do (t2, c2) <- PQ.queueFront pq+            let t = queueTime q+            t' <- readIORef t+            when (t2 < t') $ +              error "The time value is too small: processPendingEventsCore"+            when ((t2 < pointTime p) ||+                  (includingCurrentEvents && (t2 == pointTime p))) $+              do writeIORef t t2+                 PQ.dequeue pq+                 let sc = pointSpecs p+                     t0 = spcStartTime sc+                     dt = spcDT sc+                     n2 = fromIntegral $ floor ((t2 - t0) / dt)+                 c2 $ p { pointTime = t2,+                          pointIteration = n2,+                          pointPhase = -1 }+                 call q p++-- | Process the pending events synchronously, i.e. without past.+processPendingEvents :: Bool -> Dynamics ()+processPendingEvents includingCurrentEvents = Dynamics r where+  r p =+    do let q = runEventQueue $ pointRun p+           t = queueTime q+       t' <- readIORef t+       if pointTime p < t'+         then error $+              "The current time is less than " +++              "the time in the queue: processPendingEvents"+         else invokeDynamics p m+  m = processPendingEventsCore includingCurrentEvents++-- | A memoized value.+processEventsIncludingCurrent = processPendingEvents True++-- | A memoized value.+processEventsIncludingEarlier = processPendingEvents False++-- | A memoized value.+processEventsIncludingCurrentCore = processPendingEventsCore True++-- | A memoized value.+processEventsIncludingEarlierCore = processPendingEventsCore True++-- | Process the events.+processEvents :: EventProcessing -> Dynamics ()+processEvents IncludingCurrentEvents = processEventsIncludingCurrent+processEvents IncludingEarlierEvents = processEventsIncludingEarlier+processEvents IncludingCurrentEventsOrFromPast = processEventsIncludingCurrentCore+processEvents IncludingEarlierEventsOrFromPast = processEventsIncludingEarlierCore++-- | Run the 'Event' computation in the current simulation time+-- within the 'Dynamics' computation.+runEvent :: EventProcessing -> Event a -> Dynamics a+runEvent processing (Event e) =+  Dynamics $ \p ->+  do invokeDynamics p $ processEvents processing+     e p++-- | Run the 'Event' computation in the start time.+runEventInStartTime :: EventProcessing -> Event a -> Simulation a+runEventInStartTime processing e =+  runDynamicsInStartTime $ runEvent processing e++-- | Run the 'Event' computation in the stop time.+runEventInStopTime :: EventProcessing -> Event a -> Simulation a+runEventInStopTime processing e =+  runDynamicsInStopTime $ runEvent processing e++-- | Return the number of pending events that should+-- be yet actuated.+eventQueueCount :: Event Int+eventQueueCount =+  Event $ PQ.queueCount . queuePQ . runEventQueue . pointRun++-- | Actuate the event handler in the specified time points.+enqueueEventWithTimes :: [Double] -> Event () -> Event ()+enqueueEventWithTimes ts e = loop ts+  where loop []       = return ()+        loop (t : ts) = enqueueEvent t $ e >> loop ts+       +-- | Actuate the event handler in the specified time points.+enqueueEventWithPoints :: [Point] -> Event () -> Event ()+enqueueEventWithPoints xs (Event e) = loop xs+  where loop []       = return ()+        loop (x : xs) = enqueueEvent (pointTime x) $ +                        Event $ \p ->+                        do e x    -- N.B. we substitute the time point!+                           invokeEvent p $ loop xs+                           +-- | Actuate the event handler in the integration time points.+enqueueEventWithIntegTimes :: Event () -> Event ()+enqueueEventWithIntegTimes e =+  Event $ \p ->+  let points = integPoints $ pointRun p+  in invokeEvent p $ enqueueEventWithPoints points e++-- | Actuate the event handler in the start time.+enqueueEventWithStartTime :: Event () -> Event ()+enqueueEventWithStartTime e =+  Event $ \p ->+  let point = integStartPoint $ pointRun p+  in invokeEvent p $ enqueueEventWithPoints [point] e++-- | Actuate the event handler in the stop time.+enqueueEventWithStopTime :: Event () -> Event ()+enqueueEventWithStopTime e =+  Event $ \p ->+  let point = integStopPoint $ pointRun p+  in invokeEvent p $ enqueueEventWithPoints [point] e++-- | Actuate the event handler in the current time but +-- through the event queue, which allows continuing the +-- current tasks and then calling the handler after the +-- tasks are finished. The simulation time will be the same.+enqueueEventWithCurrentTime :: Event () -> Event ()+enqueueEventWithCurrentTime e =+  Event $ \p ->+  invokeEvent p $ enqueueEvent (pointTime p) e++-- | It allows cancelling the event.+data EventCancellation =+  EventCancellation { cancelEvent   :: Event (),+                      -- ^ Cancel the event.+                      eventCanceled :: Event Bool,+                      -- ^ Test whether the event was canceled.+                      eventFinished :: Event Bool+                      -- ^ Test whether the event was processed and finished.+                    }++-- | Enqueue the event with an ability to cancel it.+enqueueEventWithCancellation :: Double -> Event () -> Event EventCancellation+enqueueEventWithCancellation t e =+  Event $ \p ->+  do canceledRef <- newIORef False+     cancellableRef <- newIORef True+     finishedRef <- newIORef False+     let cancel =+           Event $ \p ->+           do x <- readIORef cancellableRef+              when x $+                writeIORef canceledRef True+         canceled =+           Event $ \p -> readIORef canceledRef+         finished =+           Event $ \p -> readIORef finishedRef+     invokeEvent p $+       enqueueEvent t $+       Event $ \p ->+       do writeIORef cancellableRef False+          x <- readIORef canceledRef+          unless x $+            do invokeEvent p e+               writeIORef finishedRef True+     return EventCancellation { cancelEvent   = cancel,+                                eventCanceled = canceled,+                                eventFinished = finished }

@@ -0,0 +1,325 @@++-- |+-- Module     : Simulation.Aivika.Internal.Process+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- A value in the 'Process' monad represents a discontinuous process that +-- can suspend in any simulation time point and then resume later in the same +-- or another time point. +-- +-- The process of this type can involve the 'Event', 'Dynamics' and 'Simulation'+-- computations. Moreover, a value in the @Process@ monad can be run within+-- the @Event@ computation.+--+-- A value of the 'ProcessId' type is just an identifier of such a process.+--+module Simulation.Aivika.Internal.Process+       (ProcessId,+        Process(..),+        invokeProcess,+        runProcess,+        runProcessInStartTime,+        runProcessInStopTime,+        enqueueProcess,+        enqueueProcessWithStartTime,+        enqueueProcessWithStopTime,+        newProcessId,+        newProcessIdWithCatch,+        holdProcess,+        interruptProcess,+        processInterrupted,+        passivateProcess,+        processPassive,+        reactivateProcess,+        processId,+        cancelProcess,+        processCanceled,+        catchProcess,+        finallyProcess,+        throwProcess) where++import Data.Maybe+import Data.IORef+import Control.Exception (IOException, throw)+import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Dynamics+import Simulation.Aivika.Internal.Event+import Simulation.Aivika.Internal.Cont++-- | Represents a process identifier.+data ProcessId = +  ProcessId { processStarted :: IORef Bool,+              processCatchFlag     :: Bool,+              processReactCont     :: IORef (Maybe (ContParams ())), +              processCancelRef     :: IORef Bool, +              processCancelToken   :: IORef Bool,+              processInterruptRef  :: IORef Bool, +              processInterruptCont :: IORef (Maybe (ContParams ())), +              processInterruptVersion :: IORef Int }++-- | Specifies a discontinuous process that can suspend at any time+-- and then resume later.+newtype Process a = Process (ProcessId -> Cont a)++-- | Invoke the process computation.+invokeProcess :: ProcessId -> Process a -> Cont a+{-# INLINE invokeProcess #-}+invokeProcess pid (Process m) = m pid++-- | Hold the process for the specified time period.+holdProcess :: Double -> Process ()+holdProcess dt =+  Process $ \pid ->+  Cont $ \c ->+  Event $ \p ->+  do let x = processInterruptCont pid+     writeIORef x $ Just c+     writeIORef (processInterruptRef pid) False+     v <- readIORef (processInterruptVersion pid)+     invokeEvent p $+       enqueueEvent (pointTime p + dt) $+       Event $ \p ->+       do v' <- readIORef (processInterruptVersion pid)+          when (v == v') $ +            do writeIORef x Nothing+               invokeEvent p $ resumeCont c ()++-- | Interrupt a process with the specified identifier if the process+-- is held by computation 'holdProcess'.+interruptProcess :: ProcessId -> Event ()+interruptProcess pid =+  Event $ \p ->+  do let x = processInterruptCont pid+     a <- readIORef x+     case a of+       Nothing -> return ()+       Just c ->+         do writeIORef x Nothing+            writeIORef (processInterruptRef pid) True+            modifyIORef (processInterruptVersion pid) $ (+) 1+            invokeEvent p $ enqueueEvent (pointTime p) $ resumeCont c ()+            +-- | Test whether the process with the specified identifier was interrupted.+processInterrupted :: ProcessId -> Event Bool+processInterrupted pid =+  Event $ \p ->+  readIORef (processInterruptRef pid)++-- | Passivate the process.+passivateProcess :: Process ()+passivateProcess =+  Process $ \pid ->+  Cont $ \c ->+  Event $ \p ->+  do let x = processReactCont 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 identifier is passivated.+processPassive :: ProcessId -> Event Bool+processPassive pid =+  Event $ \p ->+  do let x = processReactCont pid+     a <- readIORef x+     return $ isJust a++-- | Reactivate a process with the specified identifier.+reactivateProcess :: ProcessId -> Event ()+reactivateProcess pid =+  Event $ \p ->+  do let x = processReactCont pid+     a <- readIORef x+     case a of+       Nothing -> +         return ()+       Just c ->+         do writeIORef x Nothing+            invokeEvent p $ enqueueEvent (pointTime p) $ resumeCont c ()++-- | Start immediately the process with the specified identifier.+--            +-- To run the process at the specified time, you can use+-- the 'enqueueProcess' function.+runProcess :: ProcessId -> Process () -> Event ()+runProcess pid p =+  runCont m cont econt ccont (processCancelToken pid) (processCatchFlag pid)+    where cont  = return+          econt = throwEvent+          ccont = return+          m = do y <- liftIO $ readIORef (processStarted pid)+                 if y +                   then error $+                        "Another process with this identifier " +++                        "has been started already: runProcess"+                   else liftIO $ writeIORef (processStarted pid) True+                 invokeProcess pid p++-- | Start the process in the start time immediately.+runProcessInStartTime :: EventProcessing -> ProcessId -> Process () -> Simulation ()+runProcessInStartTime processing pid p =+  runEventInStartTime processing $ runProcess pid p++-- | Start the process in the stop time immediately.+runProcessInStopTime :: EventProcessing -> ProcessId -> Process () -> Simulation ()+runProcessInStopTime processing pid p =+  runEventInStopTime processing $ runProcess pid p++-- | Enqueue the process that will be then started at the specified time+-- from the event queue.+enqueueProcess :: Double -> ProcessId -> Process () -> Event ()+enqueueProcess t pid p =+  enqueueEvent t $ runProcess pid p++-- | Enqueue the process that will be then started in the start time+-- from the event queue.+enqueueProcessWithStartTime :: ProcessId -> Process () -> Event ()+enqueueProcessWithStartTime pid p =+  enqueueEventWithStartTime $ runProcess pid p++-- | Enqueue the process that will be then started in the stop time+-- from the event queue.+enqueueProcessWithStopTime :: ProcessId -> Process () -> Event ()+enqueueProcessWithStopTime pid p =+  enqueueEventWithStopTime $ runProcess pid p++-- | Return the current process identifier.+processId :: Process ProcessId+processId = Process return++-- | Create a new process identifier without exception handling.+newProcessId :: Simulation ProcessId+newProcessId =+  do x <- liftIO $ newIORef Nothing+     y <- liftIO $ newIORef False+     c <- liftIO $ newIORef False+     t <- liftIO $ newIORef False+     i <- liftIO $ newIORef False+     z <- liftIO $ newIORef Nothing+     v <- liftIO $ newIORef 0+     return ProcessId { processStarted = y,+                        processCatchFlag     = False,+                        processReactCont     = x, +                        processCancelRef     = c, +                        processCancelToken   = t,+                        processInterruptRef  = i,+                        processInterruptCont = z, +                        processInterruptVersion = v }++-- | Create a new process identifier with capabilities of catching +-- the 'IOError' exceptions and finalizing the computation. +-- The corresponded process will be slower than that one+-- which identifier is created with help of 'newProcessId'.+newProcessIdWithCatch :: Simulation ProcessId+newProcessIdWithCatch =+  do x <- liftIO $ newIORef Nothing+     y <- liftIO $ newIORef False+     c <- liftIO $ newIORef False+     t <- liftIO $ newIORef False+     i <- liftIO $ newIORef False+     z <- liftIO $ newIORef Nothing+     v <- liftIO $ newIORef 0+     return ProcessId { processStarted = y,+                        processCatchFlag     = True,+                        processReactCont     = x, +                        processCancelRef     = c, +                        processCancelToken   = t,+                        processInterruptRef  = i,+                        processInterruptCont = z, +                        processInterruptVersion = v }++-- | Cancel a process with the specified identifier.+cancelProcess :: ProcessId -> Event ()+cancelProcess pid =+  Event $ \p ->+  do z <- readIORef (processCancelRef pid) +     unless z $+       do writeIORef (processCancelRef pid) True+          writeIORef (processCancelToken pid) True++-- | Test whether the process with the specified identifier was canceled.+processCanceled :: ProcessId -> Event Bool+processCanceled pid =+  Event $ \p ->+  readIORef (processCancelRef pid)++instance Eq ProcessId where+  x == y = processReactCont x == processReactCont y    -- for the references are unique++instance Monad Process where+  return  = returnP+  m >>= k = bindP m k++instance Functor Process where+  fmap = liftM++instance SimulationLift Process where+  liftSimulation = liftSP+  +instance DynamicsLift Process where+  liftDynamics = liftDP+  +instance EventLift Process where+  liftEvent = liftEP+  +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++liftSP :: Simulation a -> Process a+{-# INLINE liftSP #-}+liftSP m = Process $ \pid -> liftSimulation m++liftDP :: Dynamics a -> Process a+{-# INLINE liftDP #-}+liftDP m = Process $ \pid -> liftDynamics m++liftEP :: Event a -> Process a+{-# INLINE liftEP #-}+liftEP m = Process $ \pid -> liftEvent m++liftIOP :: IO a -> Process a+{-# INLINE liftIOP #-}+liftIOP m = Process $ \pid -> liftIO m++-- | Exception handling within 'Process' computations.+catchProcess :: Process a -> (IOException -> Process a) -> Process a+catchProcess (Process m) h =+  Process $ \pid ->+  catchCont (m pid) $ \e ->+  let Process m' = h e in m' pid+                           +-- | A computation with finalization part.+finallyProcess :: Process a -> Process b -> Process a+finallyProcess (Process m) (Process m') =+  Process $ \pid ->+  finallyCont (m pid) (m' pid)++-- | Throw the exception with the further exception handling.+-- By some reasons, the standard 'throw' function per se is not handled +-- properly within 'Process' computations, although it will be still +-- handled if it will be hidden under the 'liftIO' function. The problem +-- arises namely with the @throw@ function, not 'IO' computations.+throwProcess :: IOException -> Process a+throwProcess = liftIO . throw+

@@ -0,0 +1,260 @@++-- |+-- Module     : Simulation.Aivika.Internal.Signal+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com>+-- License    : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability  : experimental+-- Tested with: GHC 7.6.3+--+-- This module defines the signal which we can subscribe handlers to. +-- These handlers can be disposed. The signal is triggered in the +-- current time point actuating the corresponded computations from +-- the handlers. +--++module Simulation.Aivika.Internal.Signal+       (Signal(..),+        SignalSource(..),+        newSignalSource,+        handleSignal_,+        mapSignal,+        mapSignalM,+        apSignal,+        filterSignal,+        filterSignalM,+        emptySignal,+        merge2Signals,+        merge3Signals,+        merge4Signals,+        merge5Signals) where++import Data.IORef+import Data.Monoid++import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Internal.Specs+import Simulation.Aivika.Internal.Simulation+import Simulation.Aivika.Internal.Event++-- | The signal source that can publish its signal.+data SignalSource a =+  SignalSource { publishSignal :: Signal a,+                                  -- ^ Publish the signal.+                 triggerSignal :: a -> Event ()+                                  -- ^ Trigger the signal actuating +                                  -- all its handlers at the current +                                  -- simulation time point.+               }+  +-- | The signal that can have disposable handlers.  +data Signal a =+  Signal { handleSignal :: (a -> Event ()) -> Event (Event ())+           -- ^ Subscribe the handler to the specified +           -- signal and return a nested computation +           -- that, being applied, unsubscribes the +           -- handler from this signal.+         }+  +-- | The queue of signal handlers.+data SignalHandlerQueue a =+  SignalHandlerQueue { queueStart :: IORef (Maybe (SignalHandler a)),+                       queueEnd   :: IORef (Maybe (SignalHandler a)) }+  +-- | It contains the information about the disposable queue handler.+data SignalHandler a =+  SignalHandler { handlerComp :: a -> Event (),+                  handlerPrev :: IORef (Maybe (SignalHandler a)),+                  handlerNext :: IORef (Maybe (SignalHandler a)) }++-- | Subscribe the handler to the specified signal.+-- To subscribe the disposable handlers, use function 'handleSignal'.+handleSignal_ :: Signal a -> (a -> Event ()) -> Event ()+handleSignal_ signal h = +  do x <- handleSignal signal h+     return ()+     +-- | Create a new signal source.+newSignalSource :: Simulation (SignalSource a)+newSignalSource =+  Simulation $ \r ->+  do start <- newIORef Nothing+     end <- newIORef Nothing+     let queue  = SignalHandlerQueue { queueStart = start,+                                       queueEnd   = end }+         signal = Signal { handleSignal = handle }+         source = SignalSource { publishSignal = signal, +                                 triggerSignal = trigger }+         handle h =+           Event $ \p ->+           do x <- enqueueSignalHandler queue h+              return $+                Event $ \p -> dequeueSignalHandler queue x+         trigger a =+           Event $ \p ->+           let h = queueStart queue+           in triggerSignalHandlers h a p+     return source++-- | Trigger all next signal handlers.+triggerSignalHandlers :: IORef (Maybe (SignalHandler a)) -> a -> Point -> IO ()+{-# INLINE triggerSignalHandlers #-}+triggerSignalHandlers r a p =+  do x <- readIORef r+     case x of+       Nothing -> return ()+       Just h ->+         do invokeEvent p $ handlerComp h a+            triggerSignalHandlers (handlerNext h) a p+            +-- | Enqueue the handler and return its representative in the queue.            +enqueueSignalHandler :: SignalHandlerQueue a -> (a -> Event ()) -> IO (SignalHandler a)+enqueueSignalHandler q h = +  do tail <- readIORef (queueEnd q)+     case tail of+       Nothing ->+         do prev <- newIORef Nothing+            next <- newIORef Nothing+            let handler = SignalHandler { handlerComp = h,+                                          handlerPrev = prev,+                                          handlerNext = next }+            writeIORef (queueStart q) (Just handler)+            writeIORef (queueEnd q) (Just handler)+            return handler+       Just x ->+         do prev <- newIORef tail+            next <- newIORef Nothing+            let handler = SignalHandler { handlerComp = h,+                                          handlerPrev = prev,+                                          handlerNext = next }+            writeIORef (handlerNext x) (Just handler)+            writeIORef (queueEnd q) (Just handler)+            return handler++-- | Dequeue the handler representative.+dequeueSignalHandler :: SignalHandlerQueue a -> SignalHandler a -> IO ()+dequeueSignalHandler q h = +  do prev <- readIORef (handlerPrev h)+     case prev of+       Nothing ->+         do next <- readIORef (handlerNext h)+            case next of+              Nothing ->+                do writeIORef (queueStart q) Nothing+                   writeIORef (queueEnd q) Nothing+              Just y ->+                do writeIORef (handlerPrev y) Nothing+                   writeIORef (handlerNext h) Nothing+                   writeIORef (queueStart q) next+       Just x ->+         do next <- readIORef (handlerNext h)+            case next of+              Nothing ->+                do writeIORef (handlerPrev h) Nothing+                   writeIORef (handlerNext x) Nothing+                   writeIORef (queueEnd q) prev+              Just y ->+                do writeIORef (handlerPrev h) Nothing+                   writeIORef (handlerNext h) Nothing+                   writeIORef (handlerPrev y) prev+                   writeIORef (handlerNext x) next++instance Functor Signal where+  fmap = mapSignal+  +instance Monoid (Signal a) where +  +  mempty = emptySignal+  +  mappend = merge2Signals+  +  mconcat [] = emptySignal+  mconcat [x1] = x1+  mconcat [x1, x2] = merge2Signals x1 x2+  mconcat [x1, x2, x3] = merge3Signals x1 x2 x3+  mconcat [x1, x2, x3, x4] = merge4Signals x1 x2 x3 x4+  mconcat [x1, x2, x3, x4, x5] = merge5Signals x1 x2 x3 x4 x5+  mconcat (x1 : x2 : x3 : x4 : x5 : xs) = +    mconcat $ merge5Signals x1 x2 x3 x4 x5 : xs+  +-- | Map the signal according the specified function.+mapSignal :: (a -> b) -> Signal a -> Signal b+mapSignal f m =+  Signal { handleSignal = \h -> +            handleSignal m $ h . f }++-- | Filter only those signal values that satisfy to +-- the specified predicate.+filterSignal :: (a -> Bool) -> Signal a -> Signal a+filterSignal p m =+  Signal { handleSignal = \h ->+            handleSignal m $ \a ->+            when (p a) $ h a }+  +-- | Filter only those signal values that satisfy to +-- the specified predicate.+filterSignalM :: (a -> Event Bool) -> Signal a -> Signal a+filterSignalM p m =+  Signal { handleSignal = \h ->+            handleSignal m $ \a ->+            do x <- p a+               when x $ h a }+  +-- | Merge two signals.+merge2Signals :: Signal a -> Signal a -> Signal a+merge2Signals m1 m2 =+  Signal { handleSignal = \h ->+            do x1 <- handleSignal m1 h+               x2 <- handleSignal m2 h+               return $ do { x1; x2 } }++-- | Merge three signals.+merge3Signals :: Signal a -> Signal a -> Signal a -> Signal a+merge3Signals m1 m2 m3 =+  Signal { handleSignal = \h ->+            do x1 <- handleSignal m1 h+               x2 <- handleSignal m2 h+               x3 <- handleSignal m3 h+               return $ do { x1; x2; x3 } }++-- | Merge four signals.+merge4Signals :: Signal a -> Signal a -> Signal a -> +                 Signal a -> Signal a+merge4Signals m1 m2 m3 m4 =+  Signal { handleSignal = \h ->+            do x1 <- handleSignal m1 h+               x2 <- handleSignal m2 h+               x3 <- handleSignal m3 h+               x4 <- handleSignal m4 h+               return $ do { x1; x2; x3; x4 } }+           +-- | Merge five signals.+merge5Signals :: Signal a -> Signal a -> Signal a -> +                 Signal a -> Signal a -> Signal a+merge5Signals m1 m2 m3 m4 m5 =+  Signal { handleSignal = \h ->+            do x1 <- handleSignal m1 h+               x2 <- handleSignal m2 h+               x3 <- handleSignal m3 h+               x4 <- handleSignal m4 h+               x5 <- handleSignal m5 h+               return $ do { x1; x2; x3; x4; x5 } }++-- | Compose the signal.+mapSignalM :: (a -> Event b) -> Signal a -> Signal b+mapSignalM f m =+  Signal { handleSignal = \h ->+            handleSignal m (f >=> h) }+  +-- | Transform the signal.+apSignal :: Event (a -> b) -> Signal a -> Signal b+apSignal f m =+  Signal { handleSignal = \h ->+            handleSignal m $ \a -> do { x <- f; h (x a) } }++-- | An empty signal which is never triggered.+emptySignal :: Signal a+emptySignal =+  Signal { handleSignal = \h -> return $ return () }