Error when trying to dynamically change `tstops` during `solve`

Specific Domains Modelling & Simulations
, , , ,
1

Hello,

I’m modeling a simple integrate and fire neuron that takes spike width into account by suspending the integration (and thus ignoring the input current) for a few milliseconds after a spike has been elicited. Let’s denote this “dead period” with Δ. The input current is modeled as a forcing function, and one way of providing such a forcing function is to pass it as a callback.

The issue resides in the fact that (I suppose) the integrator does not directly adapt its time step to the external forcing function, but only to the system’s dynamics. I noticed that this sometimes causes the neuron to miss variations of the input current (e.g. at the beginning) . To cope with it as it seems to be suggested by the documentation I set tstops at every time the input changes. Running this setup gives an ERROR: Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error. I cannot solve.

Here is an MWE to reproduce the error:

# Load necessary packages
using DifferentialEquations, IfElse

# Model definition

# Parameters are:
# Vth = voltage threshold
# Δ   = dead time
# u_r = reset potential
# spike_times = array where we will save spike times
# input = value of the input current (dynamically modified via a callback)
# gL,EL,C = respectively maximum conductance, "leak" nernst potential and membrane capacitance
function lif(du,u,p,t)
    Vth, Δ, u_r, spike_times, input, gL, EL, C = p 
    du[1] = ( (-1).*gL .* (u[1] .- EL) .+ input ) ./ C
end

# Model parameters and initial conditions
const u0 = -75.0
const p = [-55.0, 0.005, -85.0,Float64[], 0.0,  10.0, -75.0, 5.0]
const tspan = (0.0,40.0)
# Time step and predefined tstops (here simplified to be equally spaced across all `tspan`)
const step = 1e-2
const time_stepping = tspan[1]:step:tspan[2]
const tstops = time_stepping
# Period during which the input current is on
const active_stimulus_period  = (tspan[1]+1, tspan[2]-1)
# Array containing the current values
const ramp_current = [(210/30)*t for t in tspan[1]:step:tspan[2]]


# Utility function. It just returns the index of the closest value to `t` in `series`, in a DifferentialEquations.jl's compatible way
function get_index_of_inferior(series, t)
    i = 0
    for elem in series
        i += IfElse.ifelse(elem <= t, 1, 0)::Int64
    end
    return i
end 

# Callback defining the input current
function input_condition(u,t,integrator)
    (tspan[1]+1) < integrator.t < (tspan[2]-1)
end

function input_affect!(integrator)
    integrator.p[5] = ramp_current[get_index_of_inferior(time_stepping,integrator.t)]
end

const input_callback = DiscreteCallback(input_condition, input_affect!)

# Callback defining the reset mechanism
function condition(u,t,integrator)
    # A spike is triggered when the potential surpasses a threshold V_th
    u[1] > integrator.p[1]
end

function spike!(integrator)
    println("spike at $(integrator.t)")
    # Save the spike time
    push!(integrator.p[4], integrator.t)
    # Calculate the dead time as the minimum between Δ and the time remaining before tspan[2]
    d = min(integrator.p[2], integrator.sol.prob.tspan[2] - integrator.t)
	# Skip integration for d units of time
    integrator.t +=  d
    # Reset potential to u_r
    integrator.u[1] = integrator.p[3]
end

const lif_firing_mechanism = DiscreteCallback(condition,spike!, save_positions = (false,false))

# Define problem
prob = ODEProblem{true}(lif, [u0], tspan, deepcopy(p), callback  = CallbackSet(lif_firing_mechanism, input_callback, ),  tstops = tstops)

But if I solve it:

sol = solve(prob)

I get:

ERROR: Something went wrong. Integrator stepped past tstops but the algorithm was dtchangeable. Please report this error.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33        
  [2] handle_tstop!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, true, Vector{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, 
Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:callback, :tstops), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, 
Float64, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, 
Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\vGDGq\src\integrators\integrator_utils.jl:380
  [3] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, true, Vector{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:callback, :tstops), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, 
typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), 
Nothing}, Vector{Float64}, Vector{Any}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Vector{Any}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{1, false, GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\vGDGq\src\solve.jl:483
  [4] #__solve#502
    @ C:\Users\claud\.julia\packages\OrdinaryDiffEq\vGDGq\src\solve.jl:5 [inlined]
  [5] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:callback, :tstops), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, SciMLBase.StandardODEProblem}, ::Nothing; default_set::Bool, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:second_time, :callback, :tstops), Tuple{Bool, CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), 
