"""my-generator — a cascade data generator forked from the genesis base_generator.

A single :class:`Generator` (``cascade.interface.DataGenerator``) mixing twelve
deterministic families by configurable weight:

    ForecastPFN, SineWave, SawTooth, Step, Anomaly, Spikes, OrnsteinUhlenbeck,
    FastSpectral, NonlinearStateSpace, IntermittentCount, GarchRegime,
    MultiscaleSeasonal

Seven lightweight families are vendored from TempoPFN under ``tempo_gen/``.
Five linear-time families are implemented here to cover additional dynamics
without exceeding the trainer's generation deadline:

* ``garch_regime`` — a Markov-switching mean process with GARCH(1,1) conditional
  volatility and Student-t innovations, i.e. genuine volatility *clustering* and
  fat tails (as opposed to OrnsteinUhlenbeck's continuous stochastic-vol
  process). This is the dominant dynamic in financial, network-telemetry, and
  many operational-metric series.
* ``multiscale_seasonal`` — calendar-nested seasonality (a short period nested
  inside ~7x and ~30x multiples of itself, mirroring day-in-week-in-month
  structure), a slowly evolving amplitude envelope, an activity-level
  heteroskedastic noise floor, and occasional outage-style dips. None of the
  existing periodic families (SineWave, SawTooth) compose nested periods with a
  time-varying envelope and level shifts together.
* ``fast_spectral`` approximates smooth/composite kernel draws with random Fourier
  features; ``nonlinear_state_space`` covers interacting latent causal dynamics;
  ``intermittent_count`` covers bursty nonnegative demand and traffic.

Determinism is the load-bearing property: the emitted corpus is a pure function
of ``(seed, n_series)`` only. We seed NumPy and Python ``random`` from ``seed``,
derive every per-generator and per-series sub-seed deterministically, and use a
separate seeded RNG for length-band cropping. The upstream ``hash()``-based seed offset
(PYTHONHASHSEED-salted, not reproducible across processes) is replaced with a
stable ``zlib.crc32`` in the vendored ``abstract_classes.py``. The active
generation path is CPU-only and reproducible. Exact GP/KernelSynth/CauKer
families are intentionally inactive: their cubic covariance factorisations
cannot meet the production deadline at length 4096.
"""

from __future__ import annotations

import json
import os
import sys
import random as _py_random
from collections.abc import Iterator
from pathlib import Path

import numpy as np

from cascade.interface import DataGenerator

# The trainer imports this file by path (importlib.spec_from_file_location), so the
# vendored ``tempo_gen`` package next to it is not on sys.path by default. Add this
# file's own directory so ``import tempo_gen`` resolves however we are loaded.
# (``os``/``sys`` are not on the static-guard blocklist; only ``os.system`` is.)
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

# Vendored TempoPFN generator wrappers (Apache-2.0; see NOTICE).
from tempo_gen.synthetic_generation.anomalies.anomaly_generator_wrapper import (
    AnomalyGeneratorWrapper,
)
from tempo_gen.synthetic_generation.forecast_pfn_prior.forecast_pfn_generator_wrapper import (
    ForecastPFNGeneratorWrapper,
)
from tempo_gen.synthetic_generation.generator_params import (
    AnomalyGeneratorParams,
    ForecastPFNGeneratorParams,
    GeneratorParams,
    OrnsteinUhlenbeckProcessGeneratorParams,
    SawToothGeneratorParams,
    SineWaveGeneratorParams,
    SpikesGeneratorParams,
    StepGeneratorParams,
)
from tempo_gen.synthetic_generation.ornstein_uhlenbeck_process.ou_generator_wrapper import (
    OrnsteinUhlenbeckProcessGeneratorWrapper,
)
from tempo_gen.synthetic_generation.sawtooth.sawtooth_generator_wrapper import (
    SawToothGeneratorWrapper,
)
from tempo_gen.synthetic_generation.sine_waves.sine_wave_generator_wrapper import (
    SineWaveGeneratorWrapper,
)
from tempo_gen.synthetic_generation.spikes.spikes_generator_wrapper import (
    SpikesGeneratorWrapper,
)
from tempo_gen.synthetic_generation.steps.step_generator_wrapper import (
    StepGeneratorWrapper,
)

