aivika 5.4 → 5.5

raw patch · 4 files changed

+177/−4 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• CHANGELOG.md +7/−0
• Simulation/Aivika/Internal/Specs.hs +30/−1
• Simulation/Aivika/SystemDynamics.hs +139/−2
• aivika.cabal +1/−1

@@ -1,4 +1,11 @@ +Version 5.5+-----++* Added the 4th order Runge-Kutta 3/8-method.++* Added experimental functions delayByDT and delayIByDT.+ Version 5.4 ----- 

@@ -29,7 +29,8 @@         integStopPoint,         simulationStopPoint,         timeGrid,-        pointAt) where+        pointAt,+        delayPoint) where  import Data.IORef @@ -49,6 +50,7 @@ data Method = Euler          -- ^ Euler's method             | RungeKutta2    -- ^ the 2nd order Runge-Kutta method             | RungeKutta4    -- ^ the 4th order Runge-Kutta method+            | RungeKutta4b   -- ^ the 4th order Runge-Kutta 3/8-method             deriving (Eq, Ord, Show)  -- | It indentifies the simulation run.@@ -134,6 +136,7 @@     Euler -> [0]     RungeKutta2 -> [0, 1]     RungeKutta4 -> [0, 1, 2, 3]+    RungeKutta4b -> [0, 1, 2, 3]  -- | Returns the first and last integration phases. integPhaseBnds :: Specs -> (Int, Int)@@ -142,6 +145,7 @@     Euler -> (0, 0)     RungeKutta2 -> (0, 1)     RungeKutta4 -> (0, 3)+    RungeKutta4b -> (0, 3)  -- | Returns the first integration phase, i.e. zero. integPhaseLoBnd :: Specs -> Int@@ -154,6 +158,7 @@     Euler -> 0     RungeKutta2 -> 1     RungeKutta4 -> 3+    RungeKutta4b -> 3  -- | Returns a simulation time for the integration point specified by  -- the specs, iteration and phase.@@ -171,6 +176,10 @@             delta RungeKutta4 1 = spcDT sc / 2             delta RungeKutta4 2 = spcDT sc / 2             delta RungeKutta4 3 = spcDT sc+            delta RungeKutta4b 0 = 0+            delta RungeKutta4b 1 = spcDT sc / 3+            delta RungeKutta4b 2 = 2 * spcDT sc / 3+            delta RungeKutta4b 3 = spcDT sc  -- | Return the integration time values. integTimes :: Specs -> [Double]@@ -256,3 +265,23 @@           | i == n'   = (i, t2)           | otherwise = (i, t0 + (fromIntegral i) * dt)   in map f [0 .. n']++-- | Delay the point by the specified positive number of iterations.+delayPoint :: Point -> Int -> Point+delayPoint p dn+  | dn <= 0   = error "Expected the positive number of iterations: delayPoint"+  | otherwise =+    let sc = pointSpecs p+        n  = pointIteration p+        ph = pointPhase p+    in if ph < 0+       then let t' = pointTime p - fromIntegral dn * spcDT sc+                n' = fromIntegral $ floor $ (t' - spcStartTime sc) / spcDT sc+            in p { pointTime = t',+                   pointIteration = n',+                   pointPhase = -1 }+       else let n' = n - dn+                t' = basicTime sc n' ph+            in p { pointTime = t',+                   pointIteration = n',+                   pointPhase = ph }

