{-# LANGUAGE BangPatterns, RecursiveDo #-}

-- |
-- Module     : Simulation.Aivika.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.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 Data.Array.IO.Safe
import Data.IORef
import Control.Monad
import Control.Monad.Trans

import Simulation.Aivika.Internal.Specs
import Simulation.Aivika.Internal.Parameter
import Simulation.Aivika.Internal.Simulation
import Simulation.Aivika.Internal.Dynamics
import Simulation.Aivika.Dynamics.Extra
import Simulation.Aivika.Unboxed
import Simulation.Aivika.Table

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

--
-- Equality and Ordering
--

-- | Compare for equality.
(.==.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool
.==. :: forall a. Eq a => Dynamics a -> Dynamics a -> Dynamics Bool
(.==.) = (a -> a -> Bool) -> Dynamics a -> Dynamics a -> Dynamics 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.
(./=.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool
./=. :: forall a. Eq a => Dynamics a -> Dynamics a -> Dynamics Bool
(./=.) = (a -> a -> Bool) -> Dynamics a -> Dynamics a -> Dynamics 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.
(.<.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.<. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.<.) = (a -> a -> Bool) -> Dynamics a -> Dynamics a -> Dynamics 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.
(.>=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.>=. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.>=.) = (a -> a -> Bool) -> Dynamics a -> Dynamics a -> Dynamics 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.
(.>.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.>. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.>.) = (a -> a -> Bool) -> Dynamics a -> Dynamics a -> Dynamics 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.
(.<=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.<=. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.<=.) = (a -> a -> Bool) -> Dynamics a -> Dynamics a -> Dynamics 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 :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
maxDynamics :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics a
maxDynamics = (a -> a -> a) -> Dynamics a -> Dynamics a -> Dynamics 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 :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
minDynamics :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics a
minDynamics = (a -> a -> a) -> Dynamics a -> Dynamics a -> Dynamics 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 :: Dynamics Bool -> Dynamics a -> Dynamics a -> Dynamics a
ifDynamics :: forall a. Dynamics Bool -> Dynamics a -> Dynamics a -> Dynamics a
ifDynamics Dynamics Bool
cond Dynamics a
x Dynamics a
y =
  do Bool
a <- Dynamics Bool
cond
     if Bool
a then Dynamics a
x else Dynamics a
y

--
-- Ordinary Differential Equations
--

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

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

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

integRK4b :: Dynamics Double
             -> Dynamics Double
             -> Dynamics Double
             -> Point -> IO Double
integRK4b :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4b (Dynamics Point -> IO Double
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p =
  case Point -> Int
pointPhase Point
p of
    Int
0 -> case Point -> Int
pointIteration Point
p of
      Int
0 -> 
        Point -> IO Double
i Point
p
      Int
n -> do
        let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
            ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
            t3 :: Double
t3 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
2
            t4 :: Double
t4 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
3
            py :: Point
py = Point
p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
            p1 :: Point
p1 = Point
py
            p2 :: Point
p2 = Point
p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }
            p3 :: Point
p3 = Point
p { pointTime = t3, pointIteration = n - 1, pointPhase = 2 }
            p4 :: Point
p4 = Point
p { pointTime = t4, pointIteration = n - 1, pointPhase = 3 }
        Double
vy <- Point -> IO Double
y Point
py
        Double
k1 <- Point -> IO Double
f Point
p1
        Double
k2 <- Point -> IO Double
f Point
p2
        Double
k3 <- Point -> IO Double
f Point
p3
        Double
k4 <- Point -> IO Double
f Point
p4
        let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
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 -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point
py = Point
p { pointTime = ty, pointIteration = n, pointPhase = 0 }
          p1 :: Point
p1 = Point
py
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
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 -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
2 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
1
          py :: Point
py = Point
p { pointTime = ty, pointIteration = n, pointPhase = 0 }
          p1 :: Point
p1 = Point
py
          p2 :: Point
p2 = Point
p { pointTime = t2, pointIteration = n, pointPhase = 1 }
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      Double
k2 <- Point -> IO Double
f Point
p2
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
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 -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
3 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
1
          t3 :: Double
t3 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
2
          py :: Point
py = Point
p { pointTime = ty, pointIteration = n, pointPhase = 0 }
          p1 :: Point
p1 = Point
py
          p2 :: Point
p2 = Point
p { pointTime = t2, pointIteration = n, pointPhase = 1 }
          p3 :: Point
p3 = Point
p { pointTime = t3, pointIteration = n, pointPhase = 2 }
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      Double
k2 <- Point -> IO Double
f Point
p2
      Double
k3 <- Point -> IO Double
f Point
p3
      let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
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 -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      [Char] -> IO 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 :: Simulation [Double]
-- 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 :: Dynamics Double                  -- ^ the derivative
         -> Dynamics Double               -- ^ the initial value
         -> Simulation (Dynamics Double)  -- ^ the integral
integ :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ Dynamics Double
diff Dynamics Double
i =
  mdo Dynamics Double
y <- Dynamics Double -> Simulation (Dynamics Double)
forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memoDynamics Dynamics Double
z
      Dynamics Double
z <- (Run -> IO (Dynamics Double)) -> Simulation (Dynamics Double)
forall a. (Run -> IO a) -> Simulation a
Simulation ((Run -> IO (Dynamics Double)) -> Simulation (Dynamics Double))
-> (Run -> IO (Dynamics Double)) -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ \Run
r ->
        case Specs -> Method
spcMethod (Run -> Specs
runSpecs Run
r) of
          Method
Euler -> Dynamics Double -> IO (Dynamics Double)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> IO (Dynamics Double))
-> Dynamics Double -> IO (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEuler Dynamics Double
diff Dynamics Double
i Dynamics Double
y
          Method
RungeKutta2 -> Dynamics Double -> IO (Dynamics Double)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> IO (Dynamics Double))
-> Dynamics Double -> IO (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK2 Dynamics Double
diff Dynamics Double
i Dynamics Double
y
          Method
RungeKutta4 -> Dynamics Double -> IO (Dynamics Double)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> IO (Dynamics Double))
-> Dynamics Double -> IO (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4 Dynamics Double
diff Dynamics Double
i Dynamics Double
y
          Method
RungeKutta4b -> Dynamics Double -> IO (Dynamics Double)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> IO (Dynamics Double))
-> Dynamics Double -> IO (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4b Dynamics Double
diff Dynamics Double
i Dynamics Double
y
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

integEulerEither :: Dynamics (Either Double Double)
                    -> Dynamics Double 
                    -> Dynamics Double 
                    -> Point -> IO Double
integEulerEither :: Dynamics (Either Double Double)
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEulerEither (Dynamics Point -> IO (Either Double Double)
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p = 
  case Point -> Int
pointIteration Point
p of
    Int
0 -> 
      Point -> IO Double
i Point
p
    Int
n -> do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
          py :: Point
py = Point
p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
      Either Double Double
b <- Point -> IO (Either Double Double)
f Point
py
      case Either Double Double
b of
        Left Double
v ->
          Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
        Right Double
b -> do
          Double
a <- Point -> IO Double
y Point
py
          let !v :: Double
v = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT (Point -> Specs
pointSpecs Point
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b
          Double -> IO Double
forall a. a -> IO a
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 :: Dynamics (Either Double Double)
               -- ^ either set a new 'Left' integral value, or use a 'Right' derivative
               -> Dynamics Double
               -- ^ the initial value
               -> Simulation (Dynamics Double)
integEither :: Dynamics (Either Double Double)
-> Dynamics Double -> Simulation (Dynamics Double)
integEither Dynamics (Either Double Double)
diff Dynamics Double
i =
  mdo Dynamics Double
y <- Dynamics Double -> Simulation (Dynamics Double)
forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memoDynamics Dynamics Double
z
      Dynamics Double
z <- (Run -> IO (Dynamics Double)) -> Simulation (Dynamics Double)
forall a. (Run -> IO a) -> Simulation a
Simulation ((Run -> IO (Dynamics Double)) -> Simulation (Dynamics Double))
-> (Run -> IO (Dynamics Double)) -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ \Run
r ->
        Dynamics Double -> IO (Dynamics Double)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> IO (Dynamics Double))
-> Dynamics Double -> IO (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ Dynamics (Either Double Double)
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEulerEither Dynamics (Either Double Double)
diff Dynamics Double
i Dynamics Double
y
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics 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 :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth
smoothI :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
y) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t) Dynamics Double
i
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

-- | Return the first order exponential smooth.
--
-- This is a simplified version of the 'smoothI' function
-- without specifing the initial value.
smooth :: Dynamics Double                  -- ^ the value to smooth over time
          -> Dynamics Double               -- ^ time
          -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth
smooth :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
smooth Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI Dynamics Double
x Dynamics Double
t Dynamics 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 :: Dynamics Double                  -- ^ the value to smooth over time
            -> Dynamics Double               -- ^ time
            -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth
smooth3I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smooth3I Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y  <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
s2 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
y) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      Dynamics Double
s2 <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
s1 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
s2) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      Dynamics Double
s1 <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
s1) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      let t' :: Dynamics Double
t' = Dynamics Double
t Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
3.0
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

