Error which only appears when inside sciML_train call environment: ERROR: ArgumentError: unable to check bounds for indices of type Nothing

Hi all,

I get this type error which appears after I try to index with the output of indexin() which is of type: Array{Union{Nothing, Int64},1}. I have tried to fix this via all of the following methods:

    match_idx = convert(Array{Int64}, indexin(0:24:144, sol.t))

    match_idx = replace(indexin(0:24:144, sol.t), nothing => Int64)

    match_idx = indexin(0:24:144, sol.t)
    match_idxI = Array{Int64}(undef, length(match_data))
    for i=1:length(match_data); match_idxI[i] = match_idx[i]; end

none of which work; each produces an error about converting from nothing to Int64. However, I noticed bizarrely though, any of these options work if I run the function where this issue exists line by line in the REPL?? Its only when I call the function (which is my loss function) from sciML_train that I get the error. Moreover, I can call the function on its own without any trouble and it actually works in the first case without any fixes.

Here is a relevant excerpt of the code:

function loss(p_curr)
    sol = DifferentialEquations.solve(prob, Tsit5(), p=p_curr, saveat=tsteps, callback=cb)
    
    match_idx = indexin(0:24:144, sol.t)

    # alternatives:
    # match_idx = convert(Array{Int64}, indexin(0:24:144, sol.t))
    # match_idx = replace(indexin(0:24:144, sol.t), nothing => Int64)
    # match_idxN = indexin(0:24:144, sol.t)
    # match_idx = Array{Int64}(undef, length(match_data))
    # for i=1:length(match_data); match_idx[i] = match_idxN[i]; end

    pred = Array(sol)[:, match_idx]
    loss = sum(abs2, pred - syndata)
    return loss, sol, match_data
end

result_ode = DiffEqFlux.sciml_train(loss, p_init,
                                    ADAM(0.1),
                                    cb = callback,
                                    maxiters = 300)

Note, I do this because the callback function in my ODE problem causes the solution to include time points beyond the ‘saveat’ interval that I don’t care about. Another work around would be totally fine.

Full error:

