aivika 0.5.1 → 0.5.4

raw patch · 38 files changed

+855/−218 lines, 38 filesdep ~basedep ~mtlPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base, mtl

API changes (from Hackage documentation)

Files

• LICENSE +1/−1
• Simulation/Aivika/Dynamics.hs +5/−3
• Simulation/Aivika/Dynamics/Agent.hs +2/−2
• Simulation/Aivika/Dynamics/Base.hs +2/−2
• Simulation/Aivika/Dynamics/Buffer.hs +2/−2
• Simulation/Aivika/Dynamics/Cont.hs +2/−2
• Simulation/Aivika/Dynamics/EventQueue.hs +2/−2
• Simulation/Aivika/Dynamics/FIFO.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Cont.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Dynamics.hs +15/−2
• Simulation/Aivika/Dynamics/Internal/Fold.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Interpolate.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Memo.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Process.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Signal.hs +2/−2
• Simulation/Aivika/Dynamics/Internal/Simulation.hs +94/−5
• Simulation/Aivika/Dynamics/Internal/Time.hs +2/−2
• Simulation/Aivika/Dynamics/LIFO.hs +2/−2
• Simulation/Aivika/Dynamics/Parameter.hs +30/−6
• Simulation/Aivika/Dynamics/Process.hs +2/−2
• Simulation/Aivika/Dynamics/Random.hs +19/−13
• Simulation/Aivika/Dynamics/Ref.hs +2/−2
• Simulation/Aivika/Dynamics/Resource.hs +2/−2
• Simulation/Aivika/Dynamics/Signal.hs +2/−2
• Simulation/Aivika/Dynamics/Simulation.hs +9/−3
• Simulation/Aivika/Dynamics/SystemDynamics.hs +454/−108
• Simulation/Aivika/Dynamics/UVar.hs +2/−2
• Simulation/Aivika/Dynamics/Var.hs +2/−2
• Simulation/Aivika/PriorityQueue.hs +2/−2
• Simulation/Aivika/Queue.hs +2/−2
• Simulation/Aivika/Statistics.hs +77/−2
• Simulation/Aivika/UVector.hs +2/−2
• Simulation/Aivika/Vector.hs +2/−2
• aivika.cabal +9/−6
• examples/ChemicalReactionRec.hs +22/−0
• examples/FishBankRec.hs +53/−0
• examples/Furnace.hs +11/−13
• examples/MachRep1TimeDriven.hs +8/−8

@@ -1,4 +1,4 @@-Copyright (c) 2009, 2010, 2011, 2012 David Sorokin <david.sorokin@gmail.com>+Copyright (c) 2009, 2010, 2011, 2012, 2013 David Sorokin <david.sorokin@gmail.com>  All rights reserved. 

@@ -1,23 +1,25 @@  -- | -- Module     : Simulation.Aivika.Dynamics--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Agent--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- This module introduces an agent-based modeling. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Base--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- This module defines basic functions for the 'Dynamics' monad. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Buffer--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Cont--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.EventQueue--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.FIFO--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.4.1+-- Tested with: GHC 7.6.3 -- -- This module defines the FIFO queue. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Cont--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,13 @@ +{-# LANGUAGE RecursiveDo #-}+ -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Dynamics--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 }

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Fold--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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.

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Interpolate--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- This module defines interpolation functions. -- These functions complement the memoization, possibly except for 

@@ -3,11 +3,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Memo--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Process--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Signal--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,13 @@ +{-# LANGUAGE RecursiveDo #-}+ -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Simulation--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 }  

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Internal.Time--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- This module defines the time parameters. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.LIFO--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.4.1+-- Tested with: GHC 7.6.3 -- -- This module defines the LIFO queue. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Parameter--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Process--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Random--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Ref--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- This module defines an updatable reference that depends on the event queue. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Resource--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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'.

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Signal--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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 

@@ -1,21 +1,27 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Simulation--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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

@@ -1,21 +1,28 @@ -{-# LANGUAGE FlexibleContexts, BangPatterns #-}+{-# LANGUAGE FlexibleContexts, BangPatterns, RecursiveDo #-}  -- | -- Module     : Simulation.Aivika.Dynamics.SystemDynamics--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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

@@ -3,11 +3,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.UVar--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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.

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Dynamics.Var--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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.

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.PriorityQueue--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- An imperative heap-based priority queue. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Queue--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- An imperative double-linked queue. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Statistics--- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.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

@@ -3,11 +3,11 @@  -- | -- Module     : Simulation.Aivika.UVector--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- An imperative unboxed vector. --

@@ -1,11 +1,11 @@  -- | -- Module     : Simulation.Aivika.Vector--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>+-- Copyright  : Copyright (c) 2009-2013, David Sorokin <david.sorokin@gmail.com> -- License    : BSD3 -- Maintainer : David Sorokin <david.sorokin@gmail.com> -- Stability  : experimental--- Tested with: GHC 7.0.3+-- Tested with: GHC 7.6.3 -- -- An imperative vector. --

@@ -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 <david.sorokin@gmail.com>+copyright:       (c) 2009-2013. David Sorokin <david.sorokin@gmail.com> author:          David Sorokin maintainer:      David Sorokin <david.sorokin@gmail.com> 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

@@ -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

@@ -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

@@ -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)

@@ -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)