-- |
-- Module     : Simulation.Aivika.Internal.Specs
-- Copyright  : Copyright (c) 2009-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 8.0.1
--
-- It defines the simulation specs and related stuff.
module Simulation.Aivika.Internal.Specs
       (Specs(..),
        Method(..),
        Run(..),
        Point(..),
        EventQueue(..),
        newEventQueue,
        basicTime,
        integIterationBnds,
        integIterationHiBnd,
        integIterationLoBnd,
        integPhaseBnds,
        integPhaseHiBnd,
        integPhaseLoBnd,
        integTimes,
        integPoints,
        integPointsStartingFrom,
        integStartPoint,
        integStopPoint,
        simulationStopPoint,
        timeGrid,
        pointAt,
        delayPoint) where

import Data.IORef

import Simulation.Aivika.Generator
import qualified Simulation.Aivika.PriorityQueue as PQ

-- | It defines the simulation specs.
data Specs = Specs { Specs -> Double
spcStartTime :: Double,    -- ^ the start time
                     Specs -> Double
spcStopTime :: Double,     -- ^ the stop time
                     Specs -> Double
spcDT :: Double,           -- ^ the integration time step
                     Specs -> Method
spcMethod :: Method,       -- ^ the integration method
                     Specs -> GeneratorType
spcGeneratorType :: GeneratorType
                     -- ^ the type of the random number generator
                   }

-- | It defines the integration method.
data Method = Euler          -- ^ Euler's method
            | RungeKutta2    -- ^ the 2nd order Runge-Kutta method
            | RungeKutta4    -- ^ the 4th order Runge-Kutta method
            | RungeKutta4b   -- ^ the 4th order Runge-Kutta 3/8-method
            deriving (Method -> Method -> Bool
(Method -> Method -> Bool)
-> (Method -> Method -> Bool) -> Eq Method
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Method -> Method -> Bool
$c/= :: Method -> Method -> Bool
== :: Method -> Method -> Bool
$c== :: Method -> Method -> Bool
Eq, Eq Method
Eq Method
-> (Method -> Method -> Ordering)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Method)
-> (Method -> Method -> Method)
-> Ord Method
Method -> Method -> Bool
Method -> Method -> Ordering
Method -> Method -> Method
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Method -> Method -> Method
$cmin :: Method -> Method -> Method
max :: Method -> Method -> Method
$cmax :: Method -> Method -> Method
>= :: Method -> Method -> Bool
$c>= :: Method -> Method -> Bool
> :: Method -> Method -> Bool
$c> :: Method -> Method -> Bool
<= :: Method -> Method -> Bool
$c<= :: Method -> Method -> Bool
< :: Method -> Method -> Bool
$c< :: Method -> Method -> Bool
compare :: Method -> Method -> Ordering
$ccompare :: Method -> Method -> Ordering
$cp1Ord :: Eq Method
Ord, Int -> Method -> ShowS
[Method] -> ShowS
Method -> String
(Int -> Method -> ShowS)
-> (Method -> String) -> ([Method] -> ShowS) -> Show Method
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Method] -> ShowS
$cshowList :: [Method] -> ShowS
show :: Method -> String
$cshow :: Method -> String
showsPrec :: Int -> Method -> ShowS
$cshowsPrec :: Int -> Method -> ShowS
Show)

-- | It indentifies the simulation run.
data Run = Run { Run -> Specs
runSpecs :: Specs,  -- ^ the simulation specs
                 Run -> Int
runIndex :: Int,    -- ^ the current simulation run index
                 Run -> Int
runCount :: Int,    -- ^ the total number of runs in this experiment
                 Run -> EventQueue
runEventQueue :: EventQueue,  -- ^ the event queue
                 Run -> Generator
runGenerator :: Generator     -- ^ the random number generator
               }

-- | It defines the simulation point appended with the additional information.
data Point = Point { Point -> Specs
pointSpecs :: Specs,    -- ^ the simulation specs
                     Point -> Run
pointRun :: Run,        -- ^ the simulation run
                     Point -> Double
pointTime :: Double,    -- ^ the current time
                     Point -> Int
pointIteration :: Int,  -- ^ the current iteration
                     Point -> Int
pointPhase :: Int       -- ^ the current phase
                   }

