Silent bugs in Monte Carlo and averaging outputs

[bdestombe]

Three closely-related bugs in the Monte Carlo and averaging methods on DtsAccessor (src/dtscalibration/dts_accessor.py) that all produce silently-wrong outputs.

1. mc_remove_set_flag=True fails to remove talpha_*_mc arrays

src/dtscalibration/dts_accessor.py:1955-1956 and :2801-2802 build the cleanup list with a literal embedded double-quote:

remove_mc_set.append('talpha"_fw_mc')
remove_mc_set.append('talpha"_bw_mc')

The actual keys written into the output Dataset (e.g. :1782-1783) are talpha_fw_mc / talpha_bw_mc (no \"). The if k in out: guard at :1959 therefore always misses, and the MC arrays are never deleted. Users who set mc_remove_set_flag=True to get a lean output get a multi-GB Dataset instead.

Fix: drop the stray \" in both monte_carlo_double_ended and average_monte_carlo_double_ended.

2. ci_avg_time_flag1 is not None is always True

src/dtscalibration/dts_accessor.py:2230 (single-ended) and :2597 (double-ended):

if ci_avg_time_flag1 is not None:
    # unweighted mean
    out[\"tmpf_avg1\"] = ...

The default value of ci_avg_time_flag1 is False, not None, so the block always runs. The user's ci_avg_time_flag1=False is silently ignored, and tmpf_avg1 / tmpf_mc_avg1_var are added to the output regardless. This is inconsistent with the analogous spatial flag at :2179 which is correctly written as if ci_avg_x_flag1:.

Fix: change to if ci_avg_time_flag1: in both single-ended and double-ended average_monte_carlo_* methods.

3. np.random.normal leak breaks MC reproducibility

src/dtscalibration/dts_accessor.py:1742 uses the global RNG inside monte_carlo_double_ended:

if np.any(not_ix_sec):
    ... = np.random.normal(...)

Every other random draw in the same function uses the user-supplied state (a dask.array.random.RandomState). Whenever not_ix_sec.size > 0 (i.e. there is any x outside reference sections — the normal case), two calls to monte_carlo_double_ended with the same da_random_state produce different tmpf_mc_var results. The function's reproducibility contract is broken.

Fix: replace the global np.random.normal call with state.normal(...).

Suggested tests

• Build a synthetic double-ended dataset; run monte_carlo_double_ended twice with the same da_random_state; assert tmpf_mc_var is bit-exact equal (np.testing.assert_allclose(..., atol=0, rtol=0)).
• Run average_monte_carlo_single_ended with ci_avg_time_flag1=False, ci_avg_time_flag2=True; assert \"tmpf_avg1\" not in out and \"tmpf_mc_avg1_var\" not in out.
• Run monte_carlo_double_ended with mc_remove_set_flag=True and trans_att=[some_value]; assert \"talpha_fw_mc\" not in out and \"talpha_bw_mc\" not in out.

[bdestombe]

Two corrections / additions:

Line numbers for the RNG leak

The current revision of dts_accessor.py does not have a np.random.normal at :1742. The actual reproducibility leaks are at:

• :1343 — sst.multivariate_normal.rvs(...) in monte_carlo_single_ended
• :1736 — sst.multivariate_normal.rvs(...) in monte_carlo_double_ended
• :1758 — np.random.normal(...) for non-section alpha in monte_carlo_double_ended

da_random_state is honored only by the Stokes noise draws via state.normal(...) (lines 1410-1416, 1843-1849); all parameter draws bypass it.

The simplest fix is structural

scipy.stats.multivariate_normal.rvs accepts a random_state=... argument; pass state (or seed it from state to a np.random.Generator). For the np.random.normal call at :1758, switch to state.normal(...). If the new issue about the joint distribution of non-section alpha is fixed by extending from_i to cover all alpha indices, the line-1758 call goes away entirely.