-- | Return the third order exponential smooth.
-- 
-- This is a simplified version of the 'smooth3I' function
-- without specifying the initial value.
smooth3 :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth
smooth3 :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
smooth3 Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smooth3I Dynamics Double
x Dynamics Double
t Dynamics 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 :: Dynamics Double                  -- ^ the value to smooth over time
            -> Dynamics Double               -- ^ time
            -> Int                           -- ^ the order
            -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth
smoothNI :: Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothNI Dynamics Double
x Dynamics Double
t Int
n Dynamics Double
i =
  mdo [Dynamics Double]
s <- [Int]
-> (Int -> Simulation (Dynamics Double))
-> Simulation [Dynamics 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 (Dynamics Double))
 -> Simulation [Dynamics Double])
-> (Int -> Simulation (Dynamics Double))
-> Simulation [Dynamics 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 Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics Double
forall i e. Ix i => Array i e -> i -> e
! Int
1) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
        else Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics 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 Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics Double
forall i e. Ix i => Array i e -> i -> e
! Int
k) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      let a :: Array Int (Dynamics Double)
a  = (Int, Int) -> [Dynamics Double] -> Array Int (Dynamics Double)
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics Double]
s 
          t' :: Dynamics Double
t' = Dynamics Double
t Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Int -> Dynamics Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics 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 :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Int                           -- ^ the order
           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth
