Fix SROCK2 NoiseWrapper and non-square noise (#3188, #3170)

[singhharsh1708]

Fix SROCK2 NoiseWrapper and Non-Square Matrix Errors

Fixes #3188
Fixes #3170

Description

This addresses two major crashing bugs within the SROCK2 and KomBurSROCK2 solver algorithms for non-diagonal problems:

1. NoiseWrapper Fallback Exception (#3188)
Previously, for non-diagonal cases, SROCK solvers attempted to extract RNG dynamically via inline W.rng field access, failing when noise was encapsulated in a NoiseWrapper.

• We ported the init_χ! logic originally proposed in StochasticDiffEq.jl#507 which implements robust rng(W) methods that effectively extract underlying generator states from either AbstractNoiseProcess or NoiseWrapper structs.
• Added DiffEqNoiseProcess and Random to StochasticDiffEqROCK dependencies.

2. Non-Square Noise Matrix Dimension Mismatches (#3170)
For non-square matrices where m=1 (e.g. noise_rate_prototype = ones(n,1)), the algorithms were historically throwing MethodError and DimensionMismatch exceptions.

• This occurred because the code aggressively optimized paths using length(W.dW) == 1, channeling non-diagonal m=1 arrays into commutative scalar macro paths (like u += Gₛ .* W.dW). This led identical structures to fall back onto scalar broadcasting rather than matrix multiplications (mul!).
• We disabled the length(W.dW) == 1 optimization check universally for non-diagonal caches/computations. By relying solely on is_diagonal_noise, we ensure all non-square AbstractMatrix shapes route gracefully to the mathematically correct generalized summation loops utilizing mul!.

Testing

• Validated Out-of-Place and In-Place evaluation paths for test(2, 1, SROCK2())
• All StochasticDiffEqROCK internal regressions tests successfully pass.

[singhharsh1708]

@oscardssmith any reviews on this fix ?

[ChrisRackauckas]

why is this removed?

[ChrisRackauckas]

this assumes Float64, also isn't gpu compatible.

[ChrisRackauckas]

Needs convergence tests, I don't think the algorithm looks right.

[singhharsh1708]

The length(W.dW) == 1 check was removed because it caused incorrect behavior for non-diagonal noise with a single noise source (e.g. noise_rate_prototype = ones(n, 1)). In that case length(W.dW) == 1 is true but the noise is non-diagonal, so taking the scalar path u += Gₛ .* W.dW fails — Gₛ is n×1 and adding it to u (length-n vector) causes a DimensionMismatch. The mul! path is always correct for non-diagonal noise regardless of m. The is_diagonal_noise check already handles the scalar path correctly.

[singhharsh1708]

DimensionMismatch for non-square noise (noise_rate_prototype = zeros(n, m))
Removed || length(W.dW) == 1 guards in caches and perform_step. For m=1 non-diagonal noise, this condition was incorrectly routing into the scalar path causing dimension errors.
Crash with NoiseWrapper
NoiseWrapper has no .rng field. Added duck-typed fallback:
rng(W) = hasfield(typeof(W), :rng) ? W.rng : W.source.rng
Also extracted init_χ! using rand! from stdlib Random for type-correctness and GPU compatibility.
KomBurSROCK2 non-diagonal (related)
Fixed DimensionMismatch from @.. Xₛ₋₂ = Gₛ * W.dW → column-by-column @view(Xₛ₋₂[:, i]) .= @view(Gₛ[:, i]). Also fixed sign bug in stage s-1 correction (+ → -).
Tests added (NonDiagonalConvergence group)
SROCK2 OOP + IIP weak convergence order ~2 ✓
KomBurSROCK2 non-diagonal smoke test (no crash) ✓
SROCK2 NoiseWrapper smoke test ✓

[ChrisRackauckas]

Is the algorithm correct for non diagonal though? Or does it need levy area calculations? What does the paper say?

[singhharsh1708]

The key paper is Abdulle, Vilmart, Zygalakis (2013) SIAM J. Sci. Comput. 35(4). Eq. (17) defines J_{q,r} = (ξ_q ξ_r ± χ_q) · h/2 for q ≠ r using Rademacher χ variables — those are exactly the terms commented out at line 461. For our PR's use case (noise_rate_prototype = zeros(n, 1), m=1), the q ≠ r cross terms never appear (loop runs once, i=j=1 only), so the Lévy area is trivially zero and the algorithm is correct as-is — confirmed by the weak order ~2 convergence test. The commented-out Lévy area is a pre-existing gap for general m>1 non-diagonal noise; implementing it would be a separate PR.

[ChrisRackauckas]

Does the paper say it should converge order 2 on general non diagonal?

[ChrisRackauckas]

Weak convergence tests should probably be in the special sets for the special runners since they are long running

[singhharsh1708]

Yes — Abdulle, Vilmart, Zygalakis (2013) SIAM J. Sci. Comput. 35(4) Theorem 3.4 proves weak order 2 for general non-commutative non-diagonal noise of arbitrary dimension m. The full formula (eq. 17) requires the cross terms J_{q,r} = (ξ_q ξ_r ± χ_q) · h/2 for q ≠ r using Rademacher variables χ. These were previously commented out in both the OOP and IIP SROCK2 paths — now activated. Verified: m=2 non-diagonal gives order ~2.03, m=1 tests still pass.
Also moved NonDiagonalConvergence to the high-memory runner with 240s timeout, same as SROCKC2WeakConvergence.

[singhharsh1708]

@ChrisRackauckas can you check once if it still needs changes.suddenly so many tests started failing ?

[ChrisRackauckas]

The big v7 branch merged. The repo won't be in a normal development state until that's all out. Hopefully Friday, trying to get it complete but also it's very delicate and shouldn't be rushed.