@@ -49,6 +49,8 @@         -- * Discrete Functions         delay,         delayI,+        delayByDT,+        delayIByDT,         step,         pulse,         pulseP,@@ -237,6 +239,79 @@     _ ->        error "Incorrect phase: integRK4" +integRK4b :: Dynamics Double+             -> Dynamics Double+             -> Dynamics Double+             -> Point -> IO Double+integRK4b (Dynamics f) (Dynamics i) (Dynamics y) p =+  case pointPhase p of+    0 -> case pointIteration p of+      0 -> +        i p+      n -> do+        let sc = pointSpecs p+            ty = basicTime sc (n - 1) 0+            t1 = ty+            t2 = basicTime sc (n - 1) 1+            t3 = basicTime sc (n - 1) 2+            t4 = basicTime sc (n - 1) 3+            py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }+            p1 = py+            p2 = p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }+            p3 = p { pointTime = t3, pointIteration = n - 1, pointPhase = 2 }+            p4 = p { pointTime = t4, pointIteration = n - 1, pointPhase = 3 }+        vy <- y py+        k1 <- f p1+        k2 <- f p2+        k3 <- f p3+        k4 <- f p4+        let !v = vy + spcDT sc / 8.0 * (k1 + 3.0 * (k2 + k3) + k4)+        return v+    1 -> do+      let sc = pointSpecs p+          n  = pointIteration p+          ty = basicTime sc n 0+          t1 = ty+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+          p1 = py+      vy <- y py+      k1 <- f p1+      let !v = vy + spcDT sc / 3.0 * k1+      return v+    2 -> do+      let sc = pointSpecs p+          n  = pointIteration p+          ty = basicTime sc n 0+          t1 = ty+          t2 = basicTime sc n 1+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+          p1 = py+          p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }+      vy <- y py+      k1 <- f p1+      k2 <- f p2+      let !v = vy + spcDT sc * (- k1 / 3.0 + k2)+      return v+    3 -> do+      let sc = pointSpecs p+          n  = pointIteration p+          ty = basicTime sc n 0+          t1 = ty+          t2 = basicTime sc n 1+          t3 = basicTime sc n 2+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }+          p1 = py+          p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }+          p3 = p { pointTime = t3, pointIteration = n, pointPhase = 2 }+      vy <- y py+      k1 <- f p1+      k2 <- f p2+      k3 <- f p3+      let !v = vy + spcDT sc * (k1 - k2 + k3)+      return v+    _ -> +      error "Incorrect phase: integRK4b"+ -- | Return an integral with the specified derivative and initial value. -- -- To create a loopback, you should use the recursive do-notation.@@ -263,6 +338,7 @@           Euler -> return $ Dynamics $ integEuler diff i y           RungeKutta2 -> return $ Dynamics $ integRK2 diff i y           RungeKutta4 -> return $ Dynamics $ integRK4 diff i y+          RungeKutta4b -> return $ Dynamics $ integRK4b diff i y       return y  integEulerEither :: Dynamics (Either Double Double)@@ -632,10 +708,71 @@                                   pointIteration = n',                                   pointPhase = -1 }             | n' > n    = error $-                          "Cannot return the future data: delay. " +++                          "Cannot return the future data: delayI. " ++                           "The lag time cannot be negative."             | otherwise = error $-                          "Cannot return the current data: delay. " +++                          "Cannot return the current data: delayI. " +++                          "The lag time is too small."+      y++-- | Return the delayed value by the specified positive number of+-- integration time steps used for calculating the lag time.+delayByDT :: Dynamics a+             -- ^ the value to delay+             -> Dynamics Int+             -- ^ the delay as a multiplication of the corresponding number+             -- and the integration time step+             -> Dynamics a+             -- ^ the delayed value+delayByDT (Dynamics x) (Dynamics d) = discreteDynamics $ Dynamics r +  where+    r p = do +      let sc = pointSpecs p+          n  = pointIteration p+      a <- d p+      let p' = delayPoint p a+          n' = pointIteration p'+          y | n' < 0    = x $ p { pointTime = spcStartTime sc,+                                  pointIteration = 0, +                                  pointPhase = 0 }+            | n' < n    = x p'+            | n' > n    = error $+                          "Cannot return the future data: delayByDT. " +++                          "The lag time cannot be negative."+            | otherwise = error $+                          "Cannot return the current data: delayByDT. " +++                          "The lag time is too small."+      y+      +-- | Return the delayed value by the specified initial value and+-- a positive number of integration time steps used for calculating+-- the lag time. It allows creating a loop back.+delayIByDT :: Dynamics a+              -- ^ the value to delay+              -> Dynamics Int+              -- ^ the delay as a multiplication of the corresponding number+              -- and the integration time step+              -> Dynamics a+              -- ^ the initial value+              -> Simulation (Dynamics a)+              -- ^ the delayed value+delayIByDT (Dynamics x) (Dynamics d) (Dynamics i) = M.memoDynamics $ Dynamics r +  where+    r p = do +      let sc = pointSpecs p+          n  = pointIteration p+      a <- d p+      let p' = delayPoint p a+          n' = pointIteration p'+          y | n' < 0    = i $ p { pointTime = spcStartTime sc,+                                  pointIteration = 0, +                                  pointPhase = 0 }+            | n' < n    = x p'+            | n' > n    = error $+                          "Cannot return the future data: delayIByDT. " +++                          "The lag time cannot be negative."+            | otherwise = error $+                          "Cannot return the current data: delayIByDT. " ++                           "The lag time is too small."       y 

@@ -1,5 +1,5 @@ name:            aivika-version:         5.4+version:         5.5 synopsis:        A multi-method simulation library description:     Aivika is a discrete event simulation (DES) framework with support of activity-oriented,