Firing rate calculation sensitive to solver tolerances due to multiple peak detection

[gabrevaya]

using Neuroblox
using StochasticDiffEq
using Random

Random.seed!(1)
@named stn = STN_Adam(N_exci = 40)
sys = structural_simplify(get_system(stn))
tspan = (0.0, 1000.0)
prob = SDEProblem(sys, [], tspan, [])

sol = solve(prob, RKMil()); 
fr = firing_rate(stn, sol, threshold=0, scheduler=:dynamic)
# 60.7

sol = solve(prob, RKMil(), abstol = 1e-3, reltol = 1e-6, maxiters=1e10); 
fr = firing_rate(stn, sol, threshold=0, scheduler=:dynamic)
# 243.9

sol = solve(prob, RKMil(), abstol = 1e-4, reltol = 1e-8, maxiters=1e10);
fr = firing_rate(stn, sol, threshold=0, scheduler=:dynamic)
# 1301.275

The reason for this issue is that when reducing the tolerances, if we don't specify saveat, the timeseries tends to have more resolution, and multiple peaks are present in each spike. While discussing this issue with @harisorgn, he suggested that we might be able to fix it using peakwidths! from Peaks.jl. However, we will still have to sort out some things before implementing it since the width in that function is determined by indexes, not time, and when using any adaptive solver, time is irregularly sampled.

In any case, this is just to report the issue and be aware that this problem might happen. One workaround until this is fixed is setting the saveat kwarg:

sol = solve(prob, RKMil(), abstol = 1e-4, reltol = 1e-8, maxiters=1e10, saveat=0.1);
fr = firing_rate(stn, sol, threshold=0, scheduler=:dynamic)
# 61.75