# Active families must generate all 16,384 production series inside the trainer's
# 30-minute cap. The vendored exact GP, KernelSynth and CauKer implementations do
# an O(L^3) factorisation for every length-4096 series and cannot meet that bound
# (measured at 1.5s, 25s and 21s per series respectively). Linear-time spectral
# and nonlinear state-space priors below retain smooth-kernel and causal-dynamics
# coverage without risking disqualification before training starts.
_DEFAULT_WEIGHTS: dict[str, float] = {
    "forecast_pfn": 0.24,
    "ornstein_uhlenbeck": 0.16,
    "fast_spectral": 0.16,
    "nonlinear_state_space": 0.12,
    "multiscale_seasonal": 0.12,
    "garch_regime": 0.08,
    "intermittent_count": 0.05,
    "sine_waves": 0.025,
    "steps": 0.025,
    "sawtooth": 0.01,
    "anomalies": 0.005,
    "spikes": 0.005,
}


class _CustomBatch:
    """Mimics the vendored wrappers' ``generate_batch`` return shape (``.values``)."""

    __slots__ = ("values",)

    def __init__(self, values: np.ndarray) -> None:
        self.values = values


def _make_custom_wrapper(fn):
    """Adapt a pure-numpy ``fn(rng, length) -> np.ndarray`` into the same
    ``wrapper_cls(params).generate_batch(batch_size=..., seed=...)`` interface
    the vendored TempoPFN wrappers expose, so custom families slot into the
    existing ``_raw_stream``/``_FAMILIES`` machinery unchanged."""

    class _Wrapper:
        def __init__(self, params) -> None:
            self._length = params.length

        def generate_batch(self, *, batch_size: int, seed: int) -> _CustomBatch:
            rows = np.empty((batch_size, self._length), dtype=np.float64)
            for i in range(batch_size):
                rng = np.random.default_rng((seed + i) % _SEED_MOD)
                rows[i] = fn(rng, self._length)
            return _CustomBatch(rows)

    return _Wrapper


def _gen_garch_regime_series(rng: np.random.Generator, length: int) -> np.ndarray:
    """Markov regime-switching mean + GARCH(1,1) conditional volatility +
    Student-t innovations. Produces volatility clustering and fat tails —
    the dominant dynamic in financial/network/operational-metric series, and
    distinct from OrnsteinUhlenbeck's continuous stochastic-vol process."""
    n_regimes = int(rng.integers(2, 4))
    # Sticky transition matrix: high self-persistence, occasional switches.
    persist = rng.uniform(0.95, 0.998, size=n_regimes)
    trans = np.empty((n_regimes, n_regimes))
    for r in range(n_regimes):
        off = (1.0 - persist[r]) / max(1, n_regimes - 1)
        trans[r] = off
        trans[r, r] = persist[r]

    regime_level = rng.normal(0.0, rng.uniform(1.0, 20.0), size=n_regimes)
    regime_slope = rng.normal(0.0, 0.01, size=n_regimes)
    # GARCH(1,1) params per regime, kept stationary (alpha + beta < 1).
    omega = rng.uniform(0.01, 0.5, size=n_regimes)
    alpha = rng.uniform(0.02, 0.2, size=n_regimes)
    beta = rng.uniform(0.5, 0.9, size=n_regimes) * (1.0 - alpha) / 0.9
    df = rng.uniform(3.5, 12.0)  # Student-t degrees of freedom (fat tails)
    phi = rng.uniform(0.0, 0.5)  # short AR(1) memory on the deviation term
    scale = rng.uniform(0.2, 5.0)

    regime = int(rng.integers(0, n_regimes))
    sigma2 = omega[regime] / max(1e-6, 1.0 - alpha[regime] - beta[regime])
    dev = 0.0
    eps_prev = 0.0
    t_innov = rng.standard_t(df, size=length)
    switches = rng.random(length)

    out = np.empty(length, dtype=np.float64)
    for t in range(length):
        cum = np.cumsum(trans[regime])
        regime = int(np.searchsorted(cum, switches[t] * cum[-1]))
        sigma2 = omega[regime] + alpha[regime] * eps_prev * eps_prev + beta[regime] * sigma2
        eps = t_innov[t] * np.sqrt(max(sigma2, 1e-12))
        dev = phi * dev + eps
        out[t] = regime_level[regime] + regime_slope[regime] * t + scale * dev
        eps_prev = eps
    return out


