dde 0.1.0 → 0.2.0

raw patch · 10 files changed

+221/−27 lines, 10 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• ChangeLog.md +3/−0
• README.md +4/−0
• bench/Bench.hs +52/−2
• bench/Impl2.hs +1/−2
• bench/R1.hs +102/−0
• bench/R2.hs +37/−0
• dde.cabal +4/−2
• dde/Numeric/DDE.hs +14/−14
• dde/Numeric/DDE/Types.hs +3/−5
• examples/MackeyGlass/Main.hs +1/−2

@@ -1,5 +1,8 @@ # Changelog for dde +## 0.2.0 *March 28th 2018*+  * Improved speed for multidimensional DDEs+ ## 0.1.0 *March 21st 2018*   * Now, all dynamical variables are recorded   * Improved speed

@@ -10,3 +10,7 @@    * [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) * Pure Haskell++## Acknowledgements++* [Justus Sagemüller](https://github.com/leftaroundabout) for making this library faster

@@ -1,9 +1,12 @@ import           Criterion.Main import           Linear.V1 ( V1 (..) )+import           Linear.V2 ( V2 (..) ) import qualified Data.Vector.Storable as V  import qualified Impl1 import qualified Impl2+import qualified R1+import qualified R2  runTest1 :: Double -> Double runTest1 maxTime =@@ -19,9 +22,56 @@       (V1 s, _) = Impl2.mgModel hStep total   in s +runTest3 :: V.Vector Double -> Double+runTest3 inp =+  let hStep = recip 256.0+      delaySamples = 800+      result = R1.model hStep delaySamples inp+  in V.last result++runTest4 :: V.Vector Double -> Double+runTest4 inp =+  let hStep = recip 256.0+      delaySamples = 800+      V2 r _ = V.last $ R2.model hStep delaySamples inp+  in r+ main = defaultMain [-    bench "old version" $ whnf runTest1 maxTime-    , bench "target to optimize" $ whnf runTest2 maxTime+    bench "hard coded version" $ whnf runTest1 maxTime+    , bench "dde library version" $ whnf runTest2 maxTime++    , bench "hard coded version (2D case with external forcing)" $ whnf runTest3 inp+    , bench "dde library version (2D case with external forcing)" $ whnf runTest4 inp   ]   where     maxTime = 1000  -- 1000 time units++    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)+--+-- 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)+--+-- 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)+--+-- 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)

@@ -24,7 +24,6 @@     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'+    stepper = rk4 hStep rhs'     -- Provide the last state and the time trace     r = integ' stepper len1 totalIters totalIters (state0, hist0, inp)

