diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2016 David Sorokin + +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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/Simulation/Aivika/Branch.hs b/Simulation/Aivika/Branch.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch.hs @@ -0,0 +1,24 @@ + +-- | +-- Module : Simulation.Aivika.Branch +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 diff --git a/Simulation/Aivika/Branch/Br.hs b/Simulation/Aivika/Branch/Br.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/Br.hs @@ -0,0 +1,33 @@ + +-- | +-- Module : Simulation.Aivika.Branch.Br +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 diff --git a/Simulation/Aivika/Branch/Event.hs b/Simulation/Aivika/Branch/Event.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/Event.hs @@ -0,0 +1,211 @@ + +{-# LANGUAGE TypeFamilies #-} + +-- | +-- Module : Simulation.Aivika.Branch.Event +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 diff --git a/Simulation/Aivika/Branch/Generator.hs b/Simulation/Aivika/Branch/Generator.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/Generator.hs @@ -0,0 +1,210 @@ + +{-# LANGUAGE TypeFamilies #-} + +-- | +-- Module : Simulation.Aivika.Branch.Generator +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 diff --git a/Simulation/Aivika/Branch/Internal/Br.hs b/Simulation/Aivika/Branch/Internal/Br.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/Internal/Br.hs @@ -0,0 +1,129 @@ + +-- | +-- Module : Simulation.Aivika.Branch.Internal.Br +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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) diff --git a/Simulation/Aivika/Branch/Internal/Ref.hs b/Simulation/Aivika/Branch/Internal/Ref.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/Internal/Ref.hs @@ -0,0 +1,147 @@ + +{-# LANGUAGE BangPatterns #-} + +-- | +-- Module : Simulation.Aivika.Branch.Internal.Ref +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 () diff --git a/Simulation/Aivika/Branch/QueueStrategy.hs b/Simulation/Aivika/Branch/QueueStrategy.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/QueueStrategy.hs @@ -0,0 +1,68 @@ + +{-# LANGUAGE TypeFamilies, MultiParamTypeClasses #-} + +-- | +-- Module : Simulation.Aivika.Branch.QueueStrategy +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 diff --git a/Simulation/Aivika/Branch/Ref/Base.hs b/Simulation/Aivika/Branch/Ref/Base.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Branch/Ref/Base.hs @@ -0,0 +1,50 @@ + +{-# LANGUAGE TypeFamilies #-} + +-- | +-- Module : Simulation.Aivika.Branch.Ref.Base +-- Copyright : Copyright (c) 2016, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- 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 diff --git a/aivika-branches.cabal b/aivika-branches.cabal new file mode 100644 --- /dev/null +++ b/aivika-branches.cabal @@ -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] + . +category: Simulation +license: BSD3 +license-file: LICENSE +copyright: (c) 2016. David Sorokin +author: David Sorokin +maintainer: David Sorokin +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 diff --git a/tests/MachRep1.hs b/tests/MachRep1.hs new file mode 100644 --- /dev/null +++ b/tests/MachRep1.hs @@ -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