SciML / SciML/JumpProcesses.jl

JumpProblem(...; save_positions=...)` does not propagate to `VR_Direct`'s `ContinuousCallback

Open
#595 3 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

bug
Dominant language
Julia
Stars
150
Forks
41
Avg merge
1d 9h
Merged PRs (30d)
28

Description

When aggregating multiple VariableRateJumps with vr_aggregator = VR_Direct(), the save_positions keyword argument passed to JumpProblem(...) is silently ignored. The per-VariableRateJump default of (false, true) is used instead, inserting an extra save after every fired jump. With saveat, this produces ensemble trajectories of non-uniform length along the canonical grid a silent ensemble-average bias rather than an error.

I hit this while validating an MCWF unraveling of a Lindblad master equation. The fix turned out to be one line per jump, but the failure mode was non-obvious because the kwarg is accepted without warning.

using DifferentialEquations, JumpProcesses

const γ = (1.0, 1.0)
const H = ComplexF64[0 2 0; 2 0 2; 0 2 0]   # constant Ω = 2, Δ = 0

real_to_ψ(u) = ComplexF64[u[1]+1im*u[2], u[3]+1im*u[4], u[5]+1im*u[6]]
ψ_to_real(ψ) = Float64[real(ψ[1]),imag(ψ[1]),real(ψ[2]),imag(ψ[2]),real(ψ[3]),imag(ψ[3])]

function nojump!(du, u, p, t)
  ψ  = real_to_ψ(u)
  P2 = abs2(ψ[2]); γt = γ[1] + γ[2]
  dψ = -1im .* (H * ψ)
  dψ[1] +=  0.5 * γt * P2       * ψ[1]
  dψ[2] += -0.5 * γt * (1 - P2) * ψ[2]
  dψ[3] +=  0.5 * γt * P2       * ψ[3]
  du .= ψ_to_real(dψ)
end

rate1(u,p,t) = γ[1] * abs2(real_to_ψ(u)[2])
rate2(u,p,t) = γ[2] * abs2(real_to_ψ(u)[2])
affect1!(i)  = (fill!(i.u, 0.0); i.u[1] = 1.0)
affect2!(i)  = (fill!(i.u, 0.0); i.u[5] = 1.0)

ode = ODEProblem(nojump!, ψ_to_real(ComplexF64[1, 0, 0]), (0.0, 5.0))
ts  = collect(0:0.05:5)
N   = 500

jp_0 = JumpProblem(ode, Direct(),
      VariableRateJump(rate1, affect1!),
      VariableRateJump(rate2, affect2!);
      vr_aggregator = VR_Direct())

jp_a = JumpProblem(ode, Direct(),
      VariableRateJump(rate1, affect1!),
      VariableRateJump(rate2, affect2!);
      vr_aggregator  = VR_Direct(),
      save_positions = (false, false))   # ← silently ignored

jp_b = JumpProblem(ode, Direct(),
      VariableRateJump(rate1, affect1!; save_positions = (false, false)),
      VariableRateJump(rate2, affect2!; save_positions = (false, false));
      vr_aggregator  = VR_Direct(),
      save_positions = (false, false))

run_ens(jp) = solve(EnsembleProblem(jp), Tsit5(), EnsembleSerial();
                  trajectories = N, saveat = ts,
                  abstol = 1e-9, reltol = 1e-9)

for (label, jp) in [("(0) bare defaults", jp_0),
                  ("(A) JumpProblem kwarg only", jp_a),
                  ("(B) per-jump save_positions=(F,F)", jp_b)]
  sol   = run_ens(jp)
  n_off = count(traj -> length(traj.t) != length(ts), sol.u)
  if n_off > 0
      @warn "$label: $n_off / $N trajectories off-grid"
  else
      println("$label: clean — all $N trajectories match the saveat grid")
  end
end

Output

┌ Warning: (0) bare defaults: 261 / 500 trajectories off-grid
└ @ Main d:\Git repos\jumpAD\stage1\save_positions_repro.jl:58
┌ Warning: (A) JumpProblem kwarg only: 266 / 500 trajectories off-grid
└ @ Main d:\Git repos\jumpAD\stage1\save_positions_repro.jl:58
(B) per-jump save_positions=(F,F): clean — all 500 trajectories match the saveat grid
  • VariableRateJump's default is save_positions = (false, true) at src/jumps.jl:178, which inserts a save after each jump.
  • build_variable_integcallback (src/variable_rate.jl:380) constructs the VR_Direct ContinuousCallback by OR-reducing only the per-VariableRateJump save_positions fields.
  • JumpProblem(...; save_positions=...) is consumed by a different (constant-rate) callback path and never reaches VR_Direct.

The failure is silent and scales with jump frequency: ensemble averaging accumulates a bias that's small at low rates but dominates the signal in jump-rich regimes. There's no error or warning unless you explicitly check trajectory lengths, so it's easy to miss without a deterministic reference to compare against.


As a suggested fix I reccomend to make JumpProblem's save_positions kwarg, when explicitly passed, override the per-VariableRateJump defaults at construction or change the per-VariableRateJump default to (false, false) to match user intuition.

Happy to send a PR. Let me know which is preferred.

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start with the save_positions default in src/jumps.jl and the VR_Direct callback construction in src/variable_rate.jl around build_variable_integcallback. Reproduce the issue with the supplied EnsembleProblem example and inspect trajectory lengths with saveat. Done means an explicitly passed JumpProblem save_positions value affects VR_Direct and the ensemble trajectories stay on the requested grid.

Written by the indexing model from the issue text.

Assessment

Tech stack
julia
Domain
backend
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Quiet
Clarity
Clearly specified
Newbie friendliness
58/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.