aivika 0.4.2 → 0.4.3

raw patch · 8 files changed

+395/−166 lines, 8 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Simulation/Aivika/Dynamics/Base.hs +5/−0
• Simulation/Aivika/Dynamics/Internal/Dynamics.hs +33/−33
• Simulation/Aivika/Dynamics/Internal/Memo.hs +8/−8
• Simulation/Aivika/Dynamics/Internal/Time.hs +21/−8
• Simulation/Aivika/Dynamics/Signal.hs +30/−3
• Simulation/Aivika/Statistics.hs +281/−101
• aivika.cabal +1/−1
• examples/Furnace.hs +16/−12

@@ -17,6 +17,11 @@         dt,         time,         integTimes,+        isTimeInteg,+        integIteration,+        integIterationBnds,+        integIterationLoBnd,+        integIterationHiBnd,         -- * Interpolation and Initial Value         initDynamics,         discrete,

@@ -26,12 +26,12 @@         throwDynamics,         -- * Utilities         basicTime,-        iterationBnds,-        iterationHiBnd,-        iterationLoBnd,-        phaseBnds,-        phaseHiBnd,-        phaseLoBnd) where+        integIterationBnds,+        integIterationHiBnd,+        integIterationLoBnd,+        integPhaseBnds,+        integPhaseHiBnd,+        integPhaseLoBnd) where  import qualified Control.Exception as C import Control.Exception (IOException, throw, finally)@@ -61,50 +61,50 @@                      pointPhase :: Int       -- ^ the current phase                    } deriving (Eq, Ord, Show)            --- | Returns the iterations starting from zero.-iterations :: Specs -> [Int]-iterations sc = [i1 .. i2] where+-- | Returns the integration iterations starting from zero.+integIterations :: Specs -> [Int]+integIterations sc = [i1 .. i2] where   i1 = 0   i2 = round ((spcStopTime sc -                 spcStartTime sc) / spcDT sc) --- | Returns the first and last iterations.-iterationBnds :: Specs -> (Int, Int)-iterationBnds sc = (0, round ((spcStopTime sc - -                               spcStartTime sc) / spcDT sc))+-- | Returns the first and last integration iterations.+integIterationBnds :: Specs -> (Int, Int)+integIterationBnds sc = (0, round ((spcStopTime sc - +                                    spcStartTime sc) / spcDT sc)) --- | Returns the first iteration, i.e. zero.-iterationLoBnd :: Specs -> Int-iterationLoBnd sc = 0+-- | Returns the first integration iteration, i.e. zero.+integIterationLoBnd :: Specs -> Int+integIterationLoBnd sc = 0 --- | Returns the last iteration.-iterationHiBnd :: Specs -> Int-iterationHiBnd sc = round ((spcStopTime sc - -                            spcStartTime sc) / spcDT sc)+-- | Returns the last integration iteration.+integIterationHiBnd :: Specs -> Int+integIterationHiBnd sc = round ((spcStopTime sc - +                                 spcStartTime sc) / spcDT sc)  -- | Returns the phases for the specified simulation specs starting from zero.-phases :: Specs -> [Int]-phases sc = +integPhases :: Specs -> [Int]+integPhases sc =    case spcMethod sc of     Euler -> [0]     RungeKutta2 -> [0, 1]     RungeKutta4 -> [0, 1, 2, 3] --- | Returns the first and last phases.-phaseBnds :: Specs -> (Int, Int)-phaseBnds sc = +-- | Returns the first and last integration phases.+integPhaseBnds :: Specs -> (Int, Int)+integPhaseBnds sc =    case spcMethod sc of     Euler -> (0, 0)     RungeKutta2 -> (0, 1)     RungeKutta4 -> (0, 3) --- | Returns the first phase, i.e. zero.-phaseLoBnd :: Specs -> Int-phaseLoBnd sc = 0+-- | Returns the first integration phase, i.e. zero.+integPhaseLoBnd :: Specs -> Int+integPhaseLoBnd sc = 0                   --- | Returns the last phase, 1 for Euler's method, 2 for RK2 and 4 for RK4.-phaseHiBnd :: Specs -> Int-phaseHiBnd sc = +-- | Returns the last integration phase, 0 for Euler's method, 1 for RK2 and 3 for RK4.+integPhaseHiBnd :: Specs -> Int+integPhaseHiBnd sc =    case spcMethod sc of     Euler -> 0     RungeKutta2 -> 1@@ -159,7 +159,7 @@ runDynamicsInStopTime (Dynamics m) =   Simulation $ \r ->   do let sc = runSpecs r -         n  = iterationHiBnd sc+         n  = integIterationHiBnd sc          t  = basicTime sc n 0      m Point { pointSpecs = sc,                pointRun = r,@@ -172,7 +172,7 @@ runDynamicsInIntegTimes (Dynamics m) =   Simulation $ \r ->   do let sc = runSpecs r-         (nl, nu) = iterationBnds sc+         (nl, nu) = integIterationBnds sc          point n = Point { pointSpecs = sc,                            pointRun = r,                            pointTime = basicTime sc n 0,

@@ -43,8 +43,8 @@ memo (Dynamics m) =    Simulation $ \r ->   do let sc = runSpecs r-         (phl, phu) = phaseBnds sc-         (nl, nu)   = iterationBnds sc+         (phl, phu) = integPhaseBnds sc+         (nl, nu)   = integIterationBnds sc      arr   <- newMemoArray_ ((phl, nl), (phu, nu))      nref  <- newIORef 0      phref <- newIORef 0@@ -52,7 +52,7 @@            do let sc  = pointSpecs p                   n   = pointIteration p                   ph  = pointPhase p-                  phu = phaseHiBnd sc +                  phu = integPhaseHiBnd sc                    loop n' ph' =                      if (n' > n) || ((n' == n) && (ph' > ph))                      then @@ -80,8 +80,8 @@ umemo (Dynamics m) =    Simulation $ \r ->   do let sc = runSpecs r-         (phl, phu) = phaseBnds sc-         (nl, nu)   = iterationBnds sc+         (phl, phu) = integPhaseBnds sc+         (nl, nu)   = integIterationBnds sc      arr   <- newMemoUArray_ ((phl, nl), (phu, nu))      nref  <- newIORef 0      phref <- newIORef 0@@ -89,7 +89,7 @@            do let sc  = pointSpecs p                   n   = pointIteration p                   ph  = pointPhase p-                  phu = phaseHiBnd sc +                  phu = integPhaseHiBnd sc                    loop n' ph' =                      if (n' > n) || ((n' == n) && (ph' > ph))                      then @@ -122,7 +122,7 @@ memo0 (Dynamics m) =    Simulation $ \r ->   do let sc   = runSpecs r-         bnds = iterationBnds sc+         bnds = integIterationBnds sc      arr  <- newMemoArray_ bnds      nref <- newIORef 0      let r p =@@ -150,7 +150,7 @@ umemo0 (Dynamics m) =    Simulation $ \r ->   do let sc   = runSpecs r-         bnds = iterationBnds sc+         bnds = integIterationBnds sc      arr  <- newMemoUArray_ bnds      nref <- newIORef 0      let r p =

@@ -15,7 +15,10 @@         stoptime,          dt,          time,-        integTimes) where+        integTimes, +        isTimeInteg,+        integIteration,+        integPhase) where  import Simulation.Aivika.Dynamics.Internal.Simulation import Simulation.Aivika.Dynamics.Internal.Dynamics@@ -37,10 +40,20 @@ time = Dynamics $ return . pointTime   -- | Return the integration time points.-integTimes :: Simulation [Double]-integTimes =-  Simulation $ \r ->-  do let sc  = runSpecs r-         (nl, nu) = iterationBnds sc-         t n = basicTime sc n 0-     return $ map t [nl .. nu]+integTimes :: Specs -> [Double]+integTimes sc = map t [nl .. nu]+  where (nl, nu) = integIterationBnds sc+        t n = basicTime sc n 0+     +-- | Whether the current time is an integration time.+isTimeInteg :: Dynamics Bool+isTimeInteg = Dynamics $ \p -> return $ pointPhase p >= 0++-- | Return the integration iteration closest to the current simulation time.+integIteration :: Dynamics Int+integIteration = Dynamics $ return . pointIteration++-- | Return the integration phase for the current simulation time.+-- It is @(-1)@ for non-integration time points.+integPhase :: Dynamics Int+integPhase = Dynamics $ return . pointPhase