smoothN :: Dynamics Double
-> Dynamics Double -> Int -> Simulation (Dynamics Double)
smoothN Dynamics Double
x Dynamics Double
t Int
n = Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothNI Dynamics Double
x Dynamics Double
t Int
n Dynamics 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 :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the first order exponential delay
delay1I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay1I Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
y Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t) (Dynamics Double
i Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
t)
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ Dynamics Double
y Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t

-- | Return the first order exponential delay.
--
-- This is a simplified version of the 'delay1I' function
-- without specifying the initial value.
delay1 :: Dynamics Double                  -- ^ the value to conserve
          -> Dynamics Double               -- ^ time
          -> Simulation (Dynamics Double)  -- ^ the first order exponential delay
delay1 :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
delay1 Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay1I Dynamics Double
x Dynamics Double
t Dynamics Double
x

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

-- | Return the third order exponential delay.
--
-- This is a simplified version of the 'delay3I' function
-- without specifying the initial value.
delay3 :: Dynamics Double                  -- ^ the value to conserve
          -> Dynamics Double               -- ^ time
          -> Simulation (Dynamics Double)  -- ^ the third order exponential delay
delay3 :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
delay3 Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay3I Dynamics Double
x Dynamics Double
t Dynamics Double
x

-- | Return the n'th order exponential delay.
delayNI :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Int                           -- ^ the order
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay
delayNI :: Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
delayNI Dynamics Double
x Dynamics Double
t Int
n Dynamics Double
i =
  mdo [Dynamics Double]
