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", ) _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, } 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)), "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, ) 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") 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 == "tidal_harmonic": # Preserve the specialist's clean continuous harmonic task. # v17's positive/count observation link would otherwise turn # roughly one fifth of these scarce rows into integer counts. 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", } block = _sanitize( _domain_observation( rng, _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, ), already_nonnegative=preserve_nonnegative, positive_rate=domain_positive_rate, count_rate=domain_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) 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