@@ -40,6 +40,7 @@         SignalHistory,         signalHistorySignal,         newSignalHistory,+        newSignalHistoryThrough,         readSignalHistory) where  import Data.IORef@@ -54,6 +55,7 @@ import Simulation.Aivika.Dynamics.Internal.Dynamics import Simulation.Aivika.Dynamics.Internal.Cont import Simulation.Aivika.Dynamics.Internal.Process+import Simulation.Aivika.Dynamics.Base  import qualified Simulation.Aivika.Vector as V import qualified Simulation.Aivika.UVector as UV@@ -105,6 +107,31 @@                             signalHistoryTimes  = ts,                             signalHistoryValues = xs }        +-- | Create a history of the signal values with delay through the event queue.+-- The history will be created at the same simulation time, just the corresponded +-- handler will be subscribed to the signal after the new event will be processed +-- by the queue. +-- +-- It is very useful if we want the signal won't be triggered at the current +-- time until we complete some preparation. This is relatated to the fact that+-- the signal is updated at time of subscribing the handler. So, if we subscribe+-- to the signal which must be triggered at the current time then it will be+-- triggered. Using the event queue allows us to complete some preparation logic+-- before the signal will be triggered at the same simulation time point.+newSignalHistoryThrough :: EventQueue -> Signal a -> Dynamics (SignalHistory a)+newSignalHistoryThrough q signal =+  do ts <- liftIO UV.newVector+     xs <- liftIO V.newVector+     t  <- time+     enqueue q t $+       handleSignal_ signal $ \a ->+       Dynamics $ \p ->+       do liftIO $ UV.appendVector ts (pointTime p)+          liftIO $ V.appendVector xs a+     return SignalHistory { signalHistorySignal = signal,+                            signalHistoryTimes  = ts,+                            signalHistoryValues = xs }+        -- | Read the history of signal values. readSignalHistory :: SignalHistory a -> Dynamics (Array Int Double, Array Int a) readSignalHistory history =@@ -144,7 +171,7 @@ newSignalInIntegTimes q =   Dynamics $ \p ->   do let sc  = pointSpecs p-         (nl, nu) = iterationBnds sc+         (nl, nu) = integIterationBnds sc          point n = Point { pointSpecs = sc,                            pointRun = pointRun p,                            pointTime = basicTime sc n 0,@@ -159,7 +186,7 @@ newSignalInStartTime q =   Dynamics $ \p ->   do let sc  = pointSpecs p-         (nl, nu) = iterationBnds sc+         (nl, nu) = integIterationBnds sc          point n = Point { pointSpecs = sc,                            pointRun = pointRun p,                            pointTime = basicTime sc n 0,@@ -173,7 +200,7 @@ newSignalInStopTime q =   Dynamics $ \p ->   do let sc  = pointSpecs p-         (nl, nu) = iterationBnds sc+         (nl, nu) = integIterationBnds sc          point n = Point { pointSpecs = sc,                            pointRun = pointRun p,                            pointTime = basicTime sc n 0,

