packages feed

dde (empty) → 0.0.0

raw patch · 10 files changed

+493/−0 lines, 10 filesdep +basedep +ddedep +vectorsetup-changed

Dependencies added: base, dde, vector

Files

+ ChangeLog.md view
@@ -0,0 +1,3 @@+# Changelog for dde++## Unreleased changes
+ LICENSE view
@@ -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.
+ README.md view
@@ -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
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ dde.cabal view
@@ -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 <https://github.com/masterdezign/dde#readme>+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
+ dde/Numeric/DDE.hs view
@@ -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
+ dde/Numeric/DDE/Model.hs view
@@ -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
+ dde/Numeric/DDE/Types.hs view
@@ -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+  }
+ examples/MackeyGlass/Main.hs view
@@ -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
+ test/Spec.hs view
@@ -0,0 +1,2 @@+main :: IO ()+main = putStrLn "Test suite not yet implemented"