diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -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) diff --git a/bench/Impl2.hs b/bench/Impl2.hs --- a/bench/Impl2.hs +++ b/bench/Impl2.hs @@ -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) diff --git a/bench/R1.hs b/bench/R1.hs new file mode 100644 --- /dev/null +++ b/bench/R1.hs @@ -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 #-} + diff --git a/bench/R2.hs b/bench/R2.hs new file mode 100644 --- /dev/null +++ b/bench/R2.hs @@ -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 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: 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 category: Math @@ -93,6 +93,8 @@ other-modules: Impl1 Impl2 + R1 + R2 Numeric.DDE Numeric.DDE.Model Numeric.DDE.Types diff --git a/dde/Numeric/DDE.hs b/dde/Numeric/DDE.hs --- a/dde/Numeric/DDE.hs +++ b/dde/Numeric/DDE.hs @@ -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 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 @@ -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 diff --git a/examples/MackeyGlass/Main.hs b/examples/MackeyGlass/Main.hs --- a/examples/MackeyGlass/Main.hs +++ b/examples/MackeyGlass/Main.hs @@ -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)