@@ -1,6 +1,4 @@ -{-# LANGUAGE FlexibleContexts #-}- -- | -- Module     : Simulation.Aivika.Statistics -- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>@@ -9,111 +7,293 @@ -- Stability  : experimental -- Tested with: GHC 7.0.3 ----- Represents statistics and results.+-- Represents statistics. --+ module Simulation.Aivika.Statistics-       (Statistics, -        newStatistics,-        addStatistics,-        statisticsData,-        analyzeData,-        AnalysisResults(..),-        showResults) where +       (SamplingStats(..),+        SamplingData(..),+        samplingStatsVariance,+        samplingStatsDeviation,+        showSamplingStats,+        TimingStats(..),+        TimingData(..),+        timingStatsDeviation,+        showTimingStats) where  -import Data.Foldable-import Data.Array-import Data.Array.IO-import Control.Monad-import Control.Monad.Trans-import Control.Concurrent.MVar+-- | Defines data types that can be converted to 'Double'.+class Ord a => ConvertableToDouble a where+  +  -- | Convert the value to 'Double'.+  convertToDouble :: a -> Double+  +instance ConvertableToDouble Double where+  convertToDouble = id+  +instance ConvertableToDouble Int where+  convertToDouble = fromInteger . toInteger -import Simulation.Aivika.UVector+-- | Describes when the statistics consists of only samples +-- not bound to the simulation time.+data SamplingStats a =       +  SamplingStats { samplingStatsCount :: !Int,+                  -- ^ The total number of samples.+                  samplingStatsMin   :: !a,+                  -- ^ The minimum value among the samples.+                  samplingStatsMax   :: !a,+                  -- ^ The maximum value among the samples.+                  samplingStatsMean  :: !Double,+                  -- ^ The average value.+                  samplingStatsMean2 :: !Double +                  -- ^ The average square value.+                }+  deriving (Eq, Ord, Show)+           +-- | Specifies data type from which values we can gather the statistics.           +class SamplingData a where           +  +  -- | An empty statistics that has no samples.           +  emptySamplingStats :: SamplingStats a+           +  -- | Add a new sample to the statistics.+  addSamplingStats :: a -> SamplingStats a -> SamplingStats a --- | Represents statistics. --- --- All functions with the statistics in this module are thread-safe. Therefore --- you can use them in experiments when parallel simulations execute simultaneously.-data Statistics a = Statistics { statData :: UVector a, -                                 statLock :: MVar () }+instance SamplingData Double where --- | Create new statistics.-newStatistics :: (MArray IOUArray a IO) => IO (Statistics a)-newStatistics = -  do v <- newVector-     l <- newMVar ()-     return Statistics { statData = v, -                         statLock = l }+  emptySamplingStats =+    SamplingStats { samplingStatsCount = 0,+                    samplingStatsMin   = 1 / 0,+                    samplingStatsMax   = (-1) / 0,+                    samplingStatsMean  = 0 / 0,+                    samplingStatsMean2 = 0 / 0 }+    +  addSamplingStats = addSamplingStatsGeneric+  +instance SamplingData Int where --- | Add data to the statistics. It is thread-safe.-addStatistics :: (MArray IOUArray a IO) => Statistics a -> a -> IO ()-addStatistics s x = -  withMVar (statLock s) $ \() ->-  appendVector (statData s) x+  emptySamplingStats =+    SamplingStats { samplingStatsCount = 0,+                    samplingStatsMin   = maxBound,+                    samplingStatsMax   = minBound,+                    samplingStatsMean  = 0 / 0,+                    samplingStatsMean2 = 0 / 0 }+    +  addSamplingStats = addSamplingStatsGeneric+  +addSamplingStatsGeneric :: ConvertableToDouble a => a -> SamplingStats a -> SamplingStats a+addSamplingStatsGeneric a stats +  | isNaN x    = stats+  | count == 1 = SamplingStats { samplingStatsCount = 1,+                                 samplingStatsMin   = a,+                                 samplingStatsMax   = a,+                                 samplingStatsMean  = x,+                                 samplingStatsMean2 = x * x }+  | otherwise  = SamplingStats { samplingStatsCount = count,+                                 samplingStatsMin   = minX,+                                 samplingStatsMax   = maxX,+                                 samplingStatsMean  = meanX,+                                 samplingStatsMean2 = meanX2 }+    where count  = 1 + samplingStatsCount stats+          minX   = a `seq` min a (samplingStatsMin stats)+          maxX   = a `seq` max a (samplingStatsMax stats)+          meanX  = k1 * x + k2 * samplingStatsMean stats+          meanX2 = k1 * x * x + k2 * samplingStatsMean2 stats+          n      = fromInteger $ toInteger count+          x      = convertToDouble a+          k1     = 1.0 / n+          k2     = (n - 1.0) / n --- | Return the statistics data. It is thread-safe.-statisticsData :: (MArray IOUArray a IO) => Statistics a -> IO (Array Int a)-statisticsData s =-  withMVar (statLock s) $ \() -> freezeVector (statData s)-       --- | Represents the results of the statistic analysis.-data AnalysisResults a = -  AnalysisResults { resultsData     :: Array Int a,-                    -- ^ Statistic data.-                    resultsMean     :: Double,-                    -- ^ The average value.-                    resultsVariance :: Double,-                    -- ^ The variance.-                    resultsMin      :: a,-                    -- ^ The minimum value.-                    resultsMax      :: a -                    -- ^ The maximum value.-                  } deriving (Eq, Ord, Show)+-- | Return the variance.+samplingStatsVariance :: SamplingStats a -> Double+samplingStatsVariance stats+  | count == 1 = meanX2 - meanX * meanX+  | otherwise  = (meanX2 - meanX * meanX) * (n / (n - 1))+    where count  = samplingStatsCount stats+          meanX  = samplingStatsMean stats+          meanX2 = samplingStatsMean2 stats+          n      = fromInteger $ toInteger count+          +-- | Return the deviation.          +samplingStatsDeviation :: SamplingStats a -> Double+samplingStatsDeviation = sqrt . samplingStatsVariance --- | Analyze data.-analyzeData :: Real a => Array Int a -> AnalysisResults a-analyzeData xs =-  let (i1, i2) = bounds xs-      meanx = foldl' (\y i -> y * (1 - k i) + f i * k i) 0 [i1 .. i2]-      sqrx  = foldl' (\y i -> y * (1 - k i) + g i * k i) 0 [i1 .. i2]-      minx  = foldl' (\y i -> if i == 0 then x i else min y (x i)) 0 [i1 .. i2]-      maxx  = foldl' (\y i -> if i == 0 then x i else max y (x i)) 0 [i1 .. i2]-      x i = xs ! i-      f i = fromRational (toRational (x i))-      g i = let y = f i in y * y-      k i = 1 / fromInteger (toInteger (i - i1 + 1))-  in AnalysisResults { resultsData = xs,-                       resultsMean = meanx,-                       resultsVariance = sqrx - meanx * meanx,-                       resultsMin = minx,-                       resultsMax = maxx }-       --- | Show the results of analysis with the specified indent.       -showResults :: (Show a) => AnalysisResults a -> Int -> ShowS-showResults rs indent =-  let (i1, i2) = bounds (resultsData rs)-      tab = replicate indent ' '-  in if i1 <= i2-     then-       showString tab .-       showString "mean      = " . shows (resultsMean rs) . -       showString "\n" . -       showString tab .-       showString "deviation = " . shows (sqrt (resultsVariance rs)) . -       showString "\n" .-       showString tab .-       showString "minimum   = " . shows (resultsMin rs) . -       showString "\n" .-       showString tab .-       showString "maximum   = " . shows (resultsMax rs)-     else-       showString tab .-       showString "mean      = ---" .-       showString "\n" . -       showString tab .-       showString "deviation = ---" .-       showString "\n" . -       showString tab .-       showString "minimum   = ---" .-       showString "\n" . -       showString tab .-       showString "maximum   = ---"+-- | Show the summary of the statistics with the specified indent.       +showSamplingStats :: (Show a) => SamplingStats a -> Int -> ShowS+showSamplingStats stats indent =+  let tab = replicate indent ' '+  in showString tab .+     showString "count     = " . shows (samplingStatsCount stats) . +     showString "\n" . +     showString tab .+     showString "mean      = " . shows (samplingStatsMean stats) . +     showString "\n" . +     showString tab .+     showString "deviation = " . shows (samplingStatsDeviation stats) . +     showString "\n" .+     showString tab .+     showString "minimum   = " . shows (samplingStatsMin stats) . +     showString "\n" .+     showString tab .+     showString "maximum   = " . shows (samplingStatsMax stats)+     +-- | This is the timing statistics where data are bound to the time.+data TimingStats a =+  TimingStats { timingStatsCount     :: !Int,+                -- ^ Return the number of samples.+                timingStatsMin       :: !a,+                -- ^ Return the minimum value.+                timingStatsMax       :: !a,+                -- ^ Return the maximum value.+                timingStatsMinTime   :: !Double,+                -- ^ Return the time at which the minimum is attained.+                timingStatsMaxTime   :: !Double,+                -- ^ Return the time at which the maximum is attained.+                timingStatsStartTime :: !Double,+                -- ^ Return the start time of sampling.+                timingStatsLastTime  :: !Double,+                -- ^ Return the last time of sampling.+                timingStatsSum       :: !Double,+                -- ^ Return the sum of values.+                timingStatsSum2      :: !Double +                -- ^ Return the sum of square values.+                } deriving (Eq, Ord, Show)+                           +-- | Defines the data type from which values we can gather the timing statistics.+class TimingData a where                           +  +  -- | An empty statistics that has no samples.+  emptyTimingStats :: TimingStats a+  +  -- | Add a sample with the specified time to the statistics.+  addTimingStats :: Double -> a -> TimingStats a -> TimingStats a+  +  -- | Return the average value.+  timingStatsMean :: TimingStats a -> Double+  +  -- | Return the variance.+  timingStatsVariance :: TimingStats a -> Double+  +instance TimingData Double where+  +  emptyTimingStats = +    TimingStats { timingStatsCount     = 0,+                  timingStatsMin       = 1 / 0,+                  timingStatsMax       = (-1) / 0,+                  timingStatsMinTime   = 1 / 0,+                  timingStatsMaxTime   = (-1) / 0,+                  timingStatsStartTime = 1 / 0,+                  timingStatsLastTime  = (-1) / 0,+                  timingStatsSum       = 0,+                  timingStatsSum2      = 0 }+    +  addTimingStats      = addTimingStatsGeneric+  timingStatsMean     = timingStatsMeanGeneric+  timingStatsVariance = timingStatsVarianceGeneric++instance TimingData Int where+  +  emptyTimingStats = +    TimingStats { timingStatsCount     = 0,+                  timingStatsMin       = maxBound,+                  timingStatsMax       = minBound,+                  timingStatsMinTime   = 1 / 0,+                  timingStatsMaxTime   = (-1) / 0,+                  timingStatsStartTime = 1 / 0,+                  timingStatsLastTime  = (-1) / 0,+                  timingStatsSum       = 0,+                  timingStatsSum2      = 0 }+    +  addTimingStats      = addTimingStatsGeneric+  timingStatsMean     = timingStatsMeanGeneric+  timingStatsVariance = timingStatsVarianceGeneric++addTimingStatsGeneric :: ConvertableToDouble a => Double -> a -> TimingStats a -> TimingStats a+addTimingStatsGeneric t a stats+  | isNaN x    = stats+  | count == 1 = TimingStats { timingStatsCount     = 1,+                               timingStatsMin       = a,+                               timingStatsMax       = a,+                               timingStatsMinTime   = t,+                               timingStatsMaxTime   = t,+                               timingStatsStartTime = t,+                               timingStatsLastTime  = t,+                               timingStatsSum       = 0,+                               timingStatsSum2      = 0 }+  | t < t'     = error "The current time cannot be less than the previous one: addTimingStats"+  | otherwise  = TimingStats { timingStatsCount     = count,+                               timingStatsMin       = minX,+                               timingStatsMax       = maxX,+                               timingStatsMinTime   = minT,+                               timingStatsMaxTime   = maxT,+                               timingStatsStartTime = t0,+                               timingStatsLastTime  = t,+                               timingStatsSum       = sumX,+                               timingStatsSum2      = sumX2 }+    where count = 1 + timingStatsCount stats+          minX' = timingStatsMin stats+          maxX' = timingStatsMax stats+          minX  = a `seq` min a minX'+          maxX  = a `seq` max a maxX'+          minT | a < minX' = t+               | otherwise = timingStatsMinTime stats+          maxT | a > maxX' = t+               | otherwise = timingStatsMaxTime stats+          t0 = timingStatsStartTime stats+          t' = timingStatsLastTime stats+          x  = convertToDouble a+          sumX'  = timingStatsSum stats+          sumX   = sumX' + (t - t') * x+          sumX2' = timingStatsSum2 stats+          sumX2  = sumX2' + (t - t') * x * x+      +timingStatsMeanGeneric :: ConvertableToDouble a => TimingStats a -> Double+timingStatsMeanGeneric stats +  | t1 > t0   = sumX / (t1 - t0)+  | otherwise = minX+    where t0   = timingStatsStartTime stats+          t1   = timingStatsLastTime stats+          sumX = timingStatsSum stats+          minX = convertToDouble $ timingStatsMin stats+  +timingStatsMean2Generic :: ConvertableToDouble a => TimingStats a -> Double+timingStatsMean2Generic stats+  | t1 > t0   = sumX2 / (t1 - t0)+  | otherwise = minX * minX+    where t0    = timingStatsStartTime stats+          t1    = timingStatsLastTime stats+          sumX2 = timingStatsSum2 stats+          minX  = convertToDouble $ timingStatsMin stats++timingStatsVarianceGeneric :: ConvertableToDouble a => TimingStats a -> Double+timingStatsVarianceGeneric stats = ex2 - ex * ex+  where ex  = timingStatsMeanGeneric stats+        ex2 = timingStatsMean2Generic stats+                +-- | Return the deviation.              +timingStatsDeviation :: TimingData a => TimingStats a -> Double+timingStatsDeviation = sqrt . timingStatsVariance++-- | Show the summary of the statistics with the specified indent.       +showTimingStats :: (Show a, TimingData a) => TimingStats a -> Int -> ShowS+showTimingStats stats indent =+  let tab = replicate indent ' '+  in showString tab .+     showString "count     = " . shows (timingStatsCount stats) . +     showString "\n" . +     showString tab .+     showString "mean      = " . shows (timingStatsMean stats) . +     showString "\n" . +     showString tab .+     showString "deviation = " . shows (timingStatsDeviation stats) . +     showString "\n" .+     showString tab .+     showString "minimum   = " . shows (timingStatsMin stats) . +     showString " at t = " . shows (timingStatsMinTime stats) .+     showString "\n" .+     showString tab .+     showString "maximum   = " . shows (timingStatsMax stats) .+     showString " at t = " . shows (timingStatsMaxTime stats) .+     showString "\n" .+     showString "t in [" . shows (timingStatsStartTime stats) .+     showString ", " . shows (timingStatsLastTime stats) .+     showString "]"

@@ -1,5 +1,5 @@ name:            aivika-version:         0.4.2+version:         0.4.3 synopsis:        A multi-paradigm simulation library description:     Aivika is a small simulation library that covers many paradigms. 

@@ -46,13 +46,13 @@             -- ^ The pits for ingots.             furnacePitCount :: UVar Int,             -- ^ The count of active pits with ingots.-            furnacePitCountStats :: Statistics Int,+            furnacePitCountStats :: Ref (SamplingStats Int),             -- ^ The statistics about the active pits.             furnaceAwaitingIngots :: Q.Queue Ingot,             -- ^ The awaiting ingots in the queue.             furnaceQueueCount :: UVar Int,             -- ^ The queue count.-            furnaceQueueCountStats :: Statistics Int,+            furnaceQueueCountStats :: Ref (SamplingStats Int),             -- ^ The statistics about the queue count.             furnaceWaitCount :: Ref Int,             -- ^ The count of awaiting ingots.@@ -103,10 +103,10 @@   do normalGen <- liftIO normalGen      pits <- sequence [newPit queue | i <- [1..10]]      pitCount <- newUVar queue 0-     pitCountStats <- liftIO newStatistics+     pitCountStats <- newRef queue emptySamplingStats      awaitingIngots <- liftIO Q.newQueue      queueCount <- newUVar queue 0-     queueCountStats <- liftIO newStatistics+     queueCountStats <- newRef queue emptySamplingStats      waitCount <- newRef queue 0      waitTime <- newRef queue 0.0      heatingTime <- newRef queue 0.0@@ -198,7 +198,8 @@           t' <- time           modifyUVar (furnaceQueueCount furnace) (+ (-1))           c <- readUVar (furnaceQueueCount furnace)-          liftIO $ addStatistics (furnaceQueueCountStats furnace) c+          modifyRef (furnaceQueueCountStats furnace) $+            addSamplingStats c           loadIngot (ingot { ingotLoadTime = t',                              ingotLoadTemp = 400.0 }) pit               @@ -213,7 +214,8 @@      let furnace = ingotFurnace ingot      count <- readUVar (furnacePitCount furnace)      writeUVar (furnacePitCount furnace) (count - 1)-     liftIO $ addStatistics (furnacePitCountStats furnace) (count - 1)+     modifyRef (furnacePitCountStats furnace) $+       addSamplingStats (count - 1)            -- how long did we heat the ingot up?      t' <- time@@ -236,7 +238,8 @@      let furnace = ingotFurnace ingot      count <- readUVar (furnacePitCount furnace)      writeUVar (furnacePitCount furnace) (count + 1)-     liftIO $ addStatistics (furnacePitCountStats furnace) (count + 1)+     modifyRef (furnacePitCountStats furnace) $+       addSamplingStats (count + 1)            -- decrease the furnace temperature      h <- readRef (furnaceTemp furnace)@@ -293,7 +296,8 @@                liftIO $ Q.enqueue ingots ingot                modifyUVar (furnaceQueueCount furnace) (+ 1)                c <- readUVar (furnaceQueueCount furnace)-               liftIO $ addStatistics (furnaceQueueCountStats furnace) c+               modifyRef (furnaceQueueCountStats furnace) $+                 addSamplingStats c        else do pit:_ <- emptyPits furnace                loadIngot ingot pit        @@ -384,19 +388,19 @@             putStrLn ""                            -- the ingots in pits-          r2 <- fmap analyzeData $ liftIO $ statisticsData (furnacePitCountStats furnace)+          r2 <- readRef (furnacePitCountStats furnace)                          liftIO $ do             putStrLn "The ingots in pits: "-            putStrLn $ showResults r2 2 []+            putStrLn $ showSamplingStats r2 2 []             putStrLn ""                          -- the queue size-          r3 <- fmap analyzeData $ liftIO $ statisticsData (furnaceQueueCountStats furnace)+          r3 <- readRef (furnaceQueueCountStats furnace)                 liftIO $ do             putStrLn "The queue size: "-            putStrLn $ showResults r3 2 []+            putStrLn $ showSamplingStats r3 2 []             putStrLn ""                          -- the mean wait time in the queue