{-# LANGUAGE RecursiveDo #-}

-- |
-- Module     : Simulation.Aivika.Dynamics.Extra
-- 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
--
-- This module defines auxiliary functions such as interpolation ones
-- that complement the memoization, for example. There are scan functions too.
--

module Simulation.Aivika.Dynamics.Extra
       (-- * Interpolation
        initDynamics,
        discreteDynamics,
        interpolateDynamics,
        -- * Scans
        scanDynamics,
        scan1Dynamics) where

import Simulation.Aivika.Internal.Specs
import Simulation.Aivika.Internal.Simulation
import Simulation.Aivika.Internal.Dynamics

-- | Return the initial value.
initDynamics :: Dynamics a -> Dynamics a
{-# INLINE initDynamics #-}
initDynamics :: Dynamics a -> Dynamics a
initDynamics (Dynamics Point -> IO a
m) =
  (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO a) -> Dynamics a) -> (Point -> IO a) -> Dynamics a
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
  in Point -> IO a
m (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
0 Int
0,
             pointIteration :: Int
pointIteration = Int
0,
             pointPhase :: Int
pointPhase = Int
0 }

-- | Discretize the computation in the integration time points.
discreteDynamics :: Dynamics a -> Dynamics a
{-# INLINE discreteDynamics #-}
discreteDynamics :: Dynamics a -> Dynamics a
discreteDynamics (Dynamics Point -> IO a
m) =
  (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO a) -> Dynamics a) -> (Point -> IO a) -> Dynamics a
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  if Point -> Int
pointPhase Point
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then
    Point -> IO a
m Point
p
  else
    let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
        n :: Int
n  = Point -> Int
pointIteration Point
p
    in Point -> IO a
m (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
               pointPhase :: Int
pointPhase = Int
0 }

-- | Interpolate the computation based on the integration time points only.
-- Unlike the 'discreteDynamics' function it knows about the intermediate 
-- time points that are used in the Runge-Kutta method.
interpolateDynamics :: Dynamics a -> Dynamics a
{-# INLINE interpolateDynamics #-}
interpolateDynamics :: Dynamics a -> Dynamics a
interpolateDynamics (Dynamics Point -> IO a
m) = 
  (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO a) -> Dynamics a) -> (Point -> IO a) -> Dynamics a
forall a b. (a -> b) -> a -> b
$ \Point
p -> 
  if Point -> Int
pointPhase Point
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then 
    Point -> IO a
m Point
p
  else 
    let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
        n :: Int
n  = Point -> Int
pointIteration Point
p
    in Point -> IO a
m (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
               pointPhase :: Int
pointPhase = Int
0 }

-- | Like the standard 'scanl1' function but applied to values in 
-- the integration time points. The accumulator values are transformed
-- according to the second argument, which should be either function 
-- 'memo0Dynamics' or its unboxed version.
scan1Dynamics :: (a -> a -> a)
                 -> (Dynamics a -> Simulation (Dynamics a))
                 -> (Dynamics a -> Simulation (Dynamics a))
scan1Dynamics :: (a -> a -> a)
-> (Dynamics a -> Simulation (Dynamics a))
-> Dynamics a
-> Simulation (Dynamics a)
scan1Dynamics a -> a -> a
f Dynamics a -> Simulation (Dynamics a)
tr Dynamics a
m =
  mdo Dynamics a
y <- Dynamics a -> Simulation (Dynamics a)
tr (Dynamics a -> Simulation (Dynamics a))
-> Dynamics a -> Simulation (Dynamics a)
forall a b. (a -> b) -> a -> b
$ (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO a) -> Dynamics a) -> (Point -> IO a) -> Dynamics a
forall a b. (a -> b) -> a -> b
$ \Point
p ->
        case Point -> Int
pointIteration Point
p of
          Int
0 -> 
            Point -> Dynamics a -> IO a
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics a
m
          Int
n -> do 
            let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
                ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
            a
s <- Point -> Dynamics a -> IO a
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
py Dynamics a
y
            a
x <- Point -> Dynamics a -> IO a
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics a
m
            a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> a -> IO a
forall a b. (a -> b) -> a -> b
$! a -> a -> a
f a
s a
x
      Dynamics a -> Simulation (Dynamics a)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics a
y

-- | Like the standard 'scanl' function but applied to values in 
-- the integration time points. The accumulator values are transformed
-- according to the third argument, which should be either function
-- 'memo0Dynamics' or its unboxed version.
scanDynamics :: (a -> b -> a)
                -> a
                -> (Dynamics a -> Simulation (Dynamics a))
                -> (Dynamics b -> Simulation (Dynamics a))
scanDynamics :: (a -> b -> a)
-> a
-> (Dynamics a -> Simulation (Dynamics a))
-> Dynamics b
-> Simulation (Dynamics a)
scanDynamics a -> b -> a
f a
acc Dynamics a -> Simulation (Dynamics a)
tr Dynamics b
m =
  mdo Dynamics a
y <- Dynamics a -> Simulation (Dynamics a)
tr (Dynamics a -> Simulation (Dynamics a))
-> Dynamics a -> Simulation (Dynamics a)
forall a b. (a -> b) -> a -> b
$ (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO a) -> Dynamics a) -> (Point -> IO a) -> Dynamics a
forall a b. (a -> b) -> a -> b
$ \Point
p ->
        case Point -> Int
pointIteration Point
p of
          Int
0 -> do
            b
x <- Point -> Dynamics b -> IO b
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics b
m
            a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> a -> IO a
forall a b. (a -> b) -> a -> b
$! a -> b -> a
f a
acc b
x
          Int
n -> do 
            let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
                ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
            a
s <- Point -> Dynamics a -> IO a
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
py Dynamics a
y
            b
x <- Point -> Dynamics b -> IO b
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics b
m
            a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> a -> IO a
forall a b. (a -> b) -> a -> b
$! a -> b -> a
f a
s b
x
      Dynamics a -> Simulation (Dynamics a)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics a
y