Incorrect user-supplied dep_graph silently produces wrong results; proposal and prototype for a structural check

Whether this is filed as a bug or an enhancement is the maintainers' call. The aggregator honours the graph it is given, so nothing in the stepping loop is wrong. What is absent is any check on the graph itself, and the package holds the information needed to perform one. An error the library could catch at construction becomes a wrong published number instead.

An incorrect user-supplied dep_graph produces converged wrong results with no diagnostic, and a proposal for an opt-in structural check

Summary

For systems built from ConstantRateJumps, the documentation makes the
dependency graph the user's responsibility: "Dependency graphs are represented
as a Vector{Vector{Int}}, with the ith vector containing the indices of the
jumps for which rates must be recalculated when the ith jump occurs", and
auto-generation is available for MassActionJumps and for systems generated
from Catalyst. That obligation is discharged by hand, it is easy to discharge
incompletely, and an incomplete graph does not fail loudly. It produces a
trajectory that runs to completion and an ensemble average whose standard error
is far smaller than its own bias.

The reproducer below omits two edges from a dependency graph on a 24-jump
lattice model and obtains a coverage 2.76 per cent away from the exact
stationary value, at 36 standard errors, with no error, warning or diagnostic
from the package. A control isolates the omitted edges as the cause.

This is not a bug report against the aggregator. NRM does what the supplied
graph says. It is a request for an opt-in check on the supplied graph, since
the information needed to verify it is already available to the package, and I
have attached a working prototype of the check.

Reproducer

Full script: mwe_depgraph.jl (attached, self-contained, about one minute).

A periodic 1D lattice gas of N = 12 sites, each empty or occupied.
Adsorption fills an empty site at constant rate rA. Desorption empties an
occupied site at rD0 * exp(w * nbrOcc), so the desorption rate at site i
reads u[i-1] and u[i+1]. A jump at site i therefore changes the rates of
the desorption jumps at both of its neighbours.

Jumps are ordered 1..N adsorption, N+1..2N desorption.

# correct: a jump at site i changes u[i]; the rates reading u[i] are the two
# jumps at i plus the desorption jumps at both neighbours
function dep_correct()
    g = [Int[] for _ in 1:(2N)]
    for i in 1:N
        l = mod1(i - 1, N); r = mod1(i + 1, N)
        e = [i, N + i, N + l, N + r]
        g[i] = copy(e); g[N + i] = copy(e)
    end
    g
end

# incorrect: the two neighbour entries omitted
dep_stale() = [i <= N ? [i, N + i] : [i - N, i] for i in 1:(2N)]

Because the model satisfies detailed balance, the stationary distribution is
pi(s) ~ K^Nocc * exp(-w * pairs) with K = rA/rD0, so the reference coverage
comes from enumerating all 4096 configurations and involves no sampling and no
shared code path with the solver.

Output

JumpProcesses reproducer: an incorrect dep_graph is not detected.
w = 0 is a control. With no neighbour coupling the omitted edges
carry no information, so both graphs must agree, and they do.

w = 0.0   exact coverage = 0.50000
  correct dep_graph : 0.50027 +/- 0.00036   relative error +0.05%   1 sigma
  WRONG   dep_graph : 0.50027 +/- 0.00036   relative error +0.05%   1 sigma

w = 1.5   exact coverage = 0.31896
  correct dep_graph : 0.31918 +/- 0.00019   relative error +0.07%   1 sigma
  WRONG   dep_graph : 0.32777 +/- 0.00024   relative error +2.76%   36 sigma

Both conditions ran to completion. No error, warning or diagnostic
was emitted for the incorrect graph.

At longer runs the separation grows as expected: with 64 replications to
t = 20000 the incorrect graph gives 0.32735 +/- 0.00007 against the exact
0.31896, a bias 121 times its own standard error.

Three points about this output.

The control identifies the cause. At w = 0 the desorption rate does not
read the neighbours, the omitted edges carry no information, and the two graphs
agree to every digit, because with matched seeds they generate the same
trajectory. Whatever separates them at w = 1.5 is the omitted edges and not
the estimator or the reference.

