diff --git a/Simulation/Aivika/Experiment/DeviationChartView.hs b/Simulation/Aivika/Experiment/DeviationChartView.hs --- a/Simulation/Aivika/Experiment/DeviationChartView.hs +++ b/Simulation/Aivika/Experiment/DeviationChartView.hs @@ -1,11 +1,11 @@ -- | -- Module : Simulation.Aivika.Experiment.DeviationChartView --- Copyright : Copyright (c) 2012, David Sorokin +-- Copyright : Copyright (c) 2012-2013, David Sorokin -- License : BSD3 -- Maintainer : David Sorokin -- Stability : experimental --- Tested with: GHC 7.4.1 +-- Tested with: GHC 7.6.3 -- -- The module defines '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 () diff --git a/Simulation/Aivika/Experiment/FinalHistogramView.hs b/Simulation/Aivika/Experiment/FinalHistogramView.hs --- a/Simulation/Aivika/Experiment/FinalHistogramView.hs +++ b/Simulation/Aivika/Experiment/FinalHistogramView.hs @@ -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 $ diff --git a/Simulation/Aivika/Experiment/FinalXYChartView.hs b/Simulation/Aivika/Experiment/FinalXYChartView.hs --- a/Simulation/Aivika/Experiment/FinalXYChartView.hs +++ b/Simulation/Aivika/Experiment/FinalXYChartView.hs @@ -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 diff --git a/Simulation/Aivika/Experiment/HistogramView.hs b/Simulation/Aivika/Experiment/HistogramView.hs --- a/Simulation/Aivika/Experiment/HistogramView.hs +++ b/Simulation/Aivika/Experiment/HistogramView.hs @@ -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 $ diff --git a/Simulation/Aivika/Experiment/TimeSeriesView.hs b/Simulation/Aivika/Experiment/TimeSeriesView.hs --- a/Simulation/Aivika/Experiment/TimeSeriesView.hs +++ b/Simulation/Aivika/Experiment/TimeSeriesView.hs @@ -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 diff --git a/Simulation/Aivika/Experiment/XYChartView.hs b/Simulation/Aivika/Experiment/XYChartView.hs --- a/Simulation/Aivika/Experiment/XYChartView.hs +++ b/Simulation/Aivika/Experiment/XYChartView.hs @@ -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 diff --git a/aivika-experiment-chart.cabal b/aivika-experiment-chart.cabal --- a/aivika-experiment-chart.cabal +++ b/aivika-experiment-chart.cabal @@ -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 +copyright: (c) 2012-2013. 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 +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 diff --git a/examples/ChemicalReaction.hs b/examples/ChemicalReaction.hs --- a/examples/ChemicalReaction.hs +++ b/examples/ChemicalReaction.hs @@ -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 diff --git a/examples/DifferenceEquations.hs b/examples/DifferenceEquations.hs new file mode 100644 --- /dev/null +++ b/examples/DifferenceEquations.hs @@ -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 diff --git a/examples/Financial.hs b/examples/Financial.hs new file mode 100644 --- /dev/null +++ b/examples/Financial.hs @@ -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 diff --git a/examples/Furnace.hs b/examples/Furnace.hs new file mode 100644 --- /dev/null +++ b/examples/Furnace.hs @@ -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