{-# LANGUAGE BangPatterns, RecursiveDo, FlexibleContexts #-}

-- |
-- Module     : Simulation.Aivika.Trans.SystemDynamics
-- 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 integrals and other functions of System Dynamics.
--

module Simulation.Aivika.Trans.SystemDynamics
       (-- * Equality and Ordering
        (.==.),
        (./=.),
        (.<.),
        (.>=.),
        (.>.),
        (.<=.),
        maxDynamics,
        minDynamics,
        ifDynamics,
        -- * Ordinary Differential Equations
        integ,
        integEither,
        smoothI,
        smooth,
        smooth3I,
        smooth3,
        smoothNI,
        smoothN,
        delay1I,
        delay1,
        delay3I,
        delay3,
        delayNI,
        delayN,
        forecast,
        trend,
        -- * Difference Equations
        diffsum,
        diffsumEither,
        -- * Table Functions
        lookupDynamics,
        lookupStepwiseDynamics,
        -- * Discrete Functions
        delay,
        delayI,
        delayByDT,
        delayIByDT,
        step,
        pulse,
        pulseP,
        ramp,
        -- * Financial Functions
        npv,
        npve) where

import Data.Array

import Control.Monad
import Control.Monad.Trans
import Control.Monad.Fix

import Simulation.Aivika.Trans.Internal.Specs
import Simulation.Aivika.Trans.Internal.Parameter
import Simulation.Aivika.Trans.Internal.Simulation
import Simulation.Aivika.Trans.Internal.Dynamics
import Simulation.Aivika.Trans.Dynamics.Extra
import Simulation.Aivika.Trans.Table
import Simulation.Aivika.Trans.SD

import qualified Simulation.Aivika.Trans.Dynamics.Memo as M
import qualified Simulation.Aivika.Trans.Dynamics.Memo.Unboxed as MU

--
-- Equality and Ordering
--

-- | Compare for equality.
(.==.) :: (Monad m, Eq a) => Dynamics m a -> Dynamics m a -> Dynamics m Bool
{-# INLINE (.==.) #-}
.==. :: Dynamics m a -> Dynamics m a -> Dynamics m Bool
(.==.) = (a -> a -> Bool) -> Dynamics m a -> Dynamics m a -> Dynamics m Bool
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(==)

-- | Compare for inequality.
(./=.) :: (Monad m, Eq a) => Dynamics m a -> Dynamics m a -> Dynamics m Bool
{-# INLINE (./=.) #-}
./=. :: Dynamics m a -> Dynamics m a -> Dynamics m Bool
(./=.) = (a -> a -> Bool) -> Dynamics m a -> Dynamics m a -> Dynamics m Bool
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(/=)

-- | Compare for ordering.
(.<.) :: (Monad m, Ord a) => Dynamics m a -> Dynamics m a -> Dynamics m Bool
{-# INLINE (.<.) #-}
.<. :: Dynamics m a -> Dynamics m a -> Dynamics m Bool
(.<.) = (a -> a -> Bool) -> Dynamics m a -> Dynamics m a -> Dynamics m Bool
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(<)

-- | Compare for ordering.
(.>=.) :: (Monad m, Ord a) => Dynamics m a -> Dynamics m a -> Dynamics m Bool
{-# INLINE (.>=.) #-}
.>=. :: Dynamics m a -> Dynamics m a -> Dynamics m Bool
(.>=.) = (a -> a -> Bool) -> Dynamics m a -> Dynamics m a -> Dynamics m Bool
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(>=)

-- | Compare for ordering.
(.>.) :: (Monad m, Ord a) => Dynamics m a -> Dynamics m a -> Dynamics m Bool
{-# INLINE (.>.) #-}
.>. :: Dynamics m a -> Dynamics m a -> Dynamics m Bool
(.>.) = (a -> a -> Bool) -> Dynamics m a -> Dynamics m a -> Dynamics m Bool
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(>)

-- | Compare for ordering.
(.<=.) :: (Monad m, Ord a) => Dynamics m a -> Dynamics m a -> Dynamics m Bool
{-# INLINE (.<=.) #-}
.<=. :: Dynamics m a -> Dynamics m a -> Dynamics m Bool
(.<=.) = (a -> a -> Bool) -> Dynamics m a -> Dynamics m a -> Dynamics m Bool
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(<=)

-- | Return the maximum.
maxDynamics :: (Monad m, Ord a) => Dynamics m a -> Dynamics m a -> Dynamics m a
{-# INLINE maxDynamics #-}
maxDynamics :: Dynamics m a -> Dynamics m a -> Dynamics m a
maxDynamics = (a -> a -> a) -> Dynamics m a -> Dynamics m a -> Dynamics m a
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> a
forall a. Ord a => a -> a -> a
max

-- | Return the minimum.
minDynamics :: (Monad m, Ord a) => Dynamics m a -> Dynamics m a -> Dynamics m a
{-# INLINE minDynamics #-}
minDynamics :: Dynamics m a -> Dynamics m a -> Dynamics m a
minDynamics = (a -> a -> a) -> Dynamics m a -> Dynamics m a -> Dynamics m a
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 a -> a -> a
forall a. Ord a => a -> a -> a
min

-- | Implement the if-then-else operator.
ifDynamics :: Monad m => Dynamics m Bool -> Dynamics m a -> Dynamics m a -> Dynamics m a
{-# INLINE ifDynamics #-}
ifDynamics :: Dynamics m Bool -> Dynamics m a -> Dynamics m a -> Dynamics m a
ifDynamics Dynamics m Bool
cond Dynamics m a
x Dynamics m a
y =
  do Bool
a <- Dynamics m Bool
cond
     if Bool
a then Dynamics m a
x else Dynamics m a
y

--
-- Ordinary Differential Equations
--

integEuler :: Monad m
              => Dynamics m Double
              -> Dynamics m Double 
              -> Dynamics m Double 
              -> Point m
              -> m Double
{-# INLINABLE integEuler #-}
integEuler :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEuler (Dynamics Point m -> m Double
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p = 
  case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
    Int
0 -> 
      Point m -> m Double
i Point m
p
    Int
n -> do 
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
          py :: Point m
py = Point m
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 }
      Double
a <- Point m -> m Double
y Point m
py
      Double
b <- Point m -> m Double
f Point m
py
      let !v :: Double
v = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT (Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v

integRK2 :: Monad m
            => Dynamics m Double
            -> Dynamics m Double
            -> Dynamics m Double
            -> Point m
            -> m Double
{-# INLINABLE integRK2 #-}
integRK2 :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK2 (Dynamics Point m -> m Double
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p =
  case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointPhase Point m
p of
    Int
0 -> case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
      Int
0 ->
        Point m -> m Double
i Point m
p
      Int
n -> do
        let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
            ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
            py :: Point m
py = Point m
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 }
            p1 :: Point m
p1 = Point m
py
            p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
        Double
vy <- Point m -> m Double
y Point m
py
        Double
k1 <- Point m -> m Double
f Point m
p1
        Double
k2 <- Point m -> m Double
f Point m
p2
        let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k2)
        Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point m
p1 = Point m
py
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k1 <- Point m -> m Double
f Point m
p1
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k1
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      [Char] -> m Double
forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: integRK2"

integRK4 :: Monad m
            => Dynamics m Double
            -> Dynamics m Double
            -> Dynamics m Double
            -> Point m
            -> m Double
{-# INLINABLE integRK4 #-}
integRK4 :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4 (Dynamics Point m -> m Double
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p =
  case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointPhase Point m
p of
    Int
0 -> case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
      Int
0 -> 
        Point m -> m Double
i Point m
p
      Int
n -> do
        let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
            ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
            t3 :: Double
t3 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
2
            t4 :: Double
t4 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
3
            py :: Point m
py = Point m
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 }
            p1 :: Point m
p1 = Point m
py
            p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
            p3 :: Point m
p3 = Point m
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
2 }
            p4 :: Point m
p4 = Point m
p { pointTime :: Double
pointTime = Double
t4, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
3 }
        Double
vy <- Point m -> m Double
y Point m
py
        Double
k1 <- Point m -> m Double
f Point m
p1
        Double
k2 <- Point m -> m Double
f Point m
p2
        Double
k3 <- Point m -> m Double
f Point m
p3
        Double
k4 <- Point m -> m Double
f Point m
p4
        let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k4)
        Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point m
p1 = Point m
py
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k1 <- Point m -> m Double
f Point m
p1
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k1
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
2 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
1
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k2 <- Point m -> m Double
f Point m
p2
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k2
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
3 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t3 :: Double
t3 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
2
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p3 :: Point m
p3 = Point m
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
2 }
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k3 <- Point m -> m Double
f Point m
p3
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k3
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      [Char] -> m Double
forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: integRK4"

integRK4b :: Monad m
             => Dynamics m Double
             -> Dynamics m Double
             -> Dynamics m Double
             -> Point m
             -> m Double
{-# INLINABLE integRK4b #-}
integRK4b :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4b (Dynamics Point m -> m Double
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p =
  case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointPhase Point m
p of
    Int
0 -> case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
      Int
0 -> 
        Point m -> m Double
i Point m
p
      Int
n -> do
        let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
            ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
            t3 :: Double
t3 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
2
            t4 :: Double
t4 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
3
            py :: Point m
py = Point m
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 }
            p1 :: Point m
p1 = Point m
py
            p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
            p3 :: Point m
p3 = Point m
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
2 }
            p4 :: Point m
p4 = Point m
p { pointTime :: Double
pointTime = Double
t4, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
3 }
        Double
vy <- Point m -> m Double
y Point m
py
        Double
k1 <- Point m -> m Double
f Point m
p1
        Double
k2 <- Point m -> m Double
f Point m
p2
        Double
k3 <- Point m -> m Double
f Point m
p3
        Double
k4 <- Point m -> m Double
f Point m
p4
        let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k3) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k4)
        Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point m
p1 = Point m
py
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k1 <- Point m -> m Double
f Point m
p1
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k1
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
2 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
1
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point m
p1 = Point m
py
          p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k1 <- Point m -> m Double
f Point m
p1
      Double
k2 <- Point m -> m Double
f Point m
p2
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
* (- Double
k1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k2)
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
3 -> do
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
1
          t3 :: Double
t3 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
2
          py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point m
p1 = Point m
py
          p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
          p3 :: Point m
p3 = Point m
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
2 }
      Double
vy <- Point m -> m Double
y Point m
py
      Double
k1 <- Point m -> m Double
f Point m
p1
      Double
k2 <- Point m -> m Double
f Point m
p2
      Double
k3 <- Point m -> m Double
f Point m
p3
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k3)
      Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      [Char] -> m Double
forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: integRK4b"

-- | Return an integral with the specified derivative and initial value.
--
-- To create a loopback, you should use the recursive do-notation.
-- It allows defining the differential equations unordered as
-- in mathematics:
--
-- @
-- model = 
--   mdo a <- integ (- ka * a) 100
--       b <- integ (ka * a - kb * b) 0
--       c <- integ (kb * b) 0
--       let ka = 1
--           kb = 1
--       runDynamicsInStopTime $ sequence [a, b, c]
-- @
integ :: (MonadSD m, MonadFix m)
         => Dynamics m Double                  -- ^ the derivative
         -> Dynamics m Double                  -- ^ the initial value
         -> Simulation m (Dynamics m Double)   -- ^ the integral
{-# INLINABLE integ #-}
integ :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ Dynamics m Double
diff Dynamics m Double
i =
  mdo Dynamics m Double
y <- Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) e.
MonadMemo m e =>
Dynamics m e -> Simulation m (Dynamics m e)
MU.memoDynamics Dynamics m Double
z
      Dynamics m Double
z <- (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. (Run m -> m a) -> Simulation m a
Simulation ((Run m -> m (Dynamics m Double))
 -> Simulation m (Dynamics m Double))
-> (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ \Run m
r ->
        case Specs m -> Method
forall (m :: * -> *). Specs m -> Method
spcMethod (Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r) of
          Method
Euler -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Monad m =>
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEuler Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
          Method
RungeKutta2 -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Monad m =>
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK2 Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
          Method
RungeKutta4 -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Monad m =>
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4 Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
          Method
RungeKutta4b -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Monad m =>
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4b Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

integEulerEither :: Monad m
                    => Dynamics m (Either Double Double)
                    -> Dynamics m Double 
                    -> Dynamics m Double 
                    -> Point m
                    -> m Double
{-# INLINABLE integEulerEither #-}
integEulerEither :: Dynamics m (Either Double Double)
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEulerEither (Dynamics Point m -> m (Either Double Double)
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p = 
  case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
    Int
0 -> 
      Point m -> m Double
i Point m
p
    Int
n -> do 
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
          py :: Point m
py = Point m
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 }
      Either Double Double
b <- Point m -> m (Either Double Double)
f Point m
py
      case Either Double Double
b of
        Left Double
v ->
          Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
        Right Double
b -> do
          Double
a <- Point m -> m Double
y Point m
py
          let !v :: Double
v = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT (Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b
          Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v

-- | Like 'integ' but allows either setting a new 'Left' integral value,
-- or integrating using the 'Right' derivative directly within computation.
--
-- This function always uses Euler's method.
integEither :: (MonadSD m, MonadFix m)
               => Dynamics m (Either Double Double)
               -- ^ either set a new 'Left' integral value, or use a 'Right' derivative
               -> Dynamics m Double
               -- ^ the initial value
               -> Simulation m (Dynamics m Double)
{-# INLINABLE integEither #-}
integEither :: Dynamics m (Either Double Double)
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integEither Dynamics m (Either Double Double)
diff Dynamics m Double
i =
  mdo Dynamics m Double
y <- Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) e.
MonadMemo m e =>
Dynamics m e -> Simulation m (Dynamics m e)
MU.memoDynamics Dynamics m Double
z
      Dynamics m Double
z <- (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. (Run m -> m a) -> Simulation m a
Simulation ((Run m -> m (Dynamics m Double))
 -> Simulation m (Dynamics m Double))
-> (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ \Run m
r ->
        Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ Dynamics m (Either Double Double)
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Monad m =>
Dynamics m (Either Double Double)
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEulerEither Dynamics m (Either Double Double)
diff Dynamics m Double
i Dynamics m Double
y
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

-- | Return the first order exponential smooth.
--
-- To create a loopback, you should use the recursive do-notation
-- with help of which the function itself is defined:
--
-- @
-- smoothI x t i =
--   mdo y <- integ ((x - y) \/ t) i
--       return y
-- @     
smoothI :: (MonadSD m, MonadFix m)
           => Dynamics m Double                  -- ^ the value to smooth over time
           -> Dynamics m Double                  -- ^ time
           -> Dynamics m Double                  -- ^ the initial value
           -> Simulation m (Dynamics m Double)   -- ^ the first order exponential smooth
{-# INLINABLE smoothI #-}
smoothI :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothI Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
  mdo Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t) Dynamics m Double
i
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

-- | Return the first order exponential smooth.
--
-- This is a simplified version of the 'smoothI' function
-- without specifing the initial value.
smooth :: (MonadSD m, MonadFix m)
          => Dynamics m Double                  -- ^ the value to smooth over time
          -> Dynamics m Double                  -- ^ time
          -> Simulation m (Dynamics m Double)   -- ^ the first order exponential smooth
{-# INLINABLE smooth #-}
smooth :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
smooth Dynamics m Double
x Dynamics m Double
t = Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothI Dynamics m Double
x Dynamics m Double
t Dynamics m Double
x

-- | Return the third order exponential smooth.
--
-- To create a loopback, you should use the recursive do-notation
-- with help of which the function itself is defined:
--
-- @
-- smooth3I x t i =
--   mdo y  <- integ ((s2 - y) \/ t') i
--       s2 <- integ ((s1 - s2) \/ t') i
--       s1 <- integ ((x - s1) \/ t') i
--       let t' = t \/ 3.0
--       return y
-- @     
smooth3I :: (MonadSD m, MonadFix m)
            => Dynamics m Double                  -- ^ the value to smooth over time
            -> Dynamics m Double                  -- ^ time
            -> Dynamics m Double                  -- ^ the initial value
            -> Simulation m (Dynamics m Double)   -- ^ the third order exponential smooth
{-# INLINABLE smooth3I #-}
smooth3I :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smooth3I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
  mdo Dynamics m Double
y  <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
s2 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
      Dynamics m Double
s2 <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
s1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
s2) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
      Dynamics m Double
s1 <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
s1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
      let t' :: Dynamics m Double
t' = Dynamics m Double
t Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
3.0
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

-- | Return the third order exponential smooth.
-- 
-- This is a simplified version of the 'smooth3I' function
-- without specifying the initial value.
smooth3 :: (MonadSD m, MonadFix m)
           => Dynamics m Double                  -- ^ the value to smooth over time
           -> Dynamics m Double                  -- ^ time
           -> Simulation m (Dynamics m Double)   -- ^ the third order exponential smooth
{-# INLINABLE smooth3 #-}
smooth3 :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
smooth3 Dynamics m Double
x Dynamics m Double
t = Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smooth3I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
x

-- | Return the n'th order exponential smooth.
--
-- The result is not discrete in that sense that it may change within the integration time
-- interval depending on the integration method used. Probably, you should apply
-- the 'discreteDynamics' function to the result if you want to achieve an effect when
-- the value is not changed within the time interval, which is used sometimes.
smoothNI :: (MonadSD m, MonadFix m)
            => Dynamics m Double                  -- ^ the value to smooth over time
            -> Dynamics m Double                  -- ^ time
            -> Int                                -- ^ the order
            -> Dynamics m Double                  -- ^ the initial value
            -> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential smooth
{-# INLINABLE smoothNI #-}
smoothNI :: Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m Double
i =
  mdo [Dynamics m Double]
s <- [Int]
-> (Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
1 .. Int
n] ((Int -> Simulation m (Dynamics m Double))
 -> Simulation m [Dynamics m Double])
-> (Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double]
forall a b. (a -> b) -> a -> b
$ \Int
k ->
        if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
        then Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
        else Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
k) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
      let a :: Array Int (Dynamics m Double)
a  = (Int, Int) -> [Dynamics m Double] -> Array Int (Dynamics m Double)
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics m Double]
s 
          t' :: Dynamics m Double
t' = Dynamics m Double
t Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Int -> Dynamics m Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
n

-- | Return the n'th order exponential smooth.
--
-- This is a simplified version of the 'smoothNI' function
-- without specifying the initial value.
smoothN :: (MonadSD m, MonadFix m)
           => Dynamics m Double                  -- ^ the value to smooth over time
           -> Dynamics m Double                  -- ^ time
           -> Int                                -- ^ the order
           -> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential smooth
{-# INLINABLE smoothN #-}
smoothN :: Dynamics m Double
-> Dynamics m Double -> Int -> Simulation m (Dynamics m Double)
smoothN Dynamics m Double
x Dynamics m Double
t Int
n = Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m Double
x

-- | Return the first order exponential delay.
--
-- To create a loopback, you should use the recursive do-notation
-- with help of which the function itself is defined:
--
-- @
-- delay1I x t i =
--   mdo y <- integ (x - y \/ t) (i * t)
--       return $ y \/ t
-- @     
delay1I :: (MonadSD m, MonadFix m)
           => Dynamics m Double                  -- ^ the value to conserve
           -> Dynamics m Double                  -- ^ time
           -> Dynamics m Double                  -- ^ the initial value
           -> Simulation m (Dynamics m Double)   -- ^ the first order exponential delay
{-# INLINABLE delay1I #-}
delay1I :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delay1I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
  mdo Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t) (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t)
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t

-- | Return the first order exponential delay.
--
-- This is a simplified version of the 'delay1I' function
-- without specifying the initial value.
delay1 :: (MonadSD m, MonadFix m)
          => Dynamics m Double                  -- ^ the value to conserve
          -> Dynamics m Double                  -- ^ time
          -> Simulation m (Dynamics m Double)   -- ^ the first order exponential delay
{-# INLINABLE delay1 #-}
delay1 :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
delay1 Dynamics m Double
x Dynamics m Double
t = Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delay1I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
x

-- | Return the third order exponential delay.
delay3I :: (MonadSD m, MonadFix m)
           => Dynamics m Double                  -- ^ the value to conserve
           -> Dynamics m Double                  -- ^ time
           -> Dynamics m Double                  -- ^ the initial value
           -> Simulation m (Dynamics m Double)   -- ^ the third order exponential delay
{-# INLINABLE delay3I #-}
delay3I :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delay3I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
  mdo Dynamics m Double
y  <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
s2 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t')
      Dynamics m Double
s2 <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
s1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
s2 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t')
      Dynamics m Double
s1 <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
s1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t')
      let t' :: Dynamics m Double
t' = Dynamics m Double
t Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
3.0
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t'         

-- | Return the third order exponential delay.
--
-- This is a simplified version of the 'delay3I' function
-- without specifying the initial value.
delay3 :: (MonadSD m, MonadFix m)
          => Dynamics m Double                  -- ^ the value to conserve
          -> Dynamics m Double                  -- ^ time
          -> Simulation m (Dynamics m Double)   -- ^ the third order exponential delay
{-# INLINABLE delay3 #-}
delay3 :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
delay3 Dynamics m Double
x Dynamics m Double
t = Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delay3I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
x

-- | Return the n'th order exponential delay.
delayNI :: (MonadSD m, MonadFix m)
           => Dynamics m Double                  -- ^ the value to conserve
           -> Dynamics m Double                  -- ^ time
           -> Int                                -- ^ the order
           -> Dynamics m Double                  -- ^ the initial value
           -> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential delay
{-# INLINABLE delayNI #-}
delayNI :: Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delayNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m Double
i =
  mdo [Dynamics m Double]
s <- [Int]
-> (Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
1 .. Int
n] ((Int -> Simulation m (Dynamics m Double))
 -> Simulation m [Dynamics m Double])
-> (Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double]
forall a b. (a -> b) -> a -> b
$ \Int
k ->
        if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
        then Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- (Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t')
        else Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- (Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
k) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t')
      let a :: Array Int (Dynamics m Double)
a  = (Int, Int) -> [Dynamics m Double] -> Array Int (Dynamics m Double)
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics m Double]
s
          t' :: Dynamics m Double
t' = Dynamics m Double
t Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Int -> Dynamics m Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
n) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t'

-- | Return the n'th order exponential delay.
--
-- This is a simplified version of the 'delayNI' function
-- without specifying the initial value.
delayN :: (MonadSD m, MonadFix m)
          => Dynamics m Double                  -- ^ the value to conserve
          -> Dynamics m Double                  -- ^ time
          -> Int                                -- ^ the order
          -> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential delay
{-# INLINABLE delayN #-}
delayN :: Dynamics m Double
-> Dynamics m Double -> Int -> Simulation m (Dynamics m Double)
delayN Dynamics m Double
x Dynamics m Double
t Int
n = Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delayNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m Double
x

-- | Return the forecast.
--
-- The function has the following definition:
--
-- @
-- forecast x at hz =
--   do y <- smooth x at
--      return $ x * (1.0 + (x \/ y - 1.0) \/ at * hz)
-- @
forecast :: (MonadSD m, MonadFix m)
            => Dynamics m Double                  -- ^ the value to forecast
            -> Dynamics m Double                  -- ^ the average time
            -> Dynamics m Double                  -- ^ the time horizon
            -> Simulation m (Dynamics m Double)   -- ^ the forecast
{-# INLINABLE forecast #-}
forecast :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forecast Dynamics m Double
x Dynamics m Double
at Dynamics m Double
hz =
  do Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
smooth Dynamics m Double
x Dynamics m Double
at
     Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* (Dynamics m Double
1.0 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
1.0) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
at Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
hz)

-- | Return the trend.
--
-- The function has the following definition:
--
-- @
-- trend x at i =
--   do y <- smoothI x at (x \/ (1.0 + i * at))
--      return $ (x \/ y - 1.0) \/ at
-- @
trend :: (MonadSD m, MonadFix m)
         => Dynamics m Double                  -- ^ the value for which the trend is calculated
         -> Dynamics m Double                  -- ^ the average time
         -> Dynamics m Double                  -- ^ the initial value
         -> Simulation m (Dynamics m Double)   -- ^ the fractional change rate
{-# INLINABLE trend #-}
trend :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
trend Dynamics m Double
x Dynamics m Double
at Dynamics m Double
i =
  do Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothI Dynamics m Double
x Dynamics m Double
at (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ (Dynamics m Double
1.0 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
at))
     Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
1.0) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
at

--
-- Difference Equations
--

-- | Retun the sum for the difference equation.
-- It is like an integral returned by the 'integ' function, only now
-- the difference is used instead of derivative.
--
-- As usual, to create a loopback, you should use the recursive do-notation.
diffsum :: (MonadSD m, MonadFix m,
            MU.MonadMemo m a, Num a)
           => Dynamics m a                  -- ^ the difference
           -> Dynamics m a                  -- ^ the initial value
           -> Simulation m (Dynamics m a)   -- ^ the sum
{-# INLINABLE diffsum #-}
diffsum :: Dynamics m a -> Dynamics m a -> Simulation m (Dynamics m a)
diffsum (Dynamics Point m -> m a
diff) (Dynamics Point m -> m a
i) =
  mdo Dynamics m a
y <-
        Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
MonadMemo m e =>
Dynamics m e -> Simulation m (Dynamics m e)
MU.memo0Dynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
$
        (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m a) -> Dynamics m a)
-> (Point m -> m a) -> Dynamics m a
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
        case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
          Int
0 -> Point m -> m a
i Point m
p
          Int
n -> do 
            let Dynamics Point m -> m a
m = Dynamics m a
y
                sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
                ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point m
py = Point m
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
a <- Point m -> m a
m Point m
py
            a
b <- Point m -> m a
diff Point m
py
            let !v :: a
v = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
b
            a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return a
v
      Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m a
y

-- | Like 'diffsum' but allows either setting a new 'Left' sum value, or adding the 'Right' difference.
diffsumEither :: (MonadSD m, MonadFix m,
                  MU.MonadMemo m a, Num a)
                 => Dynamics m (Either a a)
                 -- ^ either set the 'Left' value for the sum, or add the 'Right' difference to the sum
                 -> Dynamics m a
                 -- ^ the initial value
                 -> Simulation m (Dynamics m a)
                 -- ^ the sum
{-# INLINABLE diffsumEither #-}
diffsumEither :: Dynamics m (Either a a)
-> Dynamics m a -> Simulation m (Dynamics m a)
diffsumEither (Dynamics Point m -> m (Either a a)
diff) (Dynamics Point m -> m a
i) =
  mdo Dynamics m a
y <-
        Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
MonadMemo m e =>
Dynamics m e -> Simulation m (Dynamics m e)
MU.memo0Dynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
$
        (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m a) -> Dynamics m a)
-> (Point m -> m a) -> Dynamics m a
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
        case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
          Int
0 -> Point m -> m a
i Point m
p
          Int
n -> do 
            let Dynamics Point m -> m a
m = Dynamics m a
y
                sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
                ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point m
py = Point m
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 }
            Either a a
b <- Point m -> m (Either a a)
diff Point m
py
            case Either a a
b of
              Left a
v ->
                a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return a
v
              Right a
b -> do
                a
a <- Point m -> m a
m Point m
py
                let !v :: a
v = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
b
                a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return a
v
      Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m a
y

--
-- Table Functions
--

-- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.
lookupDynamics :: Monad m => Dynamics m Double -> Array Int (Double, Double) -> Dynamics m Double
{-# INLINABLE lookupDynamics #-}
lookupDynamics :: Dynamics m Double
-> Array Int (Double, Double) -> Dynamics m Double
lookupDynamics (Dynamics Point m -> m Double
m) Array Int (Double, Double)
tbl =
  (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
  do Double
a <- Point m -> m Double
m Point m
p
     Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ Double -> Array Int (Double, Double) -> Double
tableLookup Double
a Array Int (Double, Double)
tbl

-- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function.
lookupStepwiseDynamics :: Monad m => Dynamics m Double -> Array Int (Double, Double) -> Dynamics m Double
{-# INLINABLE lookupStepwiseDynamics #-}
lookupStepwiseDynamics :: Dynamics m Double
-> Array Int (Double, Double) -> Dynamics m Double
lookupStepwiseDynamics (Dynamics Point m -> m Double
m) Array Int (Double, Double)
tbl =
  (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
  do Double
a <- Point m -> m Double
m Point m
p
     Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ Double -> Array Int (Double, Double) -> Double
tableLookupStepwise Double
a Array Int (Double, Double)
tbl

--
-- Discrete Functions
--

-- | Return the delayed value using the specified lag time.
-- This function is less accurate than 'delayByDT'.
delay :: Monad m
         => Dynamics m a          -- ^ the value to delay
         -> Dynamics m Double     -- ^ the lag time
         -> Dynamics m a          -- ^ the delayed value
{-# INLINABLE delay #-}
delay :: Dynamics m a -> Dynamics m Double -> Dynamics m a
delay (Dynamics Point m -> m a
x) (Dynamics Point m -> m Double
d) = Dynamics m a -> Dynamics m a
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m a -> Dynamics m a) -> Dynamics m a -> Dynamics m a
forall a b. (a -> b) -> a -> b
$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r 
  where
    r :: Point m -> m a
r Point m
p = do 
      let t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
          sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
      Double
a <- Point m -> m Double
d Point m
p
      let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a
          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 m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
          y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
$ Point m
p { pointTime :: Double
pointTime = Double
t',
                                  pointIteration :: Int
pointIteration = Int
n',
                                  pointPhase :: Int
pointPhase = -Int
1 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      m a
y

-- | Return the delayed value using the specified lag time and initial value.
-- Because of the latter, it allows creating a loop back.
-- This function is less accurate than 'delayIByDT'.
delayI :: MonadSD m
          => Dynamics m a                    -- ^ the value to delay
          -> Dynamics m Double               -- ^ the lag time
          -> Dynamics m a                    -- ^ the initial value
          -> Simulation m (Dynamics m a)     -- ^ the delayed value
{-# INLINABLE delayI #-}
delayI :: Dynamics m a
-> Dynamics m Double -> Dynamics m a -> Simulation m (Dynamics m a)
delayI (Dynamics Point m -> m a
x) (Dynamics Point m -> m Double
d) (Dynamics Point m -> m a
i) = Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
MonadMemo m =>
Dynamics m e -> Simulation m (Dynamics m e)
M.memo0Dynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r 
  where
    r :: Point m -> m a
r Point m
p = do 
      let t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
          sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
      Double
a <- Point m -> m Double
d Point m
p
      let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a
          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 m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
          y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
i (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
$ Point m
p { pointTime :: Double
pointTime = Double
t',
                                  pointIteration :: Int
pointIteration = Int
n',
                                  pointPhase :: Int
pointPhase = -Int
1 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      m a
y

-- | Return the delayed value by the specified positive number of
-- integration time steps used for calculating the lag time.
delayByDT :: Monad m
             => Dynamics m a
             -- ^ the value to delay
             -> Dynamics m Int
             -- ^ the delay as a multiplication of the corresponding number
             -- and the integration time step
             -> Dynamics m a
             -- ^ the delayed value
{-# INLINABLE delayByDT #-}
delayByDT :: Dynamics m a -> Dynamics m Int -> Dynamics m a
delayByDT (Dynamics Point m -> m a
x) (Dynamics Point m -> m Int
d) = Dynamics m a -> Dynamics m a
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m a -> Dynamics m a) -> Dynamics m a -> Dynamics m a
forall a b. (a -> b) -> a -> b
$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r 
  where
    r :: Point m -> m a
r Point m
p = do 
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
      Int
a <- Point m -> m Int
d Point m
p
      let p' :: Point m
p' = Point m -> Int -> Point m
forall (m :: * -> *). Point m -> Int -> Point m
delayPoint Point m
p Int
a
          n' :: Int
n' = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p'
          y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x Point m
p'
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delayByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delayByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      m a
y
      
-- | Return the delayed value by the specified initial value and
-- a positive number of integration time steps used for calculating
-- the lag time. It allows creating a loop back.
delayIByDT :: MonadSD m
              => Dynamics m a
              -- ^ the value to delay
              -> Dynamics m Int
              -- ^ the delay as a multiplication of the corresponding number
              -- and the integration time step
              -> Dynamics m a
              -- ^ the initial value
              -> Simulation m (Dynamics m a)
              -- ^ the delayed value
{-# INLINABLE delayIByDT #-}
delayIByDT :: Dynamics m a
-> Dynamics m Int -> Dynamics m a -> Simulation m (Dynamics m a)
delayIByDT (Dynamics Point m -> m a
x) (Dynamics Point m -> m Int
d) (Dynamics Point m -> m a
i) = Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
MonadMemo m =>
Dynamics m e -> Simulation m (Dynamics m e)
M.memoDynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r 
  where
    r :: Point m -> m a
r Point m
p = do 
      let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
          n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
      Int
a <- Point m -> m Int
d Point m
p
      let p' :: Point m
p' = Point m -> Int -> Point m
forall (m :: * -> *). Point m -> Int -> Point m
delayPoint Point m
p Int
a
          n' :: Int
n' = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p'
          y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
i (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x Point m
p'
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delayIByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delayIByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      m a
y

--
-- Financial Functions
--

-- | Return the Net Present Value (NPV) of the stream computed using the specified
-- discount rate, the initial value and some factor (usually 1).
--
-- It is defined in the following way:
--
-- @
-- npv stream rate init factor =
--   mdo let dt' = liftParameter dt
--       df <- integ (- df * rate) 1
--       accum <- integ (stream * df) init
--       return $ (accum + dt' * stream * df) * factor
-- @
npv :: (MonadSD m, MonadFix m)
       => Dynamics m Double                  -- ^ the stream
       -> Dynamics m Double                  -- ^ the discount rate
       -> Dynamics m Double                  -- ^ the initial value
       -> Dynamics m Double                  -- ^ factor
       -> Simulation m (Dynamics m Double)   -- ^ the Net Present Value (NPV)
{-# INLINABLE npv #-}
npv :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
npv Dynamics m Double
stream Dynamics m Double
rate Dynamics m Double
init Dynamics m Double
factor =
  mdo let dt' :: Dynamics m Double
dt' = Parameter m Double -> Dynamics m Double
forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
ParameterLift t m =>
Parameter m a -> t m a
liftParameter Parameter m Double
forall (m :: * -> *). Monad m => Parameter m Double
dt
      Dynamics m Double
df <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (- Dynamics m Double
df Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
rate) Dynamics m Double
1
      Dynamics m Double
accum <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double
init
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Dynamics m Double
accum Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
dt' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
factor

-- | Return the Net Present Value End of period (NPVE) of the stream computed
-- using the specified discount rate, the initial value and some factor.
--
-- It is defined in the following way:
--
-- @
-- npve stream rate init factor =
--   mdo let dt' = liftParameter dt
--       df <- integ (- df * rate \/ (1 + rate * dt')) (1 \/ (1 + rate * dt'))
--       accum <- integ (stream * df) init
--       return $ (accum + dt' * stream * df) * factor
-- @
npve :: (MonadSD m, MonadFix m)
        => Dynamics m Double                  -- ^ the stream
        -> Dynamics m Double                  -- ^ the discount rate
        -> Dynamics m Double                  -- ^ the initial value
        -> Dynamics m Double                  -- ^ factor
        -> Simulation m (Dynamics m Double)   -- ^ the Net Present Value End (NPVE)
{-# INLINABLE npve #-}
npve :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
npve Dynamics m Double
stream Dynamics m Double
rate Dynamics m Double
init Dynamics m Double
factor =
  mdo let dt' :: Dynamics m Double
dt' = Parameter m Double -> Dynamics m Double
forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
ParameterLift t m =>
Parameter m a -> t m a
liftParameter Parameter m Double
forall (m :: * -> *). Monad m => Parameter m Double
dt
      Dynamics m Double
df <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (- Dynamics m Double
df Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
rate Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ (Dynamics m Double
1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
rate Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
dt')) (Dynamics m Double
1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ (Dynamics m Double
1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
rate Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
dt'))
      Dynamics m Double
accum <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
(MonadSD m, MonadFix m) =>
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double
init
      Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
$ (Dynamics m Double
accum Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
dt' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
factor

-- | Computation that returns 0 until the step time and then returns the specified height.
step :: Monad m
        => Dynamics m Double
        -- ^ the height
        -> Dynamics m Double
        -- ^ the step time
        -> Dynamics m Double
{-# INLINABLE step #-}
step :: Dynamics m Double -> Dynamics m Double -> Dynamics m Double
step Dynamics m Double
h Dynamics m Double
st =
  Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$
  (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
  do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
         t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
     Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
     let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
     if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t'
       then Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
h
       else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0

-- | Computation that returns 1, starting at the time start, and lasting for the interval
-- width; 0 is returned at all other times.
pulse :: Monad m
         => Dynamics m Double
         -- ^ the time start
         -> Dynamics m Double
         -- ^ the interval width
         -> Dynamics m Double
{-# INLINABLE pulse #-}
pulse :: Dynamics m Double -> Dynamics m Double -> Dynamics m Double
pulse Dynamics m Double
st Dynamics m Double
w =
  Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$
  (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
  do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
         t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
     Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
     let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
     if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t'
       then do Double
w' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
w
               Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ if Double
t' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
st' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w' then Double
1 else Double
0
       else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0

-- | Computation that returns 1, starting at the time start, and lasting for the interval
-- width and then repeats this pattern with the specified period; 0 is returned at all
-- other times.
pulseP :: Monad m
          => Dynamics m Double
          -- ^ the time start
          -> Dynamics m Double
          -- ^ the interval width
          -> Dynamics m Double
          -- ^ the time period
          -> Dynamics m Double
{-# INLINABLE pulseP #-}
pulseP :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Dynamics m Double
pulseP Dynamics m Double
st Dynamics m Double
w Dynamics m Double
period =
  Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$
  (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
  do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
         t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
     Double
p'  <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
period
     Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
     let y' :: Double
y' = if (Double
p' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) Bool -> Bool -> Bool
&& (Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
st')
              then Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (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
- Double
st') Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
p') Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p'
              else Double
0
     let st' :: Double
st' = Double
st' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y'
     let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
     if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t'
       then do Double
w' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
w
               Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ if Double
t' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
st' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w' then Double
1 else Double
0
       else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0

-- | Computation that returns 0 until the specified time start and then
-- slopes upward until the end time and then holds constant.
ramp :: Monad m
        => Dynamics m Double
        -- ^ the slope parameter
        -> Dynamics m Double
        -- ^ the time start
        -> Dynamics m Double
        -- ^ the end time
        -> Dynamics m Double
{-# INLINABLE ramp #-}
ramp :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Dynamics m Double
ramp Dynamics m Double
slope Dynamics m Double
st Dynamics m Double
e =
  Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$
  (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
$ \Point m
p ->
  do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
         t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
     Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
     if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t
       then do Double
slope' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
slope
               Double
e' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
e
               if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
e'
                 then Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ Double
slope' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
st')
                 else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
$ Double
slope' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
e' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
st')
       else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0