diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,3 @@ +# Changelog for dde + +## Unreleased changes diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Bogdan Penkovsky (c) 2018 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Bogdan Penkovsky nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,11 @@ +# Delay Differential Equations (DDE) + +## Features + +* DDEs with multiple dynamical variables and a single delay (pull requests are welcome) +* Driven systems (with external input) +* Second and fourth order integration methods +* Example models: + * Mackey-Glass with no external input + * Driven system +* Pure Haskell diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/dde.cabal b/dde.cabal new file mode 100644 --- /dev/null +++ b/dde.cabal @@ -0,0 +1,69 @@ +-- This file has been generated from package.yaml by hpack version 0.20.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: a8f746086619c8bfab833f302e37ef3d130b22183d5388e5636b5efe2b4a1052 + +name: dde +version: 0.0.0 +synopsis: Delay differential equations +description: Please see the README on Github at +category: Math +homepage: https://github.com/masterdezign/dde#readme +bug-reports: https://github.com/masterdezign/dde/issues +author: Bogdan Penkovsky +maintainer: dev () penkovsky dot com +copyright: Bogdan Penkovsky +license: BSD3 +license-file: LICENSE +build-type: Simple +cabal-version: >= 1.10 + +extra-source-files: + ChangeLog.md + README.md + +source-repository head + type: git + location: https://github.com/masterdezign/dde + +library + hs-source-dirs: + dde + build-depends: + base >=4.7 && <5 + , vector + exposed-modules: + Numeric.DDE + Numeric.DDE.Model + Numeric.DDE.Types + other-modules: + Paths_dde + default-language: Haskell2010 + +executable mackey-glass + main-is: Main.hs + hs-source-dirs: + examples/MackeyGlass + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , dde + , vector + other-modules: + Paths_dde + default-language: Haskell2010 + +test-suite dde-test + type: exitcode-stdio-1.0 + main-is: Spec.hs + hs-source-dirs: + test + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , dde + , vector + other-modules: + Paths_dde + default-language: Haskell2010 diff --git a/dde/Numeric/DDE.hs b/dde/Numeric/DDE.hs new file mode 100644 --- /dev/null +++ b/dde/Numeric/DDE.hs @@ -0,0 +1,208 @@ +{- | + = Delay differential equations (DDE) + + + == Example: Ikeda DDE + + Below is a complete example simulating the Ikeda DDE defined as: + @tau * x(t)/dt = -x + beta * sin[x(t - tau_D)]@. + + > 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'] + > where + > -- Ikeda DDE definition + > x' = (-x + beta * (sin x_tauD)) / tau + > + > -- Constants + > tau = 0.01 + > + > -- Dynamical variable x(t) + > x = V.head xs + > + > -- Delay term x(t - tau_D) + > x_tauD = V.head hs + > + > model beta hStep len1 totalIter = DDE.integ DDE.rk4 state0 hist0 len1 hStep (ikedaRhs beta) inp + > where + > -- Initial conditions: + > -- dynamical state and delay history. + > state0 = DDE.State $ V.fromList [pi/2] + > hist0 = V.replicate len1 (pi/2) + > + > -- Input is ignored in ikedaRhs + > inp = DDE.Input $ V.replicate (totalIter + 1) 0 + > + > -- Control parameter + > beta = 2.6 + > + > main = do + > let hStep = 0.001 -- Integration step + > samplesPerDelay = round $ 1.0 / hStep + > delays = 8 + > total = delays * samplesPerDelay + > + > let (state1, trace) = model beta hStep samplesPerDelay total + > + > mapM_ print $ V.toList trace + +-} + +{-# LANGUAGE BangPatterns #-} +module Numeric.DDE ( + -- * Integrators + integ + , integ' + , integRk4_2D + , integHeun2_2D + , Input (..) + , State (..) + , HistorySnapshot (..) + + -- * Steppers + , Stepper1 (..) + , rk4 + , heun2 + ) where + +import qualified Data.Vector.Storable as V +import qualified Data.Vector.Storable.Mutable as VM +import System.IO.Unsafe ( unsafePerformIO ) + +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 + where + _rk4 hStep rhs' (State !xy) !(!x_tau1, !x_tau1') !(!u1, !u1') = State 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 + +-- | Second order Heun's stepper +heun2 :: Stepper1 +heun2 = Stepper1 _heun2 + where + _heun2 hStep rhs' (State !xy) !(!x_tau1, !x_tau1') !(!u1, !u1') = State 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) + +-- | Generic integrator for DDEs with a single delay +integ' + :: (State -> (Double, Double) -> (Double, Double) -> State) + -- ^ Iterator describing a DDE system + -> Int + -- ^ Delay length in samples + -> Int + -- ^ Number of last samples to record + -> Int + -- ^ Total number of iterations + -> (State, V.Vector Double, Input) + -- ^ Initial state vector, initial history, and external forcing + -> (State, V.Vector Double) + -- ^ 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 + 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 = do + mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..(V.length hist) - 1] + + go !v !i !xy + | i == len1 + total = do + 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 + 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' + go v (i + 1) xy' + +-- | Generic integrator that records the whole time trace $x(t)$ +integ + :: Stepper1 + -> State -- ^ Initial state x(t), y(t),... + -> V.Vector Double -- ^ Initial history for the delayed variable + -> Int -- ^ Delay length in samples + -> Double + -> RHS + -> Input -- ^ External forcing + -> (State, V.Vector Double) +integ (Stepper1 stp) state0 hist0 len1 hStep 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' + -- Record all the time trace + ! r = integ' iterator len1 totalIters totalIters (state0, hist0, inp) + +-- | RK4 integrator shortcut for 2D DDEs with zero +-- initial conditions +integRk4_2D :: Int -- ^ Delay length in samples + -> Double -- ^ Integration time step + -> RHS -- ^ DDE model + -> Input -- ^ External forcing + -> (State, V.Vector Double) +integRk4_2D len1 = integ rk4 state0 hist0 len1 + where + ! state0 = State (V.replicate 2 0.0) + ! hist0 = V.replicate len1 0 + +-- | 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 + -> Input -- ^ External forcing + -> (State, V.Vector Double) +integHeun2_2D len1 = integ heun2 state0 hist0 len1 + where + ! state0 = State (V.replicate 2 0.0) + ! hist0 = V.replicate len1 0 diff --git a/dde/Numeric/DDE/Model.hs b/dde/Numeric/DDE/Model.hs new file mode 100644 --- /dev/null +++ b/dde/Numeric/DDE/Model.hs @@ -0,0 +1,66 @@ +{- | + = 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 + ) where + +import qualified Data.Vector.Storable as V + +import Numeric.DDE.Types + + +-- | Model parameters +data Par = + -- | Mackey-Glass model (no external input) + 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 + RC { _fnl :: Double -> Double -- ^ Nonlinear transformation + , _rho :: Double -- ^ External forcing weight + , _filt :: BandpassFiltering + } + +-- | Bandpass filter (linear) parameters +data BandpassFiltering = BandpassFiltering + { _tau :: Double -- ^ System response time, s (epsilon = \tau / \tau_D) + , _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 + +-- 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 + + x = xs V.! 0 + y = xs V.! 1 + x_tau = hs V.! 0 -- Delay term diff --git a/dde/Numeric/DDE/Types.hs b/dde/Numeric/DDE/Types.hs new file mode 100644 --- /dev/null +++ b/dde/Numeric/DDE/Types.hs @@ -0,0 +1,40 @@ +module Numeric.DDE.Types ( + RHS + , HistorySnapshot (..) + , Input (..) + , InputSnapshot (..) + , State (..) + , Stepper1 (..) + ) where + +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 } + +-- | Input u(t) is one-dimensional +newtype InputSnapshot = Inp { _insnap :: Double } + +-- | 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 (as in rc-analog) +newtype HistorySnapshot = Hist { _histsnap :: V.Vector Double } + +-- | Stepper for DDEs with a single delay +-- +-- >>> _stepper stepSize rhs xyState xTau1_2 u1_2 +newtype Stepper1 = Stepper1 { + _stepper + :: Double + -> RHS + -> State + -> (Double, Double) + -> (Double, Double) + -> State + } diff --git a/examples/MackeyGlass/Main.hs b/examples/MackeyGlass/Main.hs new file mode 100644 --- /dev/null +++ b/examples/MackeyGlass/Main.hs @@ -0,0 +1,62 @@ +module Main where + +import Numeric.DDE +import Numeric.DDE.Model +import qualified Data.Vector.Storable as V + +parMG0 :: Par +parMG0 = MackeyGlass { _beta = 0.2 + , _gamma = 0.1 + } + +-- | 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 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 + -- Initial conditions + hist0 = V.replicate len1 0.2 + inp = Input (V.replicate (totalIters + 1) 0.0) + rhs' = rhs parMG0 + -- Stepper implements Runge-Kutta schema + stepper = let (Stepper1 _rk4) = rk4 + in _rk4 hStep rhs' + -- Record all the time trace + (_, r) = integ' stepper len1 totalIters totalIters (state0, hist0, inp) + +-- | Comparison with the output.dat produced by: +-- > $ xppaut -silent mg.ode +-- +-- > $ cat mg.ode +-- +-- x'=beta*delay(x,tau)/(1 + delay(x,tau)^10)-gamma*x +-- +-- par tau=17 +-- par gamma=0.1 +-- par beta=0.2 +-- +-- x(0)=0.2 +-- init x=0.2 +-- +-- @ total=1000 +-- @ bounds=10000 +-- # Maximal delay buffer (in time units) +-- @ delay=100 +-- @ dt=0.1 +-- @ trans=0 +-- @ maxstore=100000000 +-- done +runTest :: V.Vector Double +runTest = + let hStep = 0.1 + total = round(1000 / hStep) -- 1000 time units + in mgModel hStep total + +main :: IO () +main = putStrLn $ toString runTest + where + toString = unlines. map show. V.toList diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,2 @@ +main :: IO () +main = putStrLn "Test suite not yet implemented"