s <- [Int]
-> (Int -> Simulation (Dynamics Double))
-> Simulation [Dynamics 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 (Dynamics Double))
 -> Simulation [Dynamics Double])
-> (Int -> Simulation (Dynamics Double))
-> Simulation [Dynamics 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 Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- (Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics Double
forall i e. Ix i => Array i e -> i -> e
! Int
1) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
t')
        else Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics 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 Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t' Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- (Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics Double
forall i e. Ix i => Array i e -> i -> e
! Int
k) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
t')
      let a :: Array Int (Dynamics Double)
a  = (Int, Int) -> [Dynamics Double] -> Array Int (Dynamics Double)
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics Double]
s
          t' :: Dynamics Double
t' = Dynamics Double
t Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Int -> Dynamics Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Array Int (Dynamics Double)
a Array Int (Dynamics Double) -> Int -> Dynamics Double
forall i e. Ix i => Array i e -> i -> e
! Int
n) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
t'

-- | Return the n'th order exponential delay.
--
-- This is a simplified version of the 'delayNI' function
-- without specifying the initial value.
delayN :: Dynamics Double                  -- ^ the value to conserve
          -> Dynamics Double               -- ^ time
          -> Int                           -- ^ the order
          -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay
delayN :: Dynamics Double
-> Dynamics Double -> Int -> Simulation (Dynamics Double)
delayN Dynamics Double
x Dynamics Double
t Int
n = Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
delayNI Dynamics Double
x Dynamics Double
t Int
n Dynamics 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 :: Dynamics Double                  -- ^ the value to forecast
            -> Dynamics Double               -- ^ the average time
            -> Dynamics Double               -- ^ the time horizon
            -> Simulation (Dynamics Double)  -- ^ the forecast
forecast :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
forecast Dynamics Double
x Dynamics Double
at Dynamics Double
hz =
  do Dynamics Double
y <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
smooth Dynamics Double
x Dynamics Double
at
     Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* (Dynamics Double
1.0 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
+ (Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
y Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
1.0) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
at Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics 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 :: Dynamics Double                  -- ^ the value for which the trend is calculated
         -> Dynamics Double               -- ^ the average time
         -> Dynamics Double               -- ^ the initial value
         -> Simulation (Dynamics Double)  -- ^ the fractional change rate
trend :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
trend Dynamics Double
x Dynamics Double
at Dynamics Double
i =
  do Dynamics Double
y <- Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI Dynamics Double
x Dynamics Double
at (Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ (Dynamics Double
1.0 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
+ Dynamics Double
i Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
at))
     Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Dynamics Double
x Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics Double
y Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
- Dynamics Double
1.0) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ Dynamics 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 :: (Num a, Unboxed a)
           => Dynamics a               -- ^ the difference
           -> Dynamics a               -- ^ the initial value
           -> Simulation (Dynamics a)  -- ^ the sum
diffsum :: forall a.
(Num a, Unboxed a) =>
Dynamics a -> Dynamics a -> Simulation (Dynamics a)
diffsum (Dynamics Point -> IO a
diff) (Dynamics Point -> IO a
i) =
  mdo Dynamics a
y <-
        Dynamics a -> Simulation (Dynamics a)
forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memo0Dynamics (Dynamics a -> Simulation (Dynamics a))
-> Dynamics a -> Simulation (Dynamics a)
forall a b. (a -> b) -> a -> b
$
        (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO a) -> Dynamics a) -> (Point -> IO a) -> Dynamics a
