aivika-experiment-chart 0.2.1 → 0.2.2
raw patch · 11 files changed
+935/−59 lines, 11 filesdep +splitdep −MissingHdep ~aivikadep ~aivika-experimentPVP ok
version bump matches the API change (PVP)
Dependencies added: split
Dependencies removed: MissingH
Dependency ranges changed: aivika, aivika-experiment
API changes (from Hackage documentation)
Files
- Simulation/Aivika/Experiment/DeviationChartView.hs +8/−9
- Simulation/Aivika/Experiment/FinalHistogramView.hs +9/−10
- Simulation/Aivika/Experiment/FinalXYChartView.hs +1/−3
- Simulation/Aivika/Experiment/HistogramView.hs +6/−7
- Simulation/Aivika/Experiment/TimeSeriesView.hs +1/−3
- Simulation/Aivika/Experiment/XYChartView.hs +1/−3
- aivika-experiment-chart.cabal +9/−6
- examples/ChemicalReaction.hs +24/−18
- examples/DifferenceEquations.hs +71/−0
- examples/Financial.hs +246/−0
- examples/Furnace.hs +559/−0
Simulation/Aivika/Experiment/DeviationChartView.hs view
@@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Experiment.DeviationChartView--- Copyright : Copyright (c) 2012, David Sorokin <david.sorokin@gmail.com>+-- Copyright : Copyright (c) 2012-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 'DeviationChartView' that saves the deviation -- chart in the PNG file.@@ -31,14 +31,13 @@ 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.Utils (divideBy, replace) import Simulation.Aivika.Experiment.Chart (colourisePlotLines, colourisePlotFillBetween)+import Simulation.Aivika.Experiment.SamplingStatsSource import Simulation.Aivika.Dynamics import Simulation.Aivika.Dynamics.Simulation@@ -186,12 +185,12 @@ providers = flip map joinedproviders $ either id id input = flip map providers $ \provider ->- case providerToDouble provider of+ case providerToDoubleStatsSource provider of Nothing -> error $ "Cannot represent series " ++ providerName provider ++ - " as double values: simulateDeviationChart"- Just input -> input+ " as a series of double values: simulateDeviationChart"+ Just input -> samplingStatsSourceData input names = flip map joinedproviders $ \protoprovider -> case protoprovider of Left provider -> Left $ providerName provider@@ -221,7 +220,7 @@ liftIO $ withMVar lock $ \() -> forM_ (zip xs stats) $ \(x, stats) -> do y <- readArray stats i- let y' = addSamplingStats x y+ let y' = addDataToSamplingStats x y y' `seq` writeArray stats i y' return $ return ()
Simulation/Aivika/Experiment/FinalHistogramView.hs view
@@ -32,15 +32,14 @@ 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.Utils (divideBy, replace) import Simulation.Aivika.Experiment.Chart (colourisePlotBars) import Simulation.Aivika.Experiment.Histogram+import Simulation.Aivika.Experiment.ListSource import Simulation.Aivika.Dynamics import Simulation.Aivika.Dynamics.Simulation@@ -140,7 +139,7 @@ -- | The histogram item. data FinalHistogramResults = FinalHistogramResults { finalHistogramNames :: [String],- finalHistogramValues :: [IORef [Double]] }+ finalHistogramValues :: [ListRef Double] } -- | Create a new state of the view. newFinalHistogram :: FinalHistogramView -> Experiment -> FilePath -> IO FinalHistogramViewState@@ -158,7 +157,7 @@ -- | Create new histogram results. newFinalHistogramResults :: [String] -> Experiment -> IO FinalHistogramResults newFinalHistogramResults names exp =- do values <- forM names $ \_ -> liftIO $ newIORef []+ do values <- forM names $ \_ -> liftIO newListRef return FinalHistogramResults { finalHistogramNames = names, finalHistogramValues = values } @@ -171,12 +170,12 @@ providers = concat protoproviders input = flip map providers $ \provider ->- case providerToDouble provider of+ case providerToDoubleListSource provider of Nothing -> error $ "Cannot represent series " ++ providerName provider ++ - " as double values: simulateFinalHistogram"- Just input -> input+ " as a source of double values: simulateFinalHistogram"+ Just input -> listSourceData input names = map providerName providers predicate = finalHistogramPredicate $ finalHistogramView st exp = finalHistogramExperiment st@@ -203,7 +202,7 @@ do xs <- sequence input liftIO $ withMVar lock $ \() -> forM_ (zip xs values) $ \(x, values) ->- x `seq` modifyIORef values (x :)+ addDataToListRef values x return $ return () -- | Plot the histogram after the simulation is complete.@@ -224,7 +223,7 @@ Just results -> do let names = finalHistogramNames results values = finalHistogramValues results- xs <- forM values readIORef+ xs <- forM values readListRef let zs = histogramToBars . filterHistogram . histogram $ map filterData xs p = plotBars $
Simulation/Aivika/Experiment/FinalXYChartView.hs view
@@ -31,13 +31,11 @@ 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.Utils (divideBy, replace) import Simulation.Aivika.Experiment.Chart (colourisePlotLines) import Simulation.Aivika.Dynamics
Simulation/Aivika/Experiment/HistogramView.hs view
@@ -30,15 +30,14 @@ 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.Utils (divideBy, replace) import Simulation.Aivika.Experiment.Chart (colourisePlotBars) import Simulation.Aivika.Experiment.Histogram+import Simulation.Aivika.Experiment.ListSource import Simulation.Aivika.Dynamics import Simulation.Aivika.Dynamics.Simulation@@ -169,12 +168,12 @@ names = map providerName providers input = flip map providers $ \provider ->- case providerToDouble provider of+ case providerToDoubleListSource provider of Nothing -> error $ "Cannot represent series " ++ providerName provider ++ - " as double values: simulateHistogram"- Just input -> input+ " as a source of double values: simulateHistogram"+ Just input -> fmap listDataList $ listSourceData input n = experimentRunCount $ histogramExperiment st width = histogramWidth $ histogramView st height = histogramHeight $ histogramView st@@ -203,7 +202,7 @@ return $ do xs <- forM hs readSignalHistory let zs = histogramToBars . filterHistogram . build $ - map (filterData . elems . snd) xs+ map (filterData . concat . elems . snd) xs p = plotBars $ bars $ plot_bars_values ^= zs $
Simulation/Aivika/Experiment/TimeSeriesView.hs view
@@ -29,13 +29,11 @@ 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.Utils (divideBy, replace) import Simulation.Aivika.Experiment.Chart (colourisePlotLines) import Simulation.Aivika.Dynamics
Simulation/Aivika/Experiment/XYChartView.hs view
@@ -30,13 +30,11 @@ 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.Utils (divideBy, replace) import Simulation.Aivika.Experiment.Chart (colourisePlotLines) import Simulation.Aivika.Dynamics
aivika-experiment-chart.cabal view
@@ -1,5 +1,5 @@ name: aivika-experiment-chart-version: 0.2.1+version: 0.2.2 synopsis: Simulation experiments with charting for the Aivika library description: This package complements the Aivika and Aivika Experiment packages with@@ -12,17 +12,20 @@ category: Simulation license: BSD3 license-file: LICENSE-copyright: (c) 2012. David Sorokin <david.sorokin@gmail.com>+copyright: (c) 2012-2013. David Sorokin <david.sorokin@gmail.com> author: David Sorokin maintainer: David Sorokin <david.sorokin@gmail.com> homepage: http://github.com/dsorokin/aivika-experiment-chart cabal-version: >= 1.2.0 build-type: Simple-tested-with: GHC == 7.4.1+tested-with: GHC == 7.6.3 extra-source-files: examples/MachRep3.hs examples/ChemicalReaction.hs examples/BassDiffusion.hs+ examples/Financial.hs+ examples/DifferenceEquations.hs+ examples/Furnace.hs library @@ -40,10 +43,10 @@ containers >= 0.4.0.0, filepath >= 1.3.0.0, Chart >= 0.16,- MissingH >= 1.2.0.0,+ split >= 0.2.2, data-accessor >= 0.2.2.3, colour >= 2.3.3,- aivika >= 0.5.1,- aivika-experiment >= 0.2.1+ aivika >= 0.5.4,+ aivika-experiment >= 0.2.2 ghc-options: -O2
examples/ChemicalReaction.hs view
@@ -1,4 +1,6 @@ +{-# LANGUAGE RecursiveDo #-}+ import Simulation.Aivika.Dynamics import Simulation.Aivika.Dynamics.Simulation import Simulation.Aivika.Dynamics.SystemDynamics@@ -22,13 +24,17 @@ defaultExperiment { experimentSpecs = specs, experimentRunCount = 1,- experimentGenerators =+ experimentTitle = "Chemical Reaction",+ experimentDescription = "Chemical Reaction as described in " +++ "the 5-minute tutorial of Berkeley-Madonna",+ experimentGenerators = [outputView defaultExperimentSpecsView, outputView $ defaultLastValueView { lastValueSeries = ["t", "a", "b", "c"] }, outputView $ defaultTableView { tableSeries = ["t", "a", "b", "c"] }, outputView $ defaultTimeSeriesView {+ timeSeriesTitle = "Time Series", timeSeries = [Left "a", Left "b", Left "c"] }, -- outputView $ defaultTimeSeriesView { -- timeSeriesPlotTitle = "Variables a, b and c for t <= 5 or t >= 7",@@ -37,12 +43,18 @@ -- do t <- time -- return (t <= 5 || t >= 7) }, outputView $ defaultXYChartView {+ xyChartTitle = "XYChart - 1",+ xyChartPlotTitle = "b=b(a), c=c(a)", xyChartXSeries = Just "a", xyChartYSeries = [Left "b", Right "c"] }, outputView $ defaultXYChartView {+ xyChartTitle = "XYChart - 2",+ xyChartPlotTitle = "a=a(b), c=c(b)", xyChartXSeries = Just "b", xyChartYSeries = [Right "a", Right "c"] }, outputView $ defaultXYChartView {+ xyChartTitle = "XYChart - 3",+ xyChartPlotTitle = "a=a(c), b=b(c)", xyChartXSeries = Just "c", xyChartYSeries = [Right "a", Left "b"] } ] } -- outputView $ defaultXYChartView {@@ -55,22 +67,16 @@ 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)]+ mdo queue <- newQueue+ a <- integ (- ka * a) 100+ b <- integ (ka * a - kb * b) 0+ c <- integ (kb * b) 0+ let ka = 1+ kb = 1+ experimentDataInStartTime queue+ [("t", seriesEntity "time" time),+ ("a", seriesEntity "a" a),+ ("b", seriesEntity "b" b),+ ("c", seriesEntity "c" c)] main = runExperiment experiment model
+ examples/DifferenceEquations.hs view
@@ -0,0 +1,71 @@++{-# LANGUAGE RecursiveDo #-}++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.Dynamics.Random++import Simulation.Aivika.Experiment+import Simulation.Aivika.Experiment.TableView+import Simulation.Aivika.Experiment.TimeSeriesView+import Simulation.Aivika.Experiment.ExperimentSpecsView+import Simulation.Aivika.Experiment.TimingStatsView++specs = Specs { spcStartTime = 0, + spcStopTime = 10000, + spcDT = 1,+ spcMethod = RungeKutta4 }++experiment :: Experiment+experiment =+ defaultExperiment {+ experimentSpecs = specs,+ experimentRunCount = 1,+ experimentTitle = "Difference Equations",+ experimentDescription = "Difference Equations as described in " +++ "the corresponded tutorial of Berkeley-Madonna " +++ "with small modification for calculating std.",+ experimentGenerators = + [outputView defaultExperimentSpecsView,+ outputView $ defaultTableView {+ tableSeries = ["t", "x", "sumX", "sumX2", "avg", "std"] }, + outputView $ defaultTimeSeriesView {+ timeSeriesTitle = "Time Series",+ timeSeries = [Left "x", Left "avg"] },+ outputView $ defaultTimingStatsView {+ timingStatsSeries = ["x"] },+ outputView $ defaultTimeSeriesView {+ timeSeriesTitle = "Sums",+ timeSeries = [Left "sumX", Right "sumX2"] },+ outputView $ defaultTimeSeriesView {+ timeSeriesTitle = "Standard Deviation",+ timeSeries = [Left "std"] } ] }++model :: Simulation ExperimentData+model =+ mdo queue <- newQueue+ + x <- newNormal 3 0.8+ sumX <- sumDynamics x 0+ sumX2 <- sumDynamics (x * x) 0+ + -- it would be much more efficient to say:+ -- let n = fmap fromIntegral integIteration+ n <- sumDynamics 1 0++ let avg = ifDynamics (n .>. 0) (sumX / n) 0+ let std = ifDynamics (n .>. 1) (sqrt ((sumX2 - sumX * avg) / (n - 1))) 0+ + experimentDataInStartTime queue+ [("t", seriesEntity "time" time),+ ("n", seriesEntity "n" n),+ ("x", seriesEntity "x" x),+ ("sumX", seriesEntity "sumX" sumX),+ ("sumX2", seriesEntity "sumX2" sumX2),+ ("avg", seriesEntity "avg" avg),+ ("std", seriesEntity "std" std)]++main = runExperiment experiment model
+ examples/Financial.hs view
@@ -0,0 +1,246 @@++{-# LANGUAGE RecursiveDo #-}++-- This financial model is described in+-- Vensim 5 Modeling Guide, Chapter Financial Modeling and Risk.+--+-- It illustrates how you can use the Monte-Carlo simulation+-- and define external parameters. Here the system of recursive+-- diffential equations is used but the paradigm can be any+-- supported by Aivika including DES or agent-base modeling+-- or their combination.+--+-- To enable the parallel simulation, you should compile it+-- with option -threaded and then pass in other options +RTS -N2 -RTS+-- to the executable if you have a dual core processor without+-- hyper-threading. Also you can increase the number+-- of parallel threads via option -N if you have a more modern+-- processor.++import Control.Monad++-- from package aivika+import Simulation.Aivika.Dynamics+import Simulation.Aivika.Dynamics.Simulation+import Simulation.Aivika.Dynamics.Base+import Simulation.Aivika.Dynamics.SystemDynamics+import Simulation.Aivika.Dynamics.Parameter+import Simulation.Aivika.Dynamics.EventQueue++-- from package aivika-experiment+import Simulation.Aivika.Experiment+import Simulation.Aivika.Experiment.ExperimentSpecsView+import Simulation.Aivika.Experiment.TableView+import Simulation.Aivika.Experiment.FinalStatsView++-- from package aivika-experiment-chart+import Simulation.Aivika.Experiment.DeviationChartView+import Simulation.Aivika.Experiment.FinalHistogramView+import Simulation.Aivika.Experiment.TimeSeriesView++-- | The model parameters.+data Parameters =+ Parameters { paramsTaxDepreciationTime :: Simulation Double,+ paramsTaxRate :: Simulation Double,+ paramsAveragePayableDelay :: Simulation Double,+ paramsBillingProcessingTime :: Simulation Double,+ paramsBuildingTime :: Simulation Double,+ paramsDebtFinancingFraction :: Simulation Double,+ paramsDebtRetirementTime :: Simulation Double,+ paramsDiscountRate :: Simulation Double,+ paramsFractionalLossRate :: Simulation Double,+ paramsInterestRate :: Simulation Double,+ paramsPrice :: Simulation Double,+ paramsProductionCapacity :: Simulation Double,+ paramsRequiredInvestment :: Simulation Double,+ paramsVariableProductionCost :: Simulation Double }++-- | The default model parameters.+defaultParams :: Parameters+defaultParams =+ Parameters { paramsTaxDepreciationTime = 10,+ paramsTaxRate = 0.4,+ paramsAveragePayableDelay = 0.09,+ paramsBillingProcessingTime = 0.04,+ paramsBuildingTime = 1,+ paramsDebtFinancingFraction = 0.6,+ paramsDebtRetirementTime = 3,+ paramsDiscountRate = 0.12,+ paramsFractionalLossRate = 0.06,+ paramsInterestRate = 0.12,+ paramsPrice = 1,+ paramsProductionCapacity = 2400,+ paramsRequiredInvestment = 2000,+ paramsVariableProductionCost = 0.6 }++-- | Random parameters for the Monte-Carlo simulation.+randomParams :: IO Parameters+randomParams =+ do averagePayableDelay <- newRandomParameter 0.07 0.11+ billingProcessingTime <- newRandomParameter 0.03 0.05+ buildingTime <- newRandomParameter 0.8 1.2+ fractionalLossRate <- newRandomParameter 0.05 0.08+ interestRate <- newRandomParameter 0.09 0.15+ price <- newRandomParameter 0.9 1.2+ productionCapacity <- newRandomParameter 2200 2600+ requiredInvestment <- newRandomParameter 1800 2200+ variableProductionCost <- newRandomParameter 0.5 0.7+ return defaultParams { paramsAveragePayableDelay = averagePayableDelay,+ paramsBillingProcessingTime = billingProcessingTime,+ paramsBuildingTime = buildingTime,+ paramsFractionalLossRate = fractionalLossRate,+ paramsInterestRate = interestRate,+ paramsPrice = price,+ paramsProductionCapacity = productionCapacity,+ paramsRequiredInvestment = requiredInvestment,+ paramsVariableProductionCost = variableProductionCost }++-- | This is the model itself that returns experimental data.+model :: Parameters -> Simulation ExperimentData+model params =+ mdo let liftParam :: (Parameters -> Simulation a) -> Dynamics a+ liftParam f = liftSimulation $ f params++ -- the equations below are given in an arbitrary order!++ bookValue <- integ (newInvestment - taxDepreciation) 0+ let taxDepreciation = bookValue / taxDepreciationTime+ taxableIncome = grossIncome - directCosts - losses+ - interestPayments - taxDepreciation+ production = availableCapacity+ availableCapacity = ifDynamics (time .>=. buildingTime)+ productionCapacity 0+ taxDepreciationTime = liftParam paramsTaxDepreciationTime+ taxRate = liftParam paramsTaxRate+ accountsReceivable <- integ (billings - cashReceipts - losses)+ (billings / (1 / averagePayableDelay+ + fractionalLossRate))+ let averagePayableDelay =+ liftParam paramsAveragePayableDelay+ awaitingBilling <- integ (price * production - billings)+ (price * production * billingProcessingTime)+ let billingProcessingTime =+ liftParam paramsBillingProcessingTime+ billings = awaitingBilling / billingProcessingTime+ borrowing = newInvestment * debtFinancingFraction+ buildingTime = liftParam paramsBuildingTime+ cashReceipts = accountsReceivable / averagePayableDelay+ debt <- integ (borrowing - principalRepayment) 0+ let debtFinancingFraction = liftParam paramsDebtFinancingFraction+ debtRetirementTime = liftParam paramsDebtRetirementTime+ directCosts = production * variableProductionCost+ discountRate = liftParam paramsDiscountRate+ fractionalLossRate = liftParam paramsFractionalLossRate+ grossIncome = billings+ interestPayments = debt * interestRate+ interestRate = liftParam paramsInterestRate+ losses = accountsReceivable * fractionalLossRate+ netCashFlow = cashReceipts + borrowing - newInvestment+ - directCosts - interestPayments+ - principalRepayment - taxes+ netIncome = taxableIncome - taxes+ newInvestment = ifDynamics (time .>=. buildingTime)+ 0 (requiredInvestment / buildingTime)+ npvCashFlow <- npv netCashFlow discountRate 0 1+ npvIncome <- npv netIncome discountRate 0 1+ let price = liftParam paramsPrice+ principalRepayment = debt / debtRetirementTime+ productionCapacity = liftParam paramsProductionCapacity+ requiredInvestment = liftParam paramsRequiredInvestment+ taxes = taxableIncome * taxRate+ variableProductionCost = liftParam paramsVariableProductionCost++ -- we have to create an event queue to return the experimental data,+ -- although it was not used in the model itself ++ queue <- newQueue+ + experimentDataInStartTime queue+ [(netIncomeName, seriesEntity "Net income" netIncome),+ (netCashFlowName, seriesEntity "Net cash flow" netCashFlow),+ (npvIncomeName, seriesEntity "NPV income" npvIncome),+ (npvCashFlowName, seriesEntity "NPV cash flow" npvCashFlow)]++-- the names of the variables we are interested in+netIncomeName = "netIncome"+netCashFlowName = "netCashFlow"+npvIncomeName = "npvIncome"+npvCashFlowName = "npvCashFlow"++-- the simulation specs+specs = Specs 0 5 0.015625 RungeKutta4++-- | The experiment for the Monte-Carlo simulation.+monteCarloExperiment :: Experiment+monteCarloExperiment =+ defaultExperiment {+ experimentSpecs = specs,+ experimentRunCount = 1000,+ experimentTitle = "Financial Model (the Monte-Carlo simulation)",+ experimentDescription = "Financial Model (the Monte-Carlo simulation) as described in " +++ "Vensim 5 Modeling Guide, Chapter Financial Modeling and Risk.",+ experimentGenerators =+ [outputView defaultExperimentSpecsView,+ + outputView $ defaultDeviationChartView {+ deviationChartTitle = "The deviation chart for Net Income and Cash Flow",+ deviationChartSeries = [Left netIncomeName, + Left netCashFlowName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "The deviation chart for Net Present Value of Income and Cash Flow",+ deviationChartSeries = [Left npvIncomeName, + Left npvCashFlowName] },++ outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Histogram for Net Income and Cash Flow",+ finalHistogramSeries = [netIncomeName, netCashFlowName] },++ outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Histogram for Net Present Value of Income and Cash Flow",+ finalHistogramSeries = [npvIncomeName, npvCashFlowName] },++ outputView $ defaultFinalStatsView {+ finalStatsTitle = "Summary for Net Income and Cash Flow",+ finalStatsSeries = [netIncomeName, netCashFlowName] },++ outputView $ defaultFinalStatsView {+ finalStatsTitle = "Summary for Net Present Value of Income and Cash Flow",+ finalStatsSeries = [npvIncomeName, npvCashFlowName] } ] }++-- | The experiment with single simulation run.+singleExperiment :: Experiment+singleExperiment =+ defaultExperiment {+ experimentSpecs = specs,+ experimentTitle = "Financial Model",+ experimentDescription = "Financial Model as described in " +++ "Vensim 5 Modeling Guide, Chapter Financial Modeling and Risk.",+ experimentGenerators =+ [outputView defaultExperimentSpecsView,+ + outputView $ defaultTimeSeriesView {+ timeSeriesTitle = "Time series of Net Income and Cash Flow",+ timeSeries = [Left netIncomeName, + Left netCashFlowName] },+ + outputView $ defaultTimeSeriesView {+ timeSeriesTitle = "Time series of Net Present Value for Income and Cash Flow",+ timeSeries = [Left npvIncomeName, + Left npvCashFlowName] },++ outputView $ defaultTableView {+ tableTitle = "Table",+ tableSeries = [netIncomeName, netCashFlowName,+ npvIncomeName, npvCashFlowName] } ] }++main = do+ + -- run the ordinary simulation+ putStrLn "*** The simulation with default parameters..."+ runExperiment singleExperiment $ model defaultParams+ putStrLn ""++ -- run the Monte-Carlo simulation+ putStrLn "*** The Monte-Carlo simulation..."+ randomParams >>= runExperimentParallel monteCarloExperiment . model
+ examples/Furnace.hs view
@@ -0,0 +1,559 @@++-- This is a model of the Furnace. It is described in different sources [1, 2].+--+-- [1] { add a foreign source in English }+--+-- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006+--+-- This model is often used in the literature as an example of combined+-- continuous-discrete simulation but this is not a point here. It illustrates+-- how the time-driven and process-oriented simulation models can be combined+-- based on the common event queue. It still uses the differential equation but+-- it is modeled directly [3] with help of the Euler method within the time-driven+-- part of the combined model.+--+-- [3] The time bounds for such an equation are much smaller than that ones which are defined+-- by the specs. Therefore there is no sense to use the 'integ' function as it would be+-- very slow because of large allocating memory for each integral, although it is possible.+--+-- However, you can still combine the differential (and difference) equations with the DES and+-- agent-based models. The integral (as well as any 'Dynamics' computation) can be used directly+-- in the DES sub-model. But to update something from the DES sub-model that could be used aready+-- in the differential equations, you should save data with help of types 'Var' or 'UVar' as they+-- keep all the history of their past values. Also the values of these two types are managed by+-- the event queue that allows synchronizing them with the DES sub-model.+--+-- To define the external parameters for the Monte-Carlo simulation, see the Financial model.+--+-- To enable the parallel simulation, you should compile it+-- with option -threaded and then pass in other options +RTS -N2 -RTS+-- to the executable if you have a dual core processor without+-- hyper-threading. Also you can increase the number+-- of parallel threads via option -N if you have a more modern+-- processor.++import Data.Maybe+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.UVar+import Simulation.Aivika.Dynamics.Process+import Simulation.Aivika.Dynamics.Random+import Simulation.Aivika.Statistics++import Simulation.Aivika.Experiment+import Simulation.Aivika.Experiment.ExperimentSpecsView+import Simulation.Aivika.Experiment.FinalStatsView+import Simulation.Aivika.Experiment.DeviationChartView+import Simulation.Aivika.Experiment.FinalHistogramView++import qualified Simulation.Aivika.Queue as Q++-- | The simulation specs.+specs = Specs { spcStartTime = 0.0,+ -- spcStopTime = 1000.0,+ spcStopTime = 300.0,+ spcDT = 0.1,+ spcMethod = RungeKutta4 }+ +-- | Return an exponentially distributed random value with mean +-- 1 / @lambda@, where @lambda@ is a parameter of the function.+exprnd :: Double -> IO Double+exprnd lambda =+ do x <- getStdRandom random+ return (- log x / lambda)+ +-- | Return a random initial temperature of the item. +temprnd :: IO Double+temprnd =+ do x <- getStdRandom random+ return (400.0 + (600.0 - 400.0) * x)++-- | Represents the furnace.+data Furnace = + Furnace { furnaceQueue :: EventQueue,+ -- ^ The event queue.+ furnaceNormalGen :: IO Double,+ -- ^ The normal random number generator.+ furnacePits :: [Pit],+ -- ^ The pits for ingots.+ furnacePitCount :: UVar Int,+ -- ^ The count of active pits with ingots.+ 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 :: Ref (SamplingStats Int),+ -- ^ The statistics about the queue count.+ furnaceWaitCount :: Ref Int,+ -- ^ The count of awaiting ingots.+ furnaceWaitTime :: Ref Double,+ -- ^ The wait time for all loaded ingots.+ furnaceHeatingTime :: Ref Double,+ -- ^ The heating time for all unloaded ingots.+ furnaceTemp :: Ref Double,+ -- ^ The furnace temperature.+ furnaceTotalCount :: Ref Int,+ -- ^ The total count of ingots.+ furnaceLoadCount :: Ref Int,+ -- ^ The count of loaded ingots.+ furnaceUnloadCount :: Ref Int,+ -- ^ The count of unloaded ingots.+ furnaceUnloadTemps :: Ref [Double]+ -- ^ The temperatures of all unloaded ingots.+ }++-- | A pit in the furnace to place the ingots.+data Pit = + Pit { pitQueue :: EventQueue,+ -- ^ The bound dynamics queue.+ pitIngot :: Ref (Maybe Ingot),+ -- ^ The ingot in the pit.+ pitTemp :: Ref Double+ -- ^ The ingot temperature in the pit.+ }++data Ingot = + Ingot { ingotFurnace :: Furnace,+ -- ^ Return the furnace.+ ingotReceiveTime :: Double,+ -- ^ The time at which the ingot was received.+ ingotReceiveTemp :: Double,+ -- ^ The temperature with which the ingot was received.+ ingotLoadTime :: Double,+ -- ^ The time of loading in the furnace.+ ingotLoadTemp :: Double,+ -- ^ The temperature when the ingot was loaded in the furnace.+ ingotCoeff :: Double+ -- ^ The heating coefficient.+ }++-- | Create a furnace.+newFurnace :: EventQueue -> Simulation Furnace+newFurnace queue =+ do normalGen <- liftIO normalGen+ pits <- sequence [newPit queue | i <- [1..10]]+ pitCount <- newUVar queue 0+ pitCountStats <- newRef queue emptySamplingStats+ awaitingIngots <- liftIO Q.newQueue+ queueCount <- newUVar queue 0+ queueCountStats <- newRef queue emptySamplingStats+ waitCount <- newRef queue 0+ waitTime <- newRef queue 0.0+ heatingTime <- newRef queue 0.0+ h <- newRef queue 1650.0+ totalCount <- newRef queue 0+ loadCount <- newRef queue 0+ unloadCount <- newRef queue 0+ unloadTemps <- newRef queue []+ return Furnace { furnaceQueue = queue,+ furnaceNormalGen = normalGen,+ furnacePits = pits,+ furnacePitCount = pitCount,+ furnacePitCountStats = pitCountStats,+ furnaceAwaitingIngots = awaitingIngots,+ furnaceQueueCount = queueCount,+ furnaceQueueCountStats = queueCountStats,+ furnaceWaitCount = waitCount,+ furnaceWaitTime = waitTime,+ furnaceHeatingTime = heatingTime,+ furnaceTemp = h,+ furnaceTotalCount = totalCount,+ furnaceLoadCount = loadCount, + furnaceUnloadCount = unloadCount, + furnaceUnloadTemps = unloadTemps }++-- | Create a new pit.+newPit :: EventQueue -> Simulation Pit+newPit queue =+ do ingot <- newRef queue Nothing+ h' <- newRef queue 0.0+ return Pit { pitQueue = queue,+ pitIngot = ingot,+ pitTemp = h' }++-- | Create a new ingot.+newIngot :: Furnace -> Dynamics Ingot+newIngot furnace =+ do t <- time+ xi <- liftIO $ furnaceNormalGen furnace+ h' <- liftIO temprnd+ let c = 0.1 + (0.05 + xi * 0.01)+ return Ingot { ingotFurnace = furnace,+ ingotReceiveTime = t,+ ingotReceiveTemp = h',+ ingotLoadTime = t,+ ingotLoadTemp = h',+ ingotCoeff = c }++-- | Heat the ingot up in the pit if there is such an ingot.+heatPitUp :: Pit -> Dynamics ()+heatPitUp pit =+ do ingot <- readRef (pitIngot pit)+ case ingot of+ Nothing -> + return ()+ Just ingot -> do+ + -- update the temperature of the ingot.+ let furnace = ingotFurnace ingot+ dt' <- dt+ h' <- readRef (pitTemp pit)+ h <- readRef (furnaceTemp furnace)+ writeRef (pitTemp pit) $ + h' + dt' * (h - h') * ingotCoeff ingot++-- | Check whether there are ready ingots in the pits.+ingotsReady :: Furnace -> Dynamics Bool+ingotsReady furnace =+ fmap (not . null) $ + filterM (fmap (>= 2200.0) . readRef . pitTemp) $ + furnacePits furnace++-- | Try to unload the ready ingot from the specified pit.+tryUnloadPit :: Furnace -> Pit -> Dynamics ()+tryUnloadPit furnace pit =+ do h' <- readRef (pitTemp pit)+ when (h' >= 2000.0) $+ do Just ingot <- readRef (pitIngot pit) + unloadIngot ingot pit++-- | Try to load an awaiting ingot in the specified empty pit.+tryLoadPit :: Furnace -> Pit -> Dynamics () +tryLoadPit furnace pit =+ do let ingots = furnaceAwaitingIngots furnace+ flag <- liftIO $ Q.queueNull ingots+ unless flag $+ do ingot <- liftIO $ Q.queueFront ingots+ liftIO $ Q.dequeue ingots+ t' <- time+ modifyUVar (furnaceQueueCount furnace) (+ (-1))+ c <- readUVar (furnaceQueueCount furnace)+ modifyRef (furnaceQueueCountStats furnace) $+ addSamplingStats c+ loadIngot (ingot { ingotLoadTime = t',+ ingotLoadTemp = 400.0 }) pit+ +-- | Unload the ingot from the specified pit. +unloadIngot :: Ingot -> Pit -> Dynamics ()+unloadIngot ingot pit = + do h' <- readRef (pitTemp pit)+ writeRef (pitIngot pit) Nothing+ writeRef (pitTemp pit) 0.0+ + -- count the active pits+ let furnace = ingotFurnace ingot+ count <- readUVar (furnacePitCount furnace)+ writeUVar (furnacePitCount furnace) (count - 1)+ modifyRef (furnacePitCountStats furnace) $+ addSamplingStats (count - 1)+ + -- how long did we heat the ingot up?+ t' <- time+ modifyRef (furnaceHeatingTime furnace)+ (+ (t' - ingotLoadTime ingot))+ + -- what is the temperature of the unloaded ingot?+ modifyRef (furnaceUnloadTemps furnace) (h' :)+ + -- count the unloaded ingots+ modifyRef (furnaceUnloadCount furnace) (+ 1)+ +-- | Load the ingot in the specified pit+loadIngot :: Ingot -> Pit -> Dynamics ()+loadIngot ingot pit =+ do writeRef (pitIngot pit) $ Just ingot+ writeRef (pitTemp pit) $ ingotLoadTemp ingot+ + -- count the active pits+ let furnace = ingotFurnace ingot+ count <- readUVar (furnacePitCount furnace)+ writeUVar (furnacePitCount furnace) (count + 1)+ modifyRef (furnacePitCountStats furnace) $+ addSamplingStats (count + 1)+ + -- decrease the furnace temperature+ h <- readRef (furnaceTemp furnace)+ let h' = ingotLoadTemp ingot+ dh = - (h - h') / fromInteger (toInteger (count + 1))+ writeRef (furnaceTemp furnace) $ h + dh++ -- how long did we keep the ingot in the queue?+ t' <- time+ when (ingotReceiveTime ingot < t') $+ do modifyRef (furnaceWaitCount furnace) (+ 1) + modifyRef (furnaceWaitTime furnace)+ (+ (t' - ingotReceiveTime ingot))++ -- count the loaded ingots+ modifyRef (furnaceLoadCount furnace) (+ 1)+ +-- | 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+ pits' <- emptyPits furnace+ mapM_ (tryLoadPit furnace) pits'+ mapM_ heatPitUp pits+ + -- update the temperature of the furnace+ dt' <- dt+ h <- readRef (furnaceTemp furnace)+ writeRef (furnaceTemp furnace) $+ h + dt' * (2600.0 - h) * 0.2++-- | Return all empty pits.+emptyPits :: Furnace -> Dynamics [Pit]+emptyPits furnace =+ filterM (fmap isNothing . readRef . pitIngot) $+ furnacePits furnace++-- | Accept a new ingot.+acceptIngot :: Furnace -> Dynamics ()+acceptIngot furnace =+ do ingot <- newIngot furnace+ + -- counting+ modifyRef (furnaceTotalCount furnace) (+ 1)+ + -- check what to do with the new ingot+ count <- readUVar (furnacePitCount furnace)+ if count >= 10+ then do let ingots = furnaceAwaitingIngots furnace+ liftIO $ Q.enqueue ingots ingot+ modifyUVar (furnaceQueueCount furnace) (+ 1)+ c <- readUVar (furnaceQueueCount furnace)+ modifyRef (furnaceQueueCountStats furnace) $+ addSamplingStats c+ else do pit:_ <- emptyPits furnace+ loadIngot ingot pit+ +-- | Process the furnace.+processFurnace :: Furnace -> Process ()+processFurnace furnace =+ do delay <- liftIO $ exprnd (1.0 / 2.5)+ holdProcess delay+ -- we have got a new ingot+ liftDynamics $ acceptIngot furnace+ -- repeat it again+ processFurnace furnace++-- | Initialize the furnace.+initializeFurnace :: Furnace -> Dynamics ()+initializeFurnace furnace =+ do x1 <- newIngot furnace+ x2 <- newIngot furnace+ x3 <- newIngot furnace+ x4 <- newIngot furnace+ x5 <- newIngot furnace+ x6 <- newIngot furnace+ let p1 : p2 : p3 : p4 : p5 : p6 : ps = + furnacePits furnace+ loadIngot (x1 { ingotLoadTemp = 550.0 }) p1+ loadIngot (x2 { ingotLoadTemp = 600.0 }) p2+ loadIngot (x3 { ingotLoadTemp = 650.0 }) p3+ loadIngot (x4 { ingotLoadTemp = 700.0 }) p4+ loadIngot (x5 { ingotLoadTemp = 750.0 }) p5+ loadIngot (x6 { ingotLoadTemp = 800.0 }) p6+ writeRef (furnaceTotalCount furnace) 6+ writeRef (furnaceTemp furnace) 1650.0+ +-- | The simulation model that returns experimental data.+model :: Simulation ExperimentData+model =+ do queue <- newQueue+ furnace <- newFurnace queue+ pid <- newProcessID queue++ -- initialize the furnace and start its iterating in start time+ runDynamicsInStartTime $+ do initializeFurnace furnace+ startIteratingFurnace furnace+ + -- start accepting input ingots by launching the process+ runDynamicsInStartTime $+ do t0 <- starttime+ runProcess (processFurnace furnace) pid t0+ + experimentDataInStartTime queue $+ [(totalIngotCountName,+ seriesEntity "total ingot count" $+ furnaceTotalCount furnace),+ + (loadedIngotCountName,+ seriesEntity "loaded ingot count" $+ furnaceLoadCount furnace),+ + (readyIngotCountName,+ seriesEntity "ready ingot count" $+ furnaceUnloadCount furnace),++ (awaitedIngotCountName,+ seriesEntity "awaited in the queue ingot count" $+ furnaceWaitCount furnace),++ (readyIngotTempsName,+ seriesEntity "the temperature of ready ingot" $+ furnaceUnloadTemps furnace),+ + (pitCountStatsName,+ seriesEntity "the used pit count" $+ furnacePitCountStats furnace),+ + (queueCountStatsName,+ seriesEntity "the queue size" $+ furnaceQueueCountStats furnace),+ + (meanWaitTimeName,+ seriesEntity "the mean wait time" $+ readRef (furnaceWaitTime furnace) /+ fmap fromIntegral (readRef (furnaceWaitCount furnace))),++ (meanHeatingTimeName,+ seriesEntity "the mean heating time" $+ readRef (furnaceHeatingTime furnace) /+ fmap fromIntegral (readRef (furnaceUnloadCount furnace))) ]+ +totalIngotCountName = "totalIngotCount"+loadedIngotCountName = "loadedIngotCount"+readyIngotCountName = "readyIngotCount"+awaitedIngotCountName = "awaitedIngotCount"+readyIngotTempsName = "readyIngotTemps"+pitCountStatsName = "pitCountStats"+queueCountStatsName = "queueCountStats"+meanWaitTimeName = "the mean wait time in the queue"+meanHeatingTimeName = "the mean heating time"++-- | The experiment.+experiment :: Experiment+experiment =+ defaultExperiment {+ experimentSpecs = specs,+ -- experimentRunCount = 1000,+ experimentRunCount = 100,+ experimentTitle = "The Furnace model (the Monte-Carlo simulation)",+ experimentGenerators =+ [outputView defaultExperimentSpecsView,+ + outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 1.1",+ deviationChartPlotTitle = "The total, loaded and ready ingot counts",+ deviationChartSeries = [Right totalIngotCountName,+ Right loadedIngotCountName,+ Right readyIngotCountName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 1.2",+ deviationChartPlotTitle = "The awaited in the queue ingot count",+ deviationChartSeries = [Right awaitedIngotCountName] },++ outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Final Histogram - 1.1",+ finalHistogramPlotTitle = "The distribution of total, loaded and ready " +++ "ingot counts in the final time point.",+ finalHistogramSeries = [totalIngotCountName,+ loadedIngotCountName,+ readyIngotCountName] },+ + outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Final Histogram - 1.2",+ finalHistogramPlotTitle = "The distribution of the awaited in the queue " +++ "ingot count in the final time point.",+ finalHistogramSeries = [awaitedIngotCountName] },+ + outputView $ defaultFinalStatsView {+ finalStatsTitle = "Final Statistics - 1",+ finalStatsDescription = "The summary of total, loaded, ready and awaited in " +++ "the queue ingot counts in the final time point.",+ finalStatsSeries = [totalIngotCountName,+ loadedIngotCountName,+ readyIngotCountName,+ awaitedIngotCountName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 2",+ deviationChartPlotTitle = "The used pit count",+ deviationChartSeries = [Right pitCountStatsName] },+ + outputView $ defaultFinalStatsView {+ finalStatsTitle = "Final Statistics - 2",+ finalStatsDescription = "The summary of the used pit count in the final time point.",+ finalStatsSeries = [pitCountStatsName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 3",+ deviationChartPlotTitle = "The queue size",+ deviationChartSeries = [Right queueCountStatsName] },+ + outputView $ defaultFinalStatsView {+ finalStatsTitle = "Final Statistics - 3",+ finalStatsDescription = "The summary of the queue size in the final time point.",+ finalStatsSeries = [queueCountStatsName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 4",+ deviationChartPlotTitle = "The mean wait time",+ deviationChartSeries = [Right meanWaitTimeName] },++ outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Final Histogram - 4",+ finalHistogramPlotTitle = "The distribution of the mean wait time " +++ "in the final simulation time point.",+ finalHistogramSeries = [meanWaitTimeName] },+ + outputView $ defaultFinalStatsView {+ finalStatsTitle = "Final Statistics - 4",+ finalStatsDescription = "The summary of the mean wait time in " +++ "the final simulation time point.",+ finalStatsSeries = [meanWaitTimeName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 5",+ deviationChartPlotTitle = "The mean heating time",+ deviationChartSeries = [Right meanHeatingTimeName] },++ outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Final Histogram - 5",+ finalHistogramPlotTitle = "The distribution of the mean heating time " +++ "in the final simulation time point.",+ finalHistogramSeries = [meanHeatingTimeName] },+ + outputView $ defaultFinalStatsView {+ finalStatsTitle = "Final Statistics - 5",+ finalStatsDescription = "The summary of the mean heating time in " +++ "the final simulation time point.",+ finalStatsSeries = [meanHeatingTimeName] },++ outputView $ defaultDeviationChartView {+ deviationChartTitle = "Deviation Chart - 6",+ deviationChartPlotTitle = "The ready ingot temperature",+ deviationChartSeries = [Right readyIngotTempsName] },++ outputView $ defaultFinalHistogramView {+ finalHistogramTitle = "Final Histogram - 6",+ finalHistogramPlotTitle = "The distribution of the ready ingot temperature " +++ "in the final simulation time point.",+ finalHistogramSeries = [readyIngotTempsName] },+ + outputView $ defaultFinalStatsView {+ finalStatsTitle = "Final Statistics - 6",+ finalStatsDescription = "The summary of the ready ingot temperature in " +++ "the final simulation time point.",+ finalStatsSeries = [readyIngotTempsName] }+ ] }++-- | The main program that launches the simulation experiment to produce the HTML file.+main = runExperimentParallel experiment model