aivika-branches (empty) → 0.1
raw patch · 12 files changed
+1060/−0 lines, 12 filesdep +aivikadep +aivika-transformersdep +basesetup-changed
Dependencies added: aivika, aivika-transformers, base, containers, mtl, random
Files
- LICENSE +30/−0
- Setup.lhs +3/−0
- Simulation/Aivika/Branch.hs +24/−0
- Simulation/Aivika/Branch/Br.hs +33/−0
- Simulation/Aivika/Branch/Event.hs +211/−0
- Simulation/Aivika/Branch/Generator.hs +210/−0
- Simulation/Aivika/Branch/Internal/Br.hs +129/−0
- Simulation/Aivika/Branch/Internal/Ref.hs +147/−0
- Simulation/Aivika/Branch/QueueStrategy.hs +68/−0
- Simulation/Aivika/Branch/Ref/Base.hs +50/−0
- aivika-branches.cabal +51/−0
- tests/MachRep1.hs +104/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2016 David Sorokin <david.sorokin@gmail.com>++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ Simulation/Aivika/Branch.hs view
@@ -0,0 +1,24 @@++-- |+-- Module : Simulation.Aivika.Branch+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- This module re-exports the library functionality related to branching computations.+--+module Simulation.Aivika.Branch+ (-- * Modules+ module Simulation.Aivika.Branch.Br,+ module Simulation.Aivika.Branch.Event,+ module Simulation.Aivika.Branch.Generator,+ module Simulation.Aivika.Branch.QueueStrategy,+ module Simulation.Aivika.Branch.Ref.Base) where++import Simulation.Aivika.Branch.Br+import Simulation.Aivika.Branch.Event+import Simulation.Aivika.Branch.Generator+import Simulation.Aivika.Branch.QueueStrategy+import Simulation.Aivika.Branch.Ref.Base
+ Simulation/Aivika/Branch/Br.hs view
@@ -0,0 +1,33 @@++-- |+-- Module : Simulation.Aivika.Branch.Br+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- This module defines 'BrIO' as an instance of the 'MonadDES' type class.+--+module Simulation.Aivika.Branch.Br+ (BrIO,+ runBr,+ branchLevel) where++import Simulation.Aivika.Trans.Comp+import Simulation.Aivika.Trans.DES+import Simulation.Aivika.Trans.Exception+import Simulation.Aivika.Trans.Generator+import Simulation.Aivika.Trans.Event+import Simulation.Aivika.Trans.Ref.Base+import Simulation.Aivika.Trans.QueueStrategy++import Simulation.Aivika.Branch.Internal.Br+import Simulation.Aivika.Branch.Event+import Simulation.Aivika.Branch.Generator+import Simulation.Aivika.Branch.Ref.Base+import Simulation.Aivika.Branch.QueueStrategy++instance MonadDES BrIO++instance MonadComp BrIO
+ Simulation/Aivika/Branch/Event.hs view
@@ -0,0 +1,211 @@++{-# LANGUAGE TypeFamilies #-}++-- |+-- Module : Simulation.Aivika.Branch.Event+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- The module defines an event queue, where 'BrIO' is an instance of 'EventQueueing'.+-- Also it defines basic functions for branching computations.+--+module Simulation.Aivika.Branch.Event+ (branchEvent,+ futureEvent,+ futureEventWith) where++import Data.IORef++import Control.Monad+import Control.Monad.Trans++import qualified Simulation.Aivika.PriorityQueue.Pure as PQ++import Simulation.Aivika.Trans+import Simulation.Aivika.Trans.Internal.Types++import Simulation.Aivika.Branch.Internal.Br++-- | An implementation of the 'EventQueueing' type class.+instance EventQueueing BrIO where++ -- | The event queue type.+ data EventQueue BrIO =+ EventQueue { queuePQ :: IORef (PQ.PriorityQueue (Point BrIO -> BrIO ())),+ -- ^ the underlying priority queue+ queueBusy :: IORef Bool,+ -- ^ whether the queue is currently processing events+ queueTime :: IORef Double+ -- ^ the actual time of the event queue+ }++ newEventQueue specs =+ do f <- liftIO $ newIORef False+ t <- liftIO $ newIORef (spcStartTime specs)+ pq <- liftIO $ newIORef PQ.emptyQueue+ return EventQueue { queuePQ = pq,+ queueBusy = f,+ queueTime = t }++ enqueueEvent t (Event m) =+ Event $ \p ->+ Br $ \ps ->+ let pq = queuePQ $ runEventQueue $ pointRun p+ in modifyIORef pq $ \x -> PQ.enqueue x t m++ runEventWith processing (Event e) =+ Dynamics $ \p ->+ do invokeDynamics p $ processEvents processing+ e p++ eventQueueCount =+ Event $ \p ->+ Br $ \ps ->+ let pq = queuePQ $ runEventQueue $ pointRun p+ in fmap PQ.queueCount $ readIORef pq++-- | Process the pending events.+processPendingEventsCore :: Bool -> Dynamics BrIO ()+processPendingEventsCore includingCurrentEvents = Dynamics r where+ r p =+ Br $ \ps ->+ do let q = runEventQueue $ pointRun p+ f = queueBusy q+ f' <- readIORef f+ unless f' $+ do writeIORef f True+ call q p ps+ writeIORef f False+ call q p ps =+ do let pq = queuePQ q+ r = pointRun p+ f <- fmap PQ.queueNull $ readIORef pq+ unless f $+ do (t2, c2) <- fmap PQ.queueFront $ readIORef 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+ modifyIORef pq PQ.dequeue+ let sc = pointSpecs p+ t0 = spcStartTime sc+ dt = spcDT sc+ n2 = fromIntegral $ floor ((t2 - t0) / dt)+ invokeBr ps $+ c2 $ p { pointTime = t2,+ pointIteration = n2,+ pointPhase = -1 }+ call q p ps++-- | Process the pending events synchronously, i.e. without past.+processPendingEvents :: Bool -> Dynamics BrIO ()+processPendingEvents includingCurrentEvents = Dynamics r where+ r p =+ Br $ \ps ->+ 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 invokeBr ps $+ invokeDynamics p $+ processPendingEventsCore includingCurrentEvents++-- | A memoized value.+processEventsIncludingCurrent :: Dynamics BrIO ()+processEventsIncludingCurrent = processPendingEvents True++-- | A memoized value.+processEventsIncludingEarlier :: Dynamics BrIO ()+processEventsIncludingEarlier = processPendingEvents False++-- | A memoized value.+processEventsIncludingCurrentCore :: Dynamics BrIO ()+processEventsIncludingCurrentCore = processPendingEventsCore True++-- | A memoized value.+processEventsIncludingEarlierCore :: Dynamics BrIO ()+processEventsIncludingEarlierCore = processPendingEventsCore True++-- | Process the events.+processEvents :: EventProcessing -> Dynamics BrIO ()+processEvents CurrentEvents = processEventsIncludingCurrent+processEvents EarlierEvents = processEventsIncludingEarlier+processEvents CurrentEventsOrFromPast = processEventsIncludingCurrentCore+processEvents EarlierEventsOrFromPast = processEventsIncludingEarlierCore++-- | Branch a new computation and return its result leaving the current computation intact.+--+-- A new derivative branch with 'branchLevel' increased by 1 is created at the current modeling time.+-- Then the result of the specified computation for the derivative branch is returned.+--+-- The state of the current computation including its event queue and mutable references 'Ref'+-- remain intact. In some sense we copy the state of the model to the derivative branch and then+-- proceed with the derived simulation. The copying operation is relatively cheap.+branchEvent :: Event BrIO a -> Event BrIO a+branchEvent (Event m) =+ Event $ \p ->+ Br $ \ps->+ do p2 <- clonePoint p+ ps2 <- newBrParams ps+ invokeBr ps2 (m p2)++-- | Branch a new computation and return its result at the desired time+-- in the future leaving the current computation intact.+--+-- A new derivative branch with 'branchLevel' increased by 1 is created at the current modeling time.+-- All pending events are processed till the specified time for that new branch. Then the result+-- of the specified computation for the derivative branch is returned.+--+-- The state of the current computation including its event queue and mutable references 'Ref'+-- remain intact. In some sense we copy the state of the model to the derivative branch and then+-- proceed with the derived simulation. The copying operation is relatively cheap.+futureEvent :: Double -> Event BrIO a -> Event BrIO a+futureEvent = futureEventWith CurrentEvents++-- | Like 'futureEvent' but allows specifying how the pending events must be processed.+futureEventWith :: EventProcessing -> Double -> Event BrIO a -> Event BrIO a+futureEventWith processing t (Event m) =+ Event $ \p ->+ Br $ \ps ->+ do when (t < pointTime p) $+ error "The specified time is less than the current modeling time: futureEventWith"+ p2 <- clonePoint p+ ps2 <- newBrParams ps+ let sc = pointSpecs p+ t0 = spcStartTime sc+ t' = spcStopTime sc+ dt = spcDT sc+ n = fromIntegral $ floor ((t - t0) / dt)+ p' = p2 { pointTime = t,+ pointIteration = n,+ pointPhase = -1 }+ invokeBr ps2 $+ invokeDynamics p' $+ processEvents processing+ invokeBr ps2 (m p')++-- | Clone the time point.+clonePoint :: Point BrIO -> IO (Point BrIO)+clonePoint p =+ do let r = pointRun p+ q = runEventQueue r+ pq <- readIORef (queuePQ q)+ t <- readIORef (queueTime q)+ pq2 <- newIORef pq+ f2 <- newIORef False+ t2 <- newIORef t+ let q2 = EventQueue { queuePQ = pq2,+ queueBusy = f2,+ queueTime = t2 }+ r2 = r { runEventQueue = q2 }+ p2 = p { pointRun = r2 }+ return p2
+ Simulation/Aivika/Branch/Generator.hs view
@@ -0,0 +1,210 @@++{-# LANGUAGE TypeFamilies #-}++-- |+-- Module : Simulation.Aivika.Branch.Generator+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- Here is defined a random number generator,+-- where 'BrIO' is an instance of 'MonadGenerator'.+--+module Simulation.Aivika.Branch.Generator () where++import Control.Monad+import Control.Monad.Trans++import System.Random++import Data.IORef++import Simulation.Aivika.Trans+import Simulation.Aivika.Branch.Internal.Br++instance MonadGenerator BrIO where++ data Generator BrIO =+ Generator { generator01 :: BrIO Double,+ -- ^ the generator of uniform numbers from 0 to 1+ generatorNormal01 :: BrIO Double+ -- ^ the generator of normal numbers with mean 0 and variance 1+ }++ generateUniform = generateUniform01 . generator01++ generateUniformInt = generateUniformInt01 . generator01++ generateNormal = generateNormal01 . generatorNormal01++ generateExponential = generateExponential01 . generator01++ generateErlang = generateErlang01 . generator01++ generatePoisson = generatePoisson01 . generator01++ generateBinomial = generateBinomial01 . generator01++ newGenerator tp =+ case tp of+ SimpleGenerator ->+ liftIO newStdGen >>= newRandomGenerator+ SimpleGeneratorWithSeed x ->+ error "Unsupported generator type SimpleGeneratorWithSeed: newGenerator"+ CustomGenerator g ->+ g+ CustomGenerator01 g ->+ newRandomGenerator01 g++ newRandomGenerator g = + do r <- liftIO $ newIORef g+ let g01 = do g <- liftIO $ readIORef r+ let (x, g') = random g+ liftIO $ writeIORef r g'+ return x+ newRandomGenerator01 g01++ newRandomGenerator01 g01 =+ do gNormal01 <- newNormalGenerator01 g01+ return Generator { generator01 = g01,+ generatorNormal01 = gNormal01 }++-- | Generate an uniform random number with the specified minimum and maximum.+generateUniform01 :: BrIO Double+ -- ^ the generator+ -> Double+ -- ^ minimum+ -> Double+ -- ^ maximum+ -> BrIO Double+generateUniform01 g min max =+ do x <- g+ return $ min + x * (max - min)++-- | Generate an uniform random number with the specified minimum and maximum.+generateUniformInt01 :: BrIO Double+ -- ^ the generator+ -> Int+ -- ^ minimum+ -> Int+ -- ^ maximum+ -> BrIO Int+generateUniformInt01 g min max =+ do x <- g+ let min' = fromIntegral min+ max' = fromIntegral max+ return $ round (min' + x * (max' - min'))++-- | Generate a normal random number by the specified generator, mean and variance.+generateNormal01 :: BrIO Double+ -- ^ normal random numbers with mean 0 and variance 1+ -> Double+ -- ^ mean+ -> Double+ -- ^ variance+ -> BrIO Double+generateNormal01 g mu nu =+ do x <- g+ return $ mu + nu * x++-- | Create a normal random number generator with mean 0 and variance 1+-- by the specified generator of uniform random numbers from 0 to 1.+newNormalGenerator01 :: BrIO Double+ -- ^ the generator+ -> BrIO (BrIO Double)+newNormalGenerator01 g =+ do nextRef <- liftIO $ newIORef 0.0+ flagRef <- liftIO $ newIORef False+ xi1Ref <- liftIO $ newIORef 0.0+ xi2Ref <- liftIO $ newIORef 0.0+ psiRef <- liftIO $ newIORef 0.0+ let loop =+ do psi <- liftIO $ readIORef psiRef+ if (psi >= 1.0) || (psi == 0.0)+ then do g1 <- g+ g2 <- g+ let xi1 = 2.0 * g1 - 1.0+ xi2 = 2.0 * g2 - 1.0+ psi = xi1 * xi1 + xi2 * xi2+ liftIO $ writeIORef xi1Ref xi1+ liftIO $ writeIORef xi2Ref xi2+ liftIO $ writeIORef psiRef psi+ loop+ else liftIO $ writeIORef psiRef $ sqrt (- 2.0 * log psi / psi)+ return $+ do flag <- liftIO $ readIORef flagRef+ if flag+ then do liftIO $ writeIORef flagRef False+ liftIO $ readIORef nextRef+ else do liftIO $ writeIORef xi1Ref 0.0+ liftIO $ writeIORef xi2Ref 0.0+ liftIO $ writeIORef psiRef 0.0+ loop+ xi1 <- liftIO $ readIORef xi1Ref+ xi2 <- liftIO $ readIORef xi2Ref+ psi <- liftIO $ readIORef psiRef+ liftIO $ writeIORef flagRef True+ liftIO $ writeIORef nextRef $ xi2 * psi+ return $ xi1 * psi++-- | Return the exponential random number with the specified mean.+generateExponential01 :: BrIO Double+ -- ^ the generator+ -> Double+ -- ^ the mean+ -> BrIO Double+generateExponential01 g mu =+ do x <- g+ return (- log x * mu)++-- | Return the Erlang random number.+generateErlang01 :: BrIO Double+ -- ^ the generator+ -> Double+ -- ^ the scale+ -> Int+ -- ^ the shape+ -> BrIO Double+generateErlang01 g beta m =+ do x <- loop m 1+ return (- log x * beta)+ where loop m acc+ | m < 0 = error "Negative shape: generateErlang."+ | m == 0 = return acc+ | otherwise = do x <- g+ loop (m - 1) (x * acc)++-- | Generate the Poisson random number with the specified mean.+generatePoisson01 :: BrIO Double+ -- ^ the generator+ -> Double+ -- ^ the mean+ -> BrIO Int+generatePoisson01 g mu =+ do prob0 <- g+ let loop prob prod acc+ | prob <= prod = return acc+ | otherwise = loop+ (prob - prod)+ (prod * mu / fromIntegral (acc + 1))+ (acc + 1)+ loop prob0 (exp (- mu)) 0++-- | Generate a binomial random number with the specified probability and number of trials. +generateBinomial01 :: BrIO Double+ -- ^ the generator+ -> Double + -- ^ the probability+ -> Int+ -- ^ the number of trials+ -> BrIO Int+generateBinomial01 g prob trials = loop trials 0 where+ loop n acc+ | n < 0 = error "Negative number of trials: generateBinomial."+ | n == 0 = return acc+ | otherwise = do x <- g+ if x <= prob+ then loop (n - 1) (acc + 1)+ else loop (n - 1) acc
+ Simulation/Aivika/Branch/Internal/Br.hs view
@@ -0,0 +1,129 @@++-- |+-- Module : Simulation.Aivika.Branch.Internal.Br+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- This module defines a branching computation.+--+module Simulation.Aivika.Branch.Internal.Br+ (BrParams(..),+ BrIO(..),+ invokeBr,+ runBr,+ newBrParams,+ newRootBrParams,+ branchLevel) where++import Data.IORef+import Data.Maybe++import Control.Applicative+import Control.Monad+import Control.Monad.Trans+import Control.Exception (throw, catch, finally)++import Simulation.Aivika.Trans.Exception++-- | The branching computation.+newtype BrIO a = Br { unBr :: BrParams -> IO a+ -- ^ Unwrap the computation.+ }++-- | The parameters of the computation.+data BrParams =+ BrParams { brId :: !Int,+ -- ^ The branch identifier.+ brIdGenerator :: IORef Int,+ -- ^ The generator of identifiers.+ brLevel :: !Int,+ -- ^ The branch level.+ brParent :: Maybe BrParams,+ -- ^ The branch parent.+ brUniqueRef :: IORef ()+ -- ^ The unique reference to which+ -- the finalizers are attached to+ -- be garbage collected.+ }++instance Monad BrIO where++ {-# INLINE return #-}+ return = Br . const . return++ {-# INLINE (>>=) #-}+ (Br m) >>= k = Br $ \ps ->+ m ps >>= \a ->+ let m' = unBr (k a) in m' ps++instance Applicative BrIO where++ {-# INLINE pure #-}+ pure = return++ {-# INLINE (<*>) #-}+ (<*>) = ap++instance Functor BrIO where++ {-# INLINE fmap #-}+ fmap f (Br m) = Br $ fmap f . m ++instance MonadIO BrIO where++ {-# INLINE liftIO #-}+ liftIO = Br . const . liftIO++instance MonadException BrIO where++ catchComp (Br m) h = Br $ \ps ->+ catch (m ps) (\e -> unBr (h e) ps)++ finallyComp (Br m1) (Br m2) = Br $ \ps ->+ finally (m1 ps) (m2 ps)+ + throwComp e = Br $ \ps ->+ throw e++-- | Invoke the computation.+invokeBr :: BrParams -> BrIO a -> IO a+{-# INLINE invokeBr #-}+invokeBr ps (Br m) = m ps++-- | Run the branching computation.+runBr :: BrIO a -> IO a+runBr m =+ do ps <- newRootBrParams+ unBr m ps++-- | Create a new child branch.+newBrParams :: BrParams -> IO BrParams+newBrParams ps =+ do id <- atomicModifyIORef (brIdGenerator ps) $ \a ->+ let b = a + 1 in b `seq` (b, b)+ let level = 1 + brLevel ps+ uniqueRef <- newIORef ()+ return BrParams { brId = id,+ brIdGenerator = brIdGenerator ps,+ brLevel = level `seq` level,+ brParent = Just ps,+ brUniqueRef = uniqueRef }++-- | Create a root branch.+newRootBrParams :: IO BrParams+newRootBrParams =+ do genId <- newIORef 0+ uniqueRef <- newIORef ()+ return BrParams { brId = 0,+ brIdGenerator = genId,+ brLevel = 0,+ brParent = Nothing,+ brUniqueRef = uniqueRef+ }++-- | Return the current branch level starting from 0.+branchLevel :: BrIO Int+branchLevel = Br $ \ps -> return (brLevel ps)
+ Simulation/Aivika/Branch/Internal/Ref.hs view
@@ -0,0 +1,147 @@++{-# LANGUAGE BangPatterns #-}++-- |+-- Module : Simulation.Aivika.Branch.Internal.Ref+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- The implementation of mutable references.+--+module Simulation.Aivika.Branch.Internal.Ref+ (Ref,+ newEmptyRef,+ newEmptyRef0,+ newRef,+ newRef0,+ readRef,+ writeRef,+ modifyRef) where++-- import Debug.Trace++import Data.IORef+import qualified Data.IntMap as M++import System.Mem.Weak++import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Trans.Internal.Types++import Simulation.Aivika.Branch.Internal.Br++-- | A reference map.+type RefMap a = IORef (M.IntMap (IORef a))++-- | A mutable reference.+data Ref a = Ref { refMap :: RefMap a,+ -- ^ the map of actual references+ refWeakMap :: Weak (RefMap a)+ -- ^ a weak reference to the map itself+ }++instance Eq (Ref a) where+ r1 == r2 = (refMap r1) == (refMap r2)++-- | Create an empty reference.+newEmptyRef :: Simulation BrIO (Ref a)+newEmptyRef = Simulation $ const newEmptyRef0++-- | Create an empty reference.+newEmptyRef0 :: BrIO (Ref a)+newEmptyRef0 =+ Br $ \ps ->+ do rm <- newIORef M.empty+ wm <- mkWeakIORef rm $+ -- trace ("fin newEmptyRef0: " ++ show (brId ps)) $+ return ()+ return Ref { refMap = rm,+ refWeakMap = wm }++-- | Create a new reference.+newRef :: a -> Simulation BrIO (Ref a)+newRef = Simulation . const . newRef0++-- | Create a new reference.+newRef0 :: a -> BrIO (Ref a)+newRef0 a =+ Br $ \ps ->+ do r <- invokeBr ps newEmptyRef0+ ra <- newIORef a+ let !i = brId ps+ !wm = refWeakMap r+ -- mkWeakIORef (brUniqueRef ps) (trace ("fin newIORef0: " ++ show i) $ finalizeRef wm i)+ mkWeakIORef (brUniqueRef ps) (finalizeRef wm i)+ writeIORef (refMap r) $+ M.insert i ra M.empty+ return r+ +-- | Read the value of a reference.+readRef :: Ref a -> Event BrIO a+readRef r =+ Event $ \p ->+ Br $ \ps ->+ do m <- readIORef (refMap r)+ let loop ps =+ case M.lookup (brId ps) m of+ Just ra -> readIORef ra+ Nothing ->+ case brParent ps of+ Just ps' -> loop ps'+ Nothing -> error "Cannot find branch: readRef"+ loop ps++-- | Write a new value into the reference.+writeRef :: Ref a -> a -> Event BrIO ()+writeRef r a =+ Event $ \p ->+ Br $ \ps ->+ do m <- readIORef (refMap r)+ let !i = brId ps+ case M.lookup i m of+ Just ra -> a `seq` writeIORef ra a+ Nothing ->+ do ra <- a `seq` newIORef a+ let !wm = refWeakMap r+ -- mkWeakIORef (brUniqueRef ps) (trace ("fin writeRef: " ++ show i) $ finalizeRef wm i)+ mkWeakIORef (brUniqueRef ps) (finalizeRef wm i)+ atomicModifyIORef (refMap r) $ \m ->+ let m' = M.insert i ra m in (m', ())++-- | Mutate the contents of the reference.+modifyRef :: Ref a -> (a -> a) -> Event BrIO ()+modifyRef r f =+ Event $ \p ->+ Br $ \ps ->+ do m <- readIORef (refMap r)+ let !i = brId ps+ case M.lookup i m of+ Just ra ->+ do a <- readIORef ra+ let b = f a+ b `seq` writeIORef ra b+ Nothing ->+ do a <- invokeBr ps $ invokeEvent p $ readRef r+ invokeBr ps $ invokeEvent p $ writeRef r (f a)++-- | Finalize the reference cell by the specified branch identifier.+finalizeRef :: Weak (RefMap a) -> Int -> IO ()+finalizeRef wm i =+ do rm <- deRefWeak wm+ -- trace ("finalizeRef: " ++ show i) $ return ()+ case rm of+ Nothing ->+ return ()+ Just rm ->+ do m <- readIORef rm+ case M.lookup i m of+ Just ra ->+ atomicModifyIORef rm $ \m ->+ let m' = M.delete i m in (m', ())+ Nothing ->+ return ()
+ Simulation/Aivika/Branch/QueueStrategy.hs view
@@ -0,0 +1,68 @@++{-# LANGUAGE TypeFamilies, MultiParamTypeClasses #-}++-- |+-- Module : Simulation.Aivika.Branch.QueueStrategy+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- This module defines queue strategies 'FCFS' and 'LCFS' for the 'BrIO' computation.+--+module Simulation.Aivika.Branch.QueueStrategy () where++import Control.Monad.Trans++import Simulation.Aivika.Trans+import qualified Simulation.Aivika.Trans.DoubleLinkedList as LL++import Simulation.Aivika.Branch.Internal.Br+import Simulation.Aivika.Branch.Ref.Base++-- | An implementation of the 'FCFS' queue strategy.+instance QueueStrategy BrIO FCFS where++ -- | A queue used by the 'FCFS' strategy.+ newtype StrategyQueue BrIO FCFS a = FCFSQueue (LL.DoubleLinkedList BrIO a)++ newStrategyQueue s = fmap FCFSQueue LL.newList++ strategyQueueNull (FCFSQueue q) = LL.listNull q++-- | An implementation of the 'FCFS' queue strategy.+instance DequeueStrategy BrIO FCFS where++ strategyDequeue (FCFSQueue q) =+ do i <- LL.listFirst q+ LL.listRemoveFirst q+ return i++-- | An implementation of the 'FCFS' queue strategy.+instance EnqueueStrategy BrIO FCFS where++ strategyEnqueue (FCFSQueue q) i = LL.listAddLast q i++-- | An implementation of the 'LCFS' queue strategy.+instance QueueStrategy BrIO LCFS where++ -- | A queue used by the 'LCFS' strategy.+ newtype StrategyQueue BrIO LCFS a = LCFSQueue (LL.DoubleLinkedList BrIO a)++ newStrategyQueue s = fmap LCFSQueue LL.newList+ + strategyQueueNull (LCFSQueue q) = LL.listNull q++-- | An implementation of the 'LCFS' queue strategy.+instance DequeueStrategy BrIO LCFS where++ strategyDequeue (LCFSQueue q) =+ do i <- LL.listFirst q+ LL.listRemoveFirst q+ return i++-- | An implementation of the 'LCFS' queue strategy.+instance EnqueueStrategy BrIO LCFS where++ strategyEnqueue (LCFSQueue q) i = LL.listInsertFirst q i
+ Simulation/Aivika/Branch/Ref/Base.hs view
@@ -0,0 +1,50 @@++{-# LANGUAGE TypeFamilies #-}++-- |+-- Module : Simulation.Aivika.Branch.Ref.Base+-- Copyright : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>+-- License : BSD3+-- Maintainer : David Sorokin <david.sorokin@gmail.com>+-- Stability : experimental+-- Tested with: GHC 7.10.3+--+-- Here is an implementation of mutable references, where+-- 'BrIO' is an instance of 'MonadRef' and 'MonadRef0'.+--+module Simulation.Aivika.Branch.Ref.Base () where++import Simulation.Aivika.Trans.Internal.Types+import Simulation.Aivika.Trans.Comp+import Simulation.Aivika.Trans.Simulation+import Simulation.Aivika.Trans.Ref.Base++import Simulation.Aivika.Branch.Internal.Br+import qualified Simulation.Aivika.Branch.Internal.Ref as R++-- | The implementation of mutable references.+instance MonadRef BrIO where++ -- | The mutable reference.+ newtype Ref BrIO a = Ref { refValue :: R.Ref a }++ {-# INLINE newRef #-}+ newRef = fmap Ref . R.newRef ++ {-# INLINE readRef #-}+ readRef (Ref r) = R.readRef r++ {-# INLINE writeRef #-}+ writeRef (Ref r) = R.writeRef r++ {-# INLINE modifyRef #-}+ modifyRef (Ref r) = R.modifyRef r++ {-# INLINE equalRef #-}+ equalRef (Ref r1) (Ref r2) = (r1 == r2)++-- | A subtype of mutable references that can be created under more weak conditions.+instance MonadRef0 BrIO where++ {-# INLINE newRef0 #-}+ newRef0 = fmap Ref . R.newRef0
+ aivika-branches.cabal view
@@ -0,0 +1,51 @@+name: aivika-branches+version: 0.1+synopsis: Branching discrete event simulation library+description:+ This package extends the Aivika [1] library with facilities for creating branches to run+ nested simulations within simulation. For example, it can be useful for financial modeling.+ .+ \[1] <http://hackage.haskell.org/package/aivika-transformers>+ .+category: Simulation+license: BSD3+license-file: LICENSE+copyright: (c) 2016. David Sorokin <david.sorokin@gmail.com>+author: David Sorokin+maintainer: David Sorokin <david.sorokin@gmail.com>+homepage: http://www.aivikasoft.com/en/products/aivika.html+cabal-version: >= 1.6+build-type: Simple+tested-with: GHC == 7.10.3++extra-source-files: tests/MachRep1.hs++library++ exposed-modules: Simulation.Aivika.Branch+ Simulation.Aivika.Branch.Br+ Simulation.Aivika.Branch.Generator+ Simulation.Aivika.Branch.Event+ Simulation.Aivika.Branch.QueueStrategy+ Simulation.Aivika.Branch.Ref.Base++ other-modules: Simulation.Aivika.Branch.Internal.Br+ Simulation.Aivika.Branch.Internal.Ref++ build-depends: base >= 3 && < 6,+ mtl >= 2.1.1,+ containers >= 0.4.0.0,+ random >= 1.0.0.3,+ aivika >= 4.3.2,+ aivika-transformers >= 4.3.1++ extensions: TypeFamilies,+ MultiParamTypeClasses,+ BangPatterns++ ghc-options: -O2++source-repository head++ type: git+ location: https://github.com/dsorokin/aivika-branches
+ tests/MachRep1.hs view
@@ -0,0 +1,104 @@++-- It corresponds to model MachRep1 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/].+-- +-- The model description is as follows.+--+-- Two machines, which sometimes break down.+-- Up time is exponentially distributed with mean 1.0, and repair time is+-- exponentially distributed with mean 0.5. There are two repairpersons,+-- so the two machines can be repaired simultaneously if they are down+-- at the same time.+--+-- Output is long-run proportion of up time. Should get value of about+-- 0.66.++import Control.Monad+import Control.Monad.Trans++import Simulation.Aivika.Trans+import Simulation.Aivika.Branch++meanUpTime = 1.0+meanRepairTime = 0.5++specs = Specs { spcStartTime = 0.0,+ spcStopTime = 1000.0,+ spcDT = 1.0,+ spcMethod = RungeKutta4,+ spcGeneratorType = SimpleGenerator }+ +model :: Simulation BrIO (Results BrIO)+model =+ do totalUpTime <- newRef 0.0+ + let machine =+ do upTime <-+ liftParameter $+ randomExponential meanUpTime+ --+ -- r <- liftSimulation $ newRef 10+ --+ holdProcess upTime+ liftEvent $ + modifyRef totalUpTime (+ upTime)+ repairTime <-+ liftParameter $+ randomExponential meanRepairTime+ holdProcess repairTime+ machine++ runProcessInStartTime machine+ runProcessInStartTime machine+ + let maxLevel = 4++ starttime' <- liftParameter starttime+ stoptime' <- liftParameter stoptime++ let dt' = (stoptime' - starttime') / fromIntegral maxLevel+ -- dt' <- liftParameter dt+ + let forecast :: Double -> Event BrIO Double+ forecast i =+ do level <- liftComp branchLevel+ if level <= maxLevel+ then do t <- liftDynamics time+ x1 <- futureEvent (t + dt') $ forecast (i - 1)+ x2 <- futureEvent (t + dt') $ forecast (i + 1)+ let x = (x1 + x2) / 2+ x `seq` return x + else do t <- liftDynamics time+ x <- readRef totalUpTime+ return $ x / (2 * t)+ + f <- runEventInStartTime $ forecast 0+ + let upTimePropForecasted :: Event BrIO Double+ upTimePropForecasted = return f+ + let upTimeProp =+ do x <- readRef totalUpTime+ t <- liftDynamics time+ return $ x / (2 * t)++ return $+ results+ [resultSource+ "upTimeProp"+ "The long-run proportion of up time (~ 0.66)"+ upTimeProp,+ --+ resultSource+ "upTimePropForecasted"+ "The forecasted long-run proption of up time"+ upTimePropForecasted]++main :: IO ()+main =+ runBr $+ printSimulationResultsInStopTime+ printResultSourceInEnglish+ model specs