Adaptive solver ignores `step_ts` for impulse discontinuous ODE

I'm trying to integrate a jump discontinuous ODE using an adaptive solver (i.e. discontinuous at a specific timepoint).

However, the integration seems to ignore my step_ts:

import diffrax


t0, t1 = 0, 3
jump = 0.5  #  let's jump at this point
step_ts = [t0, jump, t1]


def vector_field(t, y, args):
    return -1.0 * y + 100.0 * (t == jump)  # also tried jnp.allclose but doesn't work either


term = diffrax.ODETerm(vector_field)
solver = diffrax.Tsit5()
saveat = diffrax.SaveAt(ts=step_ts)
stepsize_controller = diffrax.PIDController(
    pcoeff=0.0,
    icoeff=1.0,
    dcoeff=0.0,
    rtol=1e-4,
    atol=1e-6,
    step_ts=step_ts,
)

sol = diffrax.diffeqsolve(
    term,
    solver,
    t0=t0,
    t1=t1,
    dt0=None,
    y0=1,
    saveat=saveat,
    stepsize_controller=stepsize_controller,
)

print(sol.ts)
assert jump in sol.ts
print(sol.ys)

Output:

[0.  0.5 3. ]
[1.         0.6068427  0.04981759]  # should be something much larger since we have a large jump at t=0.5

I'm avoiding jump_ts based on https://github.com/patrick-kidger/diffrax/issues/58, which addressed the case where you have piecewise constant vector fields. However I'm not quite sure what to do when there's a delta-function in the vector field (impulse external control)

Is there a way to do this in diffrax?

Start with the reproducer using diffrax.diffeqsolve, ODETerm, and PIDController, then trace how step_ts is handled during adaptive integration and how the vector field is evaluated at the impulse time. Done means the documented discontinuous-ODE case has defined, correct behavior for an impulse at jump, with a regression test covering the result.