diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2012 David Sorokin + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +SUCH DAMAGE. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/Simulation/Aivika/Experiment/Chart.hs b/Simulation/Aivika/Experiment/Chart.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Experiment/Chart.hs @@ -0,0 +1,54 @@ + +-- | +-- Module : Simulation.Aivika.Experiment.Chart +-- Copyright : Copyright (c) 2012, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.4.1 +-- +-- The module defines some utilities used in the charting. +-- + +module Simulation.Aivika.Experiment.Chart + (colourisePlotLines, + colourisePlotFillBetween) where + +import Data.Colour +import Data.Colour.Names +import Data.Accessor + +import Graphics.Rendering.Chart + +-- | Colourise the plot lines. +colourisePlotLines :: [PlotLines x y -> PlotLines x y] +colourisePlotLines = + (plot_lines_style .> line_color ^= opaque blue) : + (plot_lines_style .> line_color ^= opaque green) : + (plot_lines_style .> line_color ^= opaque red) : + (plot_lines_style .> line_color ^= opaque black) : + (plot_lines_style .> line_color ^= opaque grey) : + (plot_lines_style .> line_color ^= opaque purple) : + (plot_lines_style .> line_color ^= opaque violet) : + (plot_lines_style .> line_color ^= opaque darkblue) : + (plot_lines_style .> line_color ^= opaque darkgreen) : + (plot_lines_style .> line_color ^= opaque darkgrey) : + (plot_lines_style .> line_color ^= opaque darkviolet) : + colourisePlotLines + +-- | Colourise the plot areas that are filling between +-- any two lines. +colourisePlotFillBetween :: [PlotFillBetween x y -> PlotFillBetween x y] +colourisePlotFillBetween = + (plot_fillbetween_style ^= solidFillStyle (withOpacity blue 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity green 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity red 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity black 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity grey 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity purple 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity violet 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity darkblue 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity darkgreen 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity darkgrey 0.4)) : + (plot_fillbetween_style ^= solidFillStyle (withOpacity darkviolet 0.4)) : + colourisePlotFillBetween diff --git a/Simulation/Aivika/Experiment/DeviationChartView.hs b/Simulation/Aivika/Experiment/DeviationChartView.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Experiment/DeviationChartView.hs @@ -0,0 +1,315 @@ + +-- | +-- Module : Simulation.Aivika.Experiment.DeviationChartView +-- Copyright : Copyright (c) 2012, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.4.1 +-- +-- The module defines 'DeviationChartView' that saves the deviation +-- chart as the PNG file. +-- + +module Simulation.Aivika.Experiment.DeviationChartView + (DeviationChartView(..), + defaultDeviationChartView) where + +import Control.Monad +import Control.Monad.Trans +import Control.Concurrent.MVar + +import qualified Data.Map as M +import Data.IORef +import Data.Maybe +import Data.Either +import Data.Array +import Data.Array.IO + +import Data.Accessor + +import System.IO +import System.FilePath + +import Data.String.Utils (replace) + +import Graphics.Rendering.Chart + +import Simulation.Aivika.Experiment +import Simulation.Aivika.Experiment.HtmlWriter +import Simulation.Aivika.Experiment.Utils (divideBy) +import Simulation.Aivika.Experiment.Chart (colourisePlotLines, colourisePlotFillBetween) + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.Signal +import Simulation.Aivika.Dynamics.EventQueue +import Simulation.Aivika.Dynamics.Base (starttime, integIterationBnds, integTimes, integIteration) +import Simulation.Aivika.Statistics + +-- | Defines the 'View' that saves the deviation chart +-- in the PNG file. +data DeviationChartView = + DeviationChartView { deviationChartTitle :: String, + -- ^ This is a title used in HTML and chart. + deviationChartDescription :: String, + -- ^ This is a description in the HTML. + deviationChartWidth :: Int, + -- ^ The width of the chart. + deviationChartHeight :: Int, + -- ^ The height of the chart. + deviationChartFileName :: FileName, + -- ^ It defines the file name for the PNG file. + -- It may include special variable @$TITLE@. + -- + -- An example is + -- + -- @ + -- deviationChartFileName = UniqueFileName \"$TITLE\", \".png\" + -- @ + deviationChartPredicate :: Dynamics Bool, + -- ^ It specifies the predicate that defines + -- when we count data when plotting the chart. + deviationChartSeries :: [Either String String], + -- ^ It contains the labels of data plotted + -- in the chart. + deviationChartPlotLines :: [PlotLines Double Double -> + PlotLines Double Double], + -- ^ Probably, an infinite sequence of plot + -- transformations based on which the plot + -- is constructed for each series. Generally, + -- it must not coincide with a sequence of + -- labels as one label may denote a whole list + -- or an array of data providers. + -- + -- Here you can define a colour or style of + -- of the plot lines. + deviationChartPlotFillBetween :: [PlotFillBetween Double Double -> + PlotFillBetween Double Double], + -- ^ Corresponds exactly to 'deviationChartPlotLines' + -- but used for plotting the deviation areas + -- by the rule of 3-sigma, while the former + -- is used for plotting the trends of the + -- random processes. + deviationChartLayout :: Layout1 Double Double -> + Layout1 Double Double + -- ^ A transformation of the plot layout, + -- where you can redefine the axes, for example. + } + +-- | The default deviation chart view. +defaultDeviationChartView :: DeviationChartView +defaultDeviationChartView = + DeviationChartView { deviationChartTitle = "Deviation Chart", + deviationChartDescription = [], + deviationChartWidth = 640, + deviationChartHeight = 480, + deviationChartFileName = UniqueFileName "$TITLE" ".png", + deviationChartPredicate = return True, + deviationChartSeries = [], + deviationChartPlotLines = colourisePlotLines, + deviationChartPlotFillBetween = colourisePlotFillBetween, + deviationChartLayout = id } + +instance View DeviationChartView where + + outputView v = + let reporter exp dir = + do st <- newDeviationChart v exp dir + return Reporter { reporterInitialise = return (), + reporterFinalise = finaliseDeviationChart st, + reporterSimulate = simulateDeviationChart st, + reporterTOCHtml = deviationChartTOCHtml st, + reporterHtml = deviationChartHtml st } + in Generator { generateReporter = reporter } + +-- | The state of the view. +data DeviationChartViewState = + DeviationChartViewState { deviationChartView :: DeviationChartView, + deviationChartExperiment :: Experiment, + deviationChartDir :: FilePath, + deviationChartFile :: IORef (Maybe FilePath), + deviationChartLock :: MVar (), + deviationChartResults :: IORef (Maybe DeviationChartResults) } + +-- | The deviation chart item. +data DeviationChartResults = + DeviationChartResults { deviationChartTimes :: IOArray Int Double, + deviationChartNames :: [Either String String], + deviationChartStats :: [IOArray Int (SamplingStats Double)] } + +-- | Create a new state of the view. +newDeviationChart :: DeviationChartView -> Experiment -> FilePath -> IO DeviationChartViewState +newDeviationChart view exp dir = + do let specs = experimentSpecs exp + bnds = integIterationBnds specs + f <- newIORef Nothing + l <- newMVar () + r <- newIORef Nothing + return DeviationChartViewState { deviationChartView = view, + deviationChartExperiment = exp, + deviationChartDir = dir, + deviationChartFile = f, + deviationChartLock = l, + deviationChartResults = r } + +-- | Create new chart results. +newDeviationChartResults :: [Either String String] -> Experiment -> IO DeviationChartResults +newDeviationChartResults names exp = + do let specs = experimentSpecs exp + bnds = integIterationBnds specs + times <- liftIO $ newListArray bnds (integTimes specs) + stats <- forM names $ \_ -> + liftIO $ newArray bnds emptySamplingStats + return DeviationChartResults { deviationChartTimes = times, + deviationChartNames = names, + deviationChartStats = stats } + +-- | Plot the time series chart during the simulation. +simulateDeviationChart :: DeviationChartViewState -> ExperimentData -> Dynamics (Dynamics ()) +simulateDeviationChart st expdata = + do let protolabels = deviationChartSeries $ deviationChartView st + protoproviders = flip map protolabels $ \protolabel -> + case protolabel of + Left label -> map Left $ experimentSeriesProviders expdata [label] + Right label -> map Right $ experimentSeriesProviders expdata [label] + joinedproviders = join protoproviders + providers = flip map joinedproviders $ either id id + input = + flip map providers $ \provider -> + case providerToDouble provider of + Nothing -> error $ + "Cannot represent series " ++ + providerName provider ++ + " as double values: simulateDeviationChart" + Just input -> input + names = flip map joinedproviders $ \protoprovider -> + case protoprovider of + Left provider -> Left $ providerName provider + Right provider -> Right $ providerName provider + predicate = deviationChartPredicate $ deviationChartView st + exp = deviationChartExperiment st + lock = deviationChartLock st + results <- liftIO $ readIORef (deviationChartResults st) + case results of + Nothing -> + liftIO $ + do results <- newDeviationChartResults names exp + writeIORef (deviationChartResults st) $ Just results + Just results -> + when (names /= deviationChartNames results) $ + error "Series with different names are returned for different runs: simulateDeviationChart" + results <- liftIO $ fmap fromJust $ readIORef (deviationChartResults st) + let stats = deviationChartStats results + t0 <- starttime + enqueue (experimentQueue expdata) t0 $ + do let h = filterSignalM (const predicate) $ + experimentSignalInIntegTimes expdata + -- we must subscribe through the event queue; + -- otherwise, we will loose a signal in the start time, + -- because the handleSignal_ function checks the event queue + handleSignal_ h $ \_ -> + do xs <- sequence input + i <- integIteration + liftIO $ withMVar lock $ \() -> + forM_ (zip xs stats) $ \(x, stats) -> + do y <- readArray stats i + writeArray stats i $ addSamplingStats x y + return $ return () + +-- | Plot the time series after the simulation is complete. +finaliseDeviationChart :: DeviationChartViewState -> IO () +finaliseDeviationChart st = + do let title = deviationChartTitle $ deviationChartView st + width = deviationChartWidth $ deviationChartView st + height = deviationChartHeight $ deviationChartView st + plotLines = deviationChartPlotLines $ deviationChartView st + plotFillBetween = deviationChartPlotFillBetween $ deviationChartView st + plotLayout = deviationChartLayout $ deviationChartView st + results <- readIORef $ deviationChartResults st + case results of + Nothing -> return () + Just results -> + do let times = deviationChartTimes results + names = deviationChartNames results + stats = deviationChartStats results + ps1 <- forM (zip3 names stats plotLines) $ \(name, stats, plotLines) -> + do xs <- getAssocs stats + zs <- forM xs $ \(i, stats) -> + do t <- readArray times i + return (t, samplingStatsMean stats) + let p = toPlot $ + plotLines $ + plot_lines_values ^= filterPlotLinesValues zs $ + plot_lines_title ^= either id id name $ + defaultPlotLines + case name of + Left _ -> return $ Left p + Right _ -> return $ Right p + ps2 <- forM (zip3 names stats plotFillBetween) $ \(name, stats, plotFillBetween) -> + do xs <- getAssocs stats + zs <- forM xs $ \(i, stats) -> + do t <- readArray times i + let mu = samplingStatsMean stats + sigma = samplingStatsDeviation stats + return (t, (mu - 3 * sigma, mu + 3 * sigma)) + let p = toPlot $ + plotFillBetween $ + plot_fillbetween_values ^= filterPlotFillBetweenValues zs $ + plot_fillbetween_title ^= either id id name $ + defaultPlotFillBetween + case name of + Left _ -> return $ Left p + Right _ -> return $ Right p + let ps = join $ flip map (zip ps1 ps2) $ \(p1, p2) -> [p2, p1] + chart = plotLayout $ + layout1_title ^= title $ + layout1_plots ^= ps $ + defaultLayout1 + file <- resolveFileName + (Just $ deviationChartDir st) + (deviationChartFileName $ deviationChartView st) $ + M.fromList [("$TITLE", title)] + renderableToPNGFile (toRenderable chart) width height file + when (experimentVerbose $ deviationChartExperiment st) $ + putStr "Generated file " >> putStrLn file + writeIORef (deviationChartFile st) $ Just file + +-- | Remove the NaN and inifity values. +filterPlotLinesValues :: [(Double, Double)] -> [[(Double, Double)]] +filterPlotLinesValues = + filter (not . null) . + divideBy (\(t, x) -> isNaN x || isInfinite x) + +-- | Remove the NaN and inifity values. +filterPlotFillBetweenValues :: [(Double, (Double, Double))] -> [(Double, (Double, Double))] +filterPlotFillBetweenValues = + join . filter (not . null) . + divideBy (\(t, (x1, x2)) -> isNaN x1 || isInfinite x1 || isNaN x2 || isInfinite x2) + +-- | Get the HTML code. +deviationChartHtml :: DeviationChartViewState -> Int -> HtmlWriter () +deviationChartHtml st index = + do header st index + file <- liftIO $ readIORef (deviationChartFile st) + case file of + Nothing -> return () + Just f -> + writeHtmlParagraph $ + writeHtmlImage (makeRelative (deviationChartDir st) f) + +header :: DeviationChartViewState -> Int -> HtmlWriter () +header st index = + do writeHtmlHeader3WithId ("id" ++ show index) $ + writeHtmlText (deviationChartTitle $ deviationChartView st) + let description = deviationChartDescription $ deviationChartView st + unless (null description) $ + writeHtmlParagraph $ + writeHtmlText description + +-- | Get the TOC item. +deviationChartTOCHtml :: DeviationChartViewState -> Int -> HtmlWriter () +deviationChartTOCHtml st index = + writeHtmlListItem $ + writeHtmlLink ("#id" ++ show index) $ + writeHtmlText (deviationChartTitle $ deviationChartView st) diff --git a/Simulation/Aivika/Experiment/TimeSeriesView.hs b/Simulation/Aivika/Experiment/TimeSeriesView.hs new file mode 100644 --- /dev/null +++ b/Simulation/Aivika/Experiment/TimeSeriesView.hs @@ -0,0 +1,251 @@ + +-- | +-- Module : Simulation.Aivika.Experiment.TimeSeriesView +-- Copyright : Copyright (c) 2012, David Sorokin +-- License : BSD3 +-- Maintainer : David Sorokin +-- Stability : experimental +-- Tested with: GHC 7.4.1 +-- +-- The module defines 'TimeSeriesView' that saves the time series +-- charts as the PNG files. +-- + +module Simulation.Aivika.Experiment.TimeSeriesView + (TimeSeriesView(..), + defaultTimeSeriesView) where + +import Control.Monad +import Control.Monad.Trans + +import qualified Data.Map as M +import Data.IORef +import Data.Maybe +import Data.Either +import Data.Array + +import Data.Accessor + +import System.IO +import System.FilePath + +import Data.String.Utils (replace) + +import Graphics.Rendering.Chart + +import Simulation.Aivika.Experiment +import Simulation.Aivika.Experiment.HtmlWriter +import Simulation.Aivika.Experiment.Utils (divideBy) +import Simulation.Aivika.Experiment.Chart (colourisePlotLines) + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.Signal +import Simulation.Aivika.Dynamics.EventQueue + +-- | Defines the 'View' that saves the time series charts +-- in the PNG files. +data TimeSeriesView = + TimeSeriesView { timeSeriesTitle :: String, + -- ^ This is a title used in HTML and chart. + timeSeriesRunTitle :: String, + -- ^ The run title for the view. It is used + -- when simulating multiple runs and it may + -- include special variables @$RUN_INDEX@, + -- @$RUN_COUNT@ and @$TITLE@. + -- + -- An example is + -- + -- @ + -- timeSeriesRunTitle = \"$TITLE / Run $RUN_INDEX of $RUN_COUNT\" + -- @ + timeSeriesDescription :: String, + -- ^ This is a description in the HTML. + timeSeriesWidth :: Int, + -- ^ The width of the chart. + timeSeriesHeight :: Int, + -- ^ The height of the chart. + timeSeriesFileName :: FileName, + -- ^ It defines the file name for each PNG file. + -- It may include special variables @$TITLE@, + -- @$RUN_INDEX@ and @$RUN_COUNT@. + -- + -- An example is + -- + -- @ + -- timeSeriesFileName = UniqueFileName \"$TITLE - $RUN_INDEX\", \".png\" + -- @ + timeSeriesPredicate :: Dynamics Bool, + -- ^ It specifies the predicate that defines + -- when we plot data in the chart. + timeSeries :: [Either String String], + -- ^ It contains the labels of data plotted + -- in the chart. + timeSeriesPlotLines :: [PlotLines Double Double -> + PlotLines Double Double], + -- ^ Probably, an infinite sequence of plot + -- transformations based on which the plot + -- is constructed for each series. Generally, + -- it must not coincide with a sequence of + -- labels as one label may denote a whole list + -- or an array of data providers. + -- + -- Here you can define a colour or style of + -- of the plot lines. + timeSeriesLayout :: Layout1 Double Double -> + Layout1 Double Double + -- ^ A transformation of the plot layout, + -- where you can redefine the axes, for example. + } + +-- | The default time series view. +defaultTimeSeriesView :: TimeSeriesView +defaultTimeSeriesView = + TimeSeriesView { timeSeriesTitle = "Time Series", + timeSeriesRunTitle = "$TITLE / Run $RUN_INDEX of $RUN_COUNT", + timeSeriesDescription = [], + timeSeriesWidth = 640, + timeSeriesHeight = 480, + timeSeriesFileName = UniqueFileName "$TITLE - $RUN_INDEX" ".png", + timeSeriesPredicate = return True, + timeSeries = [], + timeSeriesPlotLines = colourisePlotLines, + timeSeriesLayout = id } + +instance View TimeSeriesView where + + outputView v = + let reporter exp dir = + do st <- newTimeSeries v exp dir + return Reporter { reporterInitialise = return (), + reporterFinalise = return (), + reporterSimulate = simulateTimeSeries st, + reporterTOCHtml = timeSeriesTOCHtml st, + reporterHtml = timeSeriesHtml st } + in Generator { generateReporter = reporter } + +-- | The state of the view. +data TimeSeriesViewState = + TimeSeriesViewState { timeSeriesView :: TimeSeriesView, + timeSeriesExperiment :: Experiment, + timeSeriesDir :: FilePath, + timeSeriesMap :: M.Map Int FilePath } + +-- | Create a new state of the view. +newTimeSeries :: TimeSeriesView -> Experiment -> FilePath -> IO TimeSeriesViewState +newTimeSeries view exp dir = + do let n = experimentRunCount exp + fs <- forM [0..(n - 1)] $ \i -> + resolveFileName (Just dir) (timeSeriesFileName view) $ + M.fromList [("$TITLE", timeSeriesTitle view), + ("$RUN_INDEX", show $ i + 1), + ("$RUN_COUNT", show n)] + let m = M.fromList $ zip [0..(n - 1)] fs + return TimeSeriesViewState { timeSeriesView = view, + timeSeriesExperiment = exp, + timeSeriesDir = dir, + timeSeriesMap = m } + +-- | Plot the time series chart during the simulation. +simulateTimeSeries :: TimeSeriesViewState -> ExperimentData -> Dynamics (Dynamics ()) +simulateTimeSeries st expdata = + do let protolabels = timeSeries $ timeSeriesView st + labels = flip map protolabels $ either id id + providers = experimentSeriesProviders expdata labels + input = + flip map providers $ \provider -> + case providerToDouble provider of + Nothing -> error $ + "Cannot represent series " ++ + providerName provider ++ + " as double values: simulateTimeSeries" + Just input -> input + n = experimentRunCount $ timeSeriesExperiment st + width = timeSeriesWidth $ timeSeriesView st + height = timeSeriesHeight $ timeSeriesView st + predicate = timeSeriesPredicate $ timeSeriesView st + plotLines = timeSeriesPlotLines $ timeSeriesView st + plotLayout = timeSeriesLayout $ timeSeriesView st + i <- liftSimulation simulationIndex + let file = fromJust $ M.lookup (i - 1) (timeSeriesMap st) + title = + if n == 1 + then timeSeriesTitle $ timeSeriesView st + else replace "$RUN_INDEX" (show i) $ + replace "$RUN_COUNT" (show n) $ + replace "$TITLE" (timeSeriesTitle $ timeSeriesView st) + (timeSeriesRunTitle $ timeSeriesView st) + hs <- forM (zip providers input) $ \(provider, input) -> + newSignalHistoryThrough (experimentQueue expdata) $ + mapSignalM (const input) $ + filterSignalM (const predicate) $ + experimentMixedSignal expdata [provider] + return $ + do ps <- forM (zip3 hs providers plotLines) $ \(h, provider, plotLines) -> + do (ts, xs) <- readSignalHistory h + return $ + toPlot $ + plotLines $ + plot_lines_values ^= filterPlotLinesValues (zip (elems ts) (elems xs)) $ + plot_lines_title ^= providerName provider $ + defaultPlotLines + let ps' = flip map (zip ps protolabels) $ \(p, label) -> + case label of + Left _ -> Left p + Right _ -> Right p + let chart = plotLayout $ + layout1_title ^= title $ + layout1_plots ^= ps' $ + defaultLayout1 + liftIO $ + do renderableToPNGFile (toRenderable chart) width height file + when (experimentVerbose $ timeSeriesExperiment st) $ + putStr "Generated file " >> putStrLn file + +-- | Remove the NaN and inifity values. +filterPlotLinesValues :: [(Double, Double)] -> [[(Double, Double)]] +filterPlotLinesValues = + filter (not . null) . + divideBy (\(t, x) -> isNaN x || isInfinite x) + +-- | Get the HTML code. +timeSeriesHtml :: TimeSeriesViewState -> Int -> HtmlWriter () +timeSeriesHtml st index = + let n = experimentRunCount $ timeSeriesExperiment st + in if n == 1 + then timeSeriesHtmlSingle st index + else timeSeriesHtmlMultiple st index + +-- | Get the HTML code for a single run. +timeSeriesHtmlSingle :: TimeSeriesViewState -> Int -> HtmlWriter () +timeSeriesHtmlSingle st index = + do header st index + let f = fromJust $ M.lookup 0 (timeSeriesMap st) + writeHtmlParagraph $ + writeHtmlImage (makeRelative (timeSeriesDir st) f) + +-- | Get the HTML code for multiple runs. +timeSeriesHtmlMultiple :: TimeSeriesViewState -> Int -> HtmlWriter () +timeSeriesHtmlMultiple st index = + do header st index + let n = experimentRunCount $ timeSeriesExperiment st + forM_ [0..(n - 1)] $ \i -> + let f = fromJust $ M.lookup i (timeSeriesMap st) + in writeHtmlParagraph $ + writeHtmlImage (makeRelative (timeSeriesDir st) f) + +header :: TimeSeriesViewState -> Int -> HtmlWriter () +header st index = + do writeHtmlHeader3WithId ("id" ++ show index) $ + writeHtmlText (timeSeriesTitle $ timeSeriesView st) + let description = timeSeriesDescription $ timeSeriesView st + unless (null description) $ + writeHtmlParagraph $ + writeHtmlText description + +-- | Get the TOC item. +timeSeriesTOCHtml :: TimeSeriesViewState -> Int -> HtmlWriter () +timeSeriesTOCHtml st index = + writeHtmlListItem $ + writeHtmlLink ("#id" ++ show index) $ + writeHtmlText (timeSeriesTitle $ timeSeriesView st) diff --git a/aivika-experiment-chart.cabal b/aivika-experiment-chart.cabal new file mode 100644 --- /dev/null +++ b/aivika-experiment-chart.cabal @@ -0,0 +1,45 @@ +name: aivika-experiment-chart +version: 0.1 +synopsis: Simulation experiments with charting for the Aivika library +description: + This package complements the Aivika and Aivika Experiment packages with + charting capabilites. Now the simulation results can be represented + as charts. + . + It was intentionally made a separate package as it has heavy dependencies + on Haskell Charts, Cairo and GTK. + . +category: Simulation +license: BSD3 +license-file: LICENSE +copyright: (c) 2012. David Sorokin +author: David Sorokin +maintainer: David Sorokin +homepage: http://github.com/dsorokin/aivika-experiment-chart +cabal-version: >= 1.2.0 +build-type: Simple +tested-with: GHC == 7.4.1 + +extra-source-files: examples/MachRep3.hs + examples/ChemicalReaction.hs + examples/BassDiffusion.hs + +library + + exposed-modules: Simulation.Aivika.Experiment.TimeSeriesView + Simulation.Aivika.Experiment.DeviationChartView + Simulation.Aivika.Experiment.Chart + + build-depends: base >= 3 && < 6, + mtl >= 1.1.0.2, + array >= 0.3.0.0, + containers >= 0.4.0.0, + filepath >= 1.3.0.0, + Chart >= 0.16, + MissingH >= 1.2.0.0, + data-accessor >= 0.2.2.3, + colour >= 2.3.3, + aivika >= 0.4.3, + aivika-experiment >= 0.1 + + ghc-options: -O2 diff --git a/examples/BassDiffusion.hs b/examples/BassDiffusion.hs new file mode 100644 --- /dev/null +++ b/examples/BassDiffusion.hs @@ -0,0 +1,122 @@ + +import System.Random +import Data.Array +import Control.Monad +import Control.Monad.Trans + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.EventQueue +import Simulation.Aivika.Dynamics.Agent +import Simulation.Aivika.Dynamics.Ref + +import Simulation.Aivika.Experiment +import Simulation.Aivika.Experiment.DeviationChartView + +n = 100 -- the number of agents + +advertisingEffectiveness = 0.011 +contactRate = 100.0 +adoptionFraction = 0.015 + +specs = Specs { spcStartTime = 0.0, + spcStopTime = 8.0, + spcDT = 0.1, + spcMethod = RungeKutta4 } + +experiment :: Experiment +experiment = + defaultExperiment { + experimentSpecs = specs, + experimentRunCount = 20, + experimentDescription = "This is a famous Bass Diffusion model solved with help of the agent-based modelling.", + experimentGenerators = + [outputView $ defaultDeviationChartView { + deviationChartSeries = [Left "potentialAdopters", + Left "adopters"] } ] } + +exprnd :: Double -> IO Double +exprnd lambda = + do x <- getStdRandom random + return (- log x / lambda) + +boolrnd :: Double -> IO Bool +boolrnd p = + do x <- getStdRandom random + return (x <= p) + +data Person = Person { personAgent :: Agent, + personPotentialAdopter :: AgentState, + personAdopter :: AgentState } + +createPerson :: EventQueue -> Simulation Person +createPerson q = + do agent <- newAgent q + potentialAdopter <- newState agent + adopter <- newState agent + return Person { personAgent = agent, + personPotentialAdopter = potentialAdopter, + personAdopter = adopter } + +createPersons :: EventQueue -> Simulation (Array Int Person) +createPersons q = + do list <- forM [1 .. n] $ \i -> + do p <- createPerson q + return (i, p) + return $ array (1, n) list + +definePerson :: Person -> Array Int Person -> Ref Int -> Ref Int -> Simulation () +definePerson p ps potentialAdopters adopters = + do stateActivation (personPotentialAdopter p) $ + do modifyRef potentialAdopters $ \a -> a + 1 + -- add a timeout + t <- liftIO $ exprnd advertisingEffectiveness + let st = personPotentialAdopter p + st' = personAdopter p + addTimeout st t $ activateState st' + stateActivation (personAdopter p) $ + do modifyRef adopters $ \a -> a + 1 + -- add a timer that works while the state is active + let t = liftIO $ exprnd contactRate -- many times! + addTimer (personAdopter p) t $ + do i <- liftIO $ getStdRandom $ randomR (1, n) + let p' = ps ! i + st <- agentState (personAgent p') + when (st == Just (personPotentialAdopter p')) $ + do b <- liftIO $ boolrnd adoptionFraction + when b $ activateState (personAdopter p') + stateDeactivation (personPotentialAdopter p) $ + modifyRef potentialAdopters $ \a -> a - 1 + stateDeactivation (personAdopter p) $ + modifyRef adopters $ \a -> a - 1 + +definePersons :: Array Int Person -> Ref Int -> Ref Int -> Simulation () +definePersons ps potentialAdopters adopters = + forM_ (elems ps) $ \p -> + definePerson p ps potentialAdopters adopters + +activatePerson :: Person -> Dynamics () +activatePerson p = activateState (personPotentialAdopter p) + +activatePersons :: Array Int Person -> Dynamics () +activatePersons ps = + forM_ (elems ps) $ \p -> activatePerson p + +model :: Simulation ExperimentData +model = + do q <- newQueue + potentialAdopters <- newRef q 0 + adopters <- newRef q 0 + ps <- createPersons q + definePersons ps potentialAdopters adopters + runDynamicsInStartTime $ + activatePersons ps + experimentDataInStartTime q $ + [("potentialAdopters", + seriesEntity "Potential Adopters" + potentialAdopters), + ("adopters", + seriesEntity "Adopters" + adopters)] + +main = runExperiment experiment model diff --git a/examples/ChemicalReaction.hs b/examples/ChemicalReaction.hs new file mode 100644 --- /dev/null +++ b/examples/ChemicalReaction.hs @@ -0,0 +1,53 @@ + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.SystemDynamics +import Simulation.Aivika.Dynamics.EventQueue +import Simulation.Aivika.Dynamics.Base + +import Simulation.Aivika.Experiment +import Simulation.Aivika.Experiment.LastValueView +import Simulation.Aivika.Experiment.TableView +import Simulation.Aivika.Experiment.TimeSeriesView + +specs = Specs { spcStartTime = 0, + spcStopTime = 13, + spcDT = 0.01, + spcMethod = RungeKutta4 } + +experiment :: Experiment +experiment = + defaultExperiment { + experimentSpecs = specs, + experimentDescription = "Experiment Description", + experimentGenerators = + [outputView $ defaultLastValueView { + lastValueDescription = "Last Value description", + lastValueSeries = ["t", "a", "b", "c"] }, + outputView $ defaultTableView { + tableDescription = "Table description", + tableSeries = ["t", "a", "b", "c"] }, + outputView $ defaultTimeSeriesView { + timeSeries = [Left "a", Left "b", Left "c"] } ] } + +model :: Simulation ExperimentData +model = + do queue <- newQueue + integA <- newInteg 100 + integB <- newInteg 0 + integC <- newInteg 0 + let a = integValue integA + b = integValue integB + c = integValue integC + let ka = 1 + kb = 1 + integDiff integA (- ka * a) + integDiff integB (ka * a - kb * b) + integDiff integC (kb * b) + experimentDataInStartTime queue $ + [("t", seriesEntity "time" time), + ("a", seriesEntity "a" a), + ("b", seriesEntity "b" b), + ("c", seriesEntity "c" c)] + +main = runExperiment experiment model diff --git a/examples/MachRep3.hs b/examples/MachRep3.hs new file mode 100644 --- /dev/null +++ b/examples/MachRep3.hs @@ -0,0 +1,122 @@ + +-- It corresponds to model MachRep3 described in document +-- Introduction to Discrete-Event Simulation and the SimPy Language +-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. +-- SimPy is available on [http://simpy.sourceforge.net/]. +-- +-- The model description is as follows. +-- +-- Variation of models MachRep1, MachRep2. Two machines, but +-- sometimes break down. Up time is exponentially distributed with mean +-- 1.0, and repair time is exponentially distributed with mean 0.5. In +-- this example, there is only one repairperson, and she is not summoned +-- until both machines are down. We find the proportion of up time. It +-- should come out to about 0.45. + +module MachRep3Model (model) where + +import System.Random +import Control.Monad +import Control.Monad.Trans + +import Simulation.Aivika.Dynamics +import Simulation.Aivika.Dynamics.Simulation +import Simulation.Aivika.Dynamics.Base +import Simulation.Aivika.Dynamics.EventQueue +import Simulation.Aivika.Dynamics.Ref +import Simulation.Aivika.Dynamics.Resource +import Simulation.Aivika.Dynamics.Process + +import Simulation.Aivika.Experiment +import Simulation.Aivika.Experiment.LastValueView +import Simulation.Aivika.Experiment.TableView +import Simulation.Aivika.Experiment.TimeSeriesView +import Simulation.Aivika.Experiment.DeviationChartView + +specs = Specs { spcStartTime = 0.0, + spcStopTime = 1000.0, + spcDT = 1.0, + spcMethod = RungeKutta4 } + +experiment :: Experiment +experiment = + defaultExperiment { + experimentSpecs = specs, + experimentRunCount = 3, + experimentDescription = "Experiment Description", + experimentGenerators = + [outputView $ defaultLastValueView { + lastValueDescription = "Last Value description", + lastValueSeries = ["x"] }, + outputView $ defaultTableView { + tableDescription = "Table description", + tableSeries = ["x"] }, + outputView $ defaultTimeSeriesView { + timeSeries = [Left "t", Right "x"] }, + outputView $ defaultDeviationChartView { + deviationChartSeries = [Left "t", Right "x"] } ] } + +upRate = 1.0 / 1.0 -- reciprocal of mean up time +repairRate = 1.0 / 0.5 -- reciprocal of mean repair time + +exprnd :: Double -> IO Double +exprnd lambda = + do x <- getStdRandom random + return (- log x / lambda) + +model :: Simulation ExperimentData +model = + do queue <- newQueue + + -- number of machines currently up + nUp <- newRef queue 2 + + -- total up time for all machines + totalUpTime <- newRef queue 0.0 + + repairPerson <- newResource queue 1 + + pid1 <- newProcessID queue + pid2 <- newProcessID queue + + let machine :: ProcessID -> Process () + machine pid = + do startUpTime <- liftDynamics time + upTime <- liftIO $ exprnd upRate + holdProcess upTime + finishUpTime <- liftDynamics time + liftDynamics $ modifyRef totalUpTime + (+ (finishUpTime - startUpTime)) + + liftDynamics $ modifyRef nUp $ \a -> a - 1 + nUp' <- liftDynamics $ readRef nUp + if nUp' == 1 + then passivateProcess + else do n <- liftDynamics $ + resourceCount repairPerson + when (n == 1) $ + liftDynamics $ reactivateProcess pid + + requestResource repairPerson + repairTime <- liftIO $ exprnd repairRate + holdProcess repairTime + liftDynamics $ modifyRef nUp $ \a -> a + 1 + releaseResource repairPerson + + machine pid + + runDynamicsInStartTime $ + do t0 <- starttime + runProcess (machine pid2) pid1 t0 + runProcess (machine pid1) pid2 t0 + + let result = + do x <- readRef totalUpTime + y <- time + return $ x / (2 * y) + + experimentDataInStartTime queue $ + [("x", seriesEntity "The proportion of up time" result), + ("t", seriesEntity "Total up time" totalUpTime)] + +main = runExperiment experiment model