PositiveDomain callback failing on example

using DifferentialEquations
using Plots

function Scaled_Null(f,y,p,t)
    p1 = 1;
    n, g, p2, p3, p4, p5, p6, h1, h2, z = p;

    f[1] = p1-y[1]-y[1]*sum(y[3:2:end].*g);
    for i = 1:n
        f[2*i] = p2[i]-(p5[i]+p6[i]*z[i])*y[2*i];
        f[2*i+1] = z[i]*y[2*i]+p3[i]*y[2*i+1]*y[1]/(1+y[1])-p4[i]*y[2*i+1];
    end
end

function condition(u,t,integrator)
    t in integrator.p[8] || t in integrator.p[9]
 end

function affect!(integrator)
    k = findall(x->x==integrator.t,h1);
    for i = 1:length(k)
        integrator.p[10][k[i]] = rand();
    end
    k = findall(x->x==integrator.t,h2);
    for i = 1:length(k)
        integrator.p[10][k[i]] = 0.0;
    end
end

tspan = (1,3e+3); antigen_time = 1e+3;
a = 1e+2;
n = 4;
h1 = rand(1.0:antigen_time,n);
h2 = h1 .+ rand(1.0:antigen_time,n);
g = rand(n);
p2 = rand(1.0:10.0,n);
p3 = @. 0.01+(1-0.01)*$(rand(n));
p4 = @. 0.001+(0.1-0.001)*$(rand(n));
p5 = @. 0.001+(0.1-0.001)*$(rand(n));
p6 = @. 0.01+(1-0.01)*$(rand(n));
z = zeros(n);
p = [n, g, p2, p3, p4, p5, p6, h1, h2, z];

y0 = [1.0; vec(hcat(rand(1.0:a,n),rand(n))')];
c = unique(rand(1:n,rand(1:n)));
for i=1:length(c)
    y0[c[i]+1]=0.0;
end

antigen=DiscreteCallback(condition,affect!);
cb=CallbackSet(antigen,PositiveDomain());
dis=sort(unique([i for i in vcat(h1,h2) if i<=tspan[end]]));
prob=ODEProblem(Scaled_Null,y0,tspan,p);
sol=solve(prob,Rosenbrock23(),tstops=dis,callback=antigen)
sol=solve(prob,Rosenbrock23(),tstops=dis,callback=cb)
sol=solve(prob,Rosenbrock23(),tstops=dis,isoutofdomain=(y,p,t)->any(x->x<0,y),callback=antigen)

plot(sol)

Start with the supplied Julia reproducer, focusing on the PositiveDomain() callback combined with CallbackSet, Rosenbrock23(), and the tstops values. Compare the shown solves with and without PositiveDomain and with isoutofdomain; done means the failing example's behavior is understood and the callback works correctly in the affected configuration.