callbacks in differential equations are non mutable

[awage]

Describe the bug When a callback is triggered the integrator is not mutable. The solution is to use set_u! or set_state! but integrator.u should be mutable. I don't know where the problem comes from. Giving a value to the state triggers:

ERROR: LoadError: setindex!(::SVector{2, Float64}, value, ::Int) is not defined.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] setindex!(a::SVector{2, Float64}, value::Float64, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/aMe0r/src/indexing.jl:3

Minimal Working Example

using OrdinaryDiffEq
using DynamicalSystems
function forced_pendulum(u, p, t)
    @inbounds begin
    d = p[1]; F = p[2]; omega = p[3]
    du1 = u[2]
    du2 = -d*u[2] - sin(u[1])+ F*cos(omega*t)
    return SVector{2}(du1, du2)
    end
end

# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
    #@show ismutable(integrator.u)  # returns false
    x = integrator.u[1]
    if x < 0
        set_u!(integrator, [x + 2*π, integrator.u[2]])
    else
        set_u!(integrator, [x - 2*π, integrator.u[2]])
    end
end

condition(u,t,integrator) = (integrator.u[1] < -π  || integrator.u[1] > π)
cb = DiscreteCallback(condition,affect!)

 F = 1.66;  ω = 1.;  d=0.2;
df = ContinuousDynamicalSystem(forced_pendulum,rand(2), [d, F, ω])
xg = range(-pi,pi,length=150)
yg = range(-2.,4.,length=150)
# compute basin
bsn,att = basins_of_attraction((xg, yg), df;  T=2*pi/ω, diffeq=(alg=Tsit5(), callback=cb))

[Datseris]

The state here is a static vector so of course it is not mutable. What you should do in affect! is:

set_state!(integrator, new_state_that_is_also_a_static_vector)
u_modified!(integrator, true)

this should work. Can you test it and if it works close the issue?

[awage]

I completely forgot this. I used a workaround but this method is cleaner. Thanks!

using OrdinaryDiffEq
using DynamicalSystems
using Plots

# Equations of motion:
function forced_pendulum(u, p, t)
    @inbounds begin
    d = p[1]; F = p[2]; omega = p[3]
    du1 = u[2]
    du2 = -d*u[2] - sin(u[1])+ F*cos(omega*t)
    return SVector{2}(du1, du2)
    end
end

# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
    uu = integrator.u
    if integrator.u[1] < 0
        set_state!(integrator, SVector(uu[1] + 2π, uu[2]))
        u_modified!(integrator, true)
    else
        set_state!(integrator, SVector(uu[1] - 2π, uu[2]))
        u_modified!(integrator, true)
    end
end

condition(u,t,integrator) = (integrator.u[1] < -π  || integrator.u[1] > π)
cb = DiscreteCallback(condition,affect!)

 F, ω, d = 1.66, 1., 0.2
 p = [d, F, ω]
df = ContinuousDynamicalSystem(forced_pendulum,rand(2), p)

# range for forced pend
xg = range(-pi,pi,length = 300)
yg = range(-2.,4.,length = 300)

default_diffeq = (reltol = 1e-9, abstol = 1e-9,  alg = Vern9(), callback = cb)
bsn, att = basins_of_attraction((xg, yg), df; T = 2*pi/ω, diffeq = default_diffeq)