Skip to content

API

Run spectral estimation, decomposition, and interpretable detection.

Source code in src/telespec/analysis.py 
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def analyze(
    telemetry: Telemetry,
    *,
    false_alarm_probability: float = 0.01,
    bootstrap_iterations: int = 100,
    train_fraction: float = 0.3,
    random_state: int | None = 0,
) -> AnalysisResult:
    """Run spectral estimation, decomposition, and interpretable detection."""
    clean = deduplicate(telemetry)
    periodogram = generalized_lomb_scargle(clean.timestamps, clean.values)
    threshold, maxima = bootstrap_threshold(
        clean.timestamps,
        clean.values,
        periodogram.frequencies,
        false_alarm_probability=false_alarm_probability,
        iterations=bootstrap_iterations,
        random_state=random_state,
    )
    peaks = extract_peaks(periodogram, threshold, maxima)
    dominant_period = peaks[0].period if peaks else None
    segment = _longest_gap_safe_segment(clean)
    decomposition = decompose(segment, period_seconds=dominant_period)
    spectrogram = stft_spectrogram(segment, window=min(128, segment.values.size))
    feature_times, features = band_power_features(
        decomposition.timestamps,
        decomposition.remainder,
    )
    scores = robust_anomaly_scores(
        feature_times,
        features,
        train_fraction=train_fraction,
        false_positive_rate=false_alarm_probability,
    )
    findings: list[Finding] = [
        Finding(
            "periodicity",
            float(clean.timestamps[0]),
            float(clean.timestamps[-1]),
            peak.power,
            1.0 - peak.false_alarm_probability,
            (
                f"Significant {human_period(peak.period)} periodicity "
                f"(bootstrap FAP {peak.false_alarm_probability:.3g})."
            ),
            peak.period,
            peak.frequency,
            {"power": peak.power, "harmonics": list(peak.harmonics)},
        )
        for peak in peaks
    ]
    findings.extend(cusum_changes(decomposition.timestamps, decomposition.remainder))
    findings.extend(anomaly_findings(scores))
    findings.sort(key=lambda finding: (finding.start, finding.kind))
    return AnalysisResult(
        clean,
        periodogram.frequencies,
        periodogram.power,
        threshold,
        peaks,
        findings,
        decomposition.trend,
        decomposition.seasonal,
        decomposition.remainder,
        scores.times,
        scores.scores,
        {
            "false_alarm_probability": false_alarm_probability,
            "bootstrap_iterations": bootstrap_iterations,
            "train_fraction": train_fraction,
            "resampling": "linear interpolation only for STL and sliding-window features",
            "analysis_segment": {
                "start": float(segment.timestamps[0]),
                "end": float(segment.timestamps[-1]),
                "observations": int(segment.values.size),
            },
            "spectrogram": {
                "times": spectrogram.times.tolist(),
                "frequencies": spectrogram.frequencies.tolist(),
                "power": spectrogram.power.tolist(),
            },
        },
    )

Compute a floating-mean generalized Lomb-Scargle periodogram.

The direct least-squares formulation is intentionally transparent and stable for the dataset sizes targeted by telespec.

Source code in src/telespec/spectral.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def generalized_lomb_scargle(
    timestamps: FloatArray,
    values: FloatArray,
    frequencies: FloatArray | None = None,
    *,
    weights: FloatArray | None = None,
) -> Periodogram:
    """Compute a floating-mean generalized Lomb-Scargle periodogram.

    The direct least-squares formulation is intentionally transparent and
    stable for the dataset sizes targeted by telespec.
    """
    t = np.asarray(timestamps, dtype=float)
    y = np.asarray(values, dtype=float)
    if t.shape != y.shape or t.ndim != 1 or t.size < 3:
        raise ValueError("timestamps and values must be equal one-dimensional arrays")
    order = np.argsort(t)
    t = t[order] - np.min(t)
    y = y[order]
    w = np.ones(y.size, dtype=float) if weights is None else np.asarray(weights, dtype=float)[order]
    if np.any(w <= 0) or not np.all(np.isfinite(w)):
        raise ValueError("weights must be finite and positive")
    w = w / np.sum(w)
    freq = frequency_grid(t) if frequencies is None else np.asarray(frequencies, dtype=float)
    if np.any(freq <= 0):
        raise ValueError("frequencies must be positive")
    mean = float(np.dot(w, y))
    chi0 = float(np.dot(w, np.square(y - mean)))
    if chi0 <= np.finfo(float).eps:
        return Periodogram(freq, np.zeros(freq.size, dtype=float))
    root_w = np.sqrt(w)
    power = np.empty(freq.size, dtype=float)
    for index, candidate in enumerate(freq):
        phase = 2.0 * np.pi * candidate * t
        design = np.column_stack((np.cos(phase), np.sin(phase), np.ones(t.size)))
        coefficients, *_ = np.linalg.lstsq(design * root_w[:, None], y * root_w, rcond=None)
        residual = y - design @ coefficients
        chi = float(np.dot(w, np.square(residual)))
        power[index] = np.clip(1.0 - chi / chi0, 0.0, 1.0)
    return Periodogram(freq, power)
Source code in src/telespec/io.py 
136
137
138
139
140
141
142
143
144
145
146
147
def read_telemetry(
    path: str | Path,
    *,
    timestamp: str = "timestamp",
    value: str = "value",
) -> Telemetry:
    suffix = Path(path).suffix.lower()
    if suffix == ".csv":
        return read_csv(path, timestamp=timestamp, value=value)
    if suffix in {".jsonl", ".ndjson"}:
        return read_jsonl(path, timestamp=timestamp, value=value)
    raise ValueError("unsupported format; use .csv, .jsonl, or .ndjson")

Generate heavy-tailed telemetry and apply parameterized effects.

Source code in src/telespec/validation/synthetic.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def generate(
    *,
    points: int = 720,
    interval: float = 60.0,
    baseline: float = 50.0,
    daily_amplitude: float = 5.0,
    noise: float = 1.0,
    jitter: float = 0.0,
    missing: float = 0.0,
    injections: tuple[Injection, ...] = (),
    random_state: int | None = 0,
) -> SyntheticTelemetry:
    """Generate heavy-tailed telemetry and apply parameterized effects."""
    if points < 20 or interval <= 0 or not 0 <= missing < 0.9:
        raise ValueError("invalid generator parameters")
    rng = np.random.default_rng(random_state)
    times = np.arange(points, dtype=float) * interval
    if jitter:
        times += rng.normal(0.0, jitter * interval, points)
        times = np.sort(times)
    values = (
        baseline
        + daily_amplitude * np.sin(2 * np.pi * times / 86400.0)
        + rng.standard_t(df=5, size=points) * noise
    )
    for injection in injections:
        start = max(0, min(points - 1, int(points * injection.start_fraction)))
        end = max(start + 1, min(points, int(points * injection.end_fraction)))
        if injection.kind == "periodicity":
            values[start:end] += injection.magnitude * np.sin(
                2 * np.pi * times[start:end] / injection.period
            )
        elif injection.kind == "level_shift":
            values[start:end] += injection.magnitude
        elif injection.kind == "variance_burst":
            values[start:end] += rng.normal(0, injection.magnitude, end - start)
        elif injection.kind == "spikes":
            count = max(1, (end - start) // 20)
            indexes = rng.choice(np.arange(start, end), size=count, replace=False)
            values[indexes] += injection.magnitude
    keep = rng.random(points) >= missing
    telemetry = Telemetry(
        times[keep],
        values[keep],
        "latency",
        "ms",
        "synthetic",
        sampling_report(times[keep]),
    )
    return SyntheticTelemetry(telemetry, injections)