dde 0.0.1 → 0.1.0

raw patch · 10 files changed

+490/−157 lines, 10 filesdep +criteriondep +free-vector-spacesdep +lensPVP ok

version bump matches the API change (PVP)

Dependencies added: criterion, free-vector-spaces, lens, linear

API changes (from Hackage documentation)

Files

• ChangeLog.md +6/−1
• README.md +2/−3
• bench/Bench.hs +27/−0
• bench/Impl1.lhs +190/−0
• bench/Impl2.hs +30/−0
• dde.cabal +33/−2
• dde/Numeric/DDE.hs +106/−88
• dde/Numeric/DDE/Model.hs +44/−36
• dde/Numeric/DDE/Types.hs +41/−19
• examples/MackeyGlass/Main.hs +11/−8

@@ -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

@@ -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-

@@ -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

@@ -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}

@@ -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)

@@ -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 <https://github.com/masterdezign/dde#readme> 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

@@ -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

@@ -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

@@ -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))

@@ -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