generator.py

from __future__ import annotations

import json
from collections import deque
from collections.abc import Iterator
from concurrent.futures import Future, ThreadPoolExecutor
from functools import lru_cache, partial
from pathlib import Path
from queue import Full, Queue
from threading import Event, Thread

import numpy as np
from scipy.signal import lfilter

from cascade.interface import DataGenerator

_CHUNK = 2048


_STARTUP_CHUNK = 256
_RAMP_CHUNK = 1024


_SEASONAL_PERIODS = np.array(
    [4, 7, 12, 15, 24, 30, 48, 52, 60, 90, 96, 144, 168, 183, 240, 288,
     336, 365, 672, 730],
    dtype=np.float64,
)
_SEASONAL_PROBS = np.array(
    [0.01, 0.23, 0.02, 0.02, 0.07, 0.01, 0.06, 0.01, 0.08, 0.01,
     0.14, 0.07, 0.04, 0.01, 0.08, 0.08, 0.02, 0.02, 0.01, 0.01],
    dtype=np.float64,
)
_SEASONAL_PROBS /= _SEASONAL_PROBS.sum()


_SEASONAL_PAIRS = np.array(
    [[15, 60], [60, 240], [24, 168], [48, 336], [96, 672], [7, 365],
     [12, 52]],
    dtype=np.float64,
)

# Periods that dominate the revealed Cascade windows.  Counts are converted to
# a 70% empirical / 30% uniform blend so every generated chunk gets broad,
# deterministic coverage without cloning one pool snapshot's proportions.
_FORECAST_PERIODS = np.array(
    [7, 15, 24, 48, 60, 96, 144, 240, 288], dtype=np.float64
)
_FORECAST_PERIOD_COUNTS = np.array(
    [291, 18, 43, 36, 63, 147, 39, 37, 87], dtype=np.float64
)
_FORECAST_PERIOD_PROBS = (
    0.70 * _FORECAST_PERIOD_COUNTS / _FORECAST_PERIOD_COUNTS.sum()
    + 0.30 / _FORECAST_PERIODS.size
)


_FAMILIES: tuple[str, ...] = (
    "forecast_tasks",
    "tidal_harmonic",
    "trend_seasonal_ar",
    "regime_shift",
    "multiplicative",
    "ar2",
    "integrated",
    "threshold_ar",
    "chaotic",
    "spectral_gp",
    "long_memory",
    "ou_stochastic_vol",
    "physical_sensors",
    "seasonal_counts",
    "intermittent",
    "pulse_outlier",
    "conditional_stability",
    "stable_calendar_counts",
    "dock_occupancy",
    "piecewise_level",
)


# Families whose output is already a faithful observation: the artifact and
# domain-link stages would only blur the discrete structure they exist to teach.
_CLEAN_FAMILIES: frozenset[str] = frozenset(
    {"tidal_harmonic", "dock_occupancy", "piecewise_level"}
)


_DEFAULT_WEIGHTS: dict[str, float] = {
    "forecast_tasks": 0.435,
    "tidal_harmonic": 0.015,
    "trend_seasonal_ar": 0.095,
    "regime_shift": 0.095,
    "multiplicative": 0.06,
    "ar2": 0.105,
    "integrated": 0.095,
    "threshold_ar": 0.06,
    "chaotic": 0.02,
    "spectral_gp": 0.07,
    "long_memory": 0.07,
    "ou_stochastic_vol": 0.08,
    "physical_sensors": 0.08,
    "seasonal_counts": 0.045,
    "intermittent": 0.02,
    "pulse_outlier": 0.02,
    "conditional_stability": 0.07,
    "stable_calendar_counts": 0.015,
    "dock_occupancy": 0.10,
    "piecewise_level": 0.11,
}


def _validate_parameters(parameters: dict[str, float], label: str) -> None:
    if not all(np.isfinite(value) for value in parameters.values()):
        raise ValueError(f"{label} must contain only finite values")
    probability_names = {
        "forecast_task_hardness",
        "sa_clean_frac",
        "integrated_heavy_frac",
        "integrated_sv_frac",
        *(key for key in parameters if key.startswith(("observation.", "augment."))),
    }
    for name in probability_names:
        if not 0.0 <= parameters[name] <= 1.0:
            raise ValueError(f"{label}.{name} must be in [0, 1]")
    for name in (
        "tr_exc_lo",
        "tr_exc_hi",
        "gr_exc_lo",
        "gr_exc_hi",
        "sa_clean_lo",
        "sa_clean_hi",
    ):
        if parameters[name] < 0.0:
            raise ValueError(f"{label}.{name} must be non-negative")
    for lo_name, hi_name in (
        ("tr_exc_lo", "tr_exc_hi"),
        ("gr_exc_lo", "gr_exc_hi"),
        ("sa_clean_lo", "sa_clean_hi"),
        ("observation.irregular_hold_prob_lo", "observation.irregular_hold_prob_hi"),
        ("observation.shock_prob_lo", "observation.shock_prob_hi"),
    ):
        if parameters[lo_name] > parameters[hi_name]:
            raise ValueError(f"{label}.{lo_name} must be <= {hi_name}")


