diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2009, 2010, 2011, 2012 David Sorokin +Copyright (c) 2009, 2010, 2011, 2012, 2013 David Sorokin All rights reserved. diff --git a/Simulation/Aivika/Dynamics.hs b/Simulation/Aivika/Dynamics.hs --- a/Simulation/Aivika/Dynamics.hs +++ b/Simulation/Aivika/Dynamics.hs @@ -1,23 +1,25 @@ -- | -- Module : Simulation.Aivika.Dynamics --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- 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 - (Dynamics, + (-- * Dynamics + Dynamics, DynamicsLift(..), runDynamicsInStartTime, runDynamicsInStopTime, runDynamicsInIntegTimes, runDynamicsInTime, runDynamicsInTimes, + -- * Error Handling catchDynamics, finallyDynamics, throwDynamics) where diff --git a/Simulation/Aivika/Dynamics/Agent.hs b/Simulation/Aivika/Dynamics/Agent.hs --- a/Simulation/Aivika/Dynamics/Agent.hs +++ b/Simulation/Aivika/Dynamics/Agent.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Agent --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module introduces an agent-based modeling. -- diff --git a/Simulation/Aivika/Dynamics/Base.hs b/Simulation/Aivika/Dynamics/Base.hs --- a/Simulation/Aivika/Dynamics/Base.hs +++ b/Simulation/Aivika/Dynamics/Base.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Base --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module defines basic functions for the 'Dynamics' monad. -- diff --git a/Simulation/Aivika/Dynamics/Buffer.hs b/Simulation/Aivika/Dynamics/Buffer.hs --- a/Simulation/Aivika/Dynamics/Buffer.hs +++ b/Simulation/Aivika/Dynamics/Buffer.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Buffer --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Cont.hs b/Simulation/Aivika/Dynamics/Cont.hs --- a/Simulation/Aivika/Dynamics/Cont.hs +++ b/Simulation/Aivika/Dynamics/Cont.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Cont --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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 diff --git a/Simulation/Aivika/Dynamics/EventQueue.hs b/Simulation/Aivika/Dynamics/EventQueue.hs --- a/Simulation/Aivika/Dynamics/EventQueue.hs +++ b/Simulation/Aivika/Dynamics/EventQueue.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.EventQueue --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- The module introduces the event queue. Any event is the Dynamics computation, -- or, saying differently, a dynamic process that has a single purpose diff --git a/Simulation/Aivika/Dynamics/FIFO.hs b/Simulation/Aivika/Dynamics/FIFO.hs --- a/Simulation/Aivika/Dynamics/FIFO.hs +++ b/Simulation/Aivika/Dynamics/FIFO.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.FIFO --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- Tested with: GHC 7.6.3 -- -- This module defines the FIFO queue. -- diff --git a/Simulation/Aivika/Dynamics/Internal/Cont.hs b/Simulation/Aivika/Dynamics/Internal/Cont.hs --- a/Simulation/Aivika/Dynamics/Internal/Cont.hs +++ b/Simulation/Aivika/Dynamics/Internal/Cont.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Cont --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Internal/Dynamics.hs b/Simulation/Aivika/Dynamics/Internal/Dynamics.hs --- a/Simulation/Aivika/Dynamics/Internal/Dynamics.hs +++ b/Simulation/Aivika/Dynamics/Internal/Dynamics.hs @@ -1,11 +1,13 @@ +{-# LANGUAGE RecursiveDo #-} + -- | -- Module : Simulation.Aivika.Dynamics.Internal.Dynamics --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- Tested with: GHC 7.6.3 -- -- The module defines the 'Dynamics' monad representing an abstract dynamic -- process, i.e. a time varying polymorphic function. @@ -38,6 +40,7 @@ import Control.Monad import Control.Monad.Trans +import Control.Monad.Fix import Simulation.Aivika.Dynamics.Internal.Simulation @@ -295,3 +298,13 @@ -- | 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 } diff --git a/Simulation/Aivika/Dynamics/Internal/Fold.hs b/Simulation/Aivika/Dynamics/Internal/Fold.hs --- a/Simulation/Aivika/Dynamics/Internal/Fold.hs +++ b/Simulation/Aivika/Dynamics/Internal/Fold.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Fold --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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. diff --git a/Simulation/Aivika/Dynamics/Internal/Interpolate.hs b/Simulation/Aivika/Dynamics/Internal/Interpolate.hs --- a/Simulation/Aivika/Dynamics/Internal/Interpolate.hs +++ b/Simulation/Aivika/Dynamics/Internal/Interpolate.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Interpolate --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module defines interpolation functions. -- These functions complement the memoization, possibly except for diff --git a/Simulation/Aivika/Dynamics/Internal/Memo.hs b/Simulation/Aivika/Dynamics/Internal/Memo.hs --- a/Simulation/Aivika/Dynamics/Internal/Memo.hs +++ b/Simulation/Aivika/Dynamics/Internal/Memo.hs @@ -3,11 +3,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Memo --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Internal/Process.hs b/Simulation/Aivika/Dynamics/Internal/Process.hs --- a/Simulation/Aivika/Dynamics/Internal/Process.hs +++ b/Simulation/Aivika/Dynamics/Internal/Process.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Process --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Internal/Signal.hs b/Simulation/Aivika/Dynamics/Internal/Signal.hs --- a/Simulation/Aivika/Dynamics/Internal/Signal.hs +++ b/Simulation/Aivika/Dynamics/Internal/Signal.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Signal --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Internal/Simulation.hs b/Simulation/Aivika/Dynamics/Internal/Simulation.hs --- a/Simulation/Aivika/Dynamics/Internal/Simulation.hs +++ b/Simulation/Aivika/Dynamics/Internal/Simulation.hs @@ -1,11 +1,13 @@ +{-# LANGUAGE RecursiveDo #-} + -- | -- Module : Simulation.Aivika.Dynamics.Internal.Simulation --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- The module defines the 'Simulation' monad that represents a simulation run. -- @@ -18,21 +20,37 @@ 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 of the Simulation monad represents a simulation run. --- --- | A value in the 'Simulation' monad represents a simulation run. +-- | 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. @@ -98,11 +116,55 @@ 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 @@ -111,3 +173,30 @@ -- | 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 } diff --git a/Simulation/Aivika/Dynamics/Internal/Time.hs b/Simulation/Aivika/Dynamics/Internal/Time.hs --- a/Simulation/Aivika/Dynamics/Internal/Time.hs +++ b/Simulation/Aivika/Dynamics/Internal/Time.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Internal.Time --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module defines the time parameters. -- diff --git a/Simulation/Aivika/Dynamics/LIFO.hs b/Simulation/Aivika/Dynamics/LIFO.hs --- a/Simulation/Aivika/Dynamics/LIFO.hs +++ b/Simulation/Aivika/Dynamics/LIFO.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.LIFO --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- Tested with: GHC 7.6.3 -- -- This module defines the LIFO queue. -- diff --git a/Simulation/Aivika/Dynamics/Parameter.hs b/Simulation/Aivika/Dynamics/Parameter.hs --- a/Simulation/Aivika/Dynamics/Parameter.hs +++ b/Simulation/Aivika/Dynamics/Parameter.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Parameter --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- Tested with: GHC 7.6.3 -- -- This module defines the parameters of simulation experiments. -- @@ -13,22 +13,26 @@ module Simulation.Aivika.Dynamics.Parameter (newParameter, newTableParameter, - newIndexedParameter) where + 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 during the simulation run, +-- | 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 during the simulation run, +-- | 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. @@ -37,7 +41,7 @@ where (i1, i2) = bounds t n = i2 - i1 + 1 --- | Create a thread-safe parameter that returns always the same value during the simulation run, +-- | 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 = @@ -57,3 +61,23 @@ 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 diff --git a/Simulation/Aivika/Dynamics/Process.hs b/Simulation/Aivika/Dynamics/Process.hs --- a/Simulation/Aivika/Dynamics/Process.hs +++ b/Simulation/Aivika/Dynamics/Process.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Process --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Random.hs b/Simulation/Aivika/Dynamics/Random.hs --- a/Simulation/Aivika/Dynamics/Random.hs +++ b/Simulation/Aivika/Dynamics/Random.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Random --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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 @@ -23,20 +23,26 @@ import Simulation.Aivika.Dynamics.Simulation import Simulation.Aivika.Dynamics.Base --- | Return the uniform random numbers between 0.0 and 1.0 in --- the integration time points. -newRandom :: Simulation (Dynamics Double) -newRandom = - memo0 $ liftIO $ getStdRandom random +-- | 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 + x <- liftIO $ getStdRandom random + min + return x * (max - min) --- | Return the normal random numbers with mean 0.0 and variance 1.0 in --- the integration time points. -newNormal :: Simulation (Dynamics Double) -newNormal = +-- | 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 - memo0 $ liftIO g + umemo0 $ do + x <- liftIO g + mu + return x * nu --- | Normal random number generator. +-- | Normal random number generator with mean 0 and variance 1. normalGen :: IO (IO Double) normalGen = do nextRef <- newIORef 0.0 diff --git a/Simulation/Aivika/Dynamics/Ref.hs b/Simulation/Aivika/Dynamics/Ref.hs --- a/Simulation/Aivika/Dynamics/Ref.hs +++ b/Simulation/Aivika/Dynamics/Ref.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Ref --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module defines an updatable reference that depends on the event queue. -- diff --git a/Simulation/Aivika/Dynamics/Resource.hs b/Simulation/Aivika/Dynamics/Resource.hs --- a/Simulation/Aivika/Dynamics/Resource.hs +++ b/Simulation/Aivika/Dynamics/Resource.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Resource --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module defines a limited resource which can be acquired and -- then released by the discontinuous process 'Process'. diff --git a/Simulation/Aivika/Dynamics/Signal.hs b/Simulation/Aivika/Dynamics/Signal.hs --- a/Simulation/Aivika/Dynamics/Signal.hs +++ b/Simulation/Aivika/Dynamics/Signal.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Signal --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- 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 diff --git a/Simulation/Aivika/Dynamics/Simulation.hs b/Simulation/Aivika/Dynamics/Simulation.hs --- a/Simulation/Aivika/Dynamics/Simulation.hs +++ b/Simulation/Aivika/Dynamics/Simulation.hs @@ -1,21 +1,27 @@ -- | -- Module : Simulation.Aivika.Dynamics.Simulation --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- The module defines the 'Simulation' monad representing a simulation run. -- module Simulation.Aivika.Dynamics.Simulation - (Simulation, + (-- * Simulation + Simulation, SimulationLift(..), Specs(..), Method(..), runSimulation, runSimulations, + -- * Error Handling + catchSimulation, + finallySimulation, + throwSimulation, + -- * Utilities simulationIndex, simulationCount, simulationSpecs) where diff --git a/Simulation/Aivika/Dynamics/SystemDynamics.hs b/Simulation/Aivika/Dynamics/SystemDynamics.hs --- a/Simulation/Aivika/Dynamics/SystemDynamics.hs +++ b/Simulation/Aivika/Dynamics/SystemDynamics.hs @@ -1,21 +1,28 @@ -{-# LANGUAGE FlexibleContexts, BangPatterns #-} +{-# LANGUAGE FlexibleContexts, BangPatterns, RecursiveDo #-} -- | -- Module : Simulation.Aivika.Dynamics.SystemDynamics --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- This module defines integrals and other functions of System Dynamics. -- module Simulation.Aivika.Dynamics.SystemDynamics - (-- * Maximum and Minimum + (-- * Equality and Ordering + (.==.), + (./=.), + (.<.), + (.>=.), + (.>.), + (.<=.), maxDynamics, minDynamics, + ifDynamics, -- * Integrals Integ, newInteg, @@ -24,15 +31,41 @@ 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) where + lookupStepwiseD, + lookupDynamics, + lookupStepwiseDynamics, + -- * Discrete Functions + delayTrans, + delay, + delayI, + udelay, + udelayI, + -- * Financial Functions + npv, + npve) where import Data.Array import Data.Array.IO.Safe @@ -45,9 +78,33 @@ import Simulation.Aivika.Dynamics.Base -- --- Maximum and Minimum +-- 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 @@ -56,10 +113,22 @@ 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), @@ -216,93 +285,241 @@ _ -> error "Incorrect phase: integRK4" --- smoothI :: Dynamics Double -> Dynamics Double -> Dynamics Double --- -> Dynamics Double --- smoothI x t i = y where --- y = integ ((x - y) / t) i - --- smooth :: Dynamics Double -> Dynamics Double -> Dynamics Double --- smooth x t = smoothI x t x +-- | 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 --- smooth3I :: Dynamics Double -> Dynamics Double -> Dynamics Double --- -> Dynamics Double --- smooth3I x t i = y where --- y = integ ((s1 - y) / t') i --- s1 = integ ((s0 - s1) / t') i --- s0 = integ ((x - s0) / t') i --- t' = t / 3.0 +-- | 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 --- smooth3 :: Dynamics Double -> Dynamics Double -> Dynamics Double --- smooth3 x t = smooth3I x t x +-- | 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 --- smoothNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double --- -> Dynamics Double --- smoothNI x t n i = s ! n where --- s = array (1, n) [(k, f k) | k <- [1 .. n]] --- f 0 = integ ((x - s ! 0) / t') i --- f k = integ ((s ! (k - 1) - s ! k) / t') i --- t' = t / fromIntegral n +-- | 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 --- smoothN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double --- smoothN x t n = smoothNI x t n x +-- | 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 --- delay1I :: Dynamics Double -> Dynamics Double -> Dynamics Double --- -> Dynamics Double --- delay1I x t i = y where --- y = integ (x - y) (i * t) / t +-- | 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 --- delay1 :: Dynamics Double -> Dynamics Double -> Dynamics Double --- delay1 x t = delay1I x t x +-- | 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 --- delay3I :: Dynamics Double -> Dynamics Double -> Dynamics Double --- -> Dynamics Double --- delay3I x t i = y where --- y = integ (s1 - y) (i * t') / t' --- s1 = integ (s0 - s1) (i * t') / t' --- s0 = integ (x - s0) (i * t') / t' --- t' = t / 3.0 +-- | 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 --- delay3 :: Dynamics Double -> Dynamics Double -> Dynamics Double --- delay3 x t = delay3I x t x +-- | 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 --- delayNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double --- -> Dynamics Double --- delayNI x t n i = s ! n where --- s = array (1, n) [(k, f k) | k <- [1 .. n]] --- f 0 = integ (x - s ! 0) (i * t') / t' --- f k = integ (s ! (k - 1) - s ! k) (i * t') / t' --- t' = t / fromIntegral n +-- | 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' --- delayN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double --- delayN x t n = delayNI x t n x +-- | 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 --- forecast :: Dynamics Double -> Dynamics Double -> Dynamics Double --- -> Dynamics Double --- forecast x at hz = --- x * (1.0 + (x / smooth x at - 1.0) / at * hz) +-- | 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' --- trend :: Dynamics Double -> Dynamics Double -> Dynamics Double --- -> Dynamics Double --- trend x at i = --- (x / smoothI x at (x / (1.0 + i * at)) - 1.0) / at +-- | 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. -- --- Integral Functions +-- 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 an integral with the specified derivative and initial value. --- If you want to create a loopback then you should use the 'Integ' type --- directly. The 'integ' function is just a wrapper that uses this type. -integ :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double) -integ diff i = - do x <- newInteg i - integDiff x diff - return $ integValue x +-- | 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), @@ -351,13 +568,54 @@ 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 (Dynamics m) tbl = +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 = @@ -382,9 +640,10 @@ in y -- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function. -lookupStepwiseD :: Dynamics Double -> Array Int (Double, Double) - -> Dynamics Double -lookupStepwiseD (Dynamics m) tbl = +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 = @@ -404,39 +663,126 @@ | otherwise = error "Incorrect index: table" in y --- -- --- -- Discrete Functions --- -- - --- delayTrans :: Dynamics a -> Dynamics Double -> Dynamics a --- -> (Dynamics a -> Dynamics a) -> Dynamics a --- delayTrans (Dynamics x) (Dynamics d) (Dynamics i) tr = tr $ Dynamics r --- where --- r p = do --- let t = parTime p --- sc = parSpecs p --- n = parIteration p --- a <- d p --- let t' = (t - a) - spcStartTime sc --- n' = fromIntegral $ floor $ t' / spcDT sc --- y | n' < 0 = i $ p { pointTime = spcStartTime sc, --- pointIteration = 0, --- pointPhase = 0 } --- | n' < n = x $ p { pointTime = t', --- pointIteration = n', --- pointPhase = -1 } --- | n' > n = error "Cannot return the future data: delay" --- | otherwise = error "Cannot return the current data: delay" --- y +-- +-- Discrete Functions +-- --- delay :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a --- delay x d = delayTrans x d x $ memo0 discrete +-- | 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 --- delay' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a --- delay' x d = delayTrans x d x $ memo0' discrete +-- | 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 --- delayI :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a --- delayI x d i = delayTrans x d i $ memo0 discrete +-- | 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 --- delayI' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a --- delayI' x d i = delayTrans x d i $ memo0' discrete +-- | 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 diff --git a/Simulation/Aivika/Dynamics/UVar.hs b/Simulation/Aivika/Dynamics/UVar.hs --- a/Simulation/Aivika/Dynamics/UVar.hs +++ b/Simulation/Aivika/Dynamics/UVar.hs @@ -3,11 +3,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.UVar --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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. diff --git a/Simulation/Aivika/Dynamics/Var.hs b/Simulation/Aivika/Dynamics/Var.hs --- a/Simulation/Aivika/Dynamics/Var.hs +++ b/Simulation/Aivika/Dynamics/Var.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Dynamics.Var --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- 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. diff --git a/Simulation/Aivika/PriorityQueue.hs b/Simulation/Aivika/PriorityQueue.hs --- a/Simulation/Aivika/PriorityQueue.hs +++ b/Simulation/Aivika/PriorityQueue.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.PriorityQueue --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- An imperative heap-based priority queue. -- diff --git a/Simulation/Aivika/Queue.hs b/Simulation/Aivika/Queue.hs --- a/Simulation/Aivika/Queue.hs +++ b/Simulation/Aivika/Queue.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Queue --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- An imperative double-linked queue. -- diff --git a/Simulation/Aivika/Statistics.hs b/Simulation/Aivika/Statistics.hs --- a/Simulation/Aivika/Statistics.hs +++ b/Simulation/Aivika/Statistics.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Statistics --- Copyright : Copyright (c) 2009-2012, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- Represents statistics. -- @@ -15,12 +15,19 @@ SamplingData(..), samplingStatsVariance, samplingStatsDeviation, + returnSamplingStats, + listSamplingStats, + fromIntSamplingStats, showSamplingStats, TimingStats(..), TimingData(..), timingStatsDeviation, + returnTimingStats, + fromIntTimingStats, showTimingStats) where +import Data.Monoid + -- | Defines data types that can be converted to 'Double'. class Ord a => ConvertableToDouble a where @@ -58,6 +65,15 @@ -- | Add a new sample to the statistics. addSamplingStats :: a -> SamplingStats a -> SamplingStats a + -- | Combine two statistics. + combineSamplingStats :: SamplingStats a -> SamplingStats a -> SamplingStats a + +instance SamplingData a => Monoid (SamplingStats a) where + + mempty = emptySamplingStats + + mappend = combineSamplingStats + instance SamplingData Double where emptySamplingStats = @@ -69,6 +85,8 @@ addSamplingStats = addSamplingStatsGeneric + combineSamplingStats = combineSamplingStatsGeneric + instance SamplingData Int where emptySamplingStats = @@ -79,6 +97,8 @@ samplingStatsMean2 = 0 / 0 } addSamplingStats = addSamplingStatsGeneric + + combineSamplingStats = combineSamplingStatsGeneric addSamplingStatsGeneric :: ConvertableToDouble a => a -> SamplingStats a -> SamplingStats a addSamplingStatsGeneric a stats @@ -103,6 +123,37 @@ k1 = 1.0 / n k2 = (n - 1.0) / n +combineSamplingStatsGeneric :: ConvertableToDouble a => + SamplingStats a -> SamplingStats a -> SamplingStats a +combineSamplingStatsGeneric stats1 stats2 + | c1 == 0 = stats2 + | c2 == 0 = stats1 + | otherwise = SamplingStats { samplingStatsCount = c, + samplingStatsMin = minZ, + samplingStatsMax = maxZ, + samplingStatsMean = meanZ, + samplingStatsMean2 = meanZ2 } + where c1 = samplingStatsCount stats1 + c2 = samplingStatsCount stats2 + c = c1 + c2 + n1 = fromIntegral c1 + n2 = fromIntegral c2 + n = n1 + n2 + minX = samplingStatsMin stats1 + minY = samplingStatsMin stats2 + minZ = min minX minY + maxX = samplingStatsMax stats1 + maxY = samplingStatsMax stats2 + maxZ = max maxX maxY + meanX = samplingStatsMean stats1 + meanY = samplingStatsMean stats2 + meanZ = k1 * meanX + k2 * meanY + meanX2 = samplingStatsMean2 stats1 + meanY2 = samplingStatsMean2 stats2 + meanZ2 = k1 * meanX2 + k2 * meanY2 + k1 = n1 / n + k2 = n2 / n + -- | Return the variance. samplingStatsVariance :: SamplingStats a -> Double samplingStatsVariance stats @@ -117,6 +168,20 @@ samplingStatsDeviation :: SamplingStats a -> Double samplingStatsDeviation = sqrt . samplingStatsVariance +-- | Return the statistics by a single sample. +returnSamplingStats :: SamplingData a => a -> SamplingStats a +returnSamplingStats x = addSamplingStats x emptySamplingStats + +-- | Create the statistics by the specified list of data. +listSamplingStats :: SamplingData a => [a] -> SamplingStats a +listSamplingStats = foldr addSamplingStats emptySamplingStats + +-- | Convert the statistics from integer to double values. +fromIntSamplingStats :: SamplingStats Int -> SamplingStats Double +fromIntSamplingStats stats = + stats { samplingStatsMin = fromIntegral $ samplingStatsMin stats, + samplingStatsMax = fromIntegral $ samplingStatsMax stats } + -- | Show the summary of the statistics with the specified indent. showSamplingStats :: (Show a) => SamplingStats a -> Int -> ShowS showSamplingStats stats indent = @@ -272,6 +337,16 @@ -- | Return the deviation. timingStatsDeviation :: TimingData a => TimingStats a -> Double timingStatsDeviation = sqrt . timingStatsVariance + +-- | Return the statistics by single timing data. +returnTimingStats :: TimingData a => Double -> a -> TimingStats a +returnTimingStats t a = addTimingStats t a emptyTimingStats + +-- | Convert the statistics from integer to double values. +fromIntTimingStats :: TimingStats Int -> TimingStats Double +fromIntTimingStats stats = + stats { timingStatsMin = fromIntegral $ timingStatsMin stats, + timingStatsMax = fromIntegral $ timingStatsMax stats } -- | Show the summary of the statistics with the specified indent. showTimingStats :: (Show a, TimingData a) => TimingStats a -> Int -> ShowS diff --git a/Simulation/Aivika/UVector.hs b/Simulation/Aivika/UVector.hs --- a/Simulation/Aivika/UVector.hs +++ b/Simulation/Aivika/UVector.hs @@ -3,11 +3,11 @@ -- | -- Module : Simulation.Aivika.UVector --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- An imperative unboxed vector. -- diff --git a/Simulation/Aivika/Vector.hs b/Simulation/Aivika/Vector.hs --- a/Simulation/Aivika/Vector.hs +++ b/Simulation/Aivika/Vector.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Vector --- Copyright : Copyright (c) 2009-2011, David Sorokin +-- Copyright : Copyright (c) 2009-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.0.3 +-- Tested with: GHC 7.6.3 -- -- An imperative vector. -- diff --git a/aivika.cabal b/aivika.cabal --- a/aivika.cabal +++ b/aivika.cabal @@ -1,5 +1,5 @@ name: aivika -version: 0.5.1 +version: 0.5.4 synopsis: A multi-paradigm simulation library description: Aivika is a small simulation library that covers many paradigms. @@ -38,17 +38,19 @@ category: Simulation license: BSD3 license-file: LICENSE -copyright: (c) 2009-2012. David Sorokin +copyright: (c) 2009-2013. David Sorokin author: David Sorokin maintainer: David Sorokin homepage: http://github.com/dsorokin/aivika cabal-version: >= 1.2.0 build-type: Simple -tested-with: GHC == 7.4.1 +tested-with: GHC == 7.6.3 extra-source-files: examples/BassDiffusion.hs examples/ChemicalReaction.hs + examples/ChemicalReactionRec.hs examples/FishBank.hs + examples/FishBankRec.hs examples/MachRep1.hs examples/MachRep1EventDriven.hs examples/MachRep1TimeDriven.hs @@ -96,13 +98,14 @@ Simulation.Aivika.Vector Simulation.Aivika.UVector - build-depends: base >= 3 && < 6, - mtl >= 1.1.0.2, + build-depends: base >= 4.5.0.0 && < 6, + mtl >= 2.1.1, array >= 0.3.0.0, containers >= 0.4.0.0, random >= 1.0.0.3 extensions: FlexibleContexts, - BangPatterns + BangPatterns, + RecursiveDo ghc-options: -O2 diff --git a/examples/ChemicalReactionRec.hs b/examples/ChemicalReactionRec.hs new file mode 100644 --- /dev/null +++ b/examples/ChemicalReactionRec.hs @@ -0,0 +1,22 @@ + +{-# LANGUAGE RecursiveDo #-} + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.SystemDynamics + +specs = Specs { spcStartTime = 0, + spcStopTime = 13, + spcDT = 0.01, + spcMethod = RungeKutta4 } + +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] + +main = runSimulation model specs >>= print diff --git a/examples/FishBankRec.hs b/examples/FishBankRec.hs new file mode 100644 --- /dev/null +++ b/examples/FishBankRec.hs @@ -0,0 +1,53 @@ + +{-# LANGUAGE RecursiveDo #-} + +import Data.Array + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.SystemDynamics + +specs = Specs { spcStartTime = 0, + spcStopTime = 13, + spcDT = 0.01, + -- spcDT = 0.000005, + spcMethod = RungeKutta4 } + +model :: Simulation Double +model = + mdo -- integrals -- + fish <- integ (fishHatchRate - fishDeathRate - totalCatchPerYear) 1000 + ships <- integ shipBuildingRate 10 + totalProfit <- integ annualProfit 0 + -- auxiliary values -- + let annualProfit = profit + area = 100 + carryingCapacity = 1000 + catchPerShip = + lookupDynamics density $ + listArray (1, 11) [(0.0, -0.048), (1.2, 10.875), (2.4, 17.194), + (3.6, 20.548), (4.8, 22.086), (6.0, 23.344), + (7.2, 23.903), (8.4, 24.462), (9.6, 24.882), + (10.8, 25.301), (12.0, 25.86)] + deathFraction = + lookupDynamics (fish / carryingCapacity) $ + listArray (1, 11) [(0.0, 5.161), (0.1, 5.161), (0.2, 5.161), + (0.3, 5.161), (0.4, 5.161), (0.5, 5.161), + (0.6, 5.118), (0.7, 5.247), (0.8, 5.849), + (0.9, 6.151), (10.0, 6.194)] + density = fish / area + fishDeathRate = maxDynamics 0 (fish * deathFraction) + fishHatchRate = maxDynamics 0 (fish * hatchFraction) + fishPrice = 20 + fractionInvested = 0.2 + hatchFraction = 6 + operatingCost = ships * 250 + profit = revenue - operatingCost + revenue = totalCatchPerYear * fishPrice + shipBuildingRate = maxDynamics 0 (profit * fractionInvested / shipCost) + shipCost = 300 + totalCatchPerYear = maxDynamics 0 (ships * catchPerShip) + -- results -- + runDynamicsInStopTime annualProfit + +main = runSimulation model specs >>= print diff --git a/examples/Furnace.hs b/examples/Furnace.hs --- a/examples/Furnace.hs +++ b/examples/Furnace.hs @@ -257,11 +257,12 @@ -- count the loaded ingots modifyRef (furnaceLoadCount furnace) (+ 1) --- | Iterate the furnace processing. -iterateFurnace :: Furnace -> Simulation (Dynamics ()) -iterateFurnace furnace = - let pits = furnacePits furnace - in iterateDynamics $ +-- | Start iterating the furnace processing through the event queue. +startIteratingFurnace :: Furnace -> Dynamics () +startIteratingFurnace furnace = + let queue = furnaceQueue furnace + pits = furnacePits furnace + in enqueueWithIntegTimes queue $ do ready <- ingotsReady furnace when ready $ do mapM_ (tryUnloadPit furnace) pits @@ -347,12 +348,11 @@ do queue <- newQueue furnace <- newFurnace queue pid <- newProcessID queue - + + -- initialize the furnace and start its iterating in start time runDynamicsInStartTime $ - initializeFurnace furnace - - -- get the furnace iterator - iterator <- iterateFurnace furnace + do initializeFurnace furnace + startIteratingFurnace furnace -- accept input ingots runDynamicsInStartTime $ @@ -361,9 +361,7 @@ -- run the model in the final time point runDynamicsInStopTime $ - do iterator -- iterate in each time point - - -- the ingots + do -- the ingots c0 <- readRef (furnaceTotalCount furnace) c1 <- readRef (furnaceLoadCount furnace) c2 <- readRef (furnaceUnloadCount furnace) diff --git a/examples/MachRep1TimeDriven.hs b/examples/MachRep1TimeDriven.hs --- a/examples/MachRep1TimeDriven.hs +++ b/examples/MachRep1TimeDriven.hs @@ -100,15 +100,15 @@ -- create two machines with type Dynamics () m1 <- machine m2 <- machine - - -- create strictly sequential computations - c1 <- iterateDynamics m1 - c2 <- iterateDynamics m2 - + + -- start the time-driven simulation of the machines through the event queue + runDynamicsInStartTime $ + do enqueueWithIntegTimes queue m1 + enqueueWithIntegTimes queue m2 + + -- return the result in the stop time runDynamicsInStopTime $ - do c1 -- involve in the simulation - c2 -- involve in the simulation - x <- readRef totalUpTime + do x <- readRef totalUpTime y <- stoptime return $ x / (2 * y)