usdm / usdm.py
Create usdm.py
c1f5495 verified 28 days ago
26.9 kB
#!/usr/bin/env python3
"""
usdm_v2.py — Unified Subspace Diagnostic Model (corrected)
Fixes from v1:
  1. Auto-detected sentinel threshold per signal
  2. Auto-scaled window based on data length
  3. Stability classifier with MODERATE category
  4. Cross-domain comparison uses same encoder and same window
  5. Memory estimator skipped on categorical signals
  6. Failure predictor uses local-scale sentinel threshold
"""

from __future__ import annotations
import os, gzip, csv
from dataclasses import dataclass
from typing import Optional
import numpy as np

DEFAULT_K = 3
STABLE_THRESHOLD = 0.55
MODERATE_THRESHOLD = 0.42
DRIFT_THRESHOLD = 0.30
RANDOM_BASELINE = 3 / 8


# ============================================================
# AUTO-DETECTION
# ============================================================
def auto_sentinel_threshold(values: np.ndarray) -> float:
    """Detect a domain-appropriate sentinel threshold.
    Uses the ratio of the 1st percentile to the median."""
    v = values[~np.isnan(values)]
    if len(v) < 20:
        return 0.0
    med = float(np.median(np.abs(v)))
    if med <= 0:
        return 0.0
    p1 = float(np.percentile(np.abs(v), 1))
    # Sentinel is < 1% of median
    if p1 < med * 0.01:
        return med * 0.01
    # Or bottom 5% if the distribution has a natural gap
    p5 = float(np.percentile(np.abs(v), 5))
    if p1 < p5 * 0.1 and p5 > 0:
        return p1 * 2
    return 0.0