class Generator(DataGenerator):


    def __init__(self, config_dir: str, *, seed: int) -> None:
        cfg_path = Path(config_dir) / "config.json"
        cfg = json.loads(cfg_path.read_text(encoding="utf-8")) if cfg_path.is_file() else {}
        self._cfg = cfg
        self._seed = int(seed)
        self._min_len = int(cfg.get("min_length", 64))
        self._max_len = int(cfg.get("max_length", 4096))
        if self._min_len < 1 or self._max_len < self._min_len:
            raise ValueError(f"invalid length band [{self._min_len}, {self._max_len}]")
        weights = dict(_DEFAULT_WEIGHTS)
        for k, v in dict(cfg.get("family_weights", {})).items():
            if k in weights:
                weights[k] = float(v)
        w = np.asarray([weights[f] for f in _FAMILIES], dtype=np.float64)
        if not np.all(np.isfinite(w)) or w.min() < 0 or w.sum() <= 0:
            raise ValueError("family_weights must be finite, non-negative, and not all zero")
        self._weights = w / w.sum()
        curriculum = dict(cfg.get("curriculum", {}))
        self._curriculum_enabled = bool(curriculum.get("enabled", False))
        self._expected_budget_fraction = float(
            curriculum.get("expected_budget_fraction", 1.0)
        )
        if (
            not np.isfinite(self._expected_budget_fraction)
            or not 0.0 < self._expected_budget_fraction <= 1.0
        ):
            raise ValueError("curriculum.expected_budget_fraction must be in (0, 1]")
        self._curriculum_start = float(curriculum.get("start_fraction", 0.10))
        self._curriculum_end = float(curriculum.get("end_fraction", 0.70))
        if not 0.0 <= self._curriculum_start < self._curriculum_end <= 1.0:
            raise ValueError(
                "curriculum fractions must satisfy 0 <= start_fraction < end_fraction <= 1"
            )
        start_weights = dict(weights)
        for k, v in dict(curriculum.get("start_family_weights", {})).items():
            if k in start_weights:
                start_weights[k] = float(v)
        start_w = np.asarray([start_weights[f] for f in _FAMILIES], dtype=np.float64)
        if (
            not np.all(np.isfinite(start_w))
            or start_w.min() < 0
            or start_w.sum() <= 0
        ):
            raise ValueError(
                "curriculum.start_family_weights must be finite, non-negative, "
                "and not all zero"
            )
        self._start_weights = start_w / start_w.sum()

        self._tr_hi_frac = float(cfg.get("tr_hi_frac", 0.25))
        self._prefetch_depth = int(cfg.get("prefetch_depth", 2))
        if not 1 <= self._prefetch_depth <= 4:
            raise ValueError("prefetch_depth must be in [1, 4]")
        self._generation_workers = int(cfg.get("generation_workers", 1))
        if not 1 <= self._generation_workers <= 4:
            raise ValueError("generation_workers must be in [1, 4]")
        augment = dict(cfg.get("augment", {}))
        observation = dict(cfg.get("observation", {}))
        self._parameters = {
            "forecast_task_hardness": float(cfg.get("forecast_task_hardness", 1.0)),
            "tr_exc_lo": float(cfg.get("tr_exc_lo", 0.4)),
            "tr_exc_hi": float(cfg.get("tr_exc_hi", 3.0)),
            "gr_exc_lo": float(cfg.get("gr_exc_lo", 0.3)),
            "gr_exc_hi": float(cfg.get("gr_exc_hi", 2.0)),
            "sa_clean_frac": float(cfg.get("sa_clean_frac", 0.4)),
            "sa_clean_lo": float(cfg.get("sa_clean_lo", 0.02)),
            "sa_clean_hi": float(cfg.get("sa_clean_hi", 0.12)),
            "integrated_heavy_frac": float(
                cfg.get("integrated_heavy_frac", 0.25)
            ),
            "integrated_sv_frac": float(cfg.get("integrated_sv_frac", 0.30)),
            "observation.censor_rate": float(observation.get("censor_rate", 0.06)),
            "observation.quantize_rate": float(
                observation.get("quantize_rate", 0.07)
            ),
            "observation.regular_hold_rate": float(
                observation.get("regular_hold_rate", 0.04)
            ),
            "observation.irregular_hold_rate": float(
                observation.get("irregular_hold_rate", 0.0)
            ),
            "observation.irregular_hold_prob_lo": float(
                observation.get("irregular_hold_prob_lo", 0.01)
            ),
            "observation.irregular_hold_prob_hi": float(
                observation.get("irregular_hold_prob_hi", 0.10)
            ),
            "observation.shock_row_rate": float(
                observation.get("shock_row_rate", 0.0)
            ),
            "observation.shock_prob_lo": float(
                observation.get("shock_prob_lo", 0.001)
            ),
            "observation.shock_prob_hi": float(
                observation.get("shock_prob_hi", 0.015)
            ),
            "observation.positive_rate": float(
                observation.get("positive_rate", 0.0)
            ),
            "observation.count_rate": float(observation.get("count_rate", 0.0)),
            "observation.small_count_rate": float(
                observation.get("small_count_rate", 0.0)
            ),
            "augment.tsmixup": float(augment.get("tsmixup", 0.0)),
            "augment.pad_prefix": float(augment.get("pad_prefix", 0.0)),
        }
        start_parameters = dict(curriculum.get("start_parameters", {}))
        start_observation = dict(start_parameters.pop("observation", {}))
        start_augment = dict(start_parameters.pop("augment", {}))
        known_top_level = {
            key for key in self._parameters if "." not in key
        }
        unknown = set(start_parameters) - known_top_level
        unknown.update(
            f"observation.{key}"
            for key in start_observation
            if f"observation.{key}" not in self._parameters
        )
        unknown.update(
            f"augment.{key}"
            for key in start_augment
            if f"augment.{key}" not in self._parameters
        )
        if unknown:
            names = ", ".join(sorted(unknown))
            raise ValueError(f"unknown curriculum.start_parameters: {names}")
        self._start_parameters = dict(self._parameters)
        for key, value in start_parameters.items():
            self._start_parameters[key] = float(value)
        for key, value in start_observation.items():
            self._start_parameters[f"observation.{key}"] = float(value)
        for key, value in start_augment.items():
            self._start_parameters[f"augment.{key}"] = float(value)
        _validate_parameters(self._parameters, "final parameters")
        _validate_parameters(
            self._start_parameters, "curriculum.start_parameters"
        )

    @property
    def name(self) -> str:
        return str(self._cfg.get("name", "cascade-heat3-fast-learn-curriculum"))

    def _blend_at(self, token_progress: float) -> float:
        """Return the shared smoothstep blend for weights and difficulty."""
        if not self._curriculum_enabled:
            return 1.0
        position = (token_progress - self._curriculum_start) / (
            self._curriculum_end - self._curriculum_start
        )
        position = float(np.clip(position, 0.0, 1.0))
        return position * position * (3.0 - 2.0 * position)

    def _weights_at(self, token_progress: float) -> np.ndarray:
        """Blend easy-to-final family weights with a smoothstep schedule."""
        blend = self._blend_at(token_progress)
        if blend >= 1.0:
            return self._weights
        if blend <= 0.0:
            return self._start_weights
        return (1.0 - blend) * self._start_weights + blend * self._weights

    def _parameters_at(self, token_progress: float) -> dict[str, float]:
        """Blend all within-family, observation, and augmentation settings."""
        blend = self._blend_at(token_progress)
        if blend >= 1.0:
            return dict(self._parameters)
        if blend <= 0.0:
            return dict(self._start_parameters)
        return {
            key: (1.0 - blend) * self._start_parameters[key]
            + blend * final_value
            for key, final_value in self._parameters.items()
        }

    def _progress_at(self, emitted_points: float, target_points: int) -> float:
        """Calibrate nominal point progress to expected heat consumption."""
        return emitted_points / (target_points * self._expected_budget_fraction)

    def _build_chunk_with_rng(
        self,
        rng: np.random.Generator,
        lengths: np.ndarray,
        take: int,
        max_len: int,
        family_weights: np.ndarray,
        parameters: dict[str, float],
    ) -> list[np.ndarray | None]:
        """Build one chunk from the supplied deterministic random stream."""
        batch_size = int(lengths.size)
        builders = (
            partial(
                _forecast_tasks,
                hardness=parameters["forecast_task_hardness"],
            ),
            _tidal_harmonic,
            partial(
                _trend_seasonal_ar,
                hi_frac=self._tr_hi_frac,
                exc_lo=parameters["tr_exc_lo"],
                exc_hi=parameters["tr_exc_hi"],
                clean_frac=parameters["sa_clean_frac"],
                clean_lo=parameters["sa_clean_lo"],
                clean_hi=parameters["sa_clean_hi"],
            ),
            _regime_shift,
            partial(
                _multiplicative,
                hi_frac=self._tr_hi_frac,
                exc_lo=parameters["gr_exc_lo"],
                exc_hi=parameters["gr_exc_hi"],
            ),
            _ar2,
            partial(
                _integrated,
                heavy_frac=parameters["integrated_heavy_frac"],
                sv_frac=parameters["integrated_sv_frac"],
            ),
            _threshold_ar,
            _chaotic,
            _spectral_gp,
            _long_memory,
            _ou_stochastic_vol,
            _physical_sensors,
            _seasonal_counts,
            _intermittent,
            _pulse_outlier,
            _conditional_stability,
            _stable_calendar_counts,
            _dock_occupancy,
            _piecewise_level,
        )
        current_observation = {
            key.removeprefix("observation."): value
            for key, value in parameters.items()
            if key.startswith("observation.")
        }
        domain_positive_rate = current_observation.pop("positive_rate")
        domain_count_rate = current_observation.pop("count_rate")
        domain_small_count_rate = current_observation.pop("small_count_rate")
        fam_ids = rng.choice(len(_FAMILIES), size=batch_size, p=family_weights)
        chunk: list[np.ndarray | None] = [None] * batch_size
        for fam in range(len(_FAMILIES)):
            idx = np.nonzero(fam_ids == fam)[0]
            if idx.size == 0:
                continue
            block = builders[fam](rng, int(idx.size), max_len)
            family = _FAMILIES[fam]
            if family in _CLEAN_FAMILIES:
                # Preserve these specialists' exact observations. The artifact and
                # domain-link stages would turn a share of these scarce rows into
                # something else: tides into integer counts, and dock occupancy
                # into a blurred series whose runs are no longer exactly equal.
                block = _sanitize(block)
                for row, series_i in enumerate(idx):
                    length = int(lengths[series_i])
                    chunk[series_i] = np.ascontiguousarray(
                        block[row, :length], dtype=np.float64
                    )
                continue
            preserve_nonnegative = family in {
                "multiplicative",
                "physical_sensors",
                "seasonal_counts",
                "intermittent",
                "stable_calendar_counts",
            }
            artifacts = _measurement_artifacts(
                rng,
                block,
                preserve_nonnegative=preserve_nonnegative,
                preserve_integers=family in {
                    "seasonal_counts",
                    "intermittent",
                    "stable_calendar_counts",
                },
                allow_reverse=family in {
                    "trend_seasonal_ar",
                    "multiplicative",
                    "spectral_gp",
                    "long_memory",
                },
                allow_range_artifacts=family != "integrated",
                **current_observation,
            )
            block = _sanitize(
                _small_count_link(
                    rng,
                    _domain_observation(
                        rng,
                        artifacts,
                        already_nonnegative=preserve_nonnegative,
                        positive_rate=domain_positive_rate,
                        count_rate=domain_count_rate,
                    ),
                    domain_small_count_rate,
                )
            )
            for row, series_i in enumerate(idx):
                length = int(lengths[series_i])
                chunk[series_i] = np.ascontiguousarray(
                    block[row, :length], dtype=np.float64
                )

        if self._min_len == max_len:
            mix_rate = parameters["augment.tsmixup"]
            mixed = np.nonzero(rng.random(batch_size) < mix_rate)[0]
            for series_i in mixed:
                source = chunk[series_i]
                if source is None:
                    continue
                n_other = int(rng.integers(1, 3))
                others = rng.integers(0, batch_size, size=n_other)
                weights = rng.dirichlet(np.ones(n_other + 1))
                combined = weights[0] * source
                valid = True
                for j, other_i in enumerate(others):
                    other = chunk[int(other_i)]
                    if other is None:
                        valid = False
                        break
                    combined = combined + weights[j + 1] * other
                if valid:
                    chunk[series_i] = _sanitize(combined)

        pad_rate = parameters["augment.pad_prefix"]
        padded = np.nonzero(rng.random(batch_size) < pad_rate)[0]
        for series_i in padded:
            series = chunk[series_i]
            if series is None or series.size < 8:
                continue
            cut = int(rng.integers(series.size // 8, 3 * series.size // 4))
            series[:cut] = series[cut]
        return chunk

    def _build_chunk(
        self,
        seed: int,
        lengths: np.ndarray,
        take: int,
        max_len: int,
        family_weights: np.ndarray,
        parameters: dict[str, float],
    ) -> list[np.ndarray | None]:
        """Build an independently seeded chunk for ordered parallel execution."""
        return self._build_chunk_with_rng(
            np.random.default_rng(seed),
            lengths,
            take,
            max_len,
            family_weights,
            parameters,
        )

    def generate(self, n_series: int) -> Iterator[np.ndarray]:


        if n_series <= 0:
            return
        rng = np.random.default_rng(self._seed)
        max_len = self._max_len
        # stream_cpu requests token_budget // min_length + 2 rows. Recover the
        # budget so this fixed-length generator follows trainer token progress.
        target_points = max(1, max(n_series - 2, 1) * self._min_len)


        queue: Queue[object] = Queue(maxsize=self._prefetch_depth)
        stop = Event()
        done = object()

        def put(item: object) -> bool:
            while not stop.is_set():
                try:
                    queue.put(item, timeout=0.1)
                    return True
                except Full:
                    continue
            return False

        def produce() -> None:
            # Compatibility path: one worker consumes the original single RNG
            # directly, reproducing the proven merged-v12 corpus byte-for-byte.
            if self._generation_workers == 1:
                try:
                    produced = 0
                    emitted_points = 0
                    while produced < n_series and not stop.is_set():
                        if produced == 0:
                            batch_size = _STARTUP_CHUNK
                        elif produced == _STARTUP_CHUNK:
                            batch_size = _RAMP_CHUNK
                        else:
                            batch_size = _CHUNK
                        lengths = rng.integers(
                            self._min_len, max_len + 1, size=batch_size
                        )
                        take = min(batch_size, n_series - produced)
                        chunk_points = int(lengths[:take].sum())
                        midpoint_progress = self._progress_at(
                            emitted_points + 0.5 * chunk_points, target_points
                        )
                        chunk = self._build_chunk_with_rng(
                            rng,
                            lengths,
                            take,
                            max_len,
                            self._weights_at(midpoint_progress),
                            self._parameters_at(midpoint_progress),
                        )
                        if not put((chunk, take)):
                            return
                        produced += take
                        emitted_points += chunk_points
                except BaseException as exc:
                    put(exc)
                finally:
                    put(done)
                return

            executor = ThreadPoolExecutor(
                max_workers=self._generation_workers,
                thread_name_prefix="cascade-chunk",
            )
            try:
                produced = 0
                emitted_points = 0
                pending: deque[tuple[Future[list[np.ndarray | None]], int]] = deque()
                while (produced < n_series or pending) and not stop.is_set():
                    while (
                        produced < n_series
                        and len(pending) < self._generation_workers
                        and not stop.is_set()
                    ):
                        if produced == 0:
                            batch_size = _STARTUP_CHUNK
                        elif produced == _STARTUP_CHUNK:
                            batch_size = _RAMP_CHUNK
                        else:
                            batch_size = _CHUNK
                        lengths = rng.integers(
                            self._min_len, max_len + 1, size=batch_size
                        )
                        take = min(batch_size, n_series - produced)
                        chunk_points = int(lengths[:take].sum())
                        midpoint_progress = self._progress_at(
                            emitted_points + 0.5 * chunk_points, target_points
                        )
                        family_weights = self._weights_at(midpoint_progress)
                        parameters = self._parameters_at(midpoint_progress)
                        chunk_seed = int(
                            rng.integers(0, np.iinfo(np.uint64).max, dtype=np.uint64)
                        )
                        future = executor.submit(
                            self._build_chunk,
                            chunk_seed,
                            lengths,
                            take,
                            max_len,
                            family_weights,
                            parameters,
                        )
                        pending.append((future, take))
                        produced += take
                        emitted_points += chunk_points

                    if pending:
                        future, take = pending.popleft()
                        if not put((future.result(), take)):
                            return
            except BaseException as exc:
                put(exc)
            finally:
                executor.shutdown(wait=False, cancel_futures=True)
                put(done)

        producer = Thread(target=produce, name="cascade-generator", daemon=True)
        producer.start()
        try:
            while True:
                item = queue.get()
                if item is done:
                    break
                if isinstance(item, BaseException):
                    raise item
                chunk, take = item
                for arr in chunk[:take]:

                    if arr is None:
                        raise RuntimeError("internal: unfilled series slot")
                    yield arr
        finally:
            stop.set()
            producer.join(timeout=1.0)


def _ar1_batch(innov: np.ndarray, phi: np.ndarray) -> np.ndarray:


    n, L = innov.shape
    x = np.empty((n, L), dtype=np.float64)
    p = phi.reshape(n)
    for i in range(n):
        x[i] = lfilter([1.0], [1.0, -float(p[i])], innov[i])
    return x


def _ar1_grouped(
    innov: np.ndarray, group_ids: np.ndarray, phi_choices: np.ndarray
) -> np.ndarray:
    """Vectorized AR(1) filtering for a small, representative phi dictionary."""
    out = np.empty_like(innov)
    for group_id, phi in enumerate(phi_choices):
        rows = np.nonzero(group_ids == group_id)[0]
        if rows.size:
            out[rows] = lfilter(
                [1.0], [1.0, -float(phi)], innov[rows], axis=1
            )
    return out


def _ar2_batch(innov: np.ndarray, a1: np.ndarray, a2: np.ndarray) -> np.ndarray:

    n, L = innov.shape
    x = np.empty((n, L), dtype=np.float64)
    for i in range(n):
        x[i] = lfilter(
            [1.0], [1.0, -float(a1[i]), -float(a2[i])], innov[i]
        )
    return x


def _prefix_mean_std(
    x: np.ndarray, *, calibration_points: int = 512
) -> tuple[np.ndarray, np.ndarray]:

    prefix = x[:, : min(x.shape[1], calibration_points)]
    mean = prefix.mean(axis=1, keepdims=True)
    std = prefix.std(axis=1, keepdims=True)
    return mean, np.where(std < 1e-12, 1.0, std)


def _prefix_standardize(
    x: np.ndarray, *, center: bool = True, calibration_points: int = 512
) -> np.ndarray:
    mean, std = _prefix_mean_std(
        x, calibration_points=calibration_points
    )
    return (x - mean) / std if center else x / std


@lru_cache(maxsize=4)
def _seasonal_basis(L: int) -> tuple[np.ndarray, np.ndarray]:

    angle = (
        2.0
        * np.pi
        * np.arange(L, dtype=np.float64)[None, :]
        / _SEASONAL_PERIODS[:, None]
    )
    return np.sin(angle), np.cos(angle)


def _seasonal(rng: np.random.Generator, n: int, L: int, k_max: int = 3) -> np.ndarray:

    t = np.arange(L, dtype=np.float64)[None, :]
    sin_basis, cos_basis = _seasonal_basis(L)
    k = rng.integers(1, k_max + 1, size=n)
    pair = _SEASONAL_PAIRS[
        rng.integers(0, len(_SEASONAL_PAIRS), size=n)
    ]
    use_pair = rng.random(n) < 0.35
    out = np.zeros((n, L), dtype=np.float64)
    for j in range(k_max):
        active = np.nonzero(k > j)[0]
        per = rng.choice(
            _SEASONAL_PERIODS, size=n, p=_SEASONAL_PROBS
        )
        if j < 2:
            per = np.where(use_pair, pair[:, j], per)
        per = per[:, None]
        amp = rng.uniform(0.2, 2.0, size=n)[:, None]
        phase = rng.uniform(0.0, 2.0 * np.pi, size=n)[:, None]


        basis_idx = np.searchsorted(_SEASONAL_PERIODS, per[active, 0])
        component = amp[active] * (
            sin_basis[basis_idx] * np.cos(phase[active])
            + cos_basis[basis_idx] * np.sin(phase[active])
        )


        modulated = np.nonzero((k > j) & (rng.random(n) < 0.35))[0]
        if modulated.size:

            modulated_local = np.searchsorted(active, modulated)
            modulated_arg = (
                2.0 * np.pi * t / per[modulated] + phase[modulated]
            )
            m_per = np.clip(
                per[modulated] * rng.uniform(
                    4.0, 12.0, size=(modulated.size, 1)
                ),
                32.0,
                2.0 * L,
            )
            m_phase = rng.uniform(
                0.0, 2.0 * np.pi, size=(modulated.size, 1)
            )
            slow = np.sin(2.0 * np.pi * t / m_per + m_phase)
            amp_mod = 1.0 + rng.uniform(
                0.05, 0.45, size=(modulated.size, 1)
            ) * slow
            phase_mod = rng.uniform(
                0.05, 0.75, size=(modulated.size, 1)
            ) * np.sin(2.0 * np.pi * t / (1.7 * m_per) - m_phase)
            component[modulated_local] = (
                amp[modulated]
                * amp_mod
                * np.sin(modulated_arg + phase_mod)
            )
        out[active] += component
    return out


def _sparse_jumps(rng: np.random.Generator, n: int, L: int, rate: float, scale) -> np.ndarray:


    mask = rng.random((n, L)) < rate
    mask[:, 0] = False
    rows, cols = np.nonzero(mask)
    jumps = np.zeros((n, L), dtype=np.float64)
    if rows.size == 0:
        return jumps


    s = np.asarray(scale, dtype=np.float64)
    event_scale = s if s.ndim == 0 else s.reshape(n)[rows]
    jumps[rows, cols] = rng.normal(0.0, 1.0, size=rows.size) * event_scale
    return jumps


def _measurement_artifacts(
    rng: np.random.Generator,
    block: np.ndarray,
    *,
    preserve_nonnegative: bool,
    preserve_integers: bool = False,
    allow_reverse: bool = True,
    allow_range_artifacts: bool = True,
    censor_rate: float = 0.06,
    quantize_rate: float = 0.07,
    regular_hold_rate: float = 0.04,
    irregular_hold_rate: float = 0.0,
    irregular_hold_prob_lo: float = 0.01,
    irregular_hold_prob_hi: float = 0.10,
    shock_row_rate: float = 0.0,
    shock_prob_lo: float = 0.001,
    shock_prob_hi: float = 0.015,
) -> np.ndarray:


    original = np.asarray(block, dtype=np.float64)
    out = original.copy()
    n, L = out.shape

    reverse = (
        rng.random(n) < 0.06
        if allow_reverse
        else np.zeros(n, dtype=bool)
    )
    out[reverse] = out[reverse, ::-1]

    if not preserve_nonnegative:
        invert = rng.random(n) < 0.04
        out[invert] *= -1.0


    calibration_len = min(L, 512)
    shocked = np.nonzero(rng.random(n) < shock_row_rate)[0]
    if shocked.size and L > 1:
        diff = np.diff(out[shocked, :calibration_len], axis=1)
        center = np.median(diff, axis=1, keepdims=True)
        robust_scale = 1.4826 * np.median(
            np.abs(diff - center), axis=1, keepdims=True
        )
        fallback = np.maximum(np.std(diff, axis=1, keepdims=True), 1e-9)
        robust_scale = np.where(robust_scale > 1e-9, robust_scale, fallback)
        event_prob = rng.uniform(
            shock_prob_lo, shock_prob_hi, size=(shocked.size, 1)
        )
        event_rows, event_cols = np.nonzero(
            rng.random((shocked.size, L)) < event_prob
        )
        if event_rows.size:
            favored_sign = rng.choice(
                [-1.0, 1.0], size=(shocked.size, 1)
            )
            sign = np.where(
                rng.random(event_rows.size) < 0.75,
                favored_sign[event_rows, 0],
                -favored_sign[event_rows, 0],
            )
            magnitude = rng.lognormal(
                mean=np.log(4.0), sigma=0.6, size=event_rows.size
            )
            out[
                shocked[event_rows], event_cols
            ] += sign * magnitude * robust_scale[event_rows, 0]
        if preserve_nonnegative:
            np.maximum(out, 0.0, out=out)


    for row in np.nonzero(rng.random(n) < censor_rate)[0]:
        q = float(rng.uniform(0.03, 0.18))
        upper = rng.random() < 0.5
        if not allow_range_artifacts:
            continue
        calibration = out[row, :calibration_len]
        if upper:
            threshold = np.quantile(calibration, 1.0 - q)
            out[row] = np.minimum(out[row], threshold)
        else:
            threshold = np.quantile(calibration, q)
            out[row] = np.maximum(out[row], threshold)

    quantized = np.nonzero(rng.random(n) < quantize_rate)[0]
    if quantized.size:
        levels = rng.integers(16, 257, size=(quantized.size, 1))
        if allow_range_artifacts:
            x = out[quantized]
            calibration = x[:, :calibration_len]
            lo = calibration.min(axis=1, keepdims=True)
            hi = calibration.max(axis=1, keepdims=True)
            step = (hi - lo) / np.maximum(levels - 1, 1)
            safe_step = np.where(step < 1e-12, 1.0, step)
            clipped = np.clip(x, lo, hi)
            out[quantized] = (
                lo + np.rint((clipped - lo) / safe_step) * safe_step
            )


    held = np.nonzero(rng.random(n) < regular_hold_rate)[0]
    if held.size:
        factors = rng.choice([2, 4, 8], size=held.size, p=[0.55, 0.30, 0.15])
        for factor in (2, 4, 8):
            rows = held[factors == factor]
            if rows.size:
                out[rows] = np.repeat(
                    out[rows, ::factor], factor, axis=1
                )[:, :L]


    irregular = np.nonzero(rng.random(n) < irregular_hold_rate)[0]
    if irregular.size and L > 1:
        hold_prob = rng.uniform(
            irregular_hold_prob_lo,
            irregular_hold_prob_hi,
            size=(irregular.size, 1),
        )
        hold = rng.random((irregular.size, L)) < hold_prob
        hold[:, 0] = False
        source_index = np.where(
            ~hold, np.arange(L, dtype=np.int64)[None, :], 0
        )
        np.maximum.accumulate(source_index, axis=1, out=source_index)
        out[irregular] = np.take_along_axis(
            out[irregular], source_index, axis=1
        )
    if preserve_integers:
        out = np.maximum(np.rint(out), 0.0)


    degenerate = out[:, :calibration_len].std(axis=1) < 1e-9
    out[degenerate] = original[degenerate]
    return out


def _forecast_period_quota(
    rng: np.random.Generator, n: int
) -> np.ndarray:
    """Return an exact, shuffled period quota rather than iid period draws."""
    expected = n * _FORECAST_PERIOD_PROBS
    counts = np.floor(expected).astype(np.int64)
    remainder = n - int(counts.sum())
    if remainder:
        fractional = expected - counts
        order = np.argsort(-fractional, kind="stable")
        counts[order[:remainder]] += 1
    periods = np.repeat(_FORECAST_PERIODS, counts)
    rng.shuffle(periods)
    return periods


def _forecast_tasks(
    rng: np.random.Generator,
    n: int,
    L: int,
    *,
    hardness: float = 1.0,
) -> np.ndarray:
    """Compile patch-aligned, conditionally forecastable 64-step tasks.

    Every row composes a slowly changing deterministic seasonal template,
    locally changing trend, persistent residual, context-visible volatility,
    and sparse precursor/event/recovery episodes.  Event boundaries are exact
    multiples of 32 (with both 64-grid parities represented), matching the
    trainer's 32-point patches and two-patch evaluation horizon.
    """
    if n == 0 or L == 0:
        return np.empty((n, L), dtype=np.float64)

    hardness = float(np.clip(hardness, 0.0, 1.0))
    t = np.arange(L, dtype=np.float64)[None, :]
    tn = t / max(L - 1, 1)

    # Deterministic per-chunk period coverage, followed by seeded shuffling.
    period = _forecast_period_quota(rng, n)[:, None]
    phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    slow_period = period * rng.uniform(5.0, 12.0, size=(n, 1))
    slow_phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    slow_angle = 2.0 * np.pi * t / slow_period + slow_phase
    amplitude = rng.uniform(0.55, 1.8, size=(n, 1))
    amp_mod = 1.0 + rng.uniform(0.04, 0.28, size=(n, 1)) * np.sin(slow_angle)
    phase_mod = rng.uniform(0.02, 0.30, size=(n, 1)) * np.sin(
        slow_angle * rng.uniform(0.55, 0.85, size=(n, 1)) - slow_phase
    )
    angle = 2.0 * np.pi * t / period + phase + phase_mod
    seasonal = amplitude * amp_mod * np.sin(angle)
    seasonal += (
        amplitude
        * rng.uniform(0.08, 0.38, size=(n, 1))
        * np.sin(2.0 * angle + rng.uniform(0.0, 2.0 * np.pi, size=(n, 1)))
    )

    # A smooth global trend plus visible local slope corrections.  Centers sit
    # 64/96 points before a task boundary, leaving enough context to estimate
    # the new slope. ReLU ramps are continuous and substantially cheaper than
    # materialising multiple softplus fields over every 4096-point row.
    trend = (
        rng.normal(0.0, 1.1, size=(n, 1)) * (tn - 0.5)
        + rng.normal(0.0, 0.45, size=(n, 1)) * (tn - 0.5) ** 2
    )
    cuts = np.arange(256, max(257, L - 192), 32, dtype=np.int64)
    if cuts.size:
        for _ in range(2):
            center = rng.choice(cuts, size=n)[:, None] - rng.choice(
                [64, 96], size=(n, 1)
            )
            delta_slope = rng.normal(0.0, 0.0045, size=(n, 1))
            trend += delta_slope * np.maximum(t - center, 0.0)

    # Sparse stateful episodes: a 64-step visible precursor identifies sign and
    # approximate magnitude, then a 64-step event and longer smooth recovery.
    event = np.zeros((n, L), dtype=np.float64)
    event_strength = np.zeros((n, L), dtype=np.float64)
    if cuts.size:
        max_slots = 2
        slot_probability = 0.40 + 0.45 * hardness
        for _ in range(max_slots):
            active = rng.random((n, 1)) < slot_probability
            boundary = rng.choice(cuts, size=n)[:, None]
            rel = t - boundary
            sign = rng.choice([-1.0, 1.0], size=(n, 1))
            magnitude = rng.uniform(
                0.45, 1.15 + 1.10 * hardness, size=(n, 1)
            )
            precursor_fraction = rng.uniform(0.18, 0.42, size=(n, 1))

            precursor = np.where(
                (rel >= -64.0) & (rel < 0.0),
                precursor_fraction * (rel + 64.0) / 64.0,
                0.0,
            )
            onset = np.where(
                (rel >= 0.0) & (rel < 16.0),
                precursor_fraction
                + (1.0 - precursor_fraction) * rel / 16.0,
                0.0,
            )
            episode = np.where(
                (rel >= 16.0) & (rel < 64.0),
                1.0 - 0.65 * (rel - 16.0) / 48.0,
                0.0,
            )
            recovery = np.where(
                (rel >= 64.0) & (rel < 192.0),
                0.35 * (1.0 - (rel - 64.0) / 128.0),
                0.0,
            )
            profile = (precursor + onset + episode + recovery) * active
            event += sign * magnitude * profile
            event_strength += magnitude * np.abs(profile)

    # Persistent residuals with slowly varying scale.  The current volatility
    # state and event precursor make forecast width observable from context.
    vol_period = rng.uniform(384.0, 1536.0, size=(n, 1))
    vol_phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    base_sigma = rng.uniform(
        0.025 + 0.025 * hardness,
        0.10 + 0.14 * hardness,
        size=(n, 1),
    )
    sigma = base_sigma * (
        1.0
        + rng.uniform(0.10, 0.45, size=(n, 1))
        * (0.5 + 0.5 * np.sin(2.0 * np.pi * t / vol_period + vol_phase))
        + 0.18 * event_strength
    )
    innovations = rng.normal(0.0, 1.0, size=(n, L)) * sigma
    phi_choices = np.array(
        [0.50, 0.68, 0.82, 0.90 + 0.05 * hardness], dtype=np.float64
    )
    phi_ids = rng.integers(0, phi_choices.size, size=n)
    residual = _ar1_grouped(innovations, phi_ids, phi_choices)

    out = _prefix_standardize(seasonal + trend + event + residual)
    out *= np.exp(rng.uniform(np.log(0.35), np.log(6.0), size=(n, 1)))
    out += rng.normal(0.0, 2.5, size=(n, 1))
    return out


def _tidal_harmonic(
    rng: np.random.Generator, n: int, L: int
) -> np.ndarray:
    """Vectorized multi-constituent tides adapted from round-1211 chal57."""
    if n == 0 or L == 0:
        return np.empty((n, L), dtype=np.float64)

    t = np.arange(L, dtype=np.float64)[None, :]
    base_period = np.exp(rng.uniform(np.log(16.0), np.log(400.0), n))
    ratios = np.array([1.0, 0.9661, 1.0191, 1.927, 2.0787])
    n_constituents = rng.integers(2, 5, n)
    out = np.zeros((n, L), dtype=np.float64)
    for constituent, ratio in enumerate(ratios):
        active = (constituent < n_constituents).astype(np.float64)
        jitter = 1.0 + rng.uniform(-0.004, 0.004, n)
        period = base_period * ratio * jitter
        amplitude = (
            np.ones(n)
            if constituent == 0
            else rng.uniform(0.15, 0.7, n)
        ) * active
        phase = rng.uniform(0.0, 2.0 * np.pi, n)
        out += amplitude[:, None] * np.sin(
            2.0 * np.pi * t / period[:, None] + phase[:, None]
        )

    modulation_period = base_period * rng.uniform(10.0, 40.0, n)
    modulation = rng.uniform(0.1, 0.5, n)
    out *= 1.0 + modulation[:, None] * np.sin(
        2.0 * np.pi * t / modulation_period[:, None]
        + rng.uniform(0.0, 2.0 * np.pi, n)[:, None]
    )

    # Quantized phi groups preserve chal57's 0.2–0.9 AR residual range while
    # filtering whole row groups in four SciPy calls rather than one call/row.
    phi_choices = np.array([0.25, 0.48, 0.70, 0.88], dtype=np.float64)
    phi_ids = rng.integers(0, phi_choices.size, n)
    sigma = rng.uniform(0.03, 0.15, n)
    out += _ar1_grouped(
        rng.normal(0.0, 1.0, (n, L)) * sigma[:, None],
        phi_ids,
        phi_choices,
    )
    out += (
        np.cumsum(rng.normal(0.0, 1.0, (n, L)), axis=1)
        * rng.uniform(0.0, 0.02, n)[:, None]
        / np.sqrt(max(L, 1))
    )

    out = _prefix_standardize(out, calibration_points=L)
    datum = rng.random(n) < 0.4
    datum_offset = (
        np.exp(rng.uniform(np.log(2.0), np.log(50.0), n)) * datum
    )
    scale = np.exp(rng.uniform(np.log(0.02), np.log(200.0), n))
    offset = rng.normal(0.0, 3.0, n) * scale + datum_offset * scale
    return out * scale[:, None] + offset[:, None]


def _trend_seasonal_ar(rng: np.random.Generator, n: int, L: int, *,
                       hi_frac: float = 0.25, exc_lo: float = 0.4,
                       exc_hi: float = 3.0, clean_frac: float = 0.4,
                       clean_lo: float = 0.02, clean_hi: float = 0.12) -> np.ndarray:
    t = np.arange(L, dtype=np.float64)[None, :]
    level = rng.normal(0.0, 1.0, size=(n, 1))


    _hi = rng.random((n, 1)) < hi_frac
    exc = np.where(_hi, rng.normal(0.0, exc_hi, size=(n, 1)),
                   rng.normal(0.0, exc_lo, size=(n, 1)))
    tn = t / max(L - 1, 1)
    series = level + exc * tn + _seasonal(rng, n, L)
    phi = rng.uniform(0.0, 0.85, size=n)
    clean = rng.random((n, 1)) < clean_frac
    sigma = np.where(
        clean,
        rng.uniform(clean_lo, clean_hi, size=(n, 1)),
        rng.uniform(0.1, 0.6, size=(n, 1)),
    )
    # Calibrate stationary SD to sigma without changing the AR(1) ACF.
    innov = (
        rng.normal(0.0, 1.0, size=(n, L))
        * sigma
        * np.sqrt(np.maximum(1.0 - phi[:, None] ** 2, 1e-6))
    )
    return series + _ar1_batch(innov, phi)


def _regime_shift(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    level = np.cumsum(_sparse_jumps(rng, n, L, rate=3.0 / L, scale=2.0), axis=1)
    log_vol = np.cumsum(_sparse_jumps(rng, n, L, rate=3.0 / L, scale=0.5), axis=1)
    vol = np.exp(np.clip(log_vol, -3.0, 3.0)) * rng.uniform(0.1, 0.5, size=(n, 1))
    noise = rng.normal(0.0, 1.0, size=(n, L)) * vol
    seas = _seasonal(rng, n, L, k_max=2) * rng.uniform(0.0, 1.0, size=(n, 1))


    slope = rng.normal(0.0, 1.0 / L, size=(n, 1)) + np.cumsum(
        _sparse_jumps(rng, n, L, rate=2.0 / L, scale=4.0 / L), axis=1
    )
    piecewise_trend = np.cumsum(slope, axis=1)
    return level + piecewise_trend + seas + noise


def _multiplicative(rng: np.random.Generator, n: int, L: int, *,
                    hi_frac: float = 0.25, exc_lo: float = 0.3, exc_hi: float = 2.0) -> np.ndarray:
    t = np.arange(L, dtype=np.float64)[None, :]

    _hg = rng.random((n, 1)) < hi_frac
    gexc = np.where(_hg, rng.normal(0.0, exc_hi, size=(n, 1)),
                    rng.normal(0.0, exc_lo, size=(n, 1)))
    tn = t / max(L - 1, 1)
    base_level = np.exp(gexc * tn + rng.normal(0.0, 0.3, size=(n, 1)))
    amp = rng.uniform(0.1, 0.6, size=(n, 1))
    seasonal_shape = _seasonal(rng, n, L, k_max=1)
    seasonal_shape = _prefix_standardize(seasonal_shape, center=False)
    seas = 1.0 + amp * seasonal_shape
    noise = 1.0 + rng.normal(0.0, 1.0, size=(n, L)) * rng.uniform(0.02, 0.15, size=(n, 1))
    scale = rng.uniform(1.0, 50.0, size=(n, 1))
    return scale * base_level * np.clip(seas, 0.05, None) * np.clip(noise, 0.05, None)


def _ar2(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    p1 = rng.uniform(0.3, 0.98, size=n)
    p2 = rng.uniform(-0.6, 0.6, size=n)
    a2 = p2
    a1 = p1 * (1.0 - p2)
    sigma = rng.uniform(0.2, 0.8, size=(n, 1))
    burn = 512
    # Calibrate stationary SD to sigma without changing the AR(2) ACF.
    innov = (
        rng.normal(0.0, 1.0, size=(n, L + burn))
        * sigma
        * np.sqrt(
            np.maximum(
                (1.0 - p1[:, None] ** 2) * (1.0 - p2[:, None] ** 2),
                1e-6,
            )
        )
    )


    return _ar2_batch(innov, a1, a2)[:, burn:]


def _integrated(
    rng: np.random.Generator,
    n: int,
    L: int,
    *,
    heavy_frac: float = 0.25,
    sv_frac: float = 0.30,
) -> np.ndarray:


    order2 = rng.random(n) < 0.35
    drift = rng.normal(0.0, 0.02, size=(n, 1))
    sigma = rng.uniform(0.2, 1.0, size=(n, 1))
    eps = rng.normal(0.0, 1.0, size=(n, L))
    heavy = np.nonzero(rng.random(n) < heavy_frac)[0]
    if heavy.size:
        df = rng.uniform(3.0, 12.0, size=(heavy.size, 1))
        eps[heavy] = rng.standard_t(df, size=(heavy.size, L)) / np.sqrt(
            df / (df - 2.0)
        )

    stochastic = np.nonzero(rng.random(n) < sv_frac)[0]
    if stochastic.size:
        phi = 0.995
        burn = 256
        vol_innov = (
            rng.standard_normal((stochastic.size, L + burn))
            * np.sqrt(1.0 - phi * phi)
        )
        log_vol = lfilter(
            [1.0],
            [1.0, -phi],
            vol_innov,
            axis=1,
        )[:, burn:]
        log_vol *= rng.uniform(0.10, 0.55, size=(stochastic.size, 1))
        eps[stochastic] *= np.exp(np.clip(log_vol, -2.0, 2.0))

    steps = eps * sigma + drift
    walk = np.cumsum(steps, axis=1)
    walk2 = np.cumsum(walk, axis=1)
    o2 = order2[:, None]


    return np.where(o2, walk2 / max(L, 1), walk)


def _threshold_ar(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    phi_hi = rng.uniform(0.3, 0.9, size=n)
    phi_lo = rng.uniform(-0.9, 0.3, size=n)
    const_hi = rng.normal(0.0, 0.3, size=n)
    const_lo = rng.normal(0.0, 0.3, size=n)
    sigma = rng.uniform(0.2, 0.7, size=(n, 1))
    burn = 256
    total = L + burn
    innov = rng.normal(0.0, 1.0, size=(n, total)) * sigma
    x = np.empty((n, total), dtype=np.float64)
    x[:, 0] = innov[:, 0]
    for t in range(1, total):
        prev = x[:, t - 1]
        hi = prev >= 0.0
        phi = np.where(hi, phi_hi, phi_lo)
        const = np.where(hi, const_hi, const_lo)
        x[:, t] = np.clip(const + phi * prev + innov[:, t], -1e6, 1e6)
    return x[:, burn:]


def _chaotic(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    use_sine = rng.random(n) < 0.5
    r_log = rng.uniform(3.6, 4.0, size=n)
    r_sin = rng.uniform(0.85, 1.0, size=n)
    x0 = rng.uniform(0.05, 0.95, size=n)
    cur = x0.copy()
    for _ in range(64):
        nxt_log = r_log * cur * (1.0 - cur)
        nxt_sin = r_sin * np.sin(np.pi * cur)
        cur = np.clip(np.where(use_sine, nxt_sin, nxt_log), 0.0, 1.0)
    x = np.empty((n, L), dtype=np.float64)
    x[:, 0] = cur
    for t in range(1, L):
        nxt_log = r_log * cur * (1.0 - cur)
        nxt_sin = r_sin * np.sin(np.pi * cur)
        cur = np.where(use_sine, nxt_sin, nxt_log)
        cur = np.clip(cur, 0.0, 1.0)
        x[:, t] = cur
    return _prefix_standardize(x)


def _spectral_gp(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    embed = 2 * L
    f = np.fft.rfftfreq(embed)[None, :]
    lengthscale = np.exp(rng.uniform(np.log(8.0), np.log(256.0), size=(n, 1)))
    envelope = np.exp(-0.5 * (2.0 * np.pi * lengthscale * f) ** 2)
    z = rng.standard_normal((n, f.shape[1])) + 1j * rng.standard_normal((n, f.shape[1]))
    z[:, 0] = 0.0
    x = np.fft.irfft(z * np.sqrt(envelope), n=embed, axis=1)[:, :L]
    return _prefix_standardize(x)


def _long_memory(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    embed = 2 * L
    f = np.fft.rfftfreq(embed)
    safe_f = np.maximum(f, 1.0 / embed)[None, :]
    beta = rng.uniform(-0.6, 2.4, size=(n, 1))
    amp = safe_f ** (-0.5 * beta)


    multiscale = rng.random((n, 1)) < 0.4
    split_idx = rng.integers(8, max(9, f.size // 3), size=(n, 1))
    split_f = np.maximum(split_idx / embed, 1.0 / embed)
    beta_hi = rng.uniform(-0.6, 2.8, size=(n, 1))
    above = np.arange(f.size)[None, :] > split_idx
    amp_hi = split_f ** (-0.5 * beta) \
        * (safe_f / split_f) ** (-0.5 * beta_hi)
    amp = np.where(multiscale & above, amp_hi, amp)
    amp[:, 0] = 0.0
    z = rng.standard_normal((n, f.size)) + 1j * rng.standard_normal((n, f.size))
    x = np.fft.irfft(z * amp, n=embed, axis=1)[:, :L]
    integrate = rng.random(n) < 0.25
    if integrate.any():
        x[integrate] = np.cumsum(x[integrate], axis=1)
    return _prefix_standardize(x)


def _ou_stochastic_vol(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    switch_rate = np.exp(rng.uniform(np.log(0.001), np.log(0.15), size=(n, 1)))
    switches = rng.random((n, L)) < switch_rate
    switches[:, 0] = rng.random(n) < 0.5
    regime = np.bitwise_and(np.cumsum(switches, axis=1), 1).astype(np.int8)


    slow = rng.random((n, 1)) < 0.5
    phi = np.where(
        slow,
        rng.uniform(0.995, 0.9995, size=(n, 1)),
        rng.uniform(0.90, 0.99, size=(n, 1)),
    )
    mu0 = rng.normal(-2.0, 1.0, size=(n, 1))
    mu1 = rng.normal(2.0, 1.0, size=(n, 1))
    mean = np.where(regime == 0, mu0, mu1)
    seasonal_on = rng.random((n, 1)) < 0.6
    mean += seasonal_on * _seasonal(rng, n, L, k_max=3) \
        * rng.uniform(0.5, 3.0, size=(n, 1))

    log_sigma0 = rng.normal(np.log(0.3), 0.3, size=(n, 1))
    log_sigma1 = rng.normal(np.log(1.5), 0.5, size=(n, 1))
    log_sigma_mean = np.where(regime == 0, log_sigma0, log_sigma1)
    vol_rho = rng.uniform(0.951, 0.995, size=(n, 1))
    vol_eta = rng.uniform(0.03, 0.20, size=(n, 1))
    vol_eps = rng.standard_normal((n, L))
    vol_drive = (
        (1.0 - vol_rho) * log_sigma_mean
        + np.sqrt(1.0 - vol_rho * vol_rho) * vol_eta * vol_eps
    )
    log_vol = np.empty((n, L), dtype=np.float64)
    log_vol[:, 0] = log_sigma_mean[:, 0]
    for i in range(n):
        rho = float(vol_rho[i, 0])
        log_vol[i, 1:] = lfilter(
            [1.0],
            [1.0, -rho],
            vol_drive[i, 1:],
            zi=[rho * log_vol[i, 0]],
        )[0]
    vol = np.exp(np.clip(log_vol, -5.0, 5.0))

    eps = rng.standard_normal((n, L))
    heavy = np.nonzero(rng.random(n) < 0.35)[0]
    if heavy.size:


        eps[heavy] = (
            rng.standard_t(4.0, size=(heavy.size, L)) / np.sqrt(2.0)
        )
    shocks = rng.random((n, L)) < (3.0 / L)
    shock_rows, shock_cols = np.nonzero(shocks)

    eps[shock_rows, shock_cols] += rng.normal(
        0.0, 5.0, size=shock_rows.size
    )

    innovation_scale = np.sqrt(np.maximum(1.0 - phi * phi, 1e-6))
    drive = (1.0 - phi) * mean + innovation_scale * vol * eps
    out = np.empty((n, L), dtype=np.float64)
    out[:, 0] = mean[:, 0] + vol[:, 0] * eps[:, 0]
    for i in range(n):
        p = float(phi[i, 0])
        out[i, 1:] = lfilter(
            [1.0], [1.0, -p], drive[i, 1:], zi=[p * out[i, 0]]
        )[0]

    scale = np.exp(rng.uniform(np.log(0.1), np.log(50.0), size=(n, 1)))
    shift = rng.uniform(-100.0, 100.0, size=(n, 1))
    return out * scale + shift


def _physical_sensors(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    seasonal = _seasonal(rng, n, L, k_max=2)
    smooth = _spectral_gp(rng, n, L)
    fronts = np.cumsum(
        _sparse_jumps(rng, n, L, rate=5.0 / L, scale=1.0), axis=1
    )
    base = (
        seasonal * rng.uniform(0.3, 2.0, size=(n, 1))
        + smooth * rng.uniform(0.2, 1.2, size=(n, 1))
        + fronts * rng.uniform(0.2, 1.0, size=(n, 1))
    )

    kind = rng.integers(0, 4, size=n)
    out = base.copy()

    bounded = kind == 1
    if bounded.any():
        gain = rng.uniform(0.8, 3.5, size=(int(bounded.sum()), 1))
        midpoint = rng.uniform(-0.8, 0.8, size=(int(bounded.sum()), 1))
        out[bounded] = 100.0 / (1.0 + np.exp(-gain * (base[bounded] - midpoint)))

    pressure = kind == 2
    if pressure.any():
        count = int(pressure.sum())


        diffusion = np.exp(
            rng.uniform(np.log(0.03), np.log(0.20), size=(count, 1))
        )
        walk = np.cumsum(
            rng.standard_normal((count, L)) * diffusion, axis=1
        )
        level = rng.uniform(900.0, 1100.0, size=(count, 1))
        out[pressure] = (
            level + walk + 2.0 * fronts[pressure] + 0.5 * seasonal[pressure]
        )

    magnitude = kind == 3
    if magnitude.any():
        count = int(magnitude.sum())
        gusts = (rng.random((count, L)) < (8.0 / L)) \
            * rng.lognormal(0.0, 0.8, size=(count, L))
        power = rng.uniform(1.0, 1.6, size=(count, 1))
        out[magnitude] = np.abs(base[magnitude]) ** power + gusts

    return out


def _seasonal_counts(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    t = np.arange(L, dtype=np.float64)[None, :]
    period = rng.choice(
        _SEASONAL_PERIODS, size=(n, 1), p=_SEASONAL_PROBS
    )
    phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    amp = rng.uniform(0.15, 0.8, size=(n, 1))
    log_rate = amp * np.sin(2.0 * np.pi * t / period + phase)
    second = rng.random((n, 1)) < 0.55
    log_rate += second * (0.5 * amp) * np.sin(
        4.0 * np.pi * t / period + rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    )


    calendar = rng.random((n, 1)) < 0.35
    day_period = rng.choice([24, 48, 96, 144], size=(n, 1))
    day_idx = (np.floor_divide(np.arange(L)[None, :], day_period) % 7).astype(np.int64)
    day_factors = rng.normal(0.0, 0.12, size=(n, 7))
    day_factors[:, 5:] += rng.uniform(-0.8, 0.3, size=(n, 1))
    calendar_effect = np.take_along_axis(day_factors, day_idx, axis=1)
    log_rate += calendar * calendar_effect
    excursion = rng.uniform(-0.5, 0.5, size=(n, 1))
    log_rate += excursion * t / max(L - 1, 1)


    impulses = (
        (rng.random((n, L)) < (2.0 / L))
        * rng.uniform(1.0, 10.0, size=(n, L))
    )
    burst = _ar1_batch(impulses, rng.uniform(0.85, 0.995, size=(n, 1)))
    base = np.exp(rng.uniform(np.log(3.0), np.log(3000.0), size=(n, 1)))
    lam = base * np.exp(np.clip(log_rate, -5.0, 5.0)) * (1.0 + burst)
    np.clip(lam, 0.0, 1.0e7, out=lam)


    overdispersed = rng.random((n, 1)) < 0.5
    shape = rng.uniform(0.5, 4.0, size=(n, 1))
    mixed = lam * rng.gamma(shape, 1.0 / shape, size=(n, L))
    return rng.poisson(np.where(overdispersed, mixed, lam)).astype(np.float64)


def _stable_calendar_counts(
    rng: np.random.Generator, n: int, L: int
) -> np.ndarray:
    """Lower-variance calendar counts with sparse, persistent promotions."""
    t = np.arange(L, dtype=np.float64)[None, :]
    period = rng.choice(
        np.asarray([7.0, 24.0, 48.0, 96.0, 168.0, 336.0]),
        size=(n, 1),
    )
    phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    seasonal = rng.uniform(0.03, 0.20, size=(n, 1)) * np.sin(
        2.0 * np.pi * t / period + phase
    )
    calendar = rng.random((n, 1)) < 0.65
    seasonal += calendar * rng.uniform(0.01, 0.10, size=(n, 1)) * np.sin(
        2.0 * np.pi * t / (7.0 * period)
        + rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    )

    regime = np.cumsum(
        _sparse_jumps(rng, n, L, rate=1.25 / L, scale=0.15), axis=1
    )
    drift = (
        rng.normal(0.0, 0.10, size=(n, 1))
        * t
        / max(L - 1, 1)
    )
    impulses = (
        (rng.random((n, L)) < (1.0 / L))
        * rng.uniform(0.10, 0.55, size=(n, L))
    )
    promotions = _ar1_batch(
        impulses, rng.uniform(0.85, 0.98, size=n)
    )
    base = np.exp(rng.uniform(np.log(2.0), np.log(5000.0), size=(n, 1)))
    lam = base * np.exp(
        np.clip(seasonal + regime + drift + promotions, -4.0, 4.0)
    )
    np.clip(lam, 0.05, 1.0e7, out=lam)

    poisson_rows = rng.random((n, 1)) < 0.4
    poisson = rng.poisson(lam).astype(np.float64)
    smooth = np.maximum(
        np.rint(
            lam
            + rng.normal(size=(n, L))
            * np.sqrt(lam)
            * rng.uniform(0.10, 0.35, size=(n, 1))
        ),
        0.0,
    )
    return np.where(poisson_rows, poisson, smooth)


def _intermittent(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    t = np.arange(L, dtype=np.float64)[None, :]
    base_p = rng.uniform(0.03, 0.35, size=(n, 1))
    period = rng.choice([7.0, 12.0, 24.0, 48.0, 168.0], size=(n, 1))
    season = rng.uniform(0.2, 1.2, size=(n, 1)) * np.sin(
        2.0 * np.pi * t / period + rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    )
    logit = np.log(base_p / (1.0 - base_p)) + season
    p = 1.0 / (1.0 + np.exp(-logit))
    occur = (rng.random((n, L)) < p).astype(np.float64)
    magnitude = np.maximum(
        1.0,
        np.rint(
        rng.gamma(shape=2.0, scale=1.0, size=(n, L))
        * rng.uniform(1.0, 10.0, size=(n, 1))
        * np.exp(0.25 * season)
        ),
    )
    return occur * magnitude


def _pulse_outlier(rng: np.random.Generator, n: int, L: int) -> np.ndarray:


    base = _spectral_gp(rng, n, L) * rng.uniform(0.5, 2.0, size=(n, 1))
    base += _seasonal(rng, n, L, k_max=1) * rng.uniform(0.0, 1.0, size=(n, 1))
    sharp = _sparse_jumps(
        rng, n, L, rate=3.0 / L, scale=rng.uniform(3.0, 8.0, size=n)
    )
    impulses = _sparse_jumps(
        rng, n, L, rate=2.0 / L, scale=rng.uniform(2.0, 7.0, size=n)
    )
    recovery = _ar1_batch(impulses, rng.uniform(0.75, 0.995, size=n))
    series = base + sharp + recovery


    starts = rng.random((n, L)) < (2.0 / L)
    starts[:, 0] = False
    for row in range(n):
        for start in np.nonzero(starts[row])[0]:
            run = int(rng.integers(3, 65))
            end = min(int(start) + run, L)
            series[row, start:end] = series[row, start - 1]
    return series


_CS_CALM_LO = 192
_CS_CALM_HI = 1024
_CS_DYNAMIC_LO = 96
_CS_DYNAMIC_HI = 512


def _conditional_stability(
    rng: np.random.Generator, n: int, L: int
) -> np.ndarray:


    if n <= 0:
        return np.empty((0, L), dtype=np.float64)
    if L <= 0:
        return np.empty((n, 0), dtype=np.float64)

    seasonal = _seasonal(rng, n, L, k_max=2)
    ar = _ar1_batch(
        rng.normal(size=(n, L)) * rng.uniform(0.12, 0.55, size=(n, 1)),
        rng.uniform(0.35, 0.92, size=n),
    )
    smooth = _spectral_gp(rng, n, L)
    walk = np.cumsum(
        rng.normal(size=(n, L)) * rng.uniform(0.025, 0.16, size=(n, 1)),
        axis=1,
    )
    ingredients = (seasonal, ar, smooth, walk)

    kind = rng.integers(0, len(ingredients), size=n)
    calm_kind = rng.integers(0, 4, size=n)
    level = rng.normal(0.0, 2.0, size=n)
    scale = np.exp(rng.uniform(np.log(0.4), np.log(12.0), size=n))
    dynamic_amp = rng.uniform(0.6, 2.2, size=n)
    start_calm = rng.random(n) < 0.65
    out = np.empty((n, L), dtype=np.float64)

    for row in range(n):
        dynamic = ingredients[int(kind[row])][row].copy()
        calibration = dynamic[: min(L, 512)]
        calibration_mean = float(calibration.mean())
        calibration_std = float(calibration.std())
        dynamic -= calibration_mean
        if calibration_std > 1e-12:
            dynamic /= calibration_std

        current = float(level[row])
        calm = bool(start_calm[row])
        pos = 0
        segment_index = 0
        while pos < L:
            if calm:
                seg_len = int(rng.integers(_CS_CALM_LO, _CS_CALM_HI + 1))
            else:
                seg_len = int(
                    rng.integers(_CS_DYNAMIC_LO, _CS_DYNAMIC_HI + 1)
                )


            if segment_index == 0 and L >= 2 * _CS_DYNAMIC_LO:
                seg_len = min(seg_len, L - _CS_DYNAMIC_LO)
            end = min(pos + max(seg_len, 1), L)
            span = end - pos

            if calm:
                mode = int(calm_kind[row])
                if mode == 0:
                    values = np.full(span, current)
                elif mode == 1:


                    drift = rng.normal(0.0, 0.0025, size=span).cumsum()
                    drift += np.linspace(
                        0.0, float(rng.normal(0.0, 0.025)), span
                    )
                    values = current + drift
                elif mode == 2:

                    count_level = max(0.0, float(np.rint(abs(current) * 8.0)))
                    updates = rng.random(span) < 0.025
                    changes = updates * rng.choice(
                        [-1.0, 1.0], size=span
                    )
                    values = np.maximum(
                        count_level + np.cumsum(changes), 0.0
                    )
                else:


                    events = rng.random(span) < 0.012
                    values = events * rng.gamma(1.5, 0.35, size=span)
            else:
                piece = dynamic[pos:end] * float(dynamic_amp[row])
                values = piece - piece[0] + current
                if span > 1 and np.ptp(values) < 1e-10:
                    values = current + np.linspace(0.0, 1.0, span)

            out[row, pos:end] = values
            current = float(values[-1])
            pos = end
            calm = not calm
            segment_index += 1


        row_scale = 1.0 if int(calm_kind[row]) == 2 else float(scale[row])
        out[row] *= row_scale


        if L > 1 and np.ptp(out[row]) < 1e-10:
            out[row, -1] += max(1e-3, 0.01 * row_scale)
    return out


_POS_LINKS = ("offset", "softplus", "exponential")
_POS_LINK_PROBS = np.array([0.45, 0.35, 0.20], dtype=np.float64)


_DOCK_PERIODS = np.array(
    [24.0, 48.0, 96.0, 144.0, 168.0, 288.0, 336.0, 672.0], dtype=np.float64
)
_DOCK_PERIOD_PROBS = np.array(
    [0.18, 0.10, 0.16, 0.10, 0.10, 0.18, 0.09, 0.09], dtype=np.float64
)
_DOCK_PERIOD_PROBS /= _DOCK_PERIOD_PROBS.sum()
_DOCK_PHI = np.array([0.90, 0.97, 0.995, 0.9995], dtype=np.float64)


def _dock_occupancy(rng: np.random.Generator, n: int, L: int) -> np.ndarray:
    """Small-capacity integer occupancy read off a slow latent.

    Station-status feeds are the largest block of live evaluation windows and
    share four properties no other family here produces together: a handful of
    distinct integer levels, 90%+ of consecutive samples exactly equal,
    saturation at empty and full, and churn following a daily cycle. Quantizing
    a slow latent onto a small capacity reproduces all four, and latent movement
    per step measured in whole count units is what sets the flat fraction.
    """
    if n <= 0:
        return np.empty((0, L), dtype=np.float64)
    if L <= 0:
        return np.empty((n, 0), dtype=np.float64)

    t = np.arange(L, dtype=np.float64)[None, :]
    cap = rng.integers(3, 45, size=(n, 1)).astype(np.float64)
    period = rng.choice(_DOCK_PERIODS, size=(n, 1), p=_DOCK_PERIOD_PROBS)
    phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    daily = 0.5 + 0.5 * np.sin(2.0 * np.pi * t / period + phase)
    envelope = rng.uniform(0.0, 0.35, size=(n, 1)) + daily ** rng.uniform(
        1.0, 4.0, size=(n, 1)
    )
    envelope *= 0.75 + 0.25 * np.sin(
        2.0 * np.pi * t / (period * 7.0) + rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    )

    phi_ids = rng.integers(0, _DOCK_PHI.size, size=n)
    latent = _ar1_grouped(
        rng.normal(0.0, 1.0, size=(n, L)) * envelope, phi_ids, _DOCK_PHI
    )
    latent -= latent.mean(axis=1, keepdims=True)
    latent /= np.maximum(latent.std(axis=1, keepdims=True), 1e-9)

    # Swing is the span of capacity the latent covers. Small swings quantize to a
    # two- or three-level series that almost never moves; large swings press
    # against both boundaries and spend long runs exactly empty or exactly full.
    swing = np.exp(rng.uniform(np.log(0.03), np.log(0.6), size=(n, 1))) * cap
    level = rng.uniform(0.15, 0.85, size=(n, 1)) * cap
    seasonal_amp = rng.uniform(0.0, 0.30, size=(n, 1)) * cap
    occupancy = (
        level
        + swing * latent
        + seasonal_amp
        * np.sin(
            2.0 * np.pi * t / period + rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
        )
    )

    counts = np.rint(occupancy)
    span = 2.0 * cap
    reflected = cap - np.abs(np.mod(counts, span) - cap)
    out = np.where(rng.random((n, 1)) < 0.55, np.clip(counts, 0.0, cap), reflected)

    # A minority of stations read as a pure jump process: one integer level held
    # for hours, then a single dock taken or returned.
    rate = np.exp(rng.uniform(np.log(0.002), np.log(0.09), size=(n, 1))) * envelope
    moves = rng.random((n, L)) < rate
    magnitude = 1.0 + (rng.random((n, L)) < 0.18) + (rng.random((n, L)) < 0.04)
    direction = np.where(rng.random((n, L)) < 0.5, -1.0, 1.0)
    walk = np.cumsum(moves * direction * magnitude, axis=1) + np.rint(level)
    out = np.where(
        rng.random((n, 1)) < 0.35,
        cap - np.abs(np.mod(walk, span) - cap),
        out,
    )

    # An exactly constant row carries no forecastable signal and has a degenerate
    # naive-difference scale, so give the flat rows a single step.
    if L > 1:
        constant = np.ptp(out, axis=1, keepdims=True) <= 0.0
        if constant.any():
            step = np.zeros((n, L), dtype=np.float64)
            step[:, L // 2:] = 1.0
            nudged = np.where(out >= cap, out - step, out + step)
            out = np.where(constant, np.clip(nudged, 0.0, cap), out)
    return out


def _piecewise_level(rng: np.random.Generator, n: int, L: int) -> np.ndarray:
    """Levels held exactly constant for long stretches, broken by rare changes.

    The feeds this corpus loses most badly on change value on 0.1% to 8% of
    samples and are bit-identical in between, so the only good point forecast
    repeats the last value untouched and the only good interval is narrow. Their
    naive-difference scale is correspondingly tiny, which is why a forecast that
    drifts even a fraction of one unit is punished so heavily. No other family
    here sustains a level for hundreds of steps, so the model is never shown that
    persistence is the answer.
    """
    if n <= 0:
        return np.empty((0, L), dtype=np.float64)
    if L <= 0:
        return np.empty((n, 0), dtype=np.float64)

    change_rate = np.exp(rng.uniform(np.log(0.0005), np.log(0.09), size=(n, 1)))
    at = rng.random((n, L)) < change_rate
    at[:, 0] = True

    integral = rng.random((n, 1)) < 0.55
    unit_step = np.where(rng.random((n, L)) < 0.5, -1.0, 1.0) * (
        1.0 + (rng.random((n, L)) < 0.22) + (rng.random((n, L)) < 0.05)
    )
    gauss_step = rng.normal(0.0, 1.0, size=(n, L)) * rng.uniform(
        0.3, 2.5, size=(n, 1)
    )
    walk = np.cumsum(at * np.where(integral, unit_step, gauss_step), axis=1)

    cap = rng.integers(3, 40, size=(n, 1)).astype(np.float64)
    offset = walk + np.rint(rng.uniform(0.2, 0.8, size=(n, 1)) * cap)
    span = 2.0 * cap
    integer_out = np.where(
        rng.random((n, 1)) < 0.5,
        np.clip(offset, 0.0, cap),
        cap - np.abs(np.mod(offset, span) - cap),
    )

    scale = np.exp(rng.uniform(np.log(0.5), np.log(2000.0), size=(n, 1)))
    continuous_out = (walk + rng.uniform(-2.0, 2.0, size=(n, 1))) * scale
    out = np.where(integral, integer_out, continuous_out)

    # Some real feeds carry a faint reading error on top of a held level; the rest
    # are bit-exact, and those are the ones the metric punishes hardest to miss.
    noisy = (rng.random((n, 1)) < 0.25) & ~integral
    out = out + noisy * rng.normal(0.0, 1.0, size=(n, L)) * rng.uniform(
        0.002, 0.03, size=(n, 1)
    ) * scale

    if L > 1:
        constant = np.ptp(out, axis=1, keepdims=True) <= 0.0
        if constant.any():
            step = np.zeros((n, L), dtype=np.float64)
            step[:, L // 2:] = 1.0
            alternative = np.where(
                integral,
                np.clip(np.where(out >= cap, out - step, out + step), 0.0, cap),
                out + step * scale * 0.05,
            )
            out = np.where(constant, alternative, out)
    return out


_FLAT_TARGET_HI = (0.55, 0.97)
_FLAT_TARGET_LO = (0.10, 0.60)


def _small_count_link(
    rng: np.random.Generator, block: np.ndarray, rate: float
) -> np.ndarray:
    """Read selected rows off an integer grid chosen to hit a target flat share.

    Live windows are overwhelmingly integer counts, but their resolution varies
    by three orders of magnitude: a bike dock moves over seven levels while a
    national hospitalization count moves over two hundred. Picking the grid by
    capacity therefore coarsens the fine feeds badly. Picking it by target flat
    fraction instead sets the step from the row's own step-size distribution, so
    a slow row lands on a few levels and a noisy row keeps its resolution, and
    the flat fraction — the property that separates the feeds this corpus was
    losing — follows the pool's own spread. The map is monotone, so each family's
    temporal structure survives intact.
    """
    if block.size == 0 or rate <= 0.0:
        return block
    n, L = block.shape
    if L < 4:
        return block
    chosen = np.nonzero(rng.random(n) < rate)[0]
    if chosen.size == 0:
        return block

    x = block[chosen]
    m = chosen.size
    moves = np.sort(np.abs(np.diff(x, axis=1)), axis=1)
    span = np.ptp(x, axis=1, keepdims=True)

    high_band = rng.random((m, 1)) < 0.6
    target = np.where(
        high_band,
        rng.uniform(*_FLAT_TARGET_HI, size=(m, 1)),
        rng.uniform(*_FLAT_TARGET_LO, size=(m, 1)),
    )
    index = np.rint(target[:, 0] * (moves.shape[1] - 1)).astype(np.intp)
    step = 2.0 * moves[np.arange(m), index][:, None]
    # A degenerate row gives no usable step, and rescaling it would collapse it
    # to a constant. Cap the level count so no row becomes a near-continuous grid.
    usable = (step > 0.0) & (span > 0.0)
    step = np.where(usable, step, 1.0)
    step = np.maximum(step, span / 4000.0)

    counts = np.rint((x - x.min(axis=1, keepdims=True)) / step)
    np.clip(counts, 0.0, 1.0e6, out=counts)

    out = block.copy()
    out[chosen] = np.where(usable, counts, x)
    return out


def _domain_observation(
    rng: np.random.Generator,
    block: np.ndarray,
    *,
    already_nonnegative: bool,
    positive_rate: float,
    count_rate: float,
) -> np.ndarray:
    """Move selected signed rows onto positive and optionally integer domains.

    This is an observation layer, not a new family: every transformation is
    monotone, so it preserves the underlying temporal ordering while expanding
    coverage of non-negative and count-valued evaluation windows.
    """
    if already_nonnegative or block.size == 0 or positive_rate <= 0.0:
        return block

    n, _ = block.shape
    chosen = np.nonzero(rng.random(n) < positive_rate)[0]
    if chosen.size == 0:
        return block

    x = block[chosen]
    median = np.median(x, axis=1, keepdims=True)
    scale = np.median(np.abs(x - median), axis=1, keepdims=True) * 1.4826
    scale = np.where(scale < 1e-9, 1.0, scale)
    z = (x - median) / scale

    link = rng.choice(len(_POS_LINKS), size=chosen.size, p=_POS_LINK_PROBS)
    y = np.empty_like(z)

    is_offset = link == 0
    if is_offset.any():
        zo = z[is_offset]
        floor = rng.uniform(0.02, 0.8, size=(int(is_offset.sum()), 1))
        y[is_offset] = zo - zo.min(axis=1, keepdims=True) + floor

    is_softplus = link == 1
    if is_softplus.any():
        zs = z[is_softplus]
        sharp = rng.uniform(0.5, 2.0, size=(int(is_softplus.sum()), 1))
        y[is_softplus] = np.logaddexp(0.0, zs * sharp) / sharp

    is_exp = link == 2
    if is_exp.any():
        ze = z[is_exp]
        vol = rng.uniform(0.15, 0.7, size=(int(is_exp.sum()), 1))
        y[is_exp] = np.exp(np.clip(ze * vol, -12.0, 12.0))

    level = np.exp(rng.uniform(np.log(1.0), np.log(2.0e4), size=(chosen.size, 1)))
    y *= level
    np.clip(y, 0.0, 1.0e6, out=y)

    if count_rate > 0.0:
        counts = rng.random(chosen.size) < count_rate
        if counts.any():
            yc = y[counts]
            target = np.exp(
                rng.uniform(np.log(2.0), np.log(5.0e3), size=(int(counts.sum()), 1))
            )
            peak = np.maximum(yc.max(axis=1, keepdims=True), 1e-12)
            y[counts] = np.rint(yc * (target / peak))

    out = block.copy()
    out[chosen] = y
    return out


def _sanitize(block: np.ndarray) -> np.ndarray:


    x = np.asarray(block, dtype=np.float64)
    np.nan_to_num(x, copy=False, nan=0.0, posinf=1e6, neginf=-1e6)
    if x.ndim == 1:
        peak = float(np.max(np.abs(x)))
        if peak > 1e6:
            x *= 1e6 / peak
    else:
        peak = np.max(np.abs(x), axis=1, keepdims=True)
        scale = np.where(peak > 1e6, 1e6 / np.maximum(peak, 1e-12), 1.0)
        x *= scale
    return x