dde 0.2.0 → 0.3.0
raw patch · 11 files changed
+163/−73 lines, 11 filesnew-component:exe:chimera
Files
- ChangeLog.md +3/−0
- README.md +3/−2
- bench/Bench.hs +21/−22
- bench/Impl2.hs +1/−1
- bench/R2.hs +5/−2
- dde.cabal +17/−2
- dde/Numeric/DDE.hs +39/−25
- dde/Numeric/DDE/Model.hs +4/−4
- dde/Numeric/DDE/Types.hs +8/−14
- examples/Chimera/Main.hs +57/−0
- examples/MackeyGlass/Main.hs +5/−1
ChangeLog.md view
@@ -1,5 +1,8 @@ # Changelog for dde +## 0.3.0 *July 4th 2018*+ * Support DDEs with multiple delay times+ ## 0.2.0 *March 28th 2018* * Improved speed for multidimensional DDEs
README.md view
@@ -2,13 +2,14 @@ ## Features -* Autonomous DDEs with multiple dynamical variables and a single delay time (pull requests are welcome)+* Autonomous DDEs with multiple dynamical variables and multiple delay times * 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/d22c6ff82fd56c29289366a057f3d733a23844d0/dde/Numeric/DDE/Model.hs#L60)+ * [Driven system](https://github.com/masterdezign/dde/blob/80de41dd8c04c18cf439dee33bc117c96c212699/dde/Numeric/DDE/Model.hs#L60)+ * [Two delays](https://github.com/masterdezign/dde/blob/master/examples/Chimera/Main.hs) * Pure Haskell ## Acknowledgements
bench/Bench.hs view
@@ -49,29 +49,28 @@ total = round(maxTime * 256) inp = V.fromList $ map (sin. (0.001*pi*). fromIntegral) [1..total] --- benchmarking hard coded version--- time 6.546 ms (6.471 ms .. 6.627 ms)--- 0.998 R² (0.996 R² .. 0.999 R²)--- mean 6.755 ms (6.673 ms .. 6.871 ms)--- std dev 272.1 μs (193.4 μs .. 376.0 μs)--- variance introduced by outliers: 18% (moderately inflated)++-- Benchmark was run on Dell Precision T3610 --+-- time 6.377 ms (6.326 ms .. 6.430 ms)+-- 0.999 R² (0.999 R² .. 1.000 R²)+-- mean 6.433 ms (6.397 ms .. 6.481 ms)+-- std dev 125.5 μs (88.80 μs .. 170.0 μs)+ -- benchmarking dde library version--- time 6.860 ms (6.786 ms .. 6.935 ms)--- 0.998 R² (0.996 R² .. 1.000 R²)--- mean 6.843 ms (6.796 ms .. 6.917 ms)--- std dev 163.7 μs (112.2 μs .. 262.5 μs)---+-- time 6.570 ms (6.559 ms .. 6.581 ms)+-- 1.000 R² (1.000 R² .. 1.000 R²)+-- mean 6.549 ms (6.538 ms .. 6.557 ms)+-- std dev 26.70 μs (22.69 μs .. 31.91 μs)+ -- benchmarking hard coded version (2D case with external forcing)--- time 3.644 ms (3.595 ms .. 3.697 ms)--- 0.999 R² (0.998 R² .. 0.999 R²)--- mean 3.677 ms (3.645 ms .. 3.715 ms)--- std dev 106.5 μs (86.91 μs .. 130.9 μs)--- variance introduced by outliers: 12% (moderately inflated)---+-- time 3.641 ms (3.604 ms .. 3.695 ms)+-- 0.998 R² (0.996 R² .. 1.000 R²)+-- mean 3.659 ms (3.640 ms .. 3.693 ms)+-- std dev 74.47 μs (49.20 μs .. 124.5 μs)+ -- benchmarking dde library version (2D case with external forcing)--- time 6.118 ms (6.026 ms .. 6.200 ms)--- 0.999 R² (0.998 R² .. 0.999 R²)--- mean 6.241 ms (6.186 ms .. 6.336 ms)--- std dev 201.6 μs (150.1 μs .. 275.6 μs)--- variance introduced by outliers: 13% (moderately inflated)+-- time 5.964 ms (5.839 ms .. 6.067 ms)+-- 0.998 R² (0.997 R² .. 0.999 R²)+-- mean 5.958 ms (5.928 ms .. 5.989 ms)+-- std dev 96.24 μs (73.59 μs .. 144.5 μs)
bench/Impl2.hs view
@@ -26,4 +26,4 @@ -- Stepper implements Runge-Kutta schema stepper = rk4 hStep rhs' -- Provide the last state and the time trace- r = integ' stepper len1 totalIters totalIters (state0, hist0, inp)+ r = integ' stepper [len1] totalIters totalIters (state0, hist0, inp)
bench/R2.hs view
@@ -8,13 +8,16 @@ rhs :: DDE.RHS (V2 Double) rhs = DDE.RHS deriv where- deriv (V2 !x !y, DDE.Hist (V2 !x_tau _), DDE.Inp !inp) = V2 x' y'+ deriv (V2 !x !y, DDE.Hist snapshots, DDE.Inp !inp) = V2 x' y' where tau = 2.15625 delta = -3.0 f x = if x < 0 then -x else 0 + -- Delay term+ V2 !x_tau _ = head snapshots+ x' = (-x - delta * y + f (x_tau + inp)) / tau y' = x -- No input amplification, i.e. rho = 1 in rho * inp@@ -29,7 +32,7 @@ inp' = DDE.Input inp - res = DDE.integ DDE.heun2 state0 hist0 len1 dt rhs inp'+ res = DDE.integ DDE.heun2 state0 hist0 [len1] dt rhs inp' model :: Double -> Int -> V.Vector Double -> V.Vector (V2 Double) model dt delaySamples inp =
dde.cabal view
@@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 32559758517286304302f4a45a26a227e1e595de2376417a01c2b0e313ee2bbc+-- hash: 80569a250921ebee4e46574543eab800c2f9b07d4a0f2669316669a10b390ca4 name: dde-version: 0.2.0+version: 0.3.0 synopsis: Delay differential equations description: Please see the README on Github at <https://github.com/masterdezign/dde#readme> category: Math@@ -40,6 +40,21 @@ Numeric.DDE Numeric.DDE.Model Numeric.DDE.Types+ other-modules:+ Paths_dde+ default-language: Haskell2010++executable chimera+ main-is: Main.hs+ hs-source-dirs:+ examples/Chimera+ build-depends:+ base >=4.7 && <5+ , dde+ , free-vector-spaces+ , lens+ , linear+ , vector other-modules: Paths_dde default-language: Haskell2010
dde/Numeric/DDE.hs view
@@ -13,11 +13,14 @@ > > ikedaRhs beta = DDE.RHS derivative > where- > derivative ((V1 x), (DDE.Hist (V1 x_tauD)), _) = V1 x'+ > derivative ((V1 x), (DDE.Hist histSnapshots), _) = V1 x' > where > -- Ikeda DDE definition > x' = (-x + beta * (sin x_tauD)) / tau >+ > -- There is only a single delay in our model+ > V1 x_tauD = head histSnapshots+ > > -- Constants > tau = 0.01 >@@ -31,7 +34,9 @@ > -- 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+ > -- Only one delay+ > delaysInSamples = [len1]+ > (state1, trace) = DDE.integ DDE.rk4 state0 hist0 delaysInSamples hStep (ikedaRhs beta) inp > > -- Control parameter > beta = 2.6@@ -64,7 +69,7 @@ , integHeun2_2D , Input (..) , InputSnapshot (..)- , HistorySnapshot (..)+ , HistorySnapshots (..) -- * Steppers , RHS (..)@@ -92,7 +97,7 @@ 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_b = zipWith (\xp xq -> (xp ^+^ xq) ^/ 2) xy_tau1 xy_tau1' xy_tau1_c = xy_tau1_b inp1 = Inp u1 inp1_b = Inp $ (u1 + u1') / 2@@ -109,15 +114,14 @@ xy_next = xy ^+^ (hStep *^ (f1 ^+^ f2)) ^/ 2.0 {-# INLINE heun2 #-} --- | Generic integrator for DDEs (single delay time).+-- | Generic integrator for DDEs. -- Records all dynamical variables.--- integ' :: Storable state- => (state -> (HistorySnapshot state, HistorySnapshot state) -> (Double, Double) -> state)+ => (state -> (HistorySnapshots state, HistorySnapshots state) -> (Double, Double) -> state) -- ^ Iterator describing a DDE system- -> Int- -- ^ Delay length in samples+ -> [Int]+ -- ^ Delay lengths in samples -> Int -- ^ Number of last samples to record -> Int@@ -127,9 +131,13 @@ -> (state, V.Vector state) -- ^ Final state and recorded state of the first variable. -- The latter is a vector of vectors (matrix) when multiple variables are involved.-integ' iter1 len1 krecord total (xy0, hist0, Input in1) = a+integ' iter lengths krecord total (xy0, hist0, Input in1) = a where+ len1 = maximum lengths+ a = unsafePerformIO $ do+ -- The longest delay (in samples)+ v <- VM.new (len1 + total) -- Delay history -- Copy the initial history values copyHist v hist0@@ -141,6 +149,7 @@ return (xy', V.slice (len1 + total - krecord) krecord trace) -- Copy initial conditions+ -- It should be asserted that V.length hist >= len1 copyHist v hist = mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..V.length hist - 1] @@ -148,23 +157,24 @@ | i == len1 + total = return xy | otherwise = do- 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+ xy_tau1 <- mapM (\len -> VM.unsafeRead v (i - len)) lengths+ -- Note that xy_tau1 are delayed by one (discrete sample) values from xy_tau1'.+ -- Perhaps, memory access could be somehow optimized.+ xy_tau1' <- mapM (\len -> VM.unsafeRead v (i - len + 1)) lengths+ let u1 = in1 V.! (i - len1) -- Two subsequent scalar inputs u1' = in1 V.! (i - len1 + 1)- xy' = iter1 xy (Hist xy_tau1, Hist xy_tau1') (u1, u1')+ xy' = iter xy (Hist xy_tau1, Hist xy_tau1') (u1, u1') VM.unsafeWrite v i xy' go v (i + 1) xy' {-# INLINE integ' #-} --- | Generic integrator that records the whole time trace @x(t)@--- (single delay time).+-- | Generic integrator that records the whole time trace @x(t)@. integ :: (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+ -> [Int] -- ^ Delay lengths in samples -> Scalar (state Double) -- ^ Integration step -> RHS (state Double) -- ^ Derivative (DDE right-hand side) -> Input -- ^ External forcing@@ -181,48 +191,52 @@ -- | RK4 integrator shortcut for 1D DDEs with zero -- initial conditions-integRk4 :: Int -- ^ Delay length in samples+integRk4 :: [Int] -- ^ Delay lengths 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+integRk4 lengths = integ rk4 state0 hist0 lengths where state0 = V1 0.0+ len1 = maximum lengths hist0 = V.replicate len1 state0 -- | Shortcut for Heun's 2nd order 1D DDEs with zero -- initial conditions-integHeun2 :: Int -- ^ Delay length in samples+integHeun2 :: [Int] -- ^ Delay lengths 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+integHeun2 lengths = integ heun2 state0 hist0 lengths where state0 = V1 0.0+ len1 = maximum lengths hist0 = V.replicate len1 state0 -- | RK4 integrator shortcut for 2D DDEs with zero -- initial conditions-integRk4_2D :: Int -- ^ Delay length in samples+integRk4_2D :: [Int] -- ^ Delay lengths in samples -> Double -- ^ Integration time step -> RHS (V2 Double) -- ^ DDE model -> Input -- ^ External forcing -> (V2 Double, V.Vector (V2 Double))-integRk4_2D len1 = integ rk4 state0 hist0 len1+integRk4_2D lengths = integ rk4 state0 hist0 lengths where state0 = V2 0.0 0.0+ len1 = maximum lengths hist0 = V.replicate len1 state0 -- | Shortcut for Heun's 2nd order 2D DDEs with zero -- initial conditions-integHeun2_2D :: Int -- ^ Delay length in samples+integHeun2_2D :: [Int] -- ^ Delay length in samples -> Double -- ^ Integration time step -> RHS (V2 Double) -- ^ DDE model -> Input -- ^ External forcing -> (V2 Double, V.Vector (V2 Double))-integHeun2_2D len1 = integ heun2 state0 hist0 len1+integHeun2_2D lengths = integ heun2 state0 hist0 lengths where state0 = V2 0.0 0.0+ len1 = maximum lengths hist0 = V.replicate len1 state0
dde/Numeric/DDE/Model.hs view
@@ -51,11 +51,11 @@ mackeyGlassRhs :: MackeyGlass -> RHS (Linear.V1.V1 Double) mackeyGlassRhs MackeyGlass { _beta = beta, _gamma = gamma } = RHS _f where- _f (xs, Hist hs, _) = Linear.V1.V1 x'+ _f (xs, Hist snapshots, _) = Linear.V1.V1 x' where x' = beta * x_tau / (1 + x_tau^(10::Int)) - gamma * x x = xs ^._x- x_tau = hs ^._x+ x_tau = (head snapshots) ^._x -- Ikeda-like model with an integral term y(t) and external input bandpassRhs :: RC -> RHS (Linear.V2.V2 Double)@@ -64,11 +64,11 @@ _filt = BandpassFiltering { _tau = tau, _theta = theta } } = RHS _f where- _f (xs, Hist hs, Inp u) = Linear.V2.V2 x' y'+ _f (xs, Hist snapshots, Inp u) = Linear.V2.V2 x' y' where x' = (-x - y / theta + _fnl (x_tau + _rho * u)) / tau y' = x -- Integral term x = xs ^._x y = xs ^._y- x_tau = hs ^._x -- Delay term+ x_tau = (head snapshots) ^._x -- A single delay term
dde/Numeric/DDE/Types.hs view
@@ -3,7 +3,7 @@ {-# LANGUAGE FlexibleContexts #-} module Numeric.DDE.Types ( RHS (..)- , HistorySnapshot (..)+ , HistorySnapshots (..) , Input (..) , InputSnapshot (..) , Stepper (..)@@ -20,7 +20,7 @@ -- i.e. it can be a vector of any length (x(t), y(t), ...). newtype RHS state = RHS { _state- :: (state, HistorySnapshot state, InputSnapshot) -> state+ :: (state, HistorySnapshots state, InputSnapshot) -> state } -- | Input u(t) is one-dimensional@@ -29,32 +29,26 @@ -- | Vector of input data points newtype Input = Input { _input :: V.Vector Double } --- | Contains only the required snapshot of history to make steppers (e.g. Heun) work.--- There could be several delay variables-newtype HistorySnapshot state = Hist { _histsnap :: state }+-- | Contains state snapshots corresponding to each required delay length+newtype HistorySnapshots state = Hist { _histsnaps :: [state] } --- | DDE stepper (all delays are equal).+-- | DDE stepper -- -- 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 history snapshot lists -- * Two subsequent inputs -- -- The result (step) is a new state vector.-type Stepper = +type Stepper = 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))+ -> (HistorySnapshots (state Double), HistorySnapshots (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))
+ examples/Chimera/Main.hs view
@@ -0,0 +1,57 @@+import Linear ( V2 (..) )+import qualified Data.Vector.Storable as V+import qualified Numeric.DDE as DDE++rhs phi0 = DDE.RHS derivative+ where+ derivative ((V2 x y), (DDE.Hist histSnapshots), _) = V2 x' y'+ where+ -- DDE Eq. (4) from arXiv:1712.03283+ x' = (-x - delta * y + (1 - gamma) * f x_tau1 + gamma * f x_tau2) / epsilon+ y' = x++ f = airy phi0++ -- Delay terms where tau2 / tau1 = 100+ (V2 x_tau1 _):(V2 x_tau2 _):_ = histSnapshots++ -- Constants+ epsilon = 0.01+ gamma = 0.5+ delta = 0.009++airy phi0 x = beta / (1 + m * (sin (x + phi0))^2)+ where+ m = 50+ beta = 1.6++model phi0 hStep len1 len2 totalIter = (state1, V.map (\(V2 x y) -> x) trace)+ where+ -- Initial conditions:+ -- dynamical state and delay history.+ state0 = V2 0.0 0.0+ hist0 = V.fromList $ map (\n -> let x = sin(2 * pi * fromIntegral n / 1000)+ in V2 x 0.0) [1..len2]++ -- Input is ignored in `rhs`+ inp = DDE.Input $ V.replicate (totalIter + 1) 0++ delaysInSamples = [len1, len2]++ (state1, trace) = DDE.integ DDE.rk4 state0 hist0 delaysInSamples hStep (rhs phi0) inp++-- Control parameter+phi0 = -0.45++main = do+ let hStep = 0.0005 -- Integration step+ tauD1 = 1.0 -- Short delay time+ len1 = round $ tauD1 / hStep -- Samples per short delay+ len2 = 100 * len1 -- Long delay (in samples)+ delays = 10 -- 10 long delays+ total = delays * len2++ let (state1, trace) = model phi0 hStep len1 len2 total++ mapM_ print $ V.toList trace+
examples/MackeyGlass/Main.hs view
@@ -26,8 +26,12 @@ rhs' = mackeyGlassRhs parMG0 -- Stepper implements Runge-Kutta schema stepper = rk4 hStep rhs'++ -- A single delay+ delaysInSamples = [len1]+ -- Record all the time trace- (_, r) = integ' stepper len1 totalIters totalIters (state0, hist0, inp)+ (_, r) = integ' stepper delaysInSamples totalIters totalIters (state0, hist0, inp) -- | Comparison with the output.dat produced by: -- > $ xppaut -silent mg.ode