typeof(SciMLBase.FINALIZE_DEFAULT)}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}})
    @ DifferentialEquations C:\Users\claud\.julia\packages\DifferentialEquations\4jfQK\src\default_solve.jl:8
  [6] #__solve#48
    @ C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:634 [inlined]
  [7] #solve_call#28
    @ C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:384 [inlined]
  [8] solve_call
    @ C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:366 [inlined]
  [9] solve_up(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:callback, :tstops), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, StepRangeLen{Float64, 
Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, SciMLBase.StandardODEProblem}, 
::Nothing, ::Vector{Float64}, ::Vector{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:415
 [10] solve_up
    @ C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:401 [inlined]
 [11] #solve#29
    @ C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:396 [inlined]
 [12] solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, typeof(lif), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:callback, :tstops), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(condition), typeof(spike!), typeof(SciMLBase.INITIALIZE_DEFAULT), 
typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(input_condition), typeof(input_affect!), 
typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}, SciMLBase.StandardODEProblem})    
    @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\202zR\src\solve.jl:391
 [13] top-level scope
    @ REPL[40]:1

I hypothesized that the dead time caused the integrator to skip some tstops, so I modified the the spike! function to delete all the tstops that occurred between the spike time integrator.t and integrator.t + Δ:

function spike!(integrator)
    println("spike at $(integrator.t)")
    # Save the spike time
    push!(integrator.p[4], integrator.t)
    # Calculate the dead time as the minimum between Δ and the time remaining before tspan[2]
    d = min(integrator.p[2], integrator.sol.prob.tspan[2] - integrator.t)
    # Get the indexes of the tstops to be removed as all the tstops that occurred between the spike time `integrator.t` and `integrator.t + Δ`
    tstops_to_be_removed = findall(tstop -> integrator.t < tstop < (integrator.t+d), integrator.sol.prob.kwargs[:tstops])
    # Remove such tstops (notsure if it is the proper way, tried reiniti but)
    integrator.opts.tstops.valtree = [tstop for (i,tstop) in enumerate(integrator.sol.prob.kwargs[:tstops]) if i ∉ tstops_to_be_removed]
	# Skip integration for d units of time
    integrator.t +=  d
    # Reset potential to u_r
    integrator.u[1] = integrator.p[3]
end

But I get the same error.

Is there any way to use tstops together with the dead time behavior with the current given via callback? Or is there any other way to have the current passed via a callback but being sure that the integrator doesn’t miss any input variation?

Some notable behaviors:

  1. If I comment out the line integrator.t += d but leave uncommented integrator.opts.tstops.valtree = [tstop for (i,tstop) in enumerate(integrator.sol.prob.kwargs[:tstops]) if i ∉ tstops_to_be_removed] it gives the same error .
  2. If I substitute integrator.opts.tstops.valtree = ... with reinit!(integrator, integrator.u, t0 = integrator.t, erase_sol = false , tstops = [tstop for (i,tstop) in enumerate(integrator.sol.prob.kwargs[:tstops]) if i ∉ tstops_to_be_removed], reinit_callbacks = false) it gives the same error again.
  3. Commenting integrator.opts.tstops.valtree = =..., removing tstops but leaving integrator.t += d works (but I no longer have tstops).
  4. Of course if I comment both integrator.opts.tstops.valtree =... and integrator.t += d but leave the tstops it solves correctly (but with no dead time behaviour)
  5. You may confirm that if you remove tstops from the ODEProblem and solve with the same setup as in point 4., you will get a solution different from that of point 4. that misses the first spike by a few seconds.

My manifest

(dynamic_tstops) pkg> st
      Status `D:\GitHub\Academy\Resources\DifferentialEquations.jl\dynamic_tstops\Project.toml`
  [0c46a032] DifferentialEquations v7.1.0
  [615f187c] IfElse v0.1.1
  [91a5bcdd] Plots v1.29.0

Thanks in advance.

2 Likes
2

Isn’t this a duplicate of an issue?

1 Like
3

Hello @ChrisRackauckas ,

We closed that issue and opened this discourse since we thought that topic to be better suited here.

If you believe that the GitHub issue was better, we may reopen that and later post the link to the solution here.