Bug: Delay block in a feedback loop does not implement z^-1 correctly (breaks recursive digital filters)

Summary

In discrete mode (sampling_period set), a Delay block used as a z^-1 register (flip-flop)
inside a feedback loop does not behave as a unit delay. The basic
digital integrator y[k] = x[k] + y[k-1] (one Delay + one Adder)
produces wrong, solver-dependent output:

• With the default embedded adaptive solver (RKBS32), the accumulator advances at
  half rate (each value is duplicated).
• With every other tested solver (SSPRK22, RK4, EUF) and with fixed-step
  (adaptive=False), the loop fails to advance at all — only a single output
  sample is produced.

This makes it impossible to build recursive digital
structures (integrators, IIR filters, etc) from the textbook
"delay" (flip-flop) and "gain" (multiplier) primitives, even though Delay works correctly for feedforward
paths and the same recursion works correctly when realized as a single
DiscreteStateSpace/DiscreteTransferFunction block.

Environment

• pathsim 0.25.1
• Python 3.12
• numpy (any recent)

Minimal reproducible example

import numpy as np
from pathsim import Connection, Simulation
from pathsim.blocks import Adder, Delay, Scope, Source
from pathsim.solvers import RKBS32

fs = 1000.0
t_s = 1.0 / fs
N = 8
step = np.ones(N)
src = Source(lambda t: step[min(int(np.floor(t / t_s + 1e-9)), N - 1)])

# accumulator: y[k] = x[k] + y[k-1]  (one z^-1 flop in a feedback loop)
add = Adder()
z = Delay(tau=t_s, sampling_period=t_s)   # 1-sample delay
scope = Scope(sampling_period=t_s, t_wait=0.5 * t_s)

sim = Simulation(
    [src, add, z, scope],
    [Connection(src, add[0]), Connection(z, add[1]),
     Connection(add, z), Connection(add, scope)],
    dt_max=t_s, Solver=RKBS32, log=False,
)
sim.run(N * t_s)
_, data = scope.read()

print("expected:", np.cumsum(step).astype(int))
print("actual  :", np.round(np.asarray(data[0]), 2))

Expected

expected: [1 2 3 4 5 6 7 8]

Actual

actual  : [1. 2. 2. 3. 3. 4. 4. 5.]      # half rate

A fully runnable, self-contained script that reproduces this (with plots and the
sweeps below) is in pathsim_delay_feedback_bug.py.

Solver dependence

Same accumulator, sweeping solver and adaptive:

Expected in all cases: [1 2 3 4 5 6 7 8].

Ruled out (not the cause)

The half-rate output is not an artifact of the test harness:

• Scope sampling period — a finer scope does not recover the samples; it just
  oversamples the same half-rate staircase (scope=T_S/2 → [1 1 2 2 2 2 3 3 …]), so
  the accumulator genuinely updates only once per 2*T_S.
• Step size — smaller dt_max (adaptive) has no effect: T_S/1 … T_S/10 all give
  [1 2 2 3 3 4 4 5]. Fixed-step (adaptive=False) stalls at [1] for every step size.

So the behavior is independent of source, scope rate and step size, and depends only on
the solver's event handling for a Delay closing a feedback loop.

What works (for contrast — shows the bug is specific to Delay-in-a-loop)

• Feedforward Delay (no loop) is a correct z^-1: a delay line / FIR built from
  Delay + gains + Adder matches scipy.signal.lfilter to ~2e-16.
• The same integrator as a single block is exact: DiscreteTransferFunction([1, 0], [1, -1], T=t_s)
  driven by a unit step returns [1 2 3 4 5 6 7 8], and a 2nd-order IIR built as one
  DiscreteTransferFunction matches scipy.signal.lfilter to ~9e-16.

So the issue is not discrete-time filtering in general — it is specifically a Delay
block closing a feedback loop through the connection graph.

Suspected root cause

In discrete mode the Delay splits its behavior across two phases
(pathsim/blocks/delay.py):

• update(t) presents the output: self.outputs[0] = self._ring[0].
• A Schedule(t_period=sampling_period) event only sets a flag
  self._sample_next_timestep = True.
• sample(t, dt) (called at the end of a timestep) latches the input only if the
  flag is set: self._ring.append(self.inputs[0]).

For a feedforward path this is a clean unit delay. But inside a feedback loop the
"present ring[0]" and "latch inputs[0]" steps are separated by the algebraic
resolution, and the latch happens only on timesteps that coincide with the scheduled
sample event. The engine's stepping/event-localization does not place exactly one
latch per output sample when the block's newly latched input depends (through the
loop) on its own current output, so the ring advances at most once per two sampling
periods (half rate), or not at all for solvers whose stepping never aligns a latch
with the loop update.

Impact

Recursive digital blocks (integrators, biquads/IIR in direct form, delta-sigma loop
filters, accumulators) cannot be built from Delay + gain feedback the way a DSP
textbook or Simulink would draw them. The only reliable workaround is to collapse each
linear recursive loop into a single DiscreteStateSpace/DiscreteTransferFunction
block, which is not always possible (e.g. a loop containing a nonlinear element such
as a quantizer or comparator).