```
{-# LANGUAGE RecursiveDo #-}

-- |
-- Module     : Simulation.Aivika.Dynamics.Extra
-- 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
```