-- | It represents the event queue.
data EventQueue = EventQueue { EventQueue -> PriorityQueue (Point -> IO ())
queuePQ :: PQ.PriorityQueue (Point -> IO ()),
                               -- ^ the underlying priority queue
                               EventQueue -> IORef Bool
queueBusy :: IORef Bool,
                               -- ^ whether the queue is currently processing events
                               EventQueue -> IORef Double
queueTime :: IORef Double
                               -- ^ the actual time of the event queue
                             }

-- | Create a new event queue by the specified specs.
newEventQueue :: Specs -> IO EventQueue
newEventQueue :: Specs -> IO EventQueue
newEventQueue Specs
specs = 
  do IORef Bool
f <- Bool -> IO (IORef Bool)
forall a. a -> IO (IORef a)
newIORef Bool
False
     IORef Double
t <- Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef (Double -> IO (IORef Double)) -> Double -> IO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Specs -> Double
spcStartTime Specs
specs
     PriorityQueue (Point -> IO ())
pq <- IO (PriorityQueue (Point -> IO ()))
forall a. IO (PriorityQueue a)
PQ.newQueue
     EventQueue -> IO EventQueue
forall (m :: * -> *) a. Monad m => a -> m a
return EventQueue :: PriorityQueue (Point -> IO ())
-> IORef Bool -> IORef Double -> EventQueue
EventQueue { queuePQ :: PriorityQueue (Point -> IO ())
queuePQ   = PriorityQueue (Point -> IO ())
pq,
                         queueBusy :: IORef Bool
queueBusy = IORef Bool
f,
                         queueTime :: IORef Double
queueTime = IORef Double
t }

-- | Returns the integration iterations starting from zero.
integIterations :: Specs -> [Int]
integIterations :: Specs -> [Int]
integIterations Specs
sc = [Int
i1 .. Int
i2] where
  i1 :: Int
i1 = Specs -> Int
integIterationLoBnd Specs
sc
  i2 :: Int
i2 = Specs -> Int
integIterationHiBnd Specs
sc

-- | Returns the first and last integration iterations.
integIterationBnds :: Specs -> (Int, Int)
integIterationBnds :: Specs -> (Int, Int)
integIterationBnds Specs
sc = (Int
i1, Int
i2) where
  i1 :: Int
i1 = Specs -> Int
integIterationLoBnd Specs
sc
  i2 :: Int
i2 = Specs -> Int
integIterationHiBnd Specs
sc

-- | Returns the first integration iteration, i.e. zero.
integIterationLoBnd :: Specs -> Int
integIterationLoBnd :: Specs -> Int
integIterationLoBnd Specs
sc = Int
0

-- | Returns the last integration iteration.
integIterationHiBnd :: Specs -> Int
integIterationHiBnd :: Specs -> Int
integIterationHiBnd Specs
sc =
  let n :: Int
n = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round ((Specs -> Double
spcStopTime Specs
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
- 
                  Specs -> Double
spcStartTime Specs
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc)
  in if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
     then
       String -> Int
forall a. HasCallStack => String -> a
error (String -> Int) -> String -> Int
forall a b. (a -> b) -> a -> b
$
       String
"The iteration number in the stop time has a negative value. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"Either the simulation specs are incorrect, " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"or a floating point overflow occurred, " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"for example, when using a too small integration time step. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"You have to define this time step regardless of " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"whether you actually use it or not, " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"for Aivika allows combining the ordinary differential equations " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"with the discrete event simulation within one model. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"So, if you are still using the 32-bit architecture and " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"you do need a small integration time step " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"for integrating the equations " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"then you might think of using the 64-bit architecture. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"Although you could probably just forget " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"to increase the time step " String -> ShowS
forall a. [a] -> [a] -> [a]
++
       String
"after increasing the stop time: integIterationHiBnd"
     else Int
n

-- | Returns the phases for the specified simulation specs starting from zero.
integPhases :: Specs -> [Int]
integPhases :: Specs -> [Int]
integPhases Specs
sc = 
  case Specs -> Method
spcMethod Specs
sc of
    Method
Euler -> [Int
0]
    Method
RungeKutta2 -> [Int
0, Int
1]
    Method
RungeKutta4 -> [Int
0, Int
1, Int
2, Int
3]
    Method
RungeKutta4b -> [Int
0, Int
1, Int
2, Int
3]

-- | Returns the first and last integration phases.
integPhaseBnds :: Specs -> (Int, Int)
integPhaseBnds :: Specs -> (Int, Int)
integPhaseBnds Specs
sc = 
  case Specs -> Method
spcMethod Specs
sc of
    Method
Euler -> (Int
0, Int
0)
    Method
RungeKutta2 -> (Int
0, Int
1)
    Method
RungeKutta4 -> (Int
0, Int
3)
    Method
RungeKutta4b -> (Int
0, Int
3)

-- | Returns the first integration phase, i.e. zero.
integPhaseLoBnd :: Specs -> Int
integPhaseLoBnd :: Specs -> Int
integPhaseLoBnd Specs
sc = Int
0
                  
-- | Returns the last integration phase, 0 for Euler's method, 1 for RK2 and 3 for RK4.
integPhaseHiBnd :: Specs -> Int
integPhaseHiBnd :: Specs -> Int
integPhaseHiBnd Specs
sc = 
  case Specs -> Method
spcMethod Specs
sc of
    Method
Euler -> Int
0
    Method
RungeKutta2 -> Int
1
    Method
RungeKutta4 -> Int
3
    Method
RungeKutta4b -> Int
3

-- | Returns a simulation time for the integration point specified by 
-- the specs, iteration and phase.
basicTime :: Specs -> Int -> Int -> Double
basicTime :: Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
ph =
  if Int
ph Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then 
    String -> Double
forall a. HasCallStack => String -> a
error String
"Incorrect phase: basicTime"
  else
    Specs -> Double
spcStartTime Specs
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
n' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Method -> Int -> Double
forall a. (Eq a, Num a) => Method -> a -> Double
delta (Specs -> Method
spcMethod Specs
sc) Int
ph 
      where n' :: Double
n' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
            delta :: Method -> a -> Double
delta Method
Euler       a
0 = Double
0
            delta Method
RungeKutta2 a
0 = Double
0
            delta Method
RungeKutta2 a
1 = Specs -> Double
spcDT Specs
sc
            delta Method
RungeKutta4 a
0 = Double
0
            delta Method
RungeKutta4 a
1 = Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 a
2 = Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 a
3 = Specs -> Double
spcDT Specs
sc
            delta Method
RungeKutta4b a
0 = Double
0
            delta Method
RungeKutta4b a
1 = Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
            delta Method
RungeKutta4b a
2 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
            delta Method
RungeKutta4b a
3 = Specs -> Double
spcDT Specs
sc

-- | Return the integration time values.
integTimes :: Specs -> [Double]
integTimes :: Specs -> [Double]
integTimes Specs
sc = (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Double
t [Int
nl .. Int
nu]
  where (Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
        t :: Int -> Double
t Int
n = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0

-- | Return the integration time points.
integPoints :: Run -> [Point]
integPoints :: Run -> [Point]
integPoints Run
r = [Point]
points
  where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
        (Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
        points :: [Point]
points   = (Int -> Point) -> [Int] -> [Point]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Point
point [Int
nl .. Int
nu]
        point :: Int -> Point
point Int
n  = Point :: Specs -> Run -> Double -> Int -> Int -> Point
Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
                           pointRun :: Run
pointRun = Run
r,
                           pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the start time point.
integStartPoint :: Run -> Point
integStartPoint :: Run -> Point
integStartPoint Run
r = Int -> Point
point Int
nl
  where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
        (Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
        point :: Int -> Point
point Int
n  = Point :: Specs -> Run -> Double -> Int -> Int -> Point
Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
                           pointRun :: Run
pointRun = Run
r,
                           pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the integration stop time point.
integStopPoint :: Run -> Point
integStopPoint :: Run -> Point
integStopPoint Run
r = Int -> Point
point Int
nu
  where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
        (Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
        point :: Int -> Point
point Int
n  = Point :: Specs -> Run -> Double -> Int -> Int -> Point
Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
                           pointRun :: Run
pointRun = Run
r,
                           pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the simulation stop time point.
simulationStopPoint :: Run -> Point
simulationStopPoint :: Run -> Point
simulationStopPoint Run
r = Run -> Double -> Point
pointAt Run
r (Specs -> Double
spcStopTime (Specs -> Double) -> Specs -> Double
forall a b. (a -> b) -> a -> b
$ Run -> Specs
runSpecs Run
r)

-- | Return the point at the specified time.
pointAt :: Run -> Double -> Point
pointAt :: Run -> Double -> Point
pointAt Run
r Double
t = Point
p
  where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
        t0 :: Double
t0 = Specs -> Double
spcStartTime Specs
sc
        dt :: Double
dt = Specs -> Double
spcDT Specs
sc
        n :: Int
n  = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor ((Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dt)
        p :: Point
p = Point :: Specs -> Run -> Double -> Int -> Int -> Point
Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
                    pointRun :: Run
pointRun = Run
r,
                    pointTime :: Double
pointTime = Double
t,
                    pointIteration :: Int
pointIteration = Int
n,
                    pointPhase :: Int
pointPhase = -Int
1 }

-- | Return the integration time points startin from the specified iteration.
integPointsStartingFrom :: Point -> [Point]
integPointsStartingFrom :: Point -> [Point]
integPointsStartingFrom Point
p = [Point]
points
  where r :: Run
r  = Point -> Run
pointRun Point
p
        sc :: Specs
sc = Run -> Specs
runSpecs Run
r
        (Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
        n0 :: Int
n0       = if Point -> Int
pointPhase Point
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                   then Point -> Int
pointIteration Point
p
                   else Point -> Int
pointIteration Point
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        points :: [Point]
points   = (Int -> Point) -> [Int] -> [Point]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Point
point [Int
n0 .. Int
nu]
        point :: Int -> Point
point Int
n  = Point :: Specs -> Run -> Double -> Int -> Int -> Point
Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
                           pointRun :: Run
pointRun = Run
r,
                           pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the indexed time values in the grid by specified size.
timeGrid :: Specs -> Int -> [(Int, Double)]
timeGrid :: Specs -> Int -> [(Int, Double)]
timeGrid Specs
sc Int
n =
  let t0 :: Double
t0 = Specs -> Double
spcStartTime Specs
sc
      t2 :: Double
t2 = Specs -> Double
spcStopTime Specs
sc
      n' :: Int
n' = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
      dt :: Double
dt = (Double
t2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n')
      f :: Int -> (Int, Double)
f Int
i | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0    = (Int
i, Double
t0)
          | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n'   = (Int
i, Double
t2)
          | Bool
otherwise = (Int
i, Double
t0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt)
  in (Int -> (Int, Double)) -> [Int] -> [(Int, Double)]
forall a b. (a -> b) -> [a] -> [b]
map Int -> (Int, Double)
f [Int
0 .. Int
n']

-- | Delay the point by the specified positive number of iterations.
delayPoint :: Point -> Int -> Point
delayPoint :: Point -> Int -> Point
delayPoint Point
p Int
dn
  | Int
dn Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0   = String -> Point
forall a. HasCallStack => String -> a
error String
"Expected the positive number of iterations: delayPoint"
  | Bool
otherwise =
    let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
        n :: Int
n  = Point -> Int
pointIteration Point
p
        ph :: Int
ph = Point -> Int
pointPhase Point
p
    in if Int
ph Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
       then let t' :: Double
t' = Point -> Double
pointTime Point
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dn Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc
                n' :: Int
n' = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
$ (Double
t' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Specs -> Double
spcStartTime Specs
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
            in Point
p { pointTime :: Double
pointTime = Double
t',
                   pointIteration :: Int
pointIteration = Int
n',
                   pointPhase :: Int
pointPhase = -Int
1 }
       else let n' :: Int
n' = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
dn
                t' :: Double
t' = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n' Int
ph
            in Point
p { pointTime :: Double
pointTime = Double
t',
                   pointIteration :: Int
pointIteration = Int
n',
                   pointPhase :: Int
pointPhase = Int
ph }