-- |
-- 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(..),
        EventPriority(..),
        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.EventQueue 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
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
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
Ord, Int -> Method -> ShowS
[Method] -> ShowS
Method -> String
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
pointPriority :: EventPriority,    -- ^ the current priority
                     Point -> Int
pointIteration :: Int,  -- ^ the current iteration
                     Point -> Int
pointPhase :: Int       -- ^ the current phase
                   }

-- | The event priority (greater is higher).
type EventPriority = Int

-- | 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 <- forall a. a -> IO (IORef a)
newIORef Bool
False
     IORef Double
t <- forall a. a -> IO (IORef a)
newIORef forall a b. (a -> b) -> a -> b
$ Specs -> Double
spcStartTime Specs
specs
     PriorityQueue (Point -> IO ())
pq <- forall a. IO (PriorityQueue a)
PQ.newQueue
     forall (m :: * -> *) a. Monad m => a -> m a
return 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 = forall a b. (RealFrac a, Integral b) => a -> b
round ((Specs -> Double
spcStopTime Specs
sc forall a. Num a => a -> a -> a
- 
                  Specs -> Double
spcStartTime Specs
sc) forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc)
  in if Int
n forall a. Ord a => a -> a -> Bool
< Int
0
     then
       forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$
       String
"The iteration number in the stop time has a negative value. " forall a. [a] -> [a] -> [a]
++
       String
"Either the simulation specs are incorrect, " forall a. [a] -> [a] -> [a]
++
       String
"or a floating point overflow occurred, " forall a. [a] -> [a] -> [a]
++
       String
"for example, when using a too small integration time step. " forall a. [a] -> [a] -> [a]
++
       String
"You have to define this time step regardless of " forall a. [a] -> [a] -> [a]
++
       String
"whether you actually use it or not, " forall a. [a] -> [a] -> [a]
++
       String
"for Aivika allows combining the ordinary differential equations " forall a. [a] -> [a] -> [a]
++
       String
"with the discrete event simulation within one model. " forall a. [a] -> [a] -> [a]
++
       String
"So, if you are still using the 32-bit architecture and " forall a. [a] -> [a] -> [a]
++
       String
"you do need a small integration time step " forall a. [a] -> [a] -> [a]
++
       String
"for integrating the equations " forall a. [a] -> [a] -> [a]
++
       String
"then you might think of using the 64-bit architecture. " forall a. [a] -> [a] -> [a]
++
       String
"Although you could probably just forget " forall a. [a] -> [a] -> [a]
++
       String
"to increase the time step " 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 forall a. Ord a => a -> a -> Bool
< Int
0 then 
    forall a. HasCallStack => String -> a
error String
"Incorrect phase: basicTime"
  else
    Specs -> Double
spcStartTime Specs
sc forall a. Num a => a -> a -> a
+ Double
n' forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc forall a. Num a => a -> a -> a
+ forall {a}. (Eq a, Num a) => Method -> a -> Double
delta (Specs -> Method
spcMethod Specs
sc) Int
ph 
      where n' :: Double
n' = 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 forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 a
2 = Specs -> Double
spcDT Specs
sc 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 forall a. Fractional a => a -> a -> a
/ Double
3
            delta Method
RungeKutta4b a
2 = Double
2 forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc 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 = 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   = forall a b. (a -> b) -> [a] -> [b]
map Int -> Point
point [Int
nl .. Int
nu]
        point :: Int -> Point
point Int
n  = 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,
                           pointPriority :: Int
pointPriority = 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 { pointSpecs :: Specs
pointSpecs = Specs
sc,
                           pointRun :: Run
pointRun = Run
r,
                           pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
                           pointPriority :: Int
pointPriority = 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 { pointSpecs :: Specs
pointSpecs = Specs
sc,
                           pointRun :: Run
pointRun = Run
r,
                           pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
                           pointPriority :: Int
pointPriority = 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 -> Int -> Point
pointAt Run
r (Specs -> Double
spcStopTime forall a b. (a -> b) -> a -> b
$ Run -> Specs
runSpecs Run
r) Int
0

-- | Return the point at the specified time given the priority.
pointAt :: Run -> Double -> EventPriority -> Point
pointAt :: Run -> Double -> Int -> Point
pointAt Run
r Double
t Int
priority = 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  = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor ((Double
t forall a. Num a => a -> a -> a
- Double
t0) forall a. Fractional a => a -> a -> a
/ Double
dt)
        p :: Point
p = Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
                    pointRun :: Run
pointRun = Run
r,
                    pointTime :: Double
pointTime = Double
t,
                    pointPriority :: Int
pointPriority = Int
priority,
                    pointIteration :: Int
pointIteration = Int
n,
                    pointPhase :: Int
pointPhase = -Int
1 }

-- | Return the integration time points starting 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 forall a. Eq a => a -> a -> Bool
== Int
0
                   then Point -> Int
pointIteration Point
p
                   else Point -> Int
pointIteration Point
p forall a. Num a => a -> a -> a
+ Int
1
        points :: [Point]
points   = forall a b. (a -> b) -> [a] -> [b]
map Int -> Point
point [Int
n0 .. Int
nu]
        point :: Int -> Point
point Int
n  = 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,
                           pointPriority :: Int
pointPriority = 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' = forall a. Ord a => a -> a -> a
max (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
1
      dt :: Double
dt = (Double
t2 forall a. Num a => a -> a -> a
- Double
t0) forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n')
      f :: Int -> (Int, Double)
f Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
0    = (Int
i, Double
t0)
          | Int
i forall a. Eq a => a -> a -> Bool
== Int
n'   = (Int
i, Double
t2)
          | Bool
otherwise = (Int
i, Double
t0 forall a. Num a => a -> a -> a
+ (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) forall a. Num a => a -> a -> a
* Double
dt)
  in 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 forall a. Ord a => a -> a -> Bool
<= Int
0   = 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 forall a. Ord a => a -> a -> Bool
< Int
0
       then let t' :: Double
t' = Point -> Double
pointTime Point
p forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dn forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc
                n' :: Int
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor forall a b. (a -> b) -> a -> b
$ (Double
t' forall a. Num a => a -> a -> a
- Specs -> Double
spcStartTime Specs
sc) 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 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 }