@@ -0,0 +1,102 @@+{-# LANGUAGE BangPatterns #-}+module R1 where++import qualified Data.Vector.Storable as V+import qualified Data.Vector.Storable.Mutable as VM+import           System.IO.Unsafe ( unsafePerformIO )+++model :: Double -> Int -> V.Vector R -> V.Vector R+model hStep delaySamples ttrace0 = ttrace1+  where+    sample0 = Sample 0.0 0.0+    hist0 = V.replicate delaySamples 0.0++    iterator = heun2 (hStep / 2.0) (bandpassRhs (tau, delta) $ (\x -> if x < 0 then -x else 0))++    totalIters = V.length ttrace0++    -- Integrator consumes one additional input+    ttrace = ttrace0 V.++ V.singleton (V.last ttrace0)++    (_, ttrace1) = integrator iterator+                                  delaySamples+                                  totalIters+                                  (sample0, hist0, ttrace)++    (tau, delta) = (2.15625, -3.0)++type R = Double++-- Sample is a vector of any length (x, y, z, ...)+data Sample = Sample {-# UNPACK #-} !R {-# UNPACK #-} !R++getX :: Sample -> R+getX (Sample x _) = x+{-# INLINE getX #-}++-- Pair is used to define exactly a pair of values+data Pair = Pair {-# UNPACK #-} !R {-# UNPACK #-} !R++infixl 6 ..+..+(..+..) :: Sample -> Sample -> Sample+(..+..) (Sample x1 y1) (Sample x2 y2) = Sample (x1 + x2) (y1 + y2)+{-# INLINE (..+..) #-}++infixl 7 .*..+(.*..) :: R -> Sample -> Sample+(.*..) c (Sample x2 y2) = Sample (c * x2) (c * y2)+{-# INLINE (.*..) #-}++heun2+  :: R -> ((Sample, R) -> Sample) -> Sample -> Pair -> Sample+heun2 hOver2 f !x !(Pair x_h x_h2) = x_1+  where+    ! f1 = f (x, x_h)+    ! x_1' = x ..+.. 2 * hOver2 .*.. f1+    ! f2 = f (x_1', x_h2)+    ! x_1 = x ..+.. hOver2 .*.. (f1 ..+.. f2)++bandpassRhs+  :: (R, R) -> (R -> R) -> ((Sample, R) -> Sample)+bandpassRhs !(!ε, !δ) f (!(Sample x y), !x_h) = Sample x' y'+  where+    ! x' = (recip ε) * (-x - δ * y + (f x_h))+    ! y' = x++integrator+  :: (Sample -> Pair -> Sample)+    -> Int+    -> Int+    -> (Sample, V.Vector R, V.Vector R)+    -- ^ Dynamical state, delay history, external forcing+    ->  (Sample, V.Vector R)+integrator iterate1 len total (!xy0, !history0, !input) = unsafePerformIO $ do+    ! v <- VM.new total+    xy <- go v 0 xy0+    history <- V.unsafeFreeze v+    return (xy, history)+  where+    addInput !val !i = val + inputVal+      where ! inputVal = input V.! i+    {-# INLINE addInput #-}+    h i = addInput (history0 `V.unsafeIndex` i) i+    {-# INLINE h #-}+    go !v !i !xy+      | i < len - 1 = do+        let !r = iterate1 xy (Pair (h i) (h $ i + 1))+        VM.unsafeWrite v i (getX r)+        go v (i + 1) r+      | i == total = do+        return xy+      -- Iterations after the initial history has been exhausted+      | otherwise = do+        ! newX0 <- if i == len - 1+                      then return (getX xy0)+                      else VM.unsafeRead v (i - len - 1)+        ! newX <- VM.unsafeRead v (i - len)+        let !r = iterate1 xy (Pair (addInput newX0 i) (addInput newX (i + 1)))+        VM.unsafeWrite v i (getX r)+        go v (i + 1) r+{-# INLINE integrator #-}+

@@ -0,0 +1,37 @@+{-# LANGUAGE BangPatterns #-}+module R2 ( model ) where++import           Linear ( V2 (..) )+import qualified Data.Vector.Storable as V+import qualified Numeric.DDE as DDE++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'+      where+        tau = 2.15625+        delta = -3.0++        f x = if x < 0 then -x else 0++        x' = (-x - delta * y + f (x_tau + inp)) / tau+        y' = x+        -- No input amplification, i.e. rho = 1 in rho * inp++model' :: Double -> Int -> V.Vector Double -> (V2 Double, V.Vector (V2 Double))+model' dt len1 inp = res+  where+    -- Initial conditions:+    -- dynamical state and delay history.+    state0 = V2 0.0 0.0+    hist0 = V.replicate len1 state0++    inp' = DDE.Input 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 =+  let (state1, traces) = model' dt delaySamples inp+  in traces

@@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 005f02aea25cd38967b9112ce3aeb1c77f570482c9c94bdcebc344346e26de2d+-- hash: 32559758517286304302f4a45a26a227e1e595de2376417a01c2b0e313ee2bbc  name:           dde-version:        0.1.0+version:        0.2.0 synopsis:       Delay differential equations description:    Please see the README on Github at <https://github.com/masterdezign/dde#readme> category:       Math@@ -93,6 +93,8 @@   other-modules:       Impl1       Impl2+      R1+      R2       Numeric.DDE       Numeric.DDE.Model       Numeric.DDE.Types

@@ -51,6 +51,7 @@  {-# LANGUAGE BangPatterns #-} {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE RankNTypes #-} {-# LANGUAGE FlexibleContexts #-}  module Numeric.DDE (@@ -84,9 +85,7 @@  -- | Fourth order Runge-Kutta stepper rk4 :: Stepper-rk4 = Stepper _rk4-  where-    _rk4 dt (RHS rhs') !xy (Hist !xy_tau1, Hist !xy_tau1') (!u1, !u1') = xy_next+rk4 dt (RHS rhs') xy (Hist xy_tau1, Hist xy_tau1') (u1, u1') = xy_next       where         xy_next = xy ^+^ (a ^+^ 2 *^ b ^+^ 2 *^ c ^+^ d) ^/ 6         a = dt *^ rhs' (xy, Hist xy_tau1, inp1)@@ -102,14 +101,13 @@  -- | Second order Heun's stepper heun2 :: Stepper-heun2 = Stepper _heun2-  where-    _heun2 hStep (RHS rhs') !xy (!xy_tau1, !xy_tau1') (!u1, !u1') = xy_next+heun2 hStep (RHS rhs') xy (xy_tau1, xy_tau1') (u1, u1') = xy_next       where         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+{-# INLINE heun2 #-}  -- | Generic integrator for DDEs (single delay time). -- Records all dynamical variables.@@ -129,10 +127,10 @@   -> (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' iter1 len1 krecord total (xy0, hist0, Input in1) = a   where     a = unsafePerformIO $ do-      ! v <- VM.new (len1 + total)  -- Delay history+      v <- VM.new (len1 + total)  -- Delay history       -- Copy the initial history values       copyHist v hist0 @@ -143,7 +141,7 @@       return (xy', V.slice (len1 + total - krecord) krecord trace)      -- Copy initial conditions-    copyHist !v !hist =+    copyHist v hist =       mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..V.length hist - 1]      go !v !i !xy@@ -154,9 +152,10 @@         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 (Hist xy_tau1, Hist xy_tau1') (u1, u1')+            xy' = iter1 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).@@ -170,14 +169,15 @@   -> RHS (state Double)  -- ^ Derivative (DDE right-hand side)   -> Input  -- ^ External forcing   -> (state Double, V.Vector (state Double))-integ (Stepper stp) state0 hist0 len1 dt rhs' inp@(Input in1) = r+integ 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 dt rhs'+    totalIters = V.length in1 - 1+    iterator = stp dt rhs'     -- Record all the time trace-    ! r = integ' iterator len1 totalIters totalIters (state0, hist0, inp)+    r = integ' iterator len1 totalIters totalIters (state0, hist0, inp)+{-# INLINE integ #-}  -- | RK4 integrator shortcut for 1D DDEs with zero -- initial conditions

@@ -20,7 +20,7 @@ -- 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+    :: (state, HistorySnapshot state, InputSnapshot) -> state   }  -- | Input u(t) is one-dimensional@@ -44,9 +44,8 @@ -- * Two subsequent inputs -- -- The result (step) is a new state vector.-newtype Stepper = Stepper {-  _step-    :: forall state. ( Functor state, Free.VectorSpace (state Double)+type Stepper = +       forall state. ( Functor state, Free.VectorSpace (state Double)                      , Num (Free.Scalar (state Double)) )      => Free.Scalar (state Double)      -> RHS (state Double)@@ -54,7 +53,6 @@      -> (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

@@ -25,8 +25,7 @@     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'+    stepper = rk4 hStep rhs'     -- Record all the time trace     (_, r) = integ' stepper len1 totalIters totalIters (state0, hist0, inp)