diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -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) diff --git a/bench/Impl2.hs b/bench/Impl2.hs --- a/bench/Impl2.hs +++ b/bench/Impl2.hs @@ -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) diff --git a/bench/R2.hs b/bench/R2.hs --- a/bench/R2.hs +++ b/bench/R2.hs @@ -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 = 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: 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 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 diff --git a/dde/Numeric/DDE.hs b/dde/Numeric/DDE.hs --- a/dde/Numeric/DDE.hs +++ b/dde/Numeric/DDE.hs @@ -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 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 @@ -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 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 @@ -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)) diff --git a/examples/Chimera/Main.hs b/examples/Chimera/Main.hs new file mode 100644 --- /dev/null +++ b/examples/Chimera/Main.hs @@ -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 + diff --git a/examples/MackeyGlass/Main.hs b/examples/MackeyGlass/Main.hs --- a/examples/MackeyGlass/Main.hs +++ b/examples/MackeyGlass/Main.hs @@ -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