{-# LANGUAGE TypeFamilies #-}

-- |
-- Module     : Simulation.Aivika.Trans.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.Trans.Internal.Specs
       (Specs(..),
        Method(..),
        Run(..),
        Point(..),
        basicTime,
        integIterationBnds,
        integIterationHiBnd,
        integIterationLoBnd,
        integPhaseBnds,
        integPhaseHiBnd,
        integPhaseLoBnd,
        integTimes,
        integPoints,
        integPointsStartingFrom,
        integStartPoint,
        integStopPoint,
        simulationStopPoint,
        timeGrid,
        pointAt,
        delayPoint) where

import Simulation.Aivika.Trans.Internal.Types

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

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

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

-- | Returns the last integration iteration.
integIterationHiBnd :: Specs m -> Int
integIterationHiBnd :: forall (m :: * -> *). Specs m -> Int
integIterationHiBnd Specs m
sc =
  let n :: Int
n = forall a b. (RealFrac a, Integral b) => a -> b
round ((forall (m :: * -> *). Specs m -> Double
spcStopTime Specs m
sc forall a. Num a => a -> a -> a
- 
                  forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc) forall a. Fractional a => a -> a -> a
/ forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc)
  in if Int
n forall a. Ord a => a -> a -> Bool
< Int
0
     then
       forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
       [Char]
"Either the simulation specs are incorrect, " forall a. [a] -> [a] -> [a]
++
       [Char]
"or a step time is too small, because of which " forall a. [a] -> [a] -> [a]
++
       [Char]
"a floating point overflow occurred on 32-bit Haskell implementation."
     else Int
n

-- | Returns the phases for the specified simulation specs starting from zero.
integPhases :: Specs m -> [Int]
integPhases :: forall (m :: * -> *). Specs m -> [Int]
integPhases Specs m
sc = 
  case forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
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 m -> (Int, Int)
integPhaseBnds :: forall (m :: * -> *). Specs m -> (Int, Int)
integPhaseBnds Specs m
sc = 
  case forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
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 m -> Int
integPhaseLoBnd :: forall (m :: * -> *). Specs m -> Int
integPhaseLoBnd Specs m
sc = Int
0
                  
-- | Returns the last integration phase, 0 for Euler's method, 1 for RK2 and 3 for RK4.
integPhaseHiBnd :: Specs m -> Int
integPhaseHiBnd :: forall (m :: * -> *). Specs m -> Int
integPhaseHiBnd Specs m
sc = 
  case forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
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 m -> Int -> Int -> Double
{-# INLINE basicTime #-}
basicTime :: forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
ph =
  if Int
ph forall a. Ord a => a -> a -> Bool
< Int
0 then 
    forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: basicTime"
  else
    forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc forall a. Num a => a -> a -> a
+ Double
n' forall a. Num a => a -> a -> a
* forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc forall a. Num a => a -> a -> a
+ Method -> Int -> Double
delta (forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
sc) Int
ph 
      where n' :: Double
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
            delta :: Method -> Int -> Double
delta Method
Euler       Int
0 = Double
0
            delta Method
RungeKutta2 Int
0 = Double
0
            delta Method
RungeKutta2 Int
1 = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
            delta Method
RungeKutta4 Int
0 = Double
0
            delta Method
RungeKutta4 Int
1 = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 Int
2 = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 Int
3 = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
            delta Method
RungeKutta4b Int
0 = Double
0
            delta Method
RungeKutta4b Int
1 = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc forall a. Fractional a => a -> a -> a
/ Double
3
            delta Method
RungeKutta4b Int
2 = Double
2 forall a. Num a => a -> a -> a
* forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc forall a. Fractional a => a -> a -> a
/ Double
3
            delta Method
RungeKutta4b Int
3 = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc

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

-- | Return the integration time points.
integPoints :: Run m -> [Point m]
integPoints :: forall (m :: * -> *). Run m -> [Point m]
integPoints Run m
r = [Point m]
points
  where sc :: Specs m
sc = forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        points :: [Point m]
points   = forall a b. (a -> b) -> [a] -> [b]
map Int -> Point m
point [Int
nl .. Int
nu]
        point :: Int -> Point m
point Int
n  = Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
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 m -> Point m
integStartPoint :: forall (m :: * -> *). Run m -> Point m
integStartPoint Run m
r = Int -> Point m
point Int
nl
  where sc :: Specs m
sc = forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        point :: Int -> Point m
point Int
n  = Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
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 m -> Point m
integStopPoint :: forall (m :: * -> *). Run m -> Point m
integStopPoint Run m
r = Int -> Point m
point Int
nu
  where sc :: Specs m
sc = forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        point :: Int -> Point m
point Int
n  = Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
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 m -> Point m
simulationStopPoint :: forall (m :: * -> *). Run m -> Point m
simulationStopPoint Run m
r = forall (m :: * -> *). Run m -> Double -> Int -> Point m
pointAt Run m
r (forall (m :: * -> *). Specs m -> Double
spcStopTime forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r) Int
0

-- | Return the point at the specified time.
pointAt :: Run m -> Double -> EventPriority -> Point m
{-# INLINABLE pointAt #-}
pointAt :: forall (m :: * -> *). Run m -> Double -> Int -> Point m
pointAt Run m
r Double
t Int
priority = Point m
p
  where sc :: Specs m
sc = forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        t0 :: Double
t0 = forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc
        dt :: Double
dt = forall (m :: * -> *). Specs m -> Double
spcDT Specs m
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 m
p = Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                    pointRun :: Run m
pointRun = Run m
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 m -> [Point m]
integPointsStartingFrom :: forall (m :: * -> *). Point m -> [Point m]
integPointsStartingFrom Point m
p = [Point m]
points
  where r :: Run m
r  = forall (m :: * -> *). Point m -> Run m
pointRun Point m
p
        sc :: Specs m
sc = forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        n0 :: Int
n0       = if forall (m :: * -> *). Point m -> Int
pointPhase Point m
p forall a. Eq a => a -> a -> Bool
== Int
0
                   then forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
                   else forall (m :: * -> *). Point m -> Int
pointIteration Point m
p forall a. Num a => a -> a -> a
+ Int
1
        points :: [Point m]
points   = forall a b. (a -> b) -> [a] -> [b]
map Int -> Point m
point [Int
n0 .. Int
nu]
        point :: Int -> Point m
point Int
n  = Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
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 m -> Int -> [(Int, Double)]
timeGrid :: forall (m :: * -> *). Specs m -> Int -> [(Int, Double)]
timeGrid Specs m
sc Int
n =
  let t0 :: Double
t0 = forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc
      t2 :: Double
t2 = forall (m :: * -> *). Specs m -> Double
spcStopTime Specs m
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 m -> Int -> Point m
delayPoint :: forall (m :: * -> *). Point m -> Int -> Point m
delayPoint Point m
p Int
dn
  | Int
dn forall a. Ord a => a -> a -> Bool
<= Int
0   = forall a. HasCallStack => [Char] -> a
error [Char]
"Expected the positive number of iterations: delayPoint"
  | Bool
otherwise =
    let sc :: Specs m
sc = forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
        n :: Int
n  = forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
        ph :: Int
ph = forall (m :: * -> *). Point m -> Int
pointPhase Point m
p
    in if Int
ph forall a. Ord a => a -> a -> Bool
< Int
0
       then let t' :: Double
t' = forall (m :: * -> *). Point m -> Double
pointTime Point m
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
* forall (m :: * -> *). Specs m -> Double
spcDT Specs m
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
- forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc) forall a. Fractional a => a -> a -> a
/ forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
            in Point m
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' = forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n' Int
ph
            in Point m
p { pointTime :: Double
pointTime = Double
t',
                   pointIteration :: Int
pointIteration = Int
n',
                   pointPhase :: Int
pointPhase = Int
ph }