WARNING: both Flux and Plots export "scatter"; uses of it in module Main must be qualified
ERROR: ArgumentError: unable to check bounds for indices of type Nothing
Stacktrace:
 [1] checkindex(::Type{Bool}, ::Base.OneTo{Int64}, ::Nothing) at ./abstractarray.jl:561
 [2] checkindex at ./abstractarray.jl:576 [inlined]
 [3] checkbounds_indices at ./abstractarray.jl:532 [inlined] (repeats 2 times)
 [4] checkbounds at ./abstractarray.jl:485 [inlined]
 [5] checkbounds at ./abstractarray.jl:506 [inlined]
 [6] _getindex at ./multidimensional.jl:742 [inlined]
 [7] getindex at ./abstractarray.jl:1060 [inlined]
 [8] adjoint at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/lib/array.jl:31 [inlined]
 [9] _pullback(::Zygote.Context, ::typeof(getindex), ::Array{Float64,2}, ::Function, ::Array{Union{Nothing, Int64},1}) at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57
 [10] loss at ./none:11 [inlined]
 [11] _pullback(::Zygote.Context, ::typeof(loss), ::Array{Float64,1}) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface2.jl:0
 [12] #69 at /Users/willsharpless/.julia/packages/DiffEqFlux/pJFx8/src/train.jl:3 [inlined]
 [13] _pullback(::Zygote.Context, ::DiffEqFlux.var"#69#70"{typeof(loss)}, ::Array{Float64,1}, ::SciMLBase.NullParameters) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface2.jl:0
 [14] adjoint at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/lib/lib.jl:191 [inlined]
 [15] _pullback at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined]
 [16] OptimizationFunction at /Users/willsharpless/.julia/packages/SciMLBase/fypD8/src/problems/basic_problems.jl:107 [inlined]
 [17] _pullback(::Zygote.Context, ::OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::SciMLBase.NullParameters) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface2.jl:0
 [18] adjoint at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/lib/lib.jl:191 [inlined]
 [19] adjoint(::Zygote.Context, ::typeof(Core._apply_iterate), ::typeof(iterate), ::OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Tuple{Array{Float64,1},SciMLBase.NullParameters}) at ./none:0
 [20] _pullback at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined]
 [21] OptimizationFunction at /Users/willsharpless/.julia/packages/SciMLBase/fypD8/src/problems/basic_problems.jl:107 [inlined]
 [22] _pullback(::Zygote.Context, ::OptimizationFunction{false,GalacticOptim.AutoZygote,OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GalacticOptim.var"#146#156"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#149#159"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#154#164",Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::SciMLBase.NullParameters) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface2.jl:0
 [23] adjoint at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/lib/lib.jl:191 [inlined]
 [24] _pullback at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined]
 [25] #8 at /Users/willsharpless/.julia/packages/GalacticOptim/JnLwV/src/solve.jl:94 [inlined]
 [26] _pullback(::Zygote.Context, ::GalacticOptim.var"#8#13"{OptimizationProblem{false,OptimizationFunction{false,GalacticOptim.AutoZygote,OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GalacticOptim.var"#146#156"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#149#159"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#154#164",Nothing,Nothing,Nothing},Array{Float64,1},SciMLBase.NullParameters,Nothing,Nothing,Nothing,Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:cb, :maxiters),Tuple{var"#45#46",Int64}}}},Array{Float64,1},GalacticOptim.NullData}) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface2.jl:0
 [27] pullback(::Function, ::Zygote.Params) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface.jl:247
 [28] gradient(::Function, ::Zygote.Params) at /Users/willsharpless/.julia/packages/Zygote/pM10l/src/compiler/interface.jl:58
 [29] __solve(::OptimizationProblem{false,OptimizationFunction{false,GalacticOptim.AutoZygote,OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GalacticOptim.var"#146#156"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#149#159"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#154#164",Nothing,Nothing,Nothing},Array{Float64,1},SciMLBase.NullParameters,Nothing,Nothing,Nothing,Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:cb, :maxiters),Tuple{var"#45#46",Int64}}}}, ::ADAM, ::Base.Iterators.Cycle{Tuple{GalacticOptim.NullData}}; maxiters::Int64, cb::Function, progress::Bool, save_best::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/willsharpless/.julia/packages/GalacticOptim/JnLwV/src/solve.jl:93
 [30] #solve#450 at /Users/willsharpless/.julia/packages/SciMLBase/fypD8/src/solve.jl:3 [inlined]
 [31] sciml_train(::typeof(loss), ::Array{Float64,1}, ::ADAM, ::GalacticOptim.AutoZygote; lower_bounds::Nothing, upper_bounds::Nothing, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:cb, :maxiters),Tuple{var"#45#46",Int64}}}) at /Users/willsharpless/.julia/packages/DiffEqFlux/pJFx8/src/train.jl:6
 [32] top-level scope at none:1

This was an issue introduced in:

https://github.com/FluxML/NNlib.jl/pull/309

which is now fixed in NNLib.jl v0.7.18. So if you update your packages you should be fine.

Packages updated no changes in error. Solve is still occasionally returning points not on the saveat interval and attempting to add this fix with the previous method results in the original ‘unable to check bounds for indices of type Nothing’ error.

Share an example?

Yeah sure (similar to above):

function loss(p_curr)
    sol = DifferentialEquations.solve(prob, Tsit5(), p=p_curr, saveat=24, callback=cb)
    match_idx = indexin(0:24:144, sol.t)
    pred = Array(sol)[:, match_idx]
    loss = sum(abs2, pred - syndata)
    return loss, sol
end

loss(p_init)
result_ode = DiffEqFlux.sciml_train(loss, p_init, ADAM(0.1), maxiters=1000)

stack trace:

ERROR: ArgumentError: unable to check bounds for indices of type Nothing
Stacktrace:
 [1] checkindex(::Type{Bool}, ::Base.OneTo{Int64}, ::Nothing) at ./abstractarray.jl:561
 [2] checkindex at ./abstractarray.jl:576 [inlined]
 [3] checkbounds_indices at ./abstractarray.jl:532 [inlined] (repeats 2 times)
 [4] checkbounds at ./abstractarray.jl:485 [inlined]
 [5] checkbounds at ./abstractarray.jl:506 [inlined]
 [6] _getindex at ./multidimensional.jl:742 [inlined]
 [7] getindex at ./abstractarray.jl:1060 [inlined]
 [8] adjoint at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/lib/array.jl:31 [inlined]
 [9] _pullback(::Zygote.Context, ::typeof(getindex), ::Array{Float64,2}, ::Function, ::Array{Union{Nothing, Int64},1}) at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57
 [10] loss at /Users/willsharpless/Documents/Julia/example.jl:57 [inlined]
 [11] _pullback(::Zygote.Context, ::typeof(loss), ::Array{Float64,1}) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface2.jl:0
 [12] #69 at /Users/willsharpless/.julia/packages/DiffEqFlux/alPQ3/src/train.jl:3 [inlined]
 [13] _pullback(::Zygote.Context, ::DiffEqFlux.var"#69#70"{typeof(loss)}, ::Array{Float64,1}, ::SciMLBase.NullParameters) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface2.jl:0
 [14] adjoint at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/lib/lib.jl:191 [inlined]
 [15] _pullback at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined]
 [16] OptimizationFunction at /Users/willsharpless/.julia/packages/SciMLBase/yh1iq/src/problems/basic_problems.jl:107 [inlined]
 [17] _pullback(::Zygote.Context, ::OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::SciMLBase.NullParameters) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface2.jl:0
 [18] adjoint at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/lib/lib.jl:191 [inlined]
 [19] adjoint(::Zygote.Context, ::typeof(Core._apply_iterate), ::typeof(iterate), ::OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Tuple{Array{Float64,1},SciMLBase.NullParameters}) at ./none:0
 [20] _pullback at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined]
 [21] OptimizationFunction at /Users/willsharpless/.julia/packages/SciMLBase/yh1iq/src/problems/basic_problems.jl:107 [inlined]
 [22] _pullback(::Zygote.Context, ::OptimizationFunction{false,GalacticOptim.AutoZygote,OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GalacticOptim.var"#146#156"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#149#159"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#154#164",Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::SciMLBase.NullParameters) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface2.jl:0
 [23] adjoint at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/lib/lib.jl:191 [inlined]
 [24] _pullback at /Users/willsharpless/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined]
 [25] #8 at /Users/willsharpless/.julia/packages/GalacticOptim/JnLwV/src/solve.jl:94 [inlined]
 [26] _pullback(::Zygote.Context, ::GalacticOptim.var"#8#13"{OptimizationProblem{false,OptimizationFunction{false,GalacticOptim.AutoZygote,OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GalacticOptim.var"#146#156"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#149#159"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#154#164",Nothing,Nothing,Nothing},Array{Float64,1},SciMLBase.NullParameters,Nothing,Nothing,Nothing,Base.Iterators.Pairs{Symbol,Int64,Tuple{Symbol},NamedTuple{(:maxiters,),Tuple{Int64}}}},Array{Float64,1},GalacticOptim.NullData}) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface2.jl:0
 [27] pullback(::Function, ::Zygote.Params) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface.jl:247
 [28] gradient(::Function, ::Zygote.Params) at /Users/willsharpless/.julia/packages/Zygote/RxTZu/src/compiler/interface.jl:58
 [29] __solve(::OptimizationProblem{false,OptimizationFunction{false,GalacticOptim.AutoZygote,OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GalacticOptim.var"#146#156"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#149#159"{GalacticOptim.var"#145#155"{OptimizationFunction{true,GalacticOptim.AutoZygote,DiffEqFlux.var"#69#70"{typeof(loss)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing}},GalacticOptim.var"#154#164",Nothing,Nothing,Nothing},Array{Float64,1},SciMLBase.NullParameters,Nothing,Nothing,Nothing,Base.Iterators.Pairs{Symbol,Int64,Tuple{Symbol},NamedTuple{(:maxiters,),Tuple{Int64}}}}, ::ADAM, ::Base.Iterators.Cycle{Tuple{GalacticOptim.NullData}}; maxiters::Int64, cb::Function, progress::Bool, save_best::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/willsharpless/.julia/packages/GalacticOptim/JnLwV/src/solve.jl:93
 [30] #solve#468 at /Users/willsharpless/.julia/packages/SciMLBase/yh1iq/src/solve.jl:3 [inlined]
 [31] sciml_train(::typeof(loss), ::Array{Float64,1}, ::ADAM, ::GalacticOptim.AutoZygote; lower_bounds::Nothing, upper_bounds::Nothing, kwargs::Base.Iterators.Pairs{Symbol,Int64,Tuple{Symbol},NamedTuple{(:maxiters,),Tuple{Int64}}}) at /Users/willsharpless/.julia/packages/DiffEqFlux/alPQ3/src/train.jl:6
 [32] top-level scope at none:1

That example doesn’t run. Many things aren’t defined?

Apologies. Note it matters that the randomization of p_init makes the dynamics hit 0 and act on the callback.

##
using Plots, Random
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity

## Define GLV and generate synthetic data from defined parameters

tsteps = 24
n=3

function gLV!(dx,x,p,t)

    # unpack params
    n = length(x)
    mu = p[1:n]
    A = Array{eltype(x),2}(undef,3,3)
    for i=1:n; A[:,i] = p[n*i + 1: n*(i+1)]; end

    dx .= x.*(mu + A*x)

    if any(x .> 100.0) # no cc catch
        dx .= 0
    end
end

condition(u,t,integrator) = any(u .< 1e-8) #|| any(u .> 10)
function affect!(integrator)
    integrator.u[integrator.u .< 1e-8] .= 0
end
cb = ContinuousCallback(condition,affect!)

mu_true = [0.4; 0.8; 0.2].*0.1
A_true = [-1.5 0.5 -0.2; -0.7 -1.3 0.1; 0.8 -0.2 -0.9]
p_true = [mu_true;  A_true[:]]

x0 = [0.1; 0.1; 0.1].*0.1
tspan = (0.0, 144.0)

prob = ODEProblem(gLV!, x0, tspan)
sol = DifferentialEquations.solve(prob, Tsit5(), p=p_true, saveat=tsteps, callback=cb)
syndata = Array(sol) + 0.002*randn(size(Array(sol)))
plot(sol, alpha=0.3)
Plots.scatter!(sol.t, syndata')
title!("Definite")

## Random Initialize

mu_init = 0.1*rand(n)
A_init = -1*rand(n,n)
p_init = [mu_init;  A_init[:]]
print(-1*inv(A_init)*mu_init) 

What am I looking for here? Seems like that runs?

You combined the two blocks (both model definition and optimization) and ran this?

When I run all of the following code (both model definition and optimization) I receive the same error as originally listed for trying to use the correct time points with match_idx.

##
using Plots, Random
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity

## Define GLV and generate synthetic data from defined parameters

tsteps = 24
n=3

function gLV!(dx,x,p,t)

    # unpack params
    n = length(x)
    mu = p[1:n]
    A = Array{eltype(x),2}(undef,3,3)
    for i=1:n; A[:,i] = p[n*i + 1: n*(i+1)]; end

    dx .= x.*(mu + A*x)

    if any(x .> 100.0) # no cc catch
        dx .= 0
    end
end

condition(u,t,integrator) = any(u .< 1e-8) #|| any(u .> 10)
function affect!(integrator)
    integrator.u[integrator.u .< 1e-8] .= 0
end
cb = ContinuousCallback(condition,affect!)

mu_true = [0.4; 0.8; 0.2].*0.1
A_true = [-1.5 0.5 -0.2; -0.7 -1.3 0.1; 0.8 -0.2 -0.9]
p_true = [mu_true;  A_true[:]]

x0 = [0.1; 0.1; 0.1].*0.1
tspan = (0.0, 144.0)

prob = ODEProblem(gLV!, x0, tspan)
sol = DifferentialEquations.solve(prob, Tsit5(), p=p_true, saveat=tsteps, callback=cb)
syndata = Array(sol) + 0.002*randn(size(Array(sol)))
plot(sol, alpha=0.3)
Plots.scatter!(sol.t, syndata')
title!("Definite")

## Random Initialize

mu_init = 0.1*rand(n)
A_init = -1*rand(n,n)
p_init = [mu_init;  A_init[:]]
print(-1*inv(A_init)*mu_init) 

function loss(p_curr)
    sol = DifferentialEquations.solve(prob, Tsit5(), p=p_curr, saveat=24, callback=cb)
    match_idx = indexin(0:24:144, sol.t)
    pred = Array(sol)[:, match_idx]
    loss = sum(abs2, pred - syndata)
    return loss, sol
end

loss(p_init)
result_ode = DiffEqFlux.sciml_train(loss, p_init, ADAM(0.1), maxiters=1000)

The indexin(0:24:144, sol.t) is the thing that doesn’t make all that much sense. If you have something that isn’t a hit, then it gives a nothing, and then it errors when trying to index nothing. I’m not sure what you’re trying to do there.

For whatever reason, if the callback function is evoked, ie. the glv system hits a zero bound, it returns the time point that it was evoked with the saveat points in the solution. Obviously, I don’t care about these points for the fitting (and I imagine this is a temporary bug) so I use indexin() to pull out the points in the solution which correspond to the true saveat points.

It’s not a bug, it’s how your callback options are set. Do cb = ContinuousCallback(condition,affect!,save_positions = (false,false))

1 Like

Thanks Chris!