forall a b. (a -> b) -> a -> b
$ \Point
p ->
        case Point -> Int
pointIteration Point
p of
          Int
0 -> Point -> IO a
i Point
p
          Int
n -> do 
            let Dynamics Point -> IO a
m = Dynamics a
y
                sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
                ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point
py = Point
p { pointTime = ty, 
                         pointIteration = n - 1, 
                         pointPhase = 0 }
            a
a <- Point -> IO a
m Point
py
            a
b <- Point -> IO a
diff Point
py
            let !v :: a
v = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
b
            a -> IO a
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return a
v
      Dynamics a -> Simulation (Dynamics a)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics a
y

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

--
-- Table Functions
--

-- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.
lookupDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupDynamics (Dynamics Point -> IO Double
m) Array Int (Double, Double)
tbl =
  (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do Double
a <- Point -> IO Double
m Point
p
     Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO 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 :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupStepwiseDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupStepwiseDynamics (Dynamics Point -> IO Double
m) Array Int (Double, Double)
tbl =
  (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do Double
a <- Point -> IO Double
m Point
p
     Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO 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 :: Dynamics a          -- ^ the value to delay
         -> Dynamics Double  -- ^ the lag time
         -> Dynamics a       -- ^ the delayed value
delay :: forall a. Dynamics a -> Dynamics Double -> Dynamics a
delay (Dynamics Point -> IO a
x) (Dynamics Point -> IO Double
d) = Dynamics a -> Dynamics a
forall a. Dynamics a -> Dynamics a
discreteDynamics (Dynamics a -> Dynamics a) -> Dynamics a -> Dynamics a
forall a b. (a -> b) -> a -> b
$ (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let t :: Double
t  = Point -> Double
pointTime Point
p
          sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Double
a <- Point -> IO Double
d Point
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 b. Integral b => Double -> b
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 -> Double
spcStartTime Specs
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
          y :: IO a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
x (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime = spcStartTime sc,
                                  pointIteration = 0, 
                                  pointPhase = 0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime = t',
                                  pointIteration = n',
                                  pointPhase = -1 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO 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] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO 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."
      IO 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 :: Dynamics a          -- ^ the value to delay
          -> Dynamics Double  -- ^ the lag time
          -> Dynamics a       -- ^ the initial value
          -> Simulation (Dynamics a)    -- ^ the delayed value
delayI :: forall a.
Dynamics a
-> Dynamics Double -> Dynamics a -> Simulation (Dynamics a)
delayI (Dynamics Point -> IO a
x) (Dynamics Point -> IO Double
d) (Dynamics Point -> IO a
i) = Dynamics a -> Simulation (Dynamics a)
forall e. Dynamics e -> Simulation (Dynamics e)
M.memo0Dynamics (Dynamics a -> Simulation (Dynamics a))
-> Dynamics a -> Simulation (Dynamics a)
forall a b. (a -> b) -> a -> b
$ (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let t :: Double
t  = Point -> Double
pointTime Point
p
          sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Double
a <- Point -> IO Double
d Point
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 b. Integral b => Double -> b
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 -> Double
spcStartTime Specs
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
          y :: IO a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
i (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime = spcStartTime sc,
                                  pointIteration = 0, 
                                  pointPhase = 0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime = t',
                                  pointIteration = n',
                                  pointPhase = -1 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delayI. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = [Char] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO a
forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delayI. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      IO a
y

-- | Return the delayed value by the specified positive number of
-- integration time steps used for calculating the lag time.
delayByDT :: Dynamics a
             -- ^ the value to delay
             -> Dynamics Int
             -- ^ the delay as a multiplication of the corresponding number
             -- and the integration time step
             -> Dynamics a
             -- ^ the delayed value
delayByDT :: forall a. Dynamics a -> Dynamics Int -> Dynamics a
delayByDT (Dynamics Point -> IO a
x) (Dynamics Point -> IO Int
d) = Dynamics a -> Dynamics a
forall a. Dynamics a -> Dynamics a
discreteDynamics (Dynamics a -> Dynamics a) -> Dynamics a -> Dynamics a
forall a b. (a -> b) -> a -> b
$ (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Int
a <- Point -> IO Int
d Point
p
      let p' :: Point
p' = Point -> Int -> Point
delayPoint Point
p Int
a
          n' :: Int
n' = Point -> Int
pointIteration Point
p'
          y :: IO a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
x (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime = spcStartTime sc,
                                  pointIteration = 0, 
                                  pointPhase = 0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x Point
p'
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO 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] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO 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."
      IO 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 :: Dynamics a
              -- ^ the value to delay
              -> Dynamics Int
              -- ^ the delay as a multiplication of the corresponding number
              -- and the integration time step
              -> Dynamics a
              -- ^ the initial value
              -> Simulation (Dynamics a)
              -- ^ the delayed value
delayIByDT :: forall a.
Dynamics a -> Dynamics Int -> Dynamics a -> Simulation (Dynamics a)
delayIByDT (Dynamics Point -> IO a
x) (Dynamics Point -> IO Int
d) (Dynamics Point -> IO a
i) = Dynamics a -> Simulation (Dynamics a)
forall e. Dynamics e -> Simulation (Dynamics e)
M.memoDynamics (Dynamics a -> Simulation (Dynamics a))
-> Dynamics a -> Simulation (Dynamics a)
forall a b. (a -> b) -> a -> b
$ (Point -> IO a) -> Dynamics a
forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Int
a <- Point -> IO Int
d Point
p
      let p' :: Point
p' = Point -> Int -> Point
delayPoint Point
p Int
a
          n' :: Int
n' = Point -> Int
pointIteration Point
p'
          y :: IO a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
i (Point -> IO a) -> Point -> IO a
forall a b. (a -> b) -> a -> b
$ Point
p { pointTime = spcStartTime sc,
                                  pointIteration = 0, 
                                  pointPhase = 0 }
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x Point
p'
            | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO 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] -> IO a
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO a) -> [Char] -> IO 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."
      IO 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 :: Dynamics Double                  -- ^ the stream
       -> Dynamics Double               -- ^ the discount rate
       -> Dynamics Double               -- ^ the initial value
       -> Dynamics Double               -- ^ factor
       -> Simulation (Dynamics Double)  -- ^ the Net Present Value (NPV)
npv :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
npv Dynamics Double
stream Dynamics Double
rate Dynamics Double
init Dynamics Double
factor =
  mdo let dt' :: Dynamics Double
dt' = Parameter Double -> Dynamics Double
forall a. Parameter a -> Dynamics a
forall (m :: * -> *) a. ParameterLift m => Parameter a -> m a
liftParameter Parameter Double
dt
      Dynamics Double
df <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (- Dynamics Double
df Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
rate) Dynamics Double
1
      Dynamics Double
accum <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
stream Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
df) Dynamics Double
init
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Dynamics Double
accum Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
+ Dynamics Double
dt' Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
stream Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
df) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics 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 :: Dynamics Double                  -- ^ the stream
        -> Dynamics Double               -- ^ the discount rate
        -> Dynamics Double               -- ^ the initial value
        -> Dynamics Double               -- ^ factor
        -> Simulation (Dynamics Double)  -- ^ the Net Present Value End (NPVE)
npve :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
npve Dynamics Double
stream Dynamics Double
rate Dynamics Double
init Dynamics Double
factor =
  mdo let dt' :: Dynamics Double
dt' = Parameter Double -> Dynamics Double
forall a. Parameter a -> Dynamics a
forall (m :: * -> *) a. ParameterLift m => Parameter a -> m a
liftParameter Parameter Double
dt
      Dynamics Double
df <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (- Dynamics Double
df Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
rate Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ (Dynamics Double
1 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
+ Dynamics Double
rate Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
dt')) (Dynamics Double
1 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Fractional a => a -> a -> a
/ (Dynamics Double
1 Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
+ Dynamics Double
rate Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
dt'))
      Dynamics Double
accum <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
stream Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
df) Dynamics Double
init
      Dynamics Double -> Simulation (Dynamics Double)
forall a. a -> Simulation a
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics Double -> Simulation (Dynamics Double))
-> Dynamics Double -> Simulation (Dynamics Double)
forall a b. (a -> b) -> a -> b
$ (Dynamics Double
accum Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
+ Dynamics Double
dt' Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
stream Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
df) Dynamics Double -> Dynamics Double -> Dynamics Double
forall a. Num a => a -> a -> a
* Dynamics Double
factor

-- | Computation that returns 0 until the step time and then returns the specified height.
step :: Dynamics Double
        -- ^ the height
        -> Dynamics Double
        -- ^ the step time
        -> Dynamics Double
step :: Dynamics Double -> Dynamics Double -> Dynamics Double
step Dynamics Double
h Dynamics Double
st =
  Dynamics Double -> Dynamics Double
forall a. Dynamics a -> Dynamics a
discreteDynamics (Dynamics Double -> Dynamics Double)
-> Dynamics Double -> Dynamics Double
forall a b. (a -> b) -> a -> b
$
  (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
st' <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
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 -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
h
       else Double -> IO Double
forall a. a -> IO a
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 :: Dynamics Double
         -- ^ the time start
         -> Dynamics Double
         -- ^ the interval width
         -> Dynamics Double
pulse :: Dynamics Double -> Dynamics Double -> Dynamics Double
pulse Dynamics Double
st Dynamics Double
w =
  Dynamics Double -> Dynamics Double
forall a. Dynamics a -> Dynamics a
discreteDynamics (Dynamics Double -> Dynamics Double)
-> Dynamics Double -> Dynamics Double
forall a b. (a -> b) -> a -> b
$
  (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
st' <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
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 -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
w
               Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO 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 -> IO Double
forall a. a -> IO a
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 :: Dynamics Double
          -- ^ the time start
          -> Dynamics Double
          -- ^ the interval width
          -> Dynamics Double
          -- ^ the time period
          -> Dynamics Double
pulseP :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Dynamics Double
pulseP Dynamics Double
st Dynamics Double
w Dynamics Double
period =
  Dynamics Double -> Dynamics Double
forall a. Dynamics a -> Dynamics a
discreteDynamics (Dynamics Double -> Dynamics Double)
-> Dynamics Double -> Dynamics Double
forall a b. (a -> b) -> a -> b
$
  (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
p'  <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
period
     Double
st' <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics 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 b. Integral b => Double -> b
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 -> Double
spcDT Specs
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 -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
w
               Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO 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 -> IO Double
forall a. a -> IO a
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 :: Dynamics Double
        -- ^ the slope parameter
        -> Dynamics Double
        -- ^ the time start
        -> Dynamics Double
        -- ^ the end time
        -> Dynamics Double
ramp :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Dynamics Double
ramp Dynamics Double
slope Dynamics Double
st Dynamics Double
e =
  Dynamics Double -> Dynamics Double
forall a. Dynamics a -> Dynamics a
discreteDynamics (Dynamics Double -> Dynamics Double)
-> Dynamics Double -> Dynamics Double
forall a b. (a -> b) -> a -> b
$
  (Point -> IO Double) -> Dynamics Double
forall a. (Point -> IO a) -> Dynamics a
Dynamics ((Point -> IO Double) -> Dynamics Double)
-> (Point -> IO Double) -> Dynamics Double
forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
st' <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t
       then do Double
slope' <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
slope
               Double
e' <- Point -> Dynamics Double -> IO Double
forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
e
               if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
e'
                 then Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO 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 -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO 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 -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0