def _gen_multiscale_seasonal_series(rng: np.random.Generator, length: int) -> np.ndarray:
    """Calendar-nested seasonality (short period nested inside ~7x/~30x
    multiples), a slow amplitude envelope, activity-level heteroskedastic
    noise, occasional level shifts, and outage-style dips."""
    from scipy.signal import lfilter

    t = np.arange(length, dtype=np.float64)
    base_period = rng.uniform(4.0, 36.0)
    nest_ratios = (1.0, 7.0, 30.0)
    series = np.full(length, rng.normal(0.0, 5.0))

    for ratio in nest_ratios:
        period = base_period * ratio
        if period >= length * 1.5:
            continue
        amp = rng.uniform(0.3, 3.0) / ratio  # finer scales get relatively less mass
        phase = rng.uniform(0.0, 2.0 * np.pi)
        env_period = period * rng.uniform(4.0, 12.0)
        env_amp = rng.uniform(0.0, 0.6)
        env_phase = rng.uniform(0.0, 2.0 * np.pi)
        envelope = 1.0 + env_amp * np.sin(2.0 * np.pi * t / env_period + env_phase)
        series += amp * envelope * np.sin(2.0 * np.pi * t / period + phase)

    # Slow deterministic trend.
    series += rng.normal(0.0, 0.005) * t

    # A handful of persistent level shifts (structural breaks).
    n_shifts = int(rng.integers(0, 4))
    for _ in range(n_shifts):
        at = int(rng.integers(length // 10, max(length // 10 + 1, 9 * length // 10)))
        series[at:] += rng.normal(0.0, 3.0)

    # Activity-level heteroskedasticity: a slow positive AR(1) factor scales noise.
    activity_phi = rng.uniform(0.95, 0.997)
    activity_innov = rng.normal(0.0, 0.15, size=length)
    log_activity = lfilter([1.0], [1.0, -activity_phi], activity_innov)
    activity = np.exp(np.clip(log_activity - log_activity.mean(), -4.0, 4.0))
    noise_innov = rng.normal(0.0, rng.uniform(0.1, 0.6), size=length)
    ar_noise = lfilter([1.0], [1.0, -rng.uniform(0.0, 0.6)], noise_innov)
    series += activity * ar_noise

    # Occasional outage-style dips: short windows pulled toward a low floor.
    n_outages = int(rng.integers(0, 3))
    floor = series.min() - abs(rng.normal(0.0, 2.0)) - 1.0
    for _ in range(n_outages):
        width = int(rng.integers(2, max(3, length // 40 + 3)))
        start = int(rng.integers(0, max(1, length - width)))
        series[start:start + width] = floor

    return series


def _gen_fast_spectral_series(rng: np.random.Generator, length: int) -> np.ndarray:
    """Random Fourier-feature kernel draw with colored residuals.

    A finite spectral expansion approximates draws from mixtures of periodic,
    locally periodic and smooth stationary kernels in O(KL), avoiding the exact
    GP family's O(L^3) covariance factorisation.
    """
    from scipy.signal import lfilter

    t = np.arange(length, dtype=np.float64)
    n_components = int(rng.integers(3, 10))
    periods = np.exp(rng.uniform(np.log(4.0), np.log(max(8.0, length * 1.5)), n_components))
    amplitudes = rng.lognormal(mean=-0.2, sigma=0.8, size=n_components)
    phases = rng.uniform(0.0, 2.0 * np.pi, n_components)

    x = np.zeros(length, dtype=np.float64)
    for period, amplitude, phase in zip(periods, amplitudes, phases, strict=True):
        # A gentle frequency drift makes the draw locally rather than globally periodic.
        drift = rng.normal(0.0, 0.08) * t / max(1.0, length - 1)
        angle = 2.0 * np.pi * t / period + phase + drift * t / period
        x += amplitude * np.sin(angle)

    # Smooth or rough Matérn-like residual, selected per series.
    pole = rng.uniform(-0.3, 0.97)
    innovations = rng.normal(0.0, rng.uniform(0.05, 0.8), size=length)
    x += lfilter([1.0], [1.0, -pole], innovations)

    # Nonstationary envelope and low-order trend broaden the kernel family.
    envelope = np.exp(rng.uniform(-1.0, 1.0) * (t / max(1.0, length - 1) - 0.5))
    x *= envelope
    u = t / max(1.0, length - 1)
    x += rng.normal(0.0, 2.0) + rng.normal(0.0, 5.0) * u
    if rng.random() < 0.25:
        x += rng.normal(0.0, 3.0) * (u - 0.5) ** 2
    return x


def _gen_nonlinear_state_space_series(rng: np.random.Generator, length: int) -> np.ndarray:
    """Stable nonlinear latent dynamical system with exogenous forcing.

    This is a linear-time substitute for CauKer's expensive sampled SCM-GP:
    several interacting latent nodes evolve through a stable random coupling
    matrix, nonlinear activation, periodic forcing and heteroskedastic shocks.
    """
    latent_dim = int(rng.integers(3, 7))
    coupling = rng.normal(0.0, 1.0, size=(latent_dim, latent_dim))
    row_norm = np.maximum(np.sum(np.abs(coupling), axis=1, keepdims=True), 1e-12)
    coupling = coupling / row_norm * rng.uniform(0.45, 0.85)
    bias = rng.normal(0.0, 0.15, size=latent_dim)
    forcing_weights = rng.normal(0.0, 0.5, size=latent_dim)
    observation = rng.normal(0.0, 1.0, size=latent_dim)
    observation /= max(np.linalg.norm(observation), 1e-12)
    period = rng.uniform(8.0, min(512.0, max(9.0, length / 2)))
    phase = rng.uniform(0.0, 2.0 * np.pi)
    shock_prob = rng.uniform(0.001, 0.02)

    state = rng.normal(0.0, 0.5, size=latent_dim)
    out = np.empty(length, dtype=np.float64)
    for i in range(length):
        forcing = np.sin(2.0 * np.pi * i / period + phase)
        scale = rng.uniform(0.03, 0.15)
        if rng.random() < shock_prob:
            scale *= rng.uniform(4.0, 15.0)
        state = np.tanh(coupling @ state + bias + forcing_weights * forcing)
        state += rng.normal(0.0, scale, size=latent_dim)
        out[i] = observation @ state

    return rng.uniform(0.5, 20.0) * out + rng.normal(0.0, 5.0)


def _gen_intermittent_count_series(rng: np.random.Generator, length: int) -> np.ndarray:
    """Nonnegative bursty count process for demand, traffic and event series."""
    from scipy.signal import lfilter

    t = np.arange(length, dtype=np.float64)
    base = rng.lognormal(mean=rng.uniform(-1.0, 3.0), sigma=0.5)
    period = rng.uniform(4.0, min(336.0, max(5.0, length / 2)))
    seasonal = 1.0 + rng.uniform(0.1, 0.9) * np.sin(
        2.0 * np.pi * t / period + rng.uniform(0.0, 2.0 * np.pi)
    )
    latent = lfilter(
        [1.0],
        [1.0, -rng.uniform(0.7, 0.99)],
        rng.normal(0.0, rng.uniform(0.02, 0.15), size=length),
    )
    rate = np.clip(base * seasonal * np.exp(np.clip(latent, -3.0, 3.0)), 0.001, 1e4)
    counts = rng.poisson(rate).astype(np.float64)

    # Zero inflation and short promotional/incident bursts.
    counts[rng.random(length) < rng.uniform(0.0, 0.35)] = 0.0
    for _ in range(int(rng.integers(0, 5))):
        width = int(rng.integers(1, max(2, min(64, length // 20 + 1))))
        start = int(rng.integers(0, max(1, length - width)))
        counts[start:start + width] *= rng.uniform(2.0, 10.0)
    return counts


# (wrapper class, params class) per family key.
_FAMILIES: dict[str, tuple[type, type]] = {
    "forecast_pfn": (ForecastPFNGeneratorWrapper, ForecastPFNGeneratorParams),
    "ornstein_uhlenbeck": (OrnsteinUhlenbeckProcessGeneratorWrapper, OrnsteinUhlenbeckProcessGeneratorParams),
    "sine_waves": (SineWaveGeneratorWrapper, SineWaveGeneratorParams),
    "steps": (StepGeneratorWrapper, StepGeneratorParams),
    "sawtooth": (SawToothGeneratorWrapper, SawToothGeneratorParams),
    "anomalies": (AnomalyGeneratorWrapper, AnomalyGeneratorParams),
    "spikes": (SpikesGeneratorWrapper, SpikesGeneratorParams),
    "garch_regime": (_make_custom_wrapper(_gen_garch_regime_series), GeneratorParams),
    "multiscale_seasonal": (_make_custom_wrapper(_gen_multiscale_seasonal_series), GeneratorParams),
    "fast_spectral": (_make_custom_wrapper(_gen_fast_spectral_series), GeneratorParams),
    "nonlinear_state_space": (
        _make_custom_wrapper(_gen_nonlinear_state_space_series),
        GeneratorParams,
    ),
    "intermittent_count": (
        _make_custom_wrapper(_gen_intermittent_count_series),
        GeneratorParams,
    ),
}

# Keep per-series sub-seeds inside [0, 2**32) — the anomaly/spikes generators call
# np.random.seed(), which rejects seeds >= 2**32.
_SEED_MOD = 2_000_000_000


class Generator(DataGenerator):
    """Mix of vendored TempoPFN priors, emitted as a deterministic corpus."""

    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._seed = int(seed)
        self._min_len = int(cfg.get("min_length", 64))
        self._max_len = int(cfg.get("max_length", 4096))
        if not (1 <= self._min_len <= self._max_len):
            raise ValueError(f"invalid length band [{self._min_len}, {self._max_len}]")
        # Generators are drawn at this length, then random-cropped into the band.
        self._gen_len = int(cfg.get("generate_length", self._max_len))
        self._gen_len = max(self._gen_len, self._max_len)
        self._batch = max(1, int(cfg.get("batch_size", 256)))

        # Sanitisation knobs. By default we only repair non-finite values and apply
        # a generous absolute clip; we do NOT force unit scale, because varied
        # realistic scales are themselves useful signal for a from-scratch model.
        self._max_abs = float(cfg.get("max_abs_value", 1.0e6))
        self._clip_sigma = float(cfg.get("clip_sigma", 0.0))  # 0 disables sigma clip
        self._standardize = bool(cfg.get("standardize", False))

        weights = cfg.get("weights", _DEFAULT_WEIGHTS)
        # Restrict to known families with positive weight, preserve a fixed order.
        self._weights = {
            k: float(weights[k])
            for k in _FAMILIES
            if k in weights and float(weights[k]) > 0.0
        }
        if not self._weights:
            self._weights = dict(_DEFAULT_WEIGHTS)

        # The active path uses NumPy/SciPy only. Seed its global RNGs for the
        # vendored wrappers, but never mutate Torch's process-global deterministic
        # mode: local scoring imports the generator in the trainer process, and
        # strict mode makes Toto2's CUDA cumsum kernel fail before training.
        np.random.seed(self._seed % 2**31)
        _py_random.seed(self._seed)

    @property
    def name(self) -> str:
        return "tempopfn-mix-plus-garch-multiscale-v1"

    # ── allocation ──────────────────────────────────────────────────────────
    def _allocate(self, n_series: int) -> list[tuple[str, int]]:
        """Split ``n_series`` across families by weight (largest-remainder).

        Pure function of (weights, n_series) — no RNG — so the allocation is
        identical across processes.
        """
        keys = list(self._weights)
        total_w = sum(self._weights[k] for k in keys)
        raw = {k: n_series * self._weights[k] / total_w for k in keys}
        floor = {k: int(np.floor(raw[k])) for k in keys}
        assigned = sum(floor.values())
        remainder = n_series - assigned
        # Hand out the remaining slots to the largest fractional parts, breaking
        # ties by fixed key order.
        order = sorted(keys, key=lambda k: (-(raw[k] - floor[k]), keys.index(k)))
        for i in range(remainder):
            floor[order[i % len(order)]] += 1
        return [(k, floor[k]) for k in keys if floor[k] > 0]

    def _sub_seed(self, *parts: int) -> int:
        """Deterministic child seed in [0, _SEED_MOD) from the master seed."""
        ss = np.random.SeedSequence([self._seed, *parts])
        return int(ss.generate_state(1, dtype=np.uint32)[0]) % _SEED_MOD

    # ── raw draws ───────────────────────────────────────────────────────────
    def _raw_stream(self, family: str, base_seed: int, chunk: int) -> Iterator[np.ndarray]:
        """Yield an unbounded stream of raw full-length series for ``family``.

        Draws are made in contiguous batches of ``chunk`` series with
        non-overlapping per-series seeds, so the stream is a deterministic
        function of ``(base_seed, chunk)``. ``chunk`` is sized to the demand so we
        never generate a full batch to use only a handful of series.
        """
        wrapper_cls, params_cls = _FAMILIES[family]
        params = params_cls(global_seed=base_seed, length=self._gen_len)
        wrapper = wrapper_cls(params)
        chunk = max(1, chunk)
        batch_seed = base_seed
        while True:
            batch = wrapper.generate_batch(batch_size=chunk, seed=batch_seed % _SEED_MOD)
            values = np.asarray(batch.values)
            if values.ndim == 1:
                values = values[None, :]
            elif values.ndim == 3:
                # Multivariate families (CauKer) emit [batch, seq_len, channels].
                # Flatten each channel into its own univariate series so the
                # emitted corpus stays 1-D like every other family.
                values = np.moveaxis(values, 2, 1).reshape(-1, values.shape[1])
            for row in values:
                yield np.ascontiguousarray(row)
            # Advance past this batch's per-series seeds (wrapper uses seed + i).
            batch_seed += chunk

    # ── sanitisation ────────────────────────────────────────────────────────
    def _sanitize(self, arr: np.ndarray, length: int, fallback_rng: np.random.Generator) -> np.ndarray:
        """Return a finite float64 1-D array of exactly ``length`` samples."""
        x = np.asarray(arr, dtype=np.float64).ravel()
        if x.size != length:
            # Defensive: crop/pad to the requested length.
            if x.size > length:
                x = x[:length]
            else:
                x = np.concatenate([x, np.full(length - x.size, x[-1] if x.size else 0.0)])
        # Repair non-finite values, then clip to a trainer-safe magnitude.
        if not np.isfinite(x).all():
            x = np.nan_to_num(x, nan=0.0, posinf=self._max_abs, neginf=-self._max_abs)
        np.clip(x, -self._max_abs, self._max_abs, out=x)

        if self._standardize:
            std = x.std()
            if std > 1e-12:
                x = (x - x.mean()) / std
        if self._clip_sigma > 0.0:
            mu, sd = x.mean(), x.std()
            if sd > 1e-12:
                np.clip(x, mu - self._clip_sigma * sd, mu + self._clip_sigma * sd, out=x)

        if not np.isfinite(x).all():
            # Last-resort deterministic replacement (should not happen post-repair).
            x = fallback_rng.standard_normal(length)
        return np.ascontiguousarray(x, dtype=np.float64)

    # ── main entrypoint ───────────────────────────────────────────────────────
    def generate(self, n_series: int) -> Iterator[np.ndarray]:
        if n_series <= 0:
            return
        # Master RNG drives length-band cropping only — kept separate from every
        # generator's internal RNG so the crop sequence is order-deterministic.
        crop_rng = np.random.default_rng(self._sub_seed(0xC0FFEE))
        fallback_rng = np.random.default_rng(self._sub_seed(0xFA11BACC))

        allocation = self._allocate(n_series)
        emitted = 0
        for fam_idx, (family, count) in enumerate(allocation):
            base_seed = self._sub_seed(fam_idx + 1)
            main_chunk = min(self._batch, count)
            main = self._raw_stream(family, base_seed, main_chunk)
            regen = self._raw_stream(family, (base_seed + _SEED_MOD // 2) % _SEED_MOD, min(self._batch, 16))
            for _ in range(count):
                # Draw the crop window first so master-RNG state advances exactly
                # once per emitted series, regardless of any repair path.
                length = int(crop_rng.integers(self._min_len, self._max_len + 1))
                max_off = self._gen_len - length
                offset = int(crop_rng.integers(0, max_off + 1)) if max_off > 0 else 0

                raw = next(main)
                window = raw[offset:offset + length]
                series = self._sanitize(window, length, fallback_rng)
                if not np.isfinite(series).all():
                    raw2 = next(regen)
                    window = raw2[offset:offset + length]
                    series = self._sanitize(window, length, fallback_rng)
                yield series
                emitted += 1

        # Allocation sums to n_series by construction, but guard the contract.
        if emitted != n_series:  # pragma: no cover
            raise RuntimeError(f"emitted {emitted} series; expected {n_series}")