diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,8 @@ # Changelog for dde -## Unreleased changes +## 0.1.0 *March 21st 2018* + * Now, all dynamical variables are recorded + * Improved speed + +## 0.0.1 *February 23rd 2018* + * Initial release diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,12 +2,11 @@ ## Features -* Autonomous DDEs with multiple dynamical variables and a single delay (pull requests are welcome) +* Autonomous DDEs with multiple dynamical variables and a single delay time (pull requests are welcome) * Driven systems (i.e. with external input) * Non-autonomous DDEs (using driven systems with time as external input) * Second and fourth order integration methods * Example models: * [Mackey-Glass](https://github.com/masterdezign/dde/blob/master/examples/MackeyGlass/Main.hs) with no external input - * [Driven system](https://github.com/masterdezign/dde/blob/d7f636372b537b948d00097ecd09e689854b9392/dde/Numeric/DDE/Model.hs#L59) + * [Driven system](https://github.com/masterdezign/dde/blob/d22c6ff82fd56c29289366a057f3d733a23844d0/dde/Numeric/DDE/Model.hs#L60) * Pure Haskell - diff --git a/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -0,0 +1,27 @@ +import Criterion.Main +import Linear.V1 ( V1 (..) ) +import qualified Data.Vector.Storable as V + +import qualified Impl1 +import qualified Impl2 + +runTest1 :: Double -> Double +runTest1 maxTime = + let hStep = 0.1 + total = round(maxTime / hStep) + s = Impl1.mgModel hStep total + in s + +runTest2 :: Double -> Double +runTest2 maxTime = + let hStep = 0.1 + total = round(maxTime / hStep) + (V1 s, _) = Impl2.mgModel hStep total + in s + +main = defaultMain [ + bench "old version" $ whnf runTest1 maxTime + , bench "target to optimize" $ whnf runTest2 maxTime + ] + where + maxTime = 1000 -- 1000 time units diff --git a/bench/Impl1.lhs b/bench/Impl1.lhs new file mode 100644 --- /dev/null +++ b/bench/Impl1.lhs @@ -0,0 +1,190 @@ +% How to compile a pdf: +% $ cabal exec lhs2TeX -- -o MackeyGlass.tex MackeyGlass.lhs && pdflatex \ +% MackeyGlass.tex + +\documentclass{article} +%include polycode.fmt + +\title{Mackey-Glass DDE} +\author{Bogdan Penkovsky} + +\begin{document} + +\maketitle + +This document describes a numeric model of the Mackey-Glass DDE +optimized for speed. The fourth-order Runge-Kutta integration method is employed. + +> {-# LANGUAGE BangPatterns #-} +> module Impl1 ( mgModel ) where + +> import qualified Data.Vector.Storable as V +> import qualified Data.Vector.Storable.Mutable as VM +> import System.IO.Unsafe ( unsafePerformIO ) +> import Control.Monad +> import Control.Lens hiding ((<.>)) + + +Sample can be a vector of any length (x, y, z, ...). +We import V1 for single-component vector (scalar) Samples. + +> import Linear.V1 as V1 + +> type Sample = V1.V1 +> type Delay = V.Vector +> type Input = V.Vector + +Autonomous system integrator + +iterate1 is the function describing DDE system; +len1 is the number of delay elements in a delay; +krecord is the number of last samples to record; +total is the total number of iterations; + +> integrator' +> :: (VM.Storable a, Floating a) => +> (Sample a -> (a, a) -> Sample a) +> -> Parameters a +> -> Int +> -> Int +> -> Int +> -> (Sample a, Delay a, Input a) +> -> (Sample a, V.Vector a) +> integrator' iterate1 _ len1 krecord total (!xy0, !hist0, _) = a +> where +> a = unsafePerformIO $ do +> ! v <- VM.new (len1 + total) -- Delay history +> -- Copy the initial history values +> copyHist v hist0 +> +> -- Calculate the rest of the vector +> xy' <- go v len1 xy0 +> +> trace <- V.unsafeFreeze v +> return (xy', V.slice (len1 + total - krecord) krecord trace) +> +> -- Copy initial conditions +> copyHist !v !hist = +> mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..V.length hist - 1] +> +> go !v !i !xy +> | i == len1 + total = +> return xy +> | otherwise = do +> x_tau1 <- VM.unsafeRead v (i - len1) -- Delayed element +> x_tau1' <- VM.unsafeRead v (i - len1 + 1) -- The next one +> let !xy' = iterate1 xy (x_tau1, x_tau1') +> !x' = xy' ^._x +> VM.unsafeWrite v i x' +> go v (i + 1) xy' + +%if False + +> {-# SPECIALISE integrator' :: +> (Sample Float -> (Float, Float) -> Sample Float) +> -> Parameters Float +> -> Int +> -> Int +> -> Int +> -> (Sample Float, Delay Float, Input Float) +> -> (Sample Float, V.Vector Float) #-} +> {-# SPECIALISE integrator' :: +> (Sample Double -> (Double, Double) -> Sample Double) +> -> Parameters Double +> -> Int +> -> Int +> -> Int +> -> (Sample Double, Delay Double, Input Double) +> -> (Sample Double, V.Vector Double) #-} + +%endif + +> data Parameters a = Parameters { +> pBeta :: a +> , pGamma :: a +> } deriving Show + +> param1 :: Parameters Double +> param1 = Parameters { +> pBeta = 0.2 +> , pGamma = 0.1 +> } + +The Mackey-Glass model. + +> mackeyGlass +> :: Floating a => Parameters a -> ((Sample a, a) -> Sample a) +> mackeyGlass p (V1.V1 !x, !x_tau1) = V1.V1 x' +> where +> ! x' = beta * x_tau1 / (1 + x_tau1^10) - gamma * x +> beta = pBeta p +> gamma = pGamma p + +%if False + +> {-# SPECIALISE mackeyGlass :: +> Parameters Float -> (Sample Float, Float) -> Sample Float #-} +> {-# SPECIALISE mackeyGlass :: +> Parameters Double -> (Sample Double, Double) -> Sample Double #-} + +%endif + +Fourth-order Runge-Kutta for a 1D system with a single delay $\tau_1$. + +> rk4 :: Double +> -> ((Sample Double, Double) -> Sample Double) +> -> Sample Double -> (Double, Double) -> Sample Double +> rk4 hStep sys !xy (!x_tau1, !x_tau1') = xy_next +> where +> xy_next = xy + over6 * (a + x2 * b + x2 * c + d) +> over6 = V1.V1 (recip 6) +> over2 = V1.V1 (recip 2) +> x2 = V1.V1 2 +> h = V1.V1 hStep +> ! a = h * sys (xy, x_tau1) +> ! b = h * sys (xy + over2 * a, x_tau1_b) +> ! c = h * sys (xy + over2 * b, x_tau1_c) +> ! d = h * sys (xy + c, x_tau1') +> ! x_tau1_b = (x_tau1 + x_tau1') / 2 +> ! x_tau1_c = x_tau1_b + +Returns the last delay + +> fastIntegrRk4 :: Double -> Int -> Int -> Int -> (V.Vector Double, Double) +> -> Parameters Double -> (V.Vector Double, Double) +> fastIntegrRk4 hStep len1 len2 totalIters (hist0, x0) p = (data1, x1) +> where sample0 = V1.V1 x0 +> -- Iterator implements Runge-Kutta schema +> iterator = rk4 hStep (mackeyGlass p) +> -- Record only the last long delay +> (V1.V1 x1, data1) = integrator' iterator p len1 len1 totalIters (sample0, hist0, V.fromList []) + +Records the whole time trace $x(t)$ + +> fastIntegrRk4' :: Double -> Int -> Int -> (V.Vector Double, Double) +> -> Parameters Double -> (V.Vector Double, Double) +> fastIntegrRk4' hStep len1 totalIters (hist0, x0) p = (data1, x1) +> where sample0 = V1.V1 x0 +> -- Iterator implements Runge-Kutta schema +> iterator = rk4 hStep (mackeyGlass p) +> -- Record all the time trace +> (V1.V1 x1, data1) = integrator' iterator p len1 totalIters totalIters (sample0, hist0, V.fromList []) + +Constant initial conditions + +> initCondConst :: Int -> Double -> [Double] +> initCondConst = replicate + +Define a Mackey-Glass simulation, the final state is the result + +> mgModel :: Double -> Int -> Double +> mgModel hStep total = s +> where +> len1 = 17 * round (recip hStep) -- tauD = 17, delay time +> +> icond0 = V.fromList $ initCondConst len1 0.2 +> +> -- Integrate an autonomous system (no external input) +> (_, s) = fastIntegrRk4' hStep len1 total (icond0, 0.2) param1 + +\end{document} diff --git a/bench/Impl2.hs b/bench/Impl2.hs new file mode 100644 --- /dev/null +++ b/bench/Impl2.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE FlexibleContexts #-} +module Impl2 ( mgModel ) where + +import Linear.V1 +import Numeric.DDE +import Numeric.DDE.Model +import qualified Data.Vector.Storable as V + +parMG0 :: MackeyGlass +parMG0 = MackeyGlass { _beta = 0.2 + , _gamma = 0.1 + } + +-- | The Mackey-Glass model (with no external input). +-- Typical hStep = 0.1 +mgModel :: Double -> Int -> (V1 Double, V.Vector (V1 Double)) +mgModel hStep totalIters = r + where + -- Initial state x(t0) = 0.2 + state0 = V1 0.2 + len1 = 17 * round (recip hStep) -- tauD = 17, delay time + -- Initial conditions + hist0 = V.replicate len1 state0 + inp = Input (V.replicate (totalIters + 1) 0.0) + rhs' = mackeyGlassRhs parMG0 + -- Stepper implements Runge-Kutta schema + stepper = let (Stepper _rk4) = rk4 + in _rk4 hStep rhs' + -- Provide the last state and the time trace + r = integ' stepper len1 totalIters totalIters (state0, hist0, inp) diff --git a/dde.cabal b/dde.cabal --- a/dde.cabal +++ b/dde.cabal @@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: eccc053c08411c8b92867a2d841ded9b1ba3ad1c1d3ac27e70e5acc6243c84df +-- hash: 005f02aea25cd38967b9112ce3aeb1c77f570482c9c94bdcebc344346e26de2d name: dde -version: 0.0.1 +version: 0.1.0 synopsis: Delay differential equations description: Please see the README on Github at category: Math @@ -32,6 +32,9 @@ dde build-depends: base >=4.7 && <5 + , free-vector-spaces + , lens + , linear , vector exposed-modules: Numeric.DDE @@ -49,6 +52,9 @@ build-depends: base >=4.7 && <5 , dde + , free-vector-spaces + , lens + , linear , vector other-modules: Paths_dde @@ -63,7 +69,32 @@ build-depends: base >=4.7 && <5 , dde + , free-vector-spaces + , lens + , linear , vector other-modules: + Paths_dde + default-language: Haskell2010 + +benchmark dde-bench + type: exitcode-stdio-1.0 + main-is: Bench.hs + hs-source-dirs: + bench + dde + build-depends: + base >=4.7 && <5 + , criterion + , free-vector-spaces + , lens + , linear + , vector + other-modules: + Impl1 + Impl2 + Numeric.DDE + Numeric.DDE.Model + Numeric.DDE.Types Paths_dde default-language: Haskell2010 diff --git a/dde/Numeric/DDE.hs b/dde/Numeric/DDE.hs --- a/dde/Numeric/DDE.hs +++ b/dde/Numeric/DDE.hs @@ -7,33 +7,32 @@ Below is a complete example simulating the Ikeda DDE defined as: @tau * x(t)/dt = -x + beta * sin[x(t - tau_D)]@. + > import Linear ( V1 (..) ) > import qualified Data.Vector.Storable as V > import qualified Numeric.DDE as DDE > - > ikedaRhs beta ((DDE.State xs), (DDE.Hist hs), _) = DDE.State $ V.fromList [x'] + > ikedaRhs beta = DDE.RHS derivative > where - > -- Ikeda DDE definition - > x' = (-x + beta * (sin x_tauD)) / tau - > - > -- Constants - > tau = 0.01 - > - > -- Dynamical variable x(t) - > x = V.head xs + > derivative ((V1 x), (DDE.Hist (V1 x_tauD)), _) = V1 x' + > where + > -- Ikeda DDE definition + > x' = (-x + beta * (sin x_tauD)) / tau > - > -- Delay term x(t - tau_D) - > x_tauD = V.head hs + > -- Constants + > tau = 0.01 > - > model beta hStep len1 totalIter = DDE.integ DDE.rk4 state0 hist0 len1 hStep (ikedaRhs beta) inp + > model beta hStep len1 totalIter = (state1, V.map (\(V1 x) -> x) trace) > where > -- Initial conditions: > -- dynamical state and delay history. - > state0 = DDE.State $ V.fromList [pi/2] - > hist0 = V.replicate len1 (pi/2) + > state0 = V1 (pi/2) + > hist0 = V.replicate len1 state0 > > -- Input is ignored in ikedaRhs > inp = DDE.Input $ V.replicate (totalIter + 1) 0 > + > (state1, trace) = DDE.integ DDE.rk4 state0 hist0 len1 hStep (ikedaRhs beta) inp + > > -- Control parameter > beta = 2.6 > @@ -51,78 +50,73 @@ -} {-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} + module Numeric.DDE ( -- * Integrators integ , integ' + , integRk4 + , integHeun2 , integRk4_2D , integHeun2_2D , Input (..) , InputSnapshot (..) - , State (..) , HistorySnapshot (..) -- * Steppers - , Stepper1 (..) + , RHS (..) + , Stepper (..) , rk4 , heun2 ) where +import Data.VectorSpace.Free +import Linear ( (^/), V1 (..), V2 (..) ) +import Foreign.Storable ( Storable (..) ) +import System.IO.Unsafe ( unsafePerformIO ) import qualified Data.Vector.Storable as V import qualified Data.Vector.Storable.Mutable as VM -import System.IO.Unsafe ( unsafePerformIO ) -import Numeric.DDE.Types +import Numeric.DDE.Types -infixl 6 .+. -(.+.) :: V.Vector Double -> V.Vector Double -> V.Vector Double -(.+.) = V.zipWith (+) -infixl 7 *. -(*.) :: Double -> V.Vector Double -> V.Vector Double -(*.) c = V.map (* c) - -- | Fourth order Runge-Kutta stepper -rk4 :: Stepper1 -rk4 = Stepper1 _rk4 +rk4 :: Stepper +rk4 = Stepper _rk4 where - _rk4 hStep rhs' (State !xy) !(!x_tau1, !x_tau1') !(!u1, !u1') = State xy_next + _rk4 dt (RHS rhs') !xy (Hist !xy_tau1, Hist !xy_tau1') (!u1, !u1') = xy_next where - xy_next = xy .+. over6 *. (a .+. 2 *. b .+. 2 *. c .+. d) - over6 = recip 6 - over2 = recip 2 - ! (State a') = rhs' (State xy, toHist1 x_tau1, inp1) - ! a = hStep *. a' - ! (State b') = rhs' (State $ xy .+. over2 *. a, toHist1 x_tau1_b, inp1_b) - ! b = hStep *. b' - ! (State c') = rhs' (State $ xy .+. over2 *. b, toHist1 x_tau1_c, inp1_c) - ! c = hStep *. c' - ! (State d') = rhs' (State $ xy .+. c, toHist1 x_tau1', inp1') - ! d = hStep *. d' - ! x_tau1_b = (x_tau1 + x_tau1') / 2 - ! x_tau1_c = x_tau1_b - ! inp1 = Inp u1 - ! inp1_b = Inp $ (u1 + u1') / 2 - ! inp1_c = inp1_b - ! inp1' = Inp u1' - -toHist1 :: Double -> HistorySnapshot -toHist1 = Hist. V.singleton + xy_next = xy ^+^ (a ^+^ 2 *^ b ^+^ 2 *^ c ^+^ d) ^/ 6 + a = dt *^ rhs' (xy, Hist xy_tau1, inp1) + b = dt *^ rhs' (xy ^+^ a ^/ 2, Hist xy_tau1_b, inp1_b) + c = dt *^ rhs' (xy ^+^ b ^/ 2, Hist xy_tau1_c, inp1_c) + d = dt *^ rhs' (xy ^+^ c, Hist xy_tau1', inp1') + xy_tau1_b = (xy_tau1 ^+^ xy_tau1') ^/ 2 + xy_tau1_c = xy_tau1_b + inp1 = Inp u1 + inp1_b = Inp $ (u1 + u1') / 2 + inp1_c = inp1_b + inp1' = Inp u1' -- | Second order Heun's stepper -heun2 :: Stepper1 -heun2 = Stepper1 _heun2 +heun2 :: Stepper +heun2 = Stepper _heun2 where - _heun2 hStep rhs' (State !xy) !(!x_tau1, !x_tau1') !(!u1, !u1') = State xy_next + _heun2 hStep (RHS rhs') !xy (!xy_tau1, !xy_tau1') (!u1, !u1') = xy_next where - ! (State f1) = rhs' (State xy, toHist1 x_tau1, Inp u1) - ! xy' = xy .+. hStep *. f1 - ! (State f2) = rhs' (State xy', toHist1 x_tau1', Inp u1') - ! xy_next = xy .+. (hStep / 2.0) *. (f1 .+. f2) + f1 = rhs' (xy, xy_tau1, Inp u1) + xy_ = xy ^+^ hStep *^ f1 + f2 = rhs' (xy_, xy_tau1', Inp u1') + xy_next = xy ^+^ (hStep *^ (f1 ^+^ f2)) ^/ 2.0 --- | Generic integrator for DDEs with a single delay +-- | Generic integrator for DDEs (single delay time). +-- Records all dynamical variables. +-- integ' - :: (State -> (Double, Double) -> (Double, Double) -> State) + :: Storable state + => (state -> (HistorySnapshot state, HistorySnapshot state) -> (Double, Double) -> state) -- ^ Iterator describing a DDE system -> Int -- ^ Delay length in samples @@ -130,12 +124,12 @@ -- ^ Number of last samples to record -> Int -- ^ Total number of iterations - -> (State, V.Vector Double, Input) + -> (state, V.Vector state, Input) -- ^ Initial state vector, initial history, and external forcing - -> (State, V.Vector Double) + -> (state, V.Vector state) -- ^ Final state and recorded state of the first variable. - -- The latter could be a Matrix if multiple variables are needed -integ' iter1 len1 krecord total (!xy0, !hist0, !(Input in1)) = a + -- The latter is a vector of vectors (matrix) when multiple variables are involved. +integ' iter1 len1 krecord total (!xy0, !hist0, Input !in1) = a where a = unsafePerformIO $ do ! v <- VM.new (len1 + total) -- Delay history @@ -149,62 +143,86 @@ return (xy', V.slice (len1 + total - krecord) krecord trace) -- Copy initial conditions - copyHist !v !hist = do - mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..(V.length hist) - 1] + copyHist !v !hist = + mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..V.length hist - 1] go !v !i !xy - | i == len1 + total = do + | i == len1 + total = return xy | otherwise = do - x_tau1 <- VM.unsafeRead v (i - len1) -- Delayed element - x_tau1' <- VM.unsafeRead v (i - len1 + 1) -- The next one - let u1 = in1 V.! (i - len1) -- Read two subsequent inputs + xy_tau1 <- VM.unsafeRead v (i - len1) -- Two subsequent delayed states + xy_tau1' <- VM.unsafeRead v (i - len1 + 1) + let u1 = in1 V.! (i - len1) -- Two subsequent inputs u1' = in1 V.! (i - len1 + 1) - !xy' = iter1 xy (x_tau1, x_tau1') (u1, u1') - !(State xy'1) = xy' - !x' = xy'1 V.! 0 -- Read x(t) variable - VM.unsafeWrite v i x' + !xy' = iter1 xy (Hist xy_tau1, Hist xy_tau1') (u1, u1') + VM.unsafeWrite v i xy' go v (i + 1) xy' -- | Generic integrator that records the whole time trace @x(t)@ +-- (single delay time). integ - :: Stepper1 - -> State -- ^ Initial state x(t), y(t),... - -> V.Vector Double -- ^ Initial history for the delayed variable + :: (Functor state, Storable (state Double), VectorSpace (state Double), Num (Scalar (state Double))) + => Stepper + -> state Double -- ^ Initial state vector (x(t), y(t),...) + -> V.Vector (state Double) -- ^ Initial history for delayed variables -> Int -- ^ Delay length in samples - -> Double -- ^ Integration step - -> RHS -- ^ Derivative (DDE right-hand side) + -> Scalar (state Double) -- ^ Integration step + -> RHS (state Double) -- ^ Derivative (DDE right-hand side) -> Input -- ^ External forcing - -> (State, V.Vector Double) -integ (Stepper1 stp) state0 hist0 len1 hStep rhs' inp@(Input in1) = r + -> (state Double, V.Vector (state Double)) +integ (Stepper stp) state0 hist0 len1 dt rhs' inp@(Input in1) = r where -- Two subsequent inputs are needed for `rk4` and `heun2`, -- therefore subtract one ! totalIters = V.length in1 - 1 - ! iterator = stp hStep rhs' + ! iterator = stp dt rhs' -- Record all the time trace ! r = integ' iterator len1 totalIters totalIters (state0, hist0, inp) +-- | RK4 integrator shortcut for 1D DDEs with zero +-- initial conditions +integRk4 :: Int -- ^ Delay length in samples + -> Double -- ^ Integration time step + -> RHS (V1 Double) -- ^ DDE model + -> Input -- ^ External forcing + -> (V1 Double, V.Vector (V1 Double)) +integRk4 len1 = integ rk4 state0 hist0 len1 + where + state0 = V1 0.0 + hist0 = V.replicate len1 state0 + +-- | Shortcut for Heun's 2nd order 1D DDEs with zero +-- initial conditions +integHeun2 :: Int -- ^ Delay length in samples + -> Double -- ^ Integration time step + -> RHS (V1 Double) -- ^ DDE model + -> Input -- ^ External forcing + -> (V1 Double, V.Vector (V1 Double)) +integHeun2 len1 = integ heun2 state0 hist0 len1 + where + state0 = V1 0.0 + hist0 = V.replicate len1 state0 + -- | RK4 integrator shortcut for 2D DDEs with zero -- initial conditions integRk4_2D :: Int -- ^ Delay length in samples -> Double -- ^ Integration time step - -> RHS -- ^ DDE model + -> RHS (V2 Double) -- ^ DDE model -> Input -- ^ External forcing - -> (State, V.Vector Double) + -> (V2 Double, V.Vector (V2 Double)) integRk4_2D len1 = integ rk4 state0 hist0 len1 where - ! state0 = State (V.replicate 2 0.0) - ! hist0 = V.replicate len1 0 + state0 = V2 0.0 0.0 + hist0 = V.replicate len1 state0 -- | Shortcut for Heun's 2nd order 2D DDEs with zero -- initial conditions integHeun2_2D :: Int -- ^ Delay length in samples - -> Double -- ^ Integration time step - -> RHS -- ^ DDE model + -> Double -- ^ Integration time step + -> RHS (V2 Double) -- ^ DDE model -> Input -- ^ External forcing - -> (State, V.Vector Double) + -> (V2 Double, V.Vector (V2 Double)) integHeun2_2D len1 = integ heun2 state0 hist0 len1 where - ! state0 = State (V.replicate 2 0.0) - ! hist0 = V.replicate len1 0 + state0 = V2 0.0 0.0 + hist0 = V.replicate len1 state0 diff --git a/dde/Numeric/DDE/Model.hs b/dde/Numeric/DDE/Model.hs --- a/dde/Numeric/DDE/Model.hs +++ b/dde/Numeric/DDE/Model.hs @@ -1,34 +1,41 @@ {- | = Example DDE models - These example models are defined by `rhs` function. - Adding a new example model requires adjusting `rhs` and `Par`. - Switching between models is performed by changing `Par`. -} module Numeric.DDE.Model ( - Par (..) - , BandpassFiltering (..) - , rhs + BandpassFiltering (..) + , RC (..) + , MackeyGlass (..) + , mackeyGlassRhs + , bandpassRhs ) where +import Control.Lens +import Linear ( R1 (..) + , R2 (..) ) +import qualified Linear.V1 +import qualified Linear.V2 import qualified Data.Vector.Storable as V import Numeric.DDE.Types --- | Model parameters -data Par = - -- | Mackey-Glass model (no external input) +-- | Mackey-Glass model (no external input) parameters +data MackeyGlass = MackeyGlass { _beta :: Double , _gamma :: Double } - | - -- | Bandpass-filtered two-dimensional nonlinear system with - -- external input @u(t)@. - -- - -- \tau * dx(t)/dt = -x - y / theta + fnl[x(t - \tau_D) + \rho*u(t)] - -- dy(y)/dt = x +-- | Bandpass-filtered two-dimensional nonlinear system with +-- external input @u(t)@. +-- +-- \[ +-- \begin{aligned} +-- \tau dx(t)/dt &= -x - y / \theta + fnl[x(t - \tau_D) + \rho u(t)] \\ +-- dy(y)/dt &= x +-- \end{aligned} +-- \] +data RC = RC { _fnl :: Double -> Double -- ^ Nonlinear transformation , _rho :: Double -- ^ External forcing weight , _filt :: BandpassFiltering @@ -40,27 +47,28 @@ , _theta :: Double -- ^ Integration time, s (delta = \tau_D / theta) } deriving Show --- | DDE right-hand sides for example models -rhs :: Par -> RHS - --- Mackey-Glass model with no external forcing -rhs MackeyGlass { _beta = beta, _gamma = gamma } - ((State xs), (Hist hs), _) - = State (V.singleton x') - where x' = beta * x_tau / (1 + x_tau^(10::Int)) - gamma * x - x = xs V.! 0 - x_tau = hs V.! 0 +-- | Mackey-Glass model with no external forcing +mackeyGlassRhs :: MackeyGlass -> RHS (Linear.V1.V1 Double) +mackeyGlassRhs MackeyGlass { _beta = beta, _gamma = gamma } = RHS _f + where + _f (xs, Hist hs, _) = Linear.V1.V1 x' + where + x' = beta * x_tau / (1 + x_tau^(10::Int)) - gamma * x + x = xs ^._x + x_tau = hs ^._x -- Ikeda-like model with an integral term y(t) and external input -rhs RC { _fnl = _fnl, - _rho = _rho, - _filt = BandpassFiltering { _tau = _tau, _theta = _theta } - } - ((State xs), (Hist hs), (Inp u)) = State $ V.fromList [x', y'] - where - x' = (-x - (recip _theta) * y + _fnl (x_tau + _rho * u)) / _tau - y' = x -- Integral term +bandpassRhs :: RC -> RHS (Linear.V2.V2 Double) +bandpassRhs RC { _fnl = _fnl, + _rho = _rho, + _filt = BandpassFiltering { _tau = tau, _theta = theta } + } = RHS _f + where + _f (xs, Hist hs, Inp u) = Linear.V2.V2 x' y' + where + x' = (-x - y / theta + _fnl (x_tau + _rho * u)) / tau + y' = x -- Integral term - x = xs V.! 0 - y = xs V.! 1 - x_tau = hs V.! 0 -- Delay term + x = xs ^._x + y = xs ^._y + x_tau = hs ^._x -- Delay term diff --git a/dde/Numeric/DDE/Types.hs b/dde/Numeric/DDE/Types.hs --- a/dde/Numeric/DDE/Types.hs +++ b/dde/Numeric/DDE/Types.hs @@ -1,20 +1,27 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE FlexibleContexts #-} module Numeric.DDE.Types ( - RHS + RHS (..) , HistorySnapshot (..) , Input (..) , InputSnapshot (..) - , State (..) - , Stepper1 (..) + , Stepper (..) ) where +import Foreign.Storable ( Storable (..) ) +import qualified Data.VectorSpace.Free as Free import qualified Data.Vector.Storable as V --- | DDE right-hand side -type RHS = (State, HistorySnapshot, InputSnapshot) -> State - --- | State of a dynamical system, it can be a vector of any length (x(t), y(t), ...). -newtype State = State { _state :: V.Vector Double } +-- | DDE right-hand side. +-- +-- Parameter @state@ is and abstraction of a dynamical system's state, +-- i.e. it can be a vector of any length (x(t), y(t), ...). +newtype RHS state = RHS { + _state + :: Free.VectorSpace state => (state, HistorySnapshot state, InputSnapshot) -> state + } -- | Input u(t) is one-dimensional newtype InputSnapshot = Inp { _insnap :: Double } @@ -24,17 +31,32 @@ -- | Contains only the required snapshot of history to make steppers (e.g. Heun) work. -- There could be several delay variables -newtype HistorySnapshot = Hist { _histsnap :: V.Vector Double } +newtype HistorySnapshot state = Hist { _histsnap :: state } --- | Stepper for DDEs with a single delay +-- | DDE stepper (all delays are equal). -- --- >>> _stepper stepSize rhs xyState xTau1_2 u1_2 -newtype Stepper1 = Stepper1 { - _stepper - :: Double - -> RHS - -> State - -> (Double, Double) - -> (Double, Double) - -> State +-- Stepper is a function of the following arguments: +-- +-- * Integration step +-- * DDE right-hand side +-- * Current state vector @(x(t), y(t), ...)@ +-- * Two subsequent history snapshots +-- * Two subsequent inputs +-- +-- The result (step) is a new state vector. +newtype Stepper = Stepper { + _step + :: forall state. ( Functor state, Free.VectorSpace (state Double) + , Num (Free.Scalar (state Double)) ) + => Free.Scalar (state Double) + -> RHS (state Double) + -> state Double + -> (HistorySnapshot (state Double), HistorySnapshot (state Double)) + -> (Double, Double) + -> state Double } +-- NB: to allow multiple delay times, instead of +-- (HistorySnapshot state, HistorySnapshot state) +-- there should be +-- (HistorySnapshot delaystate, HistorySnapshot delaystate). +-- i.e. a vector of required delayed values (e.g. x(t-tau1), x(t-tau2), y(t-tau3)) diff --git a/examples/MackeyGlass/Main.hs b/examples/MackeyGlass/Main.hs --- a/examples/MackeyGlass/Main.hs +++ b/examples/MackeyGlass/Main.hs @@ -1,10 +1,12 @@ +{-# LANGUAGE FlexibleContexts #-} module Main where +import Linear.V1 import Numeric.DDE import Numeric.DDE.Model import qualified Data.Vector.Storable as V -parMG0 :: Par +parMG0 :: MackeyGlass parMG0 = MackeyGlass { _beta = 0.2 , _gamma = 0.1 } @@ -12,18 +14,18 @@ -- | The Mackey-Glass model (with no external input). -- This example demonstrates how to set custom initial conditions. -- Typical hStep = 0.1 -mgModel :: Double -> Int -> V.Vector Double +mgModel :: Double -> Int -> V.Vector (V1 Double) mgModel hStep totalIters = r where -- Initial state x(t0) = 0.2 - state0 = State (V.replicate 1 0.2) - len1 = 17 * (round $ recip hStep) -- tauD = 17, delay time + state0 = V1 0.2 + len1 = 17 * round (recip hStep) -- tauD = 17, delay time -- Initial conditions - hist0 = V.replicate len1 0.2 + hist0 = V.replicate len1 state0 inp = Input (V.replicate (totalIters + 1) 0.0) - rhs' = rhs parMG0 + rhs' = mackeyGlassRhs parMG0 -- Stepper implements Runge-Kutta schema - stepper = let (Stepper1 _rk4) = rk4 + stepper = let (Stepper _rk4) = rk4 in _rk4 hStep rhs' -- Record all the time trace (_, r) = integ' stepper len1 totalIters totalIters (state0, hist0, inp) @@ -54,7 +56,8 @@ runTest = let hStep = 0.1 total = round(1000 / hStep) -- 1000 time units - in mgModel hStep total + trace = mgModel hStep total + in V.map (\(V1 x) -> x) trace main :: IO () main = putStrLn $ toString runTest