Merge pull request #23 from tkoolen/periodic-callback

Add PeriodicCallback.

README.md

@@ -107,10 +107,22 @@ SavingCallback(save_func, saved_values::SavedValues;
                tdir=1)
 ```
 - `save_func(t, u, integrator)` returns the quantities which shall be saved.
-- `saved_values::SavedValues` contains vectors `t::Vector{tType}`, 
-  `saveval::Vector{savevalType}` of the saved quantities. Here, 
+- `saved_values::SavedValues` contains vectors `t::Vector{tType}`,
+  `saveval::Vector{savevalType}` of the saved quantities. Here,
   `save_func(t, u, integrator)::savevalType`.
 - `saveat` Mimicks `saveat` in `solve` for ODEs.
 - `save_everystep` Mimicks `save_everystep` in `solve` for ODEs.
 - `tdir` should be `sign(tspan[end]-tspan[1])`. It defaults to `1` and should
   be adapted if `tspan[1] > tspan[end]`.
+
+## PeriodicCallback
+
+`PeriodicCallback` can be used when a function should be called periodically in terms of integration time (as opposed to wall time), i.e. at `t = tspan[1]`, `t = tspan[1] + Δt`, `t = tspan[1] + 2Δt`, and so on. This callback can, for example, be used to model a digital controller for an analog system, running at a fixed rate.
+
+A `PeriodicCallback` can be constructed as follows:
+
+```julia
+PeriodicCallback(f, Δt::Number; kwargs...)
+```
+
+where `f` is the function to be called periodically, `Δt` is the period, and `kwargs` are keyword arguments accepted by the `DiscreteCallback` constructor.

src/DiffEqCallbacks.jl

@@ -12,5 +12,6 @@ module DiffEqCallbacks
   include("domain.jl")
   include("stepsizelimiters.jl")
   include("saving.jl")
+  include("periodic.jl")
 
 end # module

src/periodic.jl

@@ -0,0 +1,33 @@
+function PeriodicCallback(f, Δt::Number; initialize = DiffEqBase.INITIALIZE_DEFAULT, kwargs...)
+    # Value of `t` at which `f` should be called next:
+    tnext = Ref(typemax(Δt))
+    condition = (t, u, integrator) -> t == tnext[]
+
+    # Call f, update tnext, and make sure we stop at the new tnext
+    affect! = function (integrator)
+        f(integrator)
+
+        # Schedule next call to `f` using `add_tstops!`, but be careful not to keep integrating forever
+        tnew = tnext[] + Δt
+        tstops = integrator.opts.tstops
+        for i in length(tstops) : -1 : 1 # reverse iterate to encounter large elements earlier
+            if DataStructures.compare(tstops.comparer, tnew, tstops.valtree[i]) # TODO: relying on implementation details
+                tnext[] = tnew
+                add_tstop!(integrator, tnew)
+                break
+            end
+        end
+    end
+
+    # Initialization: first call to `f` should be *before* any time steps have been taken:
+    initialize_periodic = function (c, t, u, integrator)
+        @assert integrator.tdir == sign(Δt)
+        initialize(c, t, u, integrator)
+        tnext[] = t
+        affect!(integrator)
+    end
+
+    DiscreteCallback(condition, affect!; initialize = initialize_periodic, kwargs...)
+end
+
+export PeriodicCallback

test/periodic_tests.jl

@@ -0,0 +1,57 @@
+using Base.Test, OrdinaryDiffEq, DiffEqCallbacks
+
+tmin = 0.1
+tmax = 5.2
+
+for tmax_problem in [tmax; Inf]
+    # Test with both a finite and an infinite tspan for the ODEProblem.
+    #
+    # Having support for infinite tspans is one of the main reasons for implementing PeriodicCallback
+    # using add_tstop!, instead of just passing in a linspace as the tstops solve argument.
+    # (the other being that the length of the internal tstops collection would otherwise become
+    # linear in the length of the integration time interval.
+    #
+    # Testing a finite tspan is necessary because a naive implementation could add tstops after
+    # tmax and thus integrate for too long (or even indefinitely).
+
+    # Dynamics: two independent single integrators:
+    du = [0; 0]
+    u0 = [0.; 0.]
+    dynamics = (t, u) -> eltype(u).(du)
+    prob = ODEProblem(dynamics, u0, (tmin, tmax_problem))
+
+    # Callbacks periodically increase the input to the integrators:
+    Δt1 = 0.5
+    increase_du_1 = integrator -> du[1] += 1
+    periodic1_initialized = Ref(false)
+    initialize1 = (c, t, u, integrator) -> periodic1_initialized[] = true
+    periodic1 = PeriodicCallback(increase_du_1, Δt1; initialize = initialize1)
+
+    Δt2 = 1.
+    increase_du_2 = integrator -> du[2] += 1
+    periodic2 = PeriodicCallback(increase_du_2, Δt2)
+
+    # Terminate at tmax (regardless of whether the tspan of the ODE problem is infinite).
+    terminator = DiscreteCallback((t, u, integrator) -> t == tmax, terminate!)
+
+    # Solve.
+    sol = solve(prob, Tsit5(); callback = CallbackSet(terminator, periodic1, periodic2), tstops = [tmax])
+
+    # Ensure that initialize1 has been called
+    @test periodic1_initialized[]
+
+    # Make sure we haven't integrated past tmax:
+    @test sol.t[end] == tmax
+
+    # Make sure that the components of du have been incremented the appropriate number of times.
+    Δts = [Δt1, Δt2]
+    expected_num_calls = map(Δts) do Δt
+        floor(Int, (tmax - tmin) / Δt) + 1
+    end
+    @test du == expected_num_calls
+
+    # Make sure that the final state matches manual integration of the piecewise linear function
+    foreach(Δts, sol.u[end], du) do Δt, u_i, du_i
+        @test u_i ≈ Δt * sum(1 : du_i - 1) + rem(tmax - tmin, Δt) * du_i atol = 1e-5
+    end
+end

test/runtests.jl

@@ -8,4 +8,5 @@ tic()
 @time @testset "Manifold tests" begin include("manifold_tests.jl") end
 @time @testset "StepsizeLimiter tests" begin include("stepsizelimiter_tests.jl") end
 @time @testset "Saving tests" begin include("saving_tests.jl") end
+@time @testset "Periodic tests" begin include("periodic_tests.jl") end
 toc()