Delay implementation introduces artifacts in solution

[maedoc]

I had a hard time writing a descriptive title, apologies. I am testing the solvers on a non-stiff network of oscillators, checking that the coherence of the oscillators coverges to an expected analytic result and discovered that the solution appears to be wrong for non-EM solvers (why is EM unaffected?) when I set the delay to zero.

Definition of the additive noise oscillator model and history function

The functions defining the problem are as follows

struct kp
    G::Float64
    σ::Float64
    n::UInt64
    τ::Float64
end

# Kuramoto order parameter
kopz(sol) = [angle(mean(exp.(1.0im .* u_t))) for u_t in sol.u];
kop(sol) = [abs(mean(exp.(1.0im .* u_t))) for u_t in sol.u];

function kf(du,u,h,p,t)
    G = p.G/p.n
    du .= 5*2*π/1000
    for j=1:p.n
        for i=1:p.n
            if i!=j
                uⱼ = h(p, t - p.τ; idxs=j)[1]
                du[i] += G * sin(uⱼ - u[i])
            end
        end
    end
end

function kf(du,u,p,t)
    G = p.G/p.n
    du .= 5*2*π/1000
    for j=1:p.n
        for i=1:p.n
            if i!=j
                uⱼ = u[j]
                du[i] += G * sin(uⱼ - u[i])
            end
        end
    end
end

function kg(du,u,h,p,t)
    du .= sqrt(2*p.σ)
end

function kg(du,u,p,t)
    du .= sqrt(2*p.σ)
end

function kh(p,t;idxs=nothing)
    x = (1:p.n)/p.n*2*pi .+ t .* 2 .* pi
    collect(idxs==nothing ? x : x[idxs])
end

Here I use the SDDE solver but set the delay to zero

g = 0.1
p = kp(g, g/2.5, 100, 0.0); # the last zero here sets the delay
ic = kh(p, 0.0);
T = 1000
prob = SDDEProblem(kf, kg, ic, kh, (0.,T), p);
eprob = EnsembleProblem(prob)
@time sol = solve(eprob, SOSRA(), dt=1, adaptive=false, trajectories=10,saveat=0:1:T);
plot(hcat([kop(s) for s in sol]...),linealpha=0.2)

which shows the coherence of the network over time for each realization, but the result is wrong.

Suspecting the delay implementation (but not delays themselves, since they are set to zero),

prob = SDEProblem(kf, kg, ic, (0.,T), p);
eprob = EnsembleProblem(prob)
@time sol = solve(eprob, SOSRA(), dt=1, adaptive=false, trajectories=10,saveat=0:1:T);
plot(hcat([kop(s) for s in sol]...),linealpha=0.2)

which is correct. But, EM doesn't seem to mind,

prob = SDDEProblem(kf, kg, ic, kh, (0.,T), p);
eprob = EnsembleProblem(prob)
@time sol = solve(eprob, EM(), dt=1, trajectories=10,saveat=0:1:T);
plot(hcat([kop(s) for s in sol]...),linealpha=0.2)

though EM should use a shorter time step to be as accurate as SOSRA.

For reference, I include a SOSRA run with higher tolerance,

prob = SDEProblem(kf, kg, ic, (0.,T), p);
eprob = EnsembleProblem(prob)
@time sol = solve(eprob, SOSRA(), tol=1e-9, trajectories=10,saveat=0:1:T);
plot(hcat([kop(s) for s in sol]...),linealpha=0.2)

I hope this is detailed enough to explain the problem.

[ChrisRackauckas]

Thanks. @HTSykora might be interested to comment here. Something does look off. Did you run it without multithreading? It could be a thread-safety issue.

[maedoc]

Did you run it without multithreading? It could be a thread-safety issue.

This is single threaded

[ChrisRackauckas]

eprob = EnsembleProblem(prob)
@time sol = solve(eprob, SOSRA(), dt=1, adaptive=false, trajectories=10,saveat=0:1:T);

Looks like you're doing an ensemble solve with the defaults which would be EnsembleThreads.

[maedoc]

True, but I only had one thread. Anyway, here's the serial version.

prob = SDDEProblem(kf, kg, ic, kh, (0.,T), p);
eprob = EnsembleProblem(prob)
@time sol = solve(eprob, SOSRA(), EnsembleSerial(), dt=1, adaptive=false, trajectories=10,saveat=0:1:T);
plot(hcat([kop(s) for s in sol]...),linealpha=0.2)

I guess all the figures are excessive, the main point is that the steady state of the ensemble,

quantile([kop(s)[end] for s in sol],[0.05,0.95])

should be around 0.6, not 0.1.

[ChrisRackauckas]

Yes, this looks like a bug that's worth digging into, thanks. It could be something to do with aliasing maybe. @devmotion had to work through a few of these in the DelayDiffEq implementation with similar properties. Rooting them out of the stochastic version is a bit harder but now that we have a good example it should be a lot easier.

[HTSykora]

I'll look into it, but i'll need some time.

[maedoc]

I saw this paper (doi: 10.1158/2326-6066.CIR-22-0617) using SDDEs w/ the SOSRI method and made me curious if you think this bug is fixed?

edit: not trolling 👹 I swear I never want to write another SDDE solver again evar

[ChrisRackauckas]

Sorry, I was not aware of this issue.

The issue is the lack of constrained fixed point handling. Note for example in the DDE solvers, https://docs.sciml.ai/DiffEqDocs/stable/solvers/dde_solve/#Note constrained=true uses a fixed point iteration when dt is greater than the smallest lag. The reason is because any method that isn't Euler's method will use some value at a time in (t, t + dt], but if that time is greater than t + tau, then the lag h(t-tau) in (t, t + dt]. This means that "the history" isn't historical at all during this case, which means h(t-tau) is extrapolating. This is fixed through using a fixed point iteration so that the extrapolatory value converges to the interpolatory value.

However I'm seeing now that @HTSykora's code did not handle this part, as noticed by the single call to perform_step! here https://github.com/SciML/StochasticDelayDiffEq.jl/blob/master/src/solve.jl#L516. This is not good.

@devmotion did you ever look into this?

What I can at least do for now is enforce the constrained version by enforcing that dt < minimum(constant_lags) at each time step, and throw an error check there. Of course, this means that lags=0 is the worst case scenario here. Other cases should be mostly okay because adaptive SDE integrators tend to use rather small time steps and most cases do not have the lag as 0, and this would explain why the convergence tests did not pick up on this.

#67 adds the error.