def auto_window(n_total: int, target_windows: int = 20,
                min_window: int = 20, max_window: int = 200) -> int:
    """Choose a window length that produces ~target_windows windows."""
    if n_total < 2 * min_window:
        return max(n_total // 2, 10)
    w = n_total // target_windows
    return int(np.clip(w, min_window, max_window))


def is_categorical(values: np.ndarray, max_unique: int = 10) -> bool:
    """Detect categorical encodings (like A=0,C=1,G=2,T=3)."""
    v = values[~np.isnan(values)]
    if len(v) < 20:
        return False
    return len(np.unique(v)) <= max_unique


# ============================================================
# ENCODERS
# ============================================================
class ChunkSVDEncoder:
    name = "chunk_svd"

    def __init__(self, n_sub: int = 5, sub_len: int = 8, k: int = DEFAULT_K):
        self.n_sub = n_sub
        self.sub_len = sub_len
        self.k = k

    def encode(self, values, window, step):
        if len(values) < window:
            return []
        stride = max(1, (window - self.sub_len) // max(1, self.n_sub - 1))
        out = []
        for start in range(0, len(values) - window + 1, step):
            w = values[start:start + window]
            rows = []
            for i in range(self.n_sub):
                s = i * stride
                if s + self.sub_len > len(w):
                    break
                rows.append(w[s:s + self.sub_len])
            M = np.array(rows)
            if M.shape[0] < 2:
                continue
            M = M - M.mean(axis=0, keepdims=True)
            _, _, Vt = np.linalg.svd(M, full_matrices=False)
            out.append((start, Vt[:self.k].T))
        return out


class TakensEncoder:
    name = "takens"

    def __init__(self, m: int = 8, tau: int = 1, k: int = DEFAULT_K):
        self.m = m
        self.tau = tau
        self.k = k

    def encode(self, values, window, step):
        if len(values) < window:
            return []
        out = []
        for start in range(0, len(values) - window + 1, step):
            w = values[start:start + window].astype(np.float64)
            rows = []
            for t in range(len(w) - (self.m - 1) * self.tau):
                rows.append(w[t + np.arange(self.m) * self.tau])
            M = np.array(rows)
            if M.shape[0] < 2:
                continue
            M = M - M.mean(axis=0, keepdims=True)
            _, _, Vt = np.linalg.svd(M, full_matrices=False)
            out.append((start, Vt[:self.k].T))
        return out


class BinningEncoder:
    name = "binning"

    def __init__(self, n_bins: int = 16, n_blocks: int = 5, seed: int = 42):
        self.n_bins, self.n_blocks, self.seed = n_bins, n_blocks, seed

    def _hv(self, s):
        rng = np.random.default_rng(s)
        return rng.choice([-1.0, 1.0], size=8192).astype(np.float32)

    def encode(self, values, window, step):
        if len(values) < window:
            return []
        smoothed = np.convolve(values, np.ones(10)/10, mode='valid')
        pos = [self._hv(self.seed + 1000 + i) for i in range(self.n_blocks)]
        bins = [self._hv(self.seed + 5000 + b) for b in range(self.n_bins)]
        out = []
        for start in range(0, len(smoothed) - window + 1, step):
            w = smoothed[start:start + window]
            bs = len(w) // self.n_blocks
            if bs < 2:
                continue
            bm = np.array([w[i*bs:(i+1)*bs].mean()
                           for i in range(self.n_blocks)])
            lo, hi = bm.min(), bm.max()
            normed = ((bm - lo) / (hi - lo + 1e-12)
                      if hi > lo else np.full(self.n_blocks, 0.5))
            hv = np.zeros(8192, dtype=np.float32)
            for i in range(self.n_blocks):
                b = min(int(normed[i] * self.n_bins), self.n_bins - 1)
                hv += pos[i] * bins[b]
            out.append((start, np.sign(hv)))
        return out


class RotationEncoder:
    name = "rotation"

    def __init__(self, d_feat: int = 32, seed: int = 42):
        self.d_feat, self.seed = d_feat, seed

    def _hv(self, s):
        rng = np.random.default_rng(s)
        return rng.choice([-1.0, 1.0], size=8192).astype(np.float32)

    def encode(self, values, window, step):
        if len(values) < window:
            return []
        rng = np.random.default_rng(self.seed)
        freqs = rng.normal(0, 1, size=self.d_feat)
        t = np.arange(window) / window
        phi = 2 * np.pi * np.outer(freqs, t)
        out = []
        for start in range(0, len(values) - window + 1, step):
            w = values[start:start + window].astype(np.float64)
            w = w - w.mean()
            z = np.concatenate([np.cos(phi) @ w, np.sin(phi) @ w])
            z = z / (np.sqrt(np.sum(w**2)) + 1e-12)
            hv = np.zeros(8192, dtype=np.float32)
            for i in range(8):
                hv += self._hv(self.seed + 10000 + i) * np.sign(z[i % len(z)])
            out.append((start, np.sign(hv)))
        return out


ENCODERS = {
    "chunk_svd": ChunkSVDEncoder,
    "takens": TakensEncoder,
    "binning": BinningEncoder,
    "rotation": RotationEncoder,
}

POSITIVE_CONTROL = {
    "binning":   {"sine": 1.0000, "walk": 1.0000, "noise": 1.0000, "verdict": "FAIL"},
    "rotation":  {"sine": 1.0000, "walk": 1.0000, "noise": 1.0000, "verdict": "FAIL"},
    "chunk_svd": {"sine": 0.4322, "walk": 0.6057, "noise": 0.1473, "verdict": "PASS"},
    "takens":    {"sine": 0.8170, "walk": 0.9994, "noise": 0.7844, "verdict": "SENSITIVE"},
}


# ============================================================
# ALIGNMENT
# ============================================================
def subspace_align(Ua: np.ndarray, Ub: np.ndarray) -> float:
    """Requires same ambient dimension. Truncates k to min."""
    if Ua.shape[0] != Ub.shape[0]:
        return 0.0
    k = min(Ua.shape[1], Ub.shape[1])
    M = Ua[:, :k].conj().T @ Ub[:, :k]
    return float(np.sum(np.abs(M) ** 2) / k)


def trajectory(items, is_subspace=True):
    if len(items) < 2:
        return np.array([]), np.array([])
    if is_subspace:
        subs = [it[1] for it in items]
    else:
        subs = []
        for _, hv in items:
            n = len(hv) // DEFAULT_K
            M = hv[:n*DEFAULT_K].reshape(n, DEFAULT_K).astype(np.float64)
            _, _, Vt = np.linalg.svd(M, full_matrices=False)
            subs.append(Vt[:DEFAULT_K].T)
    aligns = np.array([subspace_align(subs[i-1], subs[i])
                       for i in range(1, len(subs))])
    starts = np.array([items[i][0] for i in range(1, len(items))])
    return aligns, starts


# ============================================================
# CLASSIFIERS
# ============================================================
def classify_stability(aligns: np.ndarray) -> tuple:
    """STABLE / MODERATE / DRIFT / LOST."""
    if len(aligns) < 2:
        return "UNKNOWN", 0.0
    mean = float(aligns.mean())
    slope = float(np.polyfit(np.arange(len(aligns)), aligns, 1)[0])

    if mean > STABLE_THRESHOLD and slope > -0.01:
        return "STABLE", min(1.0, mean / STABLE_THRESHOLD)
    if mean > MODERATE_THRESHOLD:
        return "MODERATE", 0.7
    if mean > DRIFT_THRESHOLD:
        if slope < -0.02:
            return "DRIFT", min(1.0, abs(slope) / 0.05)
        return "MODERATE", 0.5
    return "LOST", min(1.0, (DRIFT_THRESHOLD - mean) / DRIFT_THRESHOLD)


def detect_changepoints(aligns: np.ndarray, starts: np.ndarray,
                        threshold: float = 0.30) -> list:
    return [{"window": int(i), "start": int(starts[i]),
             "alignment": float(a)}
            for i, a in enumerate(aligns) if a < threshold]


def estimate_memory(values: np.ndarray,
                    min_period: int = 20,
                    max_period: int = 400) -> tuple:
    """Only fires on continuous signals with strong autocorrelation."""
    if is_categorical(values):
        return None, None, False
    if len(values) < 2 * max_period:
        return None, None, False
    v = values - values.mean()
    n = len(v)
    fft = np.fft.fft(v, n=2*n)
    ac = np.real(np.fft.ifft(fft * np.conj(fft)))[:n]
    ac = ac / (ac[0] + 1e-12)
    region = ac[min_period:max_period]
    if len(region) == 0:
        return None, None, False
    peak_idx = int(np.argmax(region)) + min_period
    peak_val = float(ac[peak_idx])
    noise_floor = float(np.mean(np.abs(ac[max_period:min(2*max_period, n)])))
    contrast = peak_val / (noise_floor + 1e-12)
    is_periodic = contrast > 5.0 and peak_val > 0.3
    return peak_idx, contrast, is_periodic


# ============================================================
# MAIN MODEL
# ============================================================
@dataclass
class DiagnoseResult:
    stability: str
    stability_confidence: float
    encoder_used: str
    window_used: int
    sentinel_threshold: float
    alignment_mean: float
    alignment_std: float
    alignment_min: float
    alignment_max: float
    alignment_slope: float
    n_windows: int
    failure_risk: float
    changepoints: list
    is_periodic: bool
    period_estimate: Optional[int]
    period_contrast: Optional[float]
    is_categorical: bool

    def summary(self) -> str:
        lines = [
            f"Encoder:          {self.encoder_used}",
            f"Window:           {self.window_used}",
            f"Sentinel thresh:  {self.sentinel_threshold:.4g}",
            f"Stability:        {self.stability} "
            f"(conf {self.stability_confidence:.2f})",
            f"Alignment mean:   {self.alignment_mean:.4f} "
            f"± {self.alignment_std:.4f}",
            f"Alignment range:  [{self.alignment_min:.4f}, "
            f"{self.alignment_max:.4f}]",
            f"Slope:            {self.alignment_slope:+.5f}",
            f"Windows:          {self.n_windows}",
            f"Failure risk:     {self.failure_risk:.4f}",
            f"Changepoints:     {len(self.changepoints)}",
            f"Categorical:      {self.is_categorical}",
            f"Periodic:         {self.is_periodic}",
        ]
        if self.is_periodic:
            lines.append(f"Period estimate:  {self.period_estimate}")
            lines.append(f"Period contrast:  {self.period_contrast:.2f}")
        return "\n".join(lines)


@dataclass
class CompareResult:
    alignment: float
    same_domain: bool
    interpretation: str

    def summary(self) -> str:
        return (f"Alignment:       {self.alignment:.4f}\n"
                f"Same domain:     {self.same_domain}\n"
                f"Interpretation:  {self.interpretation}")


@dataclass
class QualityReport:
    n_total: int
    n_valid: int
    n_sentinel: int
    n_dark_runs: int
    longest_dark_run: int
    sentinel_threshold: float
    mean: float
    median: float
    std: float
    cv: float

    def summary(self) -> str:
        return "\n".join([
            f"Total readings:   {self.n_total}",
            f"Valid readings:   {self.n_valid}",
            f"Sentinel values:  {self.n_sentinel} "
            f"({100*self.n_sentinel/max(self.n_total,1):.1f}%)",
            f"Sentinel thresh:  {self.sentinel_threshold:.4g}",
            f"Dark runs:        {self.n_dark_runs}",
            f"Longest run:      {self.longest_dark_run}",
            f"Mean:             {self.mean:.4f}",
            f"Median:           {self.median:.4f}",
            f"Std:              {self.std:.4f}",
            f"CV:               {self.cv:.4f}",
        ])


class UnifiedSubspaceDiagnostic:
    def __init__(self, prefer_encoder: Optional[str] = None,
                 window: Optional[int] = None,
                 step: Optional[int] = None):
        self.prefer_encoder = prefer_encoder
        self.window = window
        self.step = step

    def _pick_encoder(self, values: np.ndarray) -> str:
        if self.prefer_encoder is not None:
            return self.prefer_encoder
        if len(values) < 50:
            return "chunk_svd"
        v = values - values.mean()
        denom = np.dot(v, v) + 1e-12
        ac1 = np.dot(v[:-1], v[1:]) / denom
        fft = np.abs(np.fft.rfft(v))
        peak = np.max(fft[1:]) if len(fft) > 1 else 0
        mean_mag = np.mean(fft[1:]) + 1e-12
        periodicity = peak / mean_mag
        if periodicity > 20 and ac1 > 0.5:
            return "takens"
        return "chunk_svd"

    def diagnose(self, values) -> DiagnoseResult:
        values = np.asarray(values, dtype=np.float64)
        sentinel_thresh = auto_sentinel_threshold(values)
        is_cat = is_categorical(values)

        # Filter sentinels for the diagnostic path
        valid_mask = ~np.isnan(values) & (values >= sentinel_thresh)
        valid = values[valid_mask] if sentinel_thresh > 0 else values[~np.isnan(values)]

        # Choose window and step adaptively
        n = len(valid)
        window = self.window if self.window is not None else auto_window(n)
        step = self.step if self.step is not None else max(2, window // 8)

        # Choose encoder
        enc_name = self._pick_encoder(valid)
        encoder = ENCODERS[enc_name]()

        items = encoder.encode(valid, window, step)
        is_sub = enc_name in ("chunk_svd", "takens")
        aligns, starts = trajectory(items, is_sub)

        stability, stab_conf = classify_stability(aligns)
        cps = detect_changepoints(aligns, starts)

        # Failure prediction: use local scale
        if len(valid) == 0:
            failure = 0.0
        else:
            last_val = float(valid[-1])
            if sentinel_thresh > 0 and last_val < sentinel_thresh:
                failure = 1.0
            else:
                risk = 0.0
                if len(aligns) >= 2:
                    mean_a = float(aligns.mean())
                    slope_a = float(np.polyfit(np.arange(len(aligns)),
                                                aligns, 1)[0])
                    if mean_a < STABLE_THRESHOLD:
                        risk += 0.3 * (STABLE_THRESHOLD - mean_a)
                    if slope_a < 0:
                        risk += 0.5 * min(1.0, -slope_a / 0.05)
                    if len(aligns) > 0 and aligns.min() < 0.2:
                        risk += 0.2
                failure = float(min(1.0, risk))

        period, contrast, is_per = estimate_memory(valid)

        return DiagnoseResult(
            stability=stability,
            stability_confidence=stab_conf,
            encoder_used=enc_name,
            window_used=window,
            sentinel_threshold=sentinel_thresh,
            alignment_mean=float(aligns.mean()) if len(aligns) > 0 else 0.0,
            alignment_std=float(aligns.std()) if len(aligns) > 0 else 0.0,
            alignment_min=float(aligns.min()) if len(aligns) > 0 else 0.0,
            alignment_max=float(aligns.max()) if len(aligns) > 0 else 0.0,
            alignment_slope=(float(np.polyfit(np.arange(len(aligns)),
                                              aligns, 1)[0])
                             if len(aligns) > 1 else 0.0),
            n_windows=len(items),
            failure_risk=failure,
            changepoints=cps,
            is_periodic=is_per,
            period_estimate=period,
            period_contrast=contrast,
            is_categorical=is_cat,
        )

    def compare(self, values_a, values_b) -> CompareResult:
        """Compare two signals using the SAME encoder and SAME window."""
        a = np.asarray(values_a, dtype=np.float64)
        b = np.asarray(values_b, dtype=np.float64)
        # Filter NaN
        a = a[~np.isnan(a)]
        b = b[~np.isnan(b)]

        # Use chunk_svd with fixed params — no encoder selection, so
        # both signals land in the same ambient space.
        encoder = ChunkSVDEncoder(n_sub=5, sub_len=8, k=DEFAULT_K)
        window = auto_window(min(len(a), len(b)))
        step = max(2, window // 8)

        ia = encoder.encode(a, window, step)
        ib = encoder.encode(b, window, step)
        if not ia or not ib:
            return CompareResult(0.0, False, "insufficient data")

        # Median subspace from each
        sub_a = self._median(ia)
        sub_b = self._median(ib)
        al = subspace_align(sub_a, sub_b)

        same = al > RANDOM_BASELINE * 1.3
        if same:
            interp = f"shared structure (alignment {al:.3f} > {RANDOM_BASELINE:.3f})"
        else:
            interp = (f"distinct domains "
                      f"(alignment {al:.3f} ≤ {RANDOM_BASELINE:.3f})")
        return CompareResult(al, same, interp)

    @staticmethod
    def _median(items):
        subs = [it[1] for it in items]
        n = len(subs)
        if n <= 3:
            return subs[n // 2]
        sample = subs[::max(1, n // 30)]
        centrality = np.zeros(len(sample))
        for i, s in enumerate(sample):
            centrality[i] = sum(subspace_align(s, t) for t in sample)
        return sample[int(np.argmax(centrality))]

    def quality_report(self, values) -> QualityReport:
        values = np.asarray(values, dtype=np.float64)
        sentinel = auto_sentinel_threshold(values)
        n_total = len(values)
        if n_total == 0:
            return QualityReport(0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0)

        mask = np.isnan(values)
        if sentinel > 0:
            mask = mask | (values < sentinel)
        n_sent = int(mask.sum())
        valid = values[~mask]

        dark_runs = []
        i = 0
        while i < n_total:
            if mask[i]:
                j = i
                while j < n_total and mask[j]:
                    j += 1
                if (j - i) >= 5:
                    dark_runs.append((i, j - 1))
                i = j
            else:
                i += 1
        longest = max((e - s + 1 for s, e in dark_runs), default=0)

        if len(valid) > 0:
            mean = float(valid.mean())
            std = float(valid.std())
            median = float(np.median(valid))
            cv = std / mean if mean > 0 else 0.0
        else:
            mean = std = median = cv = 0.0

        return QualityReport(
            n_total=n_total, n_valid=int(len(valid)),
            n_sentinel=n_sent, n_dark_runs=len(dark_runs),
            longest_dark_run=longest, sentinel_threshold=sentinel,
            mean=mean, median=median, std=std, cv=cv)


# ============================================================
# LOADERS
# ============================================================
def load_ibm_t1(qubit=17, backend="ibm_fez",
                path="ibm_calibration/ibm_calibration_drift.parquet"):
    if not os.path.exists(path):
        return None
    import pandas as pd
    df = pd.read_parquet(path, columns=['backend', 'qubit_a',
                                         'property_family', 'value',
                                         'calibrated_time'])
    sub = df[(df['backend'] == backend) &
             (df['property_family'] == 'T1') &
             (df['qubit_a'].astype('Int64') == qubit)]
    sub = sub.sort_values('calibrated_time')
    return sub['value'].to_numpy(dtype=np.float64)


def load_ftir(path="amber_ftir/Spectral library_Michael Toffolo/ATR/Amber_Lithuania.CSV"):
    if not os.path.exists(path):
        return None
    rows = []
    with open(path, 'r', errors='replace') as f:
        for line in csv.reader(f):
            if len(line) >= 2:
                try:
                    rows.append((float(line[0]), float(line[1])))
                except ValueError:
                    continue
    if not rows:
        return None
    arr = np.array(rows)
    arr = arr[np.argsort(arr[:, 0])]
    return arr[:, 1]


def load_16s(fastq_path="16s_test/test/SRR13005968-test_1.fastq.gz",
             max_reads=50):
    if not os.path.exists(fastq_path):
        return None
    sequences = []
    with gzip.open(fastq_path, 'rt', errors='replace') as f:
        count = 0
        while count < max_reads:
            if not f.readline():
                break
            seq = f.readline().strip()
            f.readline()
            f.readline()
            if seq:
                sequences.append(seq)
                count += 1
    if not sequences:
        return None
    mapping = {'A': 0.0, 'C': 1.0, 'G': 2.0, 'T': 3.0, 'N': 1.5}
    all_vals = []
    for seq in sequences:
        all_vals.extend(mapping.get(c, 1.5) for c in seq[:300])
    return np.array(all_vals)


# ============================================================
# DEMO
# ============================================================
def _demo():
    print("=" * 78)
    print(" UNIFIED SUBSPACE DIAGNOSTIC MODEL v2")
    print(" Auto-adaptive windows, sentinels, and encoders")
    print("=" * 78)

    usdm = UnifiedSubspaceDiagnostic()

    # ========================================================
    # 1. IBM T1
    # ========================================================
    print("\n" + "=" * 78)
    print(" DEMO 1: IBM T1 qubit 17")
    print("=" * 78)

    t1 = load_ibm_t1(qubit=17)
    if t1 is not None:
        print(f"  Loaded {len(t1)} readings")
        qr = usdm.quality_report(t1)
        print(f"\n  QUALITY REPORT")
        for line in qr.summary().split("\n"):
            print(f"    {line}")

        result = usdm.diagnose(t1)
        print(f"\n  DIAGNOSTIC (auto-windowed, auto-sentinel)")
        for line in result.summary().split("\n"):
            print(f"    {line}")

        # Compare against a known-stable qubit
        t1_stable = load_ibm_t1(qubit=0)
        if t1_stable is not None:
            result2 = usdm.diagnose(t1_stable)
            print(f"\n  REFERENCE: qubit 0 (known stable)")
            print(f"    Stability: {result2.stability}")
            print(f"    Alignment: {result2.alignment_mean:.4f} ± "
                  f"{result2.alignment_std:.4f}")
            print(f"    Window:    {result2.window_used}")

    # ========================================================
    # 2. FTIR
    # ========================================================
    print("\n" + "=" * 78)
    print(" DEMO 2: FTIR amber")
    print("=" * 78)

    ftir = load_ftir()
    if ftir is not None:
        print(f"  Loaded {len(ftir)} spectral points")
        result = usdm.diagnose(ftir)
        for line in result.summary().split("\n"):
            print(f"    {line}")

    # ========================================================
    # 3. 16S
    # ========================================================
    print("\n" + "=" * 78)
    print(" DEMO 3: 16S rRNA")
    print("=" * 78)

    seq = load_16s()
    if seq is not None:
        print(f"  Loaded {len(seq)} nucleotide positions")
        result = usdm.diagnose(seq)
        for line in result.summary().split("\n"):
            print(f"    {line}")

    # ========================================================
    # 4. Cross-domain (using same encoder + same window)
    # ========================================================
    print("\n" + "=" * 78)
    print(" DEMO 4: Cross-domain comparisons")
    print("=" * 78)

    if t1 is not None and ftir is not None:
        cmp = usdm.compare(t1, ftir)
        print(f"\n  IBM T1 vs FTIR")
        for line in cmp.summary().split("\n"):
            print(f"    {line}")

    if ftir is not None and seq is not None:
        cmp = usdm.compare(ftir, seq)
        print(f"\n  FTIR vs 16S")
        for line in cmp.summary().split("\n"):
            print(f"    {line}")

    if t1 is not None and seq is not None:
        cmp = usdm.compare(t1, seq)
        print(f"\n  IBM T1 vs 16S")
        for line in cmp.summary().split("\n"):
            print(f"    {line}")

    # ========================================================
    # 5. Synthetic
    # ========================================================
    print("\n" + "=" * 78)
    print(" DEMO 5: Synthetic signals")
    print("=" * 78)

    n = 2400
    t = np.linspace(0, 8 * np.pi, n)
    rng = np.random.default_rng(42)
    signals = {
        "sine": np.sin(t) + 0.05 * rng.normal(0, 1, n),
        "walk": np.cumsum(rng.normal(0, 1, n)) / np.sqrt(n),
        "noise": rng.normal(0, 1, n),
    }

    print(f"\n  {'signal':>8s}  {'encoder':>10s}  {'stability':>10s}  "
          f"{'mean':>8s}  {'std':>8s}  {'window':>7s}")
    print(f"  {'-'*8}  {'-'*10}  {'-'*10}  {'-'*8}  {'-'*8}  {'-'*7}")
    for name, sig in signals.items():
        r = usdm.diagnose(sig)
        print(f"  {name:>8s}  {r.encoder_used:>10s}  {r.stability:>10s}  "
              f"{r.alignment_mean:>8.4f}  {r.alignment_std:>8.4f}  "
              f"{r.window_used:>7d}")

    # ========================================================
    # 6. Periodic
    # ========================================================
    print("\n" + "=" * 78)
    print(" DEMO 6: Periodic signals")
    print("=" * 78)

    print(f"\n  {'period':>7s}  {'detected':>9s}  {'contrast':>9s}  "
          f"{'verdict':>10s}")
    print(f"  {'-'*7}  {'-'*9}  {'-'*9}  {'-'*10}")
    for period in [40, 80, 160, 320]:
        t2 = np.arange(n)
        sig = np.sin(2 * np.pi * t2 / period) + 0.05 * rng.normal(0, 1, n)
        p, c, is_p = estimate_memory(sig, max_period=500)
        detected = p if p else "---"
        contrast = f"{c:.2f}" if c else "---"
        verdict = "PERIODIC" if is_p else "none"
        print(f"  {period:>7d}  {detected!s:>9s}  {contrast:>9s}  "
              f"{verdict:>10s}")

    print("\n  Done.")


if __name__ == "__main__":
    _demo()