module Simulation.Aivika.Dynamics.SystemDynamics
(
(.==.),
(./=.),
(.<.),
(.>=.),
(.>.),
(.<=.),
maxDynamics,
minDynamics,
ifDynamics,
Integ,
newInteg,
integInit,
integValue,
integDiff,
integ,
smoothI,
smooth,
smooth3I,
smooth3,
smoothNI,
smoothN,
delay1I,
delay1,
delay3I,
delay3,
delayNI,
delayN,
forecast,
trend,
Sum,
newSum,
sumInit,
sumValue,
sumDiff,
sumDynamics,
lookupD,
lookupStepwiseD,
lookupDynamics,
lookupStepwiseDynamics,
delayTrans,
delay,
delayI,
udelay,
udelayI,
npv,
npve) where
import Data.Array
import Data.Array.IO.Safe
import Data.IORef
import Control.Monad
import Control.Monad.Trans
import Simulation.Aivika.Dynamics.Internal.Simulation
import Simulation.Aivika.Dynamics.Internal.Dynamics
import Simulation.Aivika.Dynamics.Base
(.==.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool
(.==.) = liftM2 (==)
(./=.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool
(./=.) = liftM2 (/=)
(.<.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
(.<.) = liftM2 (<)
(.>=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
(.>=.) = liftM2 (>=)
(.>.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
(.>.) = liftM2 (>)
(.<=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
(.<=.) = liftM2 (<=)
maxDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
maxDynamics = liftM2 max
minDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
minDynamics = liftM2 min
ifDynamics :: Dynamics Bool -> Dynamics a -> Dynamics a -> Dynamics a
ifDynamics cond x y =
do a <- cond
if a then x else y
data Integ = Integ { integInit :: Dynamics Double,
integExternal :: IORef (Dynamics Double),
integInternal :: IORef (Dynamics Double) }
newInteg :: Dynamics Double -> Simulation Integ
newInteg i =
do r1 <- liftIO $ newIORef $ initDynamics i
r2 <- liftIO $ newIORef $ initDynamics i
let integ = Integ { integInit = i,
integExternal = r1,
integInternal = r2 }
z = Dynamics $ \p ->
do (Dynamics m) <- readIORef (integInternal integ)
m p
y <- umemo z
liftIO $ writeIORef (integExternal integ) y
return integ
integValue :: Integ -> Dynamics Double
integValue integ =
Dynamics $ \p ->
do (Dynamics m) <- readIORef (integExternal integ)
m p
integDiff :: Integ -> Dynamics Double -> Simulation ()
integDiff integ diff =
do let z = Dynamics $ \p ->
do y <- readIORef (integExternal integ)
let i = integInit integ
case spcMethod (pointSpecs p) of
Euler -> integEuler diff i y p
RungeKutta2 -> integRK2 diff i y p
RungeKutta4 -> integRK4 diff i y p
liftIO $ writeIORef (integInternal integ) z
integEuler :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Point -> IO Double
integEuler (Dynamics f) (Dynamics i) (Dynamics y) p =
case pointIteration p of
0 ->
i p
n -> do
let sc = pointSpecs p
ty = basicTime sc (n 1) 0
py = p { pointTime = ty, pointIteration = n 1, pointPhase = 0 }
a <- y py
b <- f py
let !v = a + spcDT (pointSpecs p) * b
return v
integRK2 :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Point -> IO Double
integRK2 (Dynamics f) (Dynamics i) (Dynamics y) p =
case pointPhase p of
0 -> case pointIteration p of
0 ->
i p
n -> do
let sc = pointSpecs p
ty = basicTime sc (n 1) 0
t1 = ty
t2 = basicTime sc (n 1) 1
py = p { pointTime = ty, pointIteration = n 1, pointPhase = 0 }
p1 = py
p2 = p { pointTime = t2, pointIteration = n 1, pointPhase = 1 }
vy <- y py
k1 <- f p1
k2 <- f p2
let !v = vy + spcDT sc / 2.0 * (k1 + k2)
return v
1 -> do
let sc = pointSpecs p
n = pointIteration p
ty = basicTime sc n 0
t1 = ty
py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
p1 = py
vy <- y py
k1 <- f p1
let !v = vy + spcDT sc * k1
return v
_ ->
error "Incorrect phase: integRK2"
integRK4 :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Point -> IO Double
integRK4 (Dynamics f) (Dynamics i) (Dynamics y) p =
case pointPhase p of
0 -> case pointIteration p of
0 ->
i p
n -> do
let sc = pointSpecs p
ty = basicTime sc (n 1) 0
t1 = ty
t2 = basicTime sc (n 1) 1
t3 = basicTime sc (n 1) 2
t4 = basicTime sc (n 1) 3
py = p { pointTime = ty, pointIteration = n 1, pointPhase = 0 }
p1 = py
p2 = p { pointTime = t2, pointIteration = n 1, pointPhase = 1 }
p3 = p { pointTime = t3, pointIteration = n 1, pointPhase = 2 }
p4 = p { pointTime = t4, pointIteration = n 1, pointPhase = 3 }
vy <- y py
k1 <- f p1
k2 <- f p2
k3 <- f p3
k4 <- f p4
let !v = vy + spcDT sc / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
return v
1 -> do
let sc = pointSpecs p
n = pointIteration p
ty = basicTime sc n 0
t1 = ty
py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
p1 = py
vy <- y py
k1 <- f p1
let !v = vy + spcDT sc / 2.0 * k1
return v
2 -> do
let sc = pointSpecs p
n = pointIteration p
ty = basicTime sc n 0
t2 = basicTime sc n 1
py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }
vy <- y py
k2 <- f p2
let !v = vy + spcDT sc / 2.0 * k2
return v
3 -> do
let sc = pointSpecs p
n = pointIteration p
ty = basicTime sc n 0
t3 = basicTime sc n 2
py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
p3 = p { pointTime = t3, pointIteration = n, pointPhase = 2 }
vy <- y py
k3 <- f p3
let !v = vy + spcDT sc * k3
return v
_ ->
error "Incorrect phase: integRK4"
integ :: Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
integ diff i =
mdo y <- umemo z
z <- Simulation $ \r ->
case spcMethod (runSpecs r) of
Euler -> return $ Dynamics $ integEuler diff i y
RungeKutta2 -> return $ Dynamics $ integRK2 diff i y
RungeKutta4 -> return $ Dynamics $ integRK4 diff i y
return y
smoothI :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI x t i =
mdo y <- integ ((x y) / t) i
return y
smooth :: Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smooth x t = smoothI x t x
smooth3I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
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
smooth3 :: Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smooth3 x t = smooth3I x t x
smoothNI :: Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothNI x t n i =
mdo s <- forM [1 .. n] $ \k ->
if k == 1
then integ ((x a ! 1) / t') i
else integ ((a ! (k 1) a ! k) / t') i
let a = listArray (1, n) s
t' = t / fromIntegral n
return $ a ! n
smoothN :: Dynamics Double
-> Dynamics Double
-> Int
-> Simulation (Dynamics Double)
smoothN x t n = smoothNI x t n x
delay1I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay1I x t i =
mdo y <- integ (x y / t) (i * t)
return $ y / t
delay1 :: Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay1 x t = delay1I x t x
delay3I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay3I x t i =
mdo y <- integ (s2 / t' y / t') (i * t')
s2 <- integ (s1 / t' s2 / t') (i * t')
s1 <- integ (x s1 / t') (i * t')
let t' = t / 3.0
return $ y / t'
delay3 :: Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay3 x t = delay3I x t x
delayNI :: Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
delayNI x t n i =
mdo s <- forM [1 .. n] $ \k ->
if k == 1
then integ (x (a ! 1) / t') (i * t')
else integ ((a ! (k 1)) / t' (a ! k) / t') (i * t')
let a = listArray (1, n) s
t' = t / fromIntegral n
return $ (a ! n) / t'
delayN :: Dynamics Double
-> Dynamics Double
-> Int
-> Simulation (Dynamics Double)
delayN x t n = delayNI x t n x
forecast :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
forecast x at hz =
do y <- smooth x at
return $ x * (1.0 + (x / y 1.0) / at * hz)
trend :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
trend x at i =
do y <- smoothI x at (x / (1.0 + i * at))
return $ (x / y 1.0) / at
data Sum a = Sum { sumInit :: Dynamics a,
sumExternal :: IORef (Dynamics a),
sumInternal :: IORef (Dynamics a) }
newSum :: (MArray IOUArray a IO, Num a) => Dynamics a -> Simulation (Sum a)
newSum i =
do r1 <- liftIO $ newIORef $ initDynamics i
r2 <- liftIO $ newIORef $ initDynamics i
let sum = Sum { sumInit = i,
sumExternal = r1,
sumInternal = r2 }
z = Dynamics $ \p ->
do (Dynamics m) <- readIORef (sumInternal sum)
m p
y <- umemo0 z
liftIO $ writeIORef (sumExternal sum) y
return sum
sumValue :: Sum a -> Dynamics a
sumValue sum =
Dynamics $ \p ->
do (Dynamics m) <- readIORef (sumExternal sum)
m p
sumDiff :: (MArray IOUArray a IO, Num a) => Sum a -> Dynamics a -> Simulation ()
sumDiff sum (Dynamics diff) =
do let z = Dynamics $ \p ->
case pointIteration p of
0 -> do
let Dynamics i = sumInit sum
i p
n -> do
Dynamics y <- readIORef (sumExternal sum)
let sc = pointSpecs p
ty = basicTime sc (n 1) 0
py = p { pointTime = ty,
pointIteration = n 1,
pointPhase = 0 }
a <- y py
b <- diff py
let !v = a + b
return v
liftIO $ writeIORef (sumInternal sum) z
sumDynamics :: (MArray IOUArray a IO, Num a)
=> Dynamics a
-> Dynamics a
-> Simulation (Dynamics a)
sumDynamics (Dynamics diff) (Dynamics i) =
mdo y <- umemo z
z <- Simulation $ \r ->
return $ Dynamics $ \p ->
case pointIteration p of
0 -> i p
n -> do
let Dynamics m = y
sc = pointSpecs p
ty = basicTime sc (n 1) 0
py = p { pointTime = ty,
pointIteration = n 1,
pointPhase = 0 }
a <- m py
b <- diff py
let !v = a + b
return v
return y
lookupD :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupD = lookupDynamics
lookupStepwiseD :: Dynamics Double
-> Array Int (Double, Double)
-> Dynamics Double
lookupStepwiseD = lookupStepwiseDynamics
lookupDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupDynamics (Dynamics m) tbl =
Dynamics (\p -> do a <- m p; return $ find first last a) where
(first, last) = bounds tbl
find left right x =
if left > right then
error "Incorrect index: table"
else
let index = (left + 1 + right) `div` 2
x1 = fst $ tbl ! index
in if x1 <= x then
let y | index < right = find index right x
| right == last = snd $ tbl ! right
| otherwise =
let x2 = fst $ tbl ! (index + 1)
y1 = snd $ tbl ! index
y2 = snd $ tbl ! (index + 1)
in y1 + (y2 y1) * (x x1) / (x2 x1)
in y
else
let y | left < index = find left (index 1) x
| left == first = snd $ tbl ! left
| otherwise = error "Incorrect index: table"
in y
lookupStepwiseDynamics :: Dynamics Double
-> Array Int (Double, Double)
-> Dynamics Double
lookupStepwiseDynamics (Dynamics m) tbl =
Dynamics (\p -> do a <- m p; return $ find first last a) where
(first, last) = bounds tbl
find left right x =
if left > right then
error "Incorrect index: table"
else
let index = (left + 1 + right) `div` 2
x1 = fst $ tbl ! index
in if x1 <= x then
let y | index < right = find index right x
| right == last = snd $ tbl ! right
| otherwise = snd $ tbl ! right
in y
else
let y | left < index = find left (index 1) x
| left == first = snd $ tbl ! left
| otherwise = error "Incorrect index: table"
in y
delayTrans :: Dynamics a
-> Dynamics Double
-> Dynamics a
-> (Dynamics a -> Simulation (Dynamics a))
-> Simulation (Dynamics a)
delayTrans (Dynamics x) (Dynamics d) (Dynamics i) tr = tr $ Dynamics r
where
r p = do
let t = pointTime p
sc = pointSpecs p
n = pointIteration p
a <- d p
let t' = t a
n' = fromIntegral $ floor $ (t' spcStartTime sc) / spcDT sc
y | n' < 0 = i $ p { pointTime = spcStartTime sc,
pointIteration = 0,
pointPhase = 0 }
| n' < n = x $ p { pointTime = t',
pointIteration = n',
pointPhase = 1 }
| n' > n = error $
"Cannot return the future data: delayTrans. " ++
"The lag time cannot be negative."
| otherwise = error $
"Cannot return the current data: delayTrans. " ++
"The lag time is too small."
y
delay :: Dynamics a
-> Dynamics Double
-> Simulation (Dynamics a)
delay x d = delayTrans x d x memo0
delayI :: Dynamics a
-> Dynamics Double
-> Dynamics a
-> Simulation (Dynamics a)
delayI x d i = delayTrans x d i memo0
udelay :: (MArray IOUArray a IO, Num a)
=> Dynamics a
-> Dynamics Double
-> Simulation (Dynamics a)
udelay x d = delayTrans x d x umemo0
udelayI :: (MArray IOUArray a IO, Num a)
=> Dynamics a
-> Dynamics Double
-> Dynamics a
-> Simulation (Dynamics a)
udelayI x d i = delayTrans x d i umemo0
npv :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
npv stream rate init factor =
mdo df <- integ ( df * rate) 1
accum <- integ (stream * df) init
return $ (accum + dt * stream * df) * factor
npve :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
npve stream rate init factor =
mdo df <- integ ( df * rate / (1 + rate * dt)) (1 / (1 + rate * dt))
accum <- integ (stream * df) init
return $ (accum + dt * stream * df) * factor