The failure is invisible from the output. The incorrect run converges. Its
standard error is a faithful description of its own sampling variability. A user
without an independently computed answer sees a converged result and has no
signal that anything is wrong.

The error is not monotone in the coupling. At w = 3.0 the bias falls to
1.16 per cent, because under strong repulsion occupied sites are largely
isolated and stale and current rates coincide more often. A dependency graph
exercised in a weakly coupled regime can therefore pass and still be wrong where
the coupling matters.

Where this shows up

Lattice and surface models in which a rate at one site depends on the occupancy
of neighbouring sites. The spatial interface does not cover this case as far as
I can tell: SpatialMassActionJump with NSM or DirectCRDirect handles
mass-action reactions within a subvolume plus hopping between subvolumes, where
a reaction rate reads only its own subvolume's species counts. A rate that reads
a neighbour's state falls back to hand-written ConstantRateJumps and a
hand-written dependency graph, which is the configuration above.

Proposal

An opt-in structural check, run once before any trajectory, with no cost to the
stepping loop. The needed quantities are already in the JumpProblem.

reads[j]  = state indices whose perturbation changes the rate of jump j
writes[i] = state indices modified by the affect! of jump i
required  = { (i, j) : writes[i] intersects reads[j] }

Every required pair must appear in the declared graph: a jump that writes a
variable that some rate reads is exactly a jump after which that rate must be
recomputed. Anything required and not declared is reported, naming the pair.

Two shapes, either of which would have caught this:

1. validate_dep_graph(jprob; states = ..., nsamples = ...), a test-time
   utility a user calls in their own test suite.
2. JumpProblem(...; check_dep_graph = true), running the same check at
   construction over the initial condition plus sampled states, warning rather
   than erroring.

A prototype is attached as validate_depgraph.jl. It discovers reads by
perturbing each state component and writes by applying each affect! to a
mock integrator, then reports the required-but-undeclared pairs. On the
reproducer, run over 40 sampled states and taking under a second:

Structural check of the declared dependency graph, no trajectory run.

correct dep_graph
  required edges          : 96
  declared edges          : 96
  required but NOT declared: 0
  VERDICT: consistent

incorrect dep_graph
  required edges          : 96
  declared edges          : 48
  required but NOT declared: 48
      jump 1 writes a variable that the rate of jump 14 reads
      jump 1 writes a variable that the rate of jump 24 reads
      jump 2 writes a variable that the rate of jump 13 reads
      jump 2 writes a variable that the rate of jump 15 reads
      jump 3 writes a variable that the rate of jump 14 reads
      jump 3 writes a variable that the rate of jump 16 reads
      ... and 42 more
  VERDICT: INCOMPLETE, listed rates are never recomputed

The 48 undeclared pairs are the 24 sites times the two neighbour edges each,
which is exactly what was removed. The check names them rather than reporting a
boolean, so the user is told which rates are never recomputed and after which
jumps.

Limits of the proposal, which are real

The check samples the state space, so it falsifies rather than proves: a
dependency that never manifests at any sampled state is not found. It assumes
rates are deterministic functions of (u, p, t), so a rate reading external
mutable state is outside its reach. Perturbing a state component requires
knowing what a valid neighbouring state is, which is trivial for occupancy
models and needs a user-supplied sampler in general. It is a cheap filter that
catches the common hand-editing error, not a correctness proof for the graph.

The prototype uses a mock integrator to observe writes; inside the package the
real integrator would be used, and the sampled states could default to those
reachable from u0.

I am willing to turn this into a PR if the maintainers think either shape is
worth having, and to take direction on which.

Environment

Julia 1.12.7, Windows x86-64
JumpProcesses v9.32.0

Attached: mwe_depgraph.jl (reproducer), validate_depgraph.jl (prototype
check).

Start with the attached mwe_dep_graph.jl reproducer and validate_depgraph.jl prototype, then inspect the JumpProblem construction path and existing dependency-graph tests. Compare the proposed validate_dep_graph utility with the check_dep_graph construction option. Done means an agreed API and tests detect the omitted edges while preserving the stepping loop behavior.