Skills › Scientific Agent Skills › pathogen-variant-surveillance › lapis_client.py
lapis_client.py
787 lines · skills/pathogen-variant-surveillance/scripts/lapis_client.py at 5f2bc03
1#!/usr/bin/env python3
2# /// script
3# requires-python = ">=3.11"
4# dependencies = []
5# ///
6"""Minimal LAPIS client plus the pure helpers the four CLIs share.
8Standard library only. Network access to the GenSpectrum family of LAPIS
9instances (cov-spectrum.org, genspectrum.org, pathoplexus.org) is required for
10the request functions; every helper below the ``--- pure helpers ---`` mark is
11offline and independently testable.
13The registry and live contracts were reviewed on 2026-10-01. Field names and
14lineage indexes are per-instance. Numeric/date fields support ranges, but date
15analyses require date-typed columns. Actual response versions are checked before
16combining counts. See references/lapis-api.md for deployment-specific behavior.
17"""
18from __future__ import annotations
20import json
21import hashlib
22import math
23import re
24import time
25import urllib.error
26import urllib.parse
27import urllib.request
28from datetime import date, timedelta
29from typing import Any, Iterable, Sequence
31# Instances verified reachable on 2026-10-01. The registry is a convenience,
32# not an authority: any LAPIS deployment works via --base-url, and every script
33# introspects the schema at runtime rather than trusting this table.
34INSTANCES: dict[str, str] = {
35 "sars-cov-2": "https://lapis.cov-spectrum.org/open/v2",
36 "influenza-a": "https://lapis.genspectrum.org/influenza-a",
37 "h1n1pdm": "https://lapis.genspectrum.org/h1n1pdm",
38 "h3n2": "https://lapis.genspectrum.org/h3n2",
39 "h5n1": "https://lapis.genspectrum.org/h5n1",
40 "rsv-a": "https://lapis.pathoplexus.org/rsv-a",
41 "rsv-b": "https://lapis.pathoplexus.org/rsv-b",
42 "hmpv": "https://lapis.pathoplexus.org/hmpv",
43 "measles": "https://lapis.pathoplexus.org/measles",
44 "mpox": "https://lapis.pathoplexus.org/mpox",
45 "west-nile": "https://lapis.pathoplexus.org/west-nile",
46 "dengue": "https://lapis.pathoplexus.org/dengue",
47 "ebola-zaire": "https://lapis.pathoplexus.org/ebola-zaire",
48 "ebola-sudan": "https://lapis.pathoplexus.org/ebola-sudan",
49 "cchf": "https://lapis.pathoplexus.org/cchf",
50}
52# pango-designation is the authority for SARS-CoV-2 lineage names. Fetched at
53# run time, never cached in this repository: withdrawals and redesignations
54# land continuously and a stale copy is worse than no copy.
55PANGO_ALIAS_URL = (
56 "https://raw.githubusercontent.com/cov-lineages/pango-designation/master/"
57 "pango_designation/alias_key.json"
58)
59PANGO_NOTES_URL = (
60 "https://raw.githubusercontent.com/cov-lineages/pango-designation/master/"
61 "lineage_notes.txt"
62)
64USER_AGENT = "scientific-agent-skills-pathogen-variant-surveillance/1.3"
65TIMEOUT = 60
66MAX_ATTEMPTS = 3
67RETRY_STATUS = {429, 500, 502, 503, 504}
69# Collection- and submission-date fields, in the order they should be preferred
70# when an instance declares more than one. Anything not listed is still found by
71# the substring fallback in pick_date_field().
72COLLECTION_DATE_FIELDS = (
73 "date",
74 "sampleCollectionDateRangeLower",
75 "sampleCollectionDate",
76 "collectionDate",
77 "dateCollected",
78)
79SUBMISSION_DATE_FIELDS = (
80 "dateSubmitted",
81 "earliestReleaseDate",
82 "ncbiReleaseDate",
83 "releasedDate",
84 "submittedDate",
85 "dateReleased",
86)
87# Lineage-like columns for instances that declare no lineage index at all.
88LINEAGE_FIELDS = (
89 "pangoLineage",
90 "nextcladePangoLineage",
91 "lineage",
92 "clade",
93 "nextstrainClade",
94 "subtype",
95 "genotype",
96 "serotype",
97)
98RANGE_TYPES = {"date", "int", "float"}
100# Two-sided 95% normal quantile, shared by the interval and the slope CI.
101Z_95 = 1.959964
104class LapisError(RuntimeError):
105 """A LAPIS request failed in a way the caller cannot paper over."""
108class SnapshotChanged(LapisError):
109 """Responses used in one analysis came from different database versions."""
112_DATA_VERSIONS: dict[str, str] = {}
115def record_version(base_url: str, payload: Any, headers: Any, path: str) -> None:
116 """Compare actual response versions; /sample/info alone cannot pin a query."""
117 version = headers.get("lapis-data-version")
118 if isinstance(payload, dict):
119 info = payload.get("info") or {}
120 version = version or info.get("dataVersion") or payload.get("dataVersion")
121 if version is None:
122 if path in {"sample/aggregated", "sample/aminoAcidMutations", "sample/nucleotideMutations"}:
123 raise LapisError("data response has no dataVersion; cannot verify snapshot consistency")
124 return
125 version = str(version)
126 previous = _DATA_VERSIONS.setdefault(base_url, version)
127 if version != previous:
128 raise SnapshotChanged(f"database changed from {previous} to {version}; rerun the entire analysis")
131# --- network ----------------------------------------------------------------
134def _encode(params: dict[str, Any]) -> str:
135 """Encode query parameters, repeating keys for list values."""
136 pairs: list[tuple[str, str]] = []
137 for key, value in params.items():
138 if value is None:
139 continue
140 if isinstance(value, (list, tuple, set)):
141 pairs.extend((key, str(item)) for item in value if item is not None)
142 else:
143 pairs.append((key, str(value)))
144 return urllib.parse.urlencode(pairs)
147def _error_detail(body: bytes) -> str:
148 """Pull the human-readable reason out of a LAPIS error body.
150 Two envelopes are in use: ``{"error": {...}, "info": {...}}`` and a bare
151 RFC 7807 ``{"type": ..., "detail": ...}``. Both carry ``detail``, and on a
152 bad filter key that string enumerates every key the instance accepts.
153 """
154 try:
155 payload = json.loads(body.decode("utf-8", "replace"))
156 except (ValueError, AttributeError):
157 return sanitize(body.decode("utf-8", "replace"), limit=400) if body else ""
158 if isinstance(payload, dict):
159 inner = payload.get("error") if isinstance(payload.get("error"), dict) else payload
160 detail = inner.get("detail") or inner.get("title")
161 if detail:
162 return sanitize(detail, limit=1200)
163 return sanitize(json.dumps(payload), limit=1200)
166def request(base_url: str, path: str, params: dict[str, Any] | None = None) -> Any:
167 """GET a JSON document from a LAPIS instance, retrying transient failures."""
168 url = f"{base_url.rstrip('/')}/{path.lstrip('/')}"
169 if params:
170 query = _encode(params)
171 if query:
172 url = f"{url}?{query}"
173 req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT})
175 last: str = ""
176 for attempt in range(MAX_ATTEMPTS):
177 try:
178 with urllib.request.urlopen(req, timeout=TIMEOUT) as response:
179 payload = json.load(response)
180 record_version(base_url, payload, response.headers, path)
181 return payload
182 except urllib.error.HTTPError as exc:
183 detail = _error_detail(exc.read())
184 last = f"HTTP {exc.code}: {detail}"
185 if exc.code not in RETRY_STATUS:
186 raise LapisError(f"{url}\n {last}") from exc
187 except (urllib.error.URLError, TimeoutError, json.JSONDecodeError) as exc:
188 last = str(exc)
189 if attempt < MAX_ATTEMPTS - 1:
190 time.sleep(1.5 * (attempt + 1))
191 raise LapisError(f"request failed after {MAX_ATTEMPTS} attempts: {url}\n {last}")
194def resolve_base_url(instance: str | None, base_url: str | None) -> str:
195 """Turn ``--instance`` / ``--base-url`` into a URL, or explain the options."""
196 # Each CLI invocation starts a fresh consistency check, even in one interpreter.
197 _DATA_VERSIONS.clear()
198 if base_url:
199 return base_url.rstrip("/")
200 if not instance:
201 raise LapisError("give --instance NAME or --base-url URL")
202 try:
203 return INSTANCES[instance]
204 except KeyError:
205 known = ", ".join(sorted(INSTANCES))
206 raise LapisError(
207 f"unknown instance {instance!r}. Known: {known}. "
208 "Any other LAPIS deployment works via --base-url."
209 ) from None
212_SCHEMA_CACHE: dict[str, dict] = {}
215def describe_instance(base_url: str) -> dict:
216 """Fetch and cache ``/sample/databaseConfig`` for an instance.
218 Returns ``{"name", "openness", "primary_key", "types", "lineage_indexed",
219 "features"}`` where ``types`` maps every declared metadata field to its
220 declared type. Everything downstream chooses field names from this rather
221 than assuming the SARS-CoV-2 vocabulary.
222 """
223 if base_url in _SCHEMA_CACHE:
224 return _SCHEMA_CACHE[base_url]
225 schema = request(base_url, "sample/databaseConfig").get("schema", {})
226 metadata = schema.get("metadata", [])
227 described = {
228 "name": schema.get("instanceName", ""),
229 "openness": schema.get("opennessLevel", ""),
230 "primary_key": schema.get("primaryKey", ""),
231 "types": {m["name"]: m.get("type", "string") for m in metadata},
232 "lineage_indexed": [m["name"] for m in metadata if m.get("generateLineageIndex")],
233 "features": [f.get("name") for f in schema.get("features", [])],
234 }
235 _SCHEMA_CACHE[base_url] = described
236 return described
239def data_version(base_url: str) -> str:
240 """The instance's current data version, for stamping any result you keep."""
241 if base_url in _DATA_VERSIONS:
242 return _DATA_VERSIONS[base_url]
243 return str(request(base_url, "sample/info").get("dataVersion", ""))
246def surveillance_filters(schema: dict, filters: dict) -> dict:
247 """Avoid counting multiple versions or revoked records when those fields exist."""
248 defaults = {}
249 if "versionStatus" in schema["types"]:
250 defaults["versionStatus"] = "LATEST_VERSION"
251 if "isRevocation" in schema["types"]:
252 defaults["isRevocation"] = "false"
253 if "dataUseTerms" in schema["types"]:
254 defaults["dataUseTerms"] = "OPEN"
255 return {**defaults, **filters}
258def aggregated(
259 base_url: str, filters: dict[str, Any], fields: Sequence[str] = ()
260) -> list[dict]:
261 """Grouped counts. ``fields`` is the group-by, not a projection.
263 These helpers retrieve every group and sort client-side for compatibility.
264 Modern deployments also support ordered limit/offset queries.
265 """
266 params = dict(filters)
267 if fields:
268 params["fields"] = ",".join(fields)
269 return request(base_url, "sample/aggregated", params).get("data", [])
272def count(base_url: str, filters: dict[str, Any]) -> int:
273 """Total sequences matching a filter."""
274 rows = aggregated(base_url, filters)
275 return int(rows[0]["count"]) if rows else 0
278def mutations(
279 base_url: str,
280 filters: dict[str, Any],
281 *,
282 amino_acid: bool = True,
283 min_proportion: float = 0.05,
284) -> list[dict]:
285 """Mutations carried by the matching sequences, with per-site proportions.
287 Rows carry ``mutation``, ``count``, ``coverage``, ``proportion``,
288 ``sequenceName`` (the gene or segment), ``position``, ``mutationFrom`` and
289 ``mutationTo``. ``proportion`` is over ``coverage`` — sequences that
290 actually resolved that site — not over every matching sequence.
291 """
292 endpoint = "sample/aminoAcidMutations" if amino_acid else "sample/nucleotideMutations"
293 params = dict(filters)
294 params["minProportion"] = min_proportion
295 return request(base_url, endpoint, params).get("data", [])
298def lineage_definition(base_url: str, column: str) -> dict:
299 """The instance's live lineage tree for a column: parents and aliases.
301 Only served for columns carrying a lineage index; a 400 means the column
302 has none, and there is no ancestry to query.
303 """
304 return request(base_url, f"sample/lineageDefinition/{column}")
307# SHA-256 digests of the fetched bytes, keyed by URL. HTTP ETags are opaque
308# validators, not Git object IDs; retaining the files is required for replay.
309PANGO_BLOBS: dict[str, str] = {}
312def _fetch_text(url: str) -> str:
313 req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT})
314 try:
315 with urllib.request.urlopen(req, timeout=TIMEOUT) as response:
316 content = response.read()
317 PANGO_BLOBS[url] = hashlib.sha256(content).hexdigest()
318 return content.decode("utf-8", "replace")
319 except (urllib.error.URLError, TimeoutError) as exc:
320 raise LapisError(f"could not fetch {url}: {exc}") from exc
323def pango_provenance() -> str:
324 """Content digests of the pango-designation files fetched this run."""
325 return " ".join(
326 f"{url.rsplit('/', 1)[-1]}@sha256:{sha}" for url, sha in sorted(PANGO_BLOBS.items())
327 )
330def fetch_pango_aliases() -> dict[str, Any]:
331 """The live ``alias_key.json`` from pango-designation.
333 Values are a string (the unaliased parent path) for ordinary lineages and a
334 **list of parents** for recombinants. LAPIS does not carry the recombinant
335 parentage — its lineage tree roots every X* lineage — so this file is the
336 only place ``XFG -> [LF.7, LP.8.1.2]`` is recorded.
337 """
338 try:
339 return json.loads(_fetch_text(PANGO_ALIAS_URL))
340 except (ValueError, TypeError) as exc:
341 raise LapisError("invalid Pango alias JSON") from exc
344def fetch_lineage_notes() -> dict[str, dict[str, str]]:
345 """Parse ``lineage_notes.txt`` into ``{name: {status, note}}``.
347 Names prefixed ``*`` have been withdrawn or redesignated. Those entries are
348 the reason a model's remembered lineage facts are not merely stale but
349 wrong, so they are surfaced rather than filtered out.
350 """
351 parsed: dict[str, dict[str, str]] = {}
352 for line in _fetch_text(PANGO_NOTES_URL).splitlines():
353 if not line.strip() or line.startswith("Lineage\t"):
354 continue
355 name, _, note = line.partition("\t")
356 name = name.strip()
357 if not name:
358 continue
359 withdrawn = name.startswith("*")
360 parsed[name.lstrip("*")] = {
361 "status": "withdrawn" if withdrawn else "designated",
362 "note": note.strip(),
363 }
364 return parsed
367# --- pure helpers -----------------------------------------------------------
370def supports_range(schema: dict, field: str) -> bool:
371 """True when ``<field>From`` / ``<field>To`` exist for this field.
373 LAPIS derives range filters from the declared type, so a date recorded as a
374 string (H5N1's ``sampleCollectionDate``) has none.
375 """
376 return schema.get("types", {}).get(field) in RANGE_TYPES
379def range_keys(field: str) -> tuple[str, str]:
380 """The inclusive lower/upper filter keys for a range-capable field."""
381 return f"{field}From", f"{field}To"
384def pick_date_field(schema: dict, role: str = "collection", preferred: str | None = None) -> str:
385 """Choose a range-capable date field for ``collection`` or ``submission``.
387 Raises rather than guessing: a silently wrong date column produces a
388 plausible time series of the wrong thing.
389 """
390 types = schema.get("types", {})
391 if preferred:
392 if preferred not in types:
393 raise LapisError(f"{preferred!r} is not a field on this instance")
394 if types[preferred] != "date":
395 raise LapisError(
396 f"{preferred!r} is typed {types[preferred]!r} on this instance, so LAPIS "
397 f"cannot use it as a date range filter"
398 )
399 return preferred
401 ordered = COLLECTION_DATE_FIELDS if role == "collection" else SUBMISSION_DATE_FIELDS
402 for name in ordered:
403 if types.get(name) == "date":
404 return name
405 needles = ("collect",) if role == "collection" else ("submit", "release")
406 for name, kind in sorted(types.items()):
407 low = name.lower()
408 if kind == "date" and any(n in low for n in needles) and "_seg" not in low:
409 return name
410 raise LapisError(
411 f"no range-capable {role} date field on this instance; "
412 f"date-typed fields are: {sorted(n for n, t in types.items() if t == 'date')}"
413 )
416def looks_like_lineage(name: str) -> bool:
417 """True for field names that plausibly hold a lineage, clade, or genotype.
419 Deliberately loose, because the naming is not standardised across
420 deployments: seasonal influenza splits the call per segment (``cladeHA``,
421 ``cladeNA``), CCHF names it after the segment (``lineage_S``), mpox carries
422 both ``clade`` and ``outbreakLineage``. Per-segment quality columns
423 (``completeness_seg3``) are excluded.
424 """
425 low = name.lower()
426 # "nextcladeQcOverallScore" contains "clade" and is not a lineage call.
427 if any(token in low for token in ("_seg", "qc", "score", "coverage", "version")):
428 return False
429 return any(
430 token in low for token in ("clade", "lineage", "genotype", "serotype", "subtype")
431 )
434def lineage_field_candidates(schema: dict) -> list[tuple[str, bool]]:
435 """Every plausible lineage column as ``(name, indexed)``, best first.
437 Indexed columns come first because only they support descendant queries.
438 Within a tier, the names this skill knows outrank the ones it guessed, and
439 an HA-derived clade outranks NA — for influenza the HA clade is the one
440 that antigenic and vaccine-strain discussion refers to.
441 """
442 types = schema.get("types", {})
443 indexed = set(schema.get("lineage_indexed", []))
444 known = {name: rank for rank, name in enumerate(LINEAGE_FIELDS)}
446 def sort_key(name: str) -> tuple:
447 return (
448 0 if name in indexed else 1,
449 known.get(name, len(LINEAGE_FIELDS)),
450 0 if name.upper().endswith("HA") else 1,
451 name,
452 )
454 names = [n for n in types if looks_like_lineage(n)]
455 return [(n, n in indexed) for n in sorted(names, key=sort_key)]
458def pick_lineage_field(schema: dict, preferred: str | None = None) -> tuple[str, bool]:
459 """Choose the lineage/clade column. Returns ``(field, has_lineage_index)``.
461 The index flag decides whether a trailing ``*`` expands to descendants or
462 matches nothing, so it travels with the field name everywhere.
463 """
464 types = schema.get("types", {})
465 if preferred:
466 if preferred not in types:
467 raise LapisError(f"{preferred!r} is not a field on this instance")
468 return preferred, preferred in set(schema.get("lineage_indexed", []))
469 candidates = lineage_field_candidates(schema)
470 if not candidates:
471 raise LapisError(
472 "no lineage-like column on this instance; pass --lineage-field with one of: "
473 + ", ".join(sorted(types))
474 )
475 return candidates[0]
478def lineage_filter(value: str, has_index: bool, include_sublineages: bool) -> str:
479 """Build the lineage filter value, refusing the query that lies.
481 ``XFG*`` on an indexed column means XFG and every descendant. The same
482 string on an unindexed column is a literal match against a name no sequence
483 carries, and LAPIS answers 0 without complaint.
484 """
485 stripped = value.rstrip("*")
486 wants_children = include_sublineages or value.endswith("*")
487 if not wants_children:
488 return stripped
489 if not has_index:
490 raise LapisError(
491 f"{value!r} asks for sublineages, but this column carries no lineage index, "
492 "so a trailing '*' would match nothing and report 0. Query the exact name, "
493 "or enumerate descendants yourself."
494 )
495 return f"{stripped}*"
498def wilson_interval(successes: int, total: int, z: float = Z_95) -> tuple[float, float]:
499 """Wilson score interval for a binomial proportion.
501 Used instead of the normal approximation because surveillance weeks are
502 routinely small or zero-count, where Wald intervals leave the unit interval
503 and report zero width at p=0.
504 """
505 if total <= 0:
506 return (0.0, 1.0)
507 k, n = float(successes), float(total)
508 denom = n + z * z
509 centre = (k + z * z / 2.0) / denom
510 half = (z / denom) * math.sqrt(k * (n - k) / n + z * z / 4.0)
511 return (max(0.0, centre - half), min(1.0, centre + half))
514def iso_week_start(value: str) -> str | None:
515 """Monday of the ISO week containing an ISO date, or None if unparseable.
517 LAPIS returns nulls and partial dates for sequences whose collection date
518 was never reported; those are counted separately, never silently binned.
519 """
520 if not value or len(value) < 10:
521 return None
522 try:
523 parsed = date.fromisoformat(value[:10])
524 except ValueError:
525 return None
526 return (parsed - timedelta(days=parsed.weekday())).isoformat()
529def bin_weekly(rows: Iterable[dict], date_key: str) -> tuple[dict[str, int], int]:
530 """Sum ``count`` into ISO weeks. Returns ``(weeks, undated)``."""
531 weeks: dict[str, int] = {}
532 undated = 0
533 for row in rows:
534 week = iso_week_start(str(row.get(date_key) or ""))
535 n = int(row.get("count") or 0)
536 if week is None:
537 undated += n
538 continue
539 weeks[week] = weeks.get(week, 0) + n
540 return weeks, undated
543def week_range(start: str, end: str) -> list[str]:
544 """Every ISO week start from ``start`` to ``end`` inclusive, gaps included."""
545 first, last = date.fromisoformat(start), date.fromisoformat(end)
546 out: list[str] = []
547 cursor = first - timedelta(days=first.weekday())
548 while cursor <= last:
549 out.append(cursor.isoformat())
550 cursor += timedelta(days=7)
551 return out
554def flag_low_coverage(weekly_totals: dict[str, int], fraction: float = 0.4) -> dict[str, bool]:
555 """Flag low counts relative to the older half; this is not completeness."""
556 if not weekly_totals:
557 return {}
558 ordered = sorted(weekly_totals)
559 settled = ordered[: max(1, len(ordered) // 2)]
560 counts = sorted(weekly_totals[w] for w in settled)
561 mid = len(counts) // 2
562 median = counts[mid] if len(counts) % 2 else (counts[mid - 1] + counts[mid]) / 2
563 threshold = median * fraction
564 return {week: weekly_totals[week] <= 0 or weekly_totals[week] < threshold for week in ordered}
567def logit_slope(
568 points: Sequence[tuple[float, int, int]],
569 *,
570 min_successes: int = 5,
571 min_nonzero_weeks: int = 3,
572) -> dict[str, float] | None:
573 """Weighted least-squares slope of log-odds against time.
575 ``points`` are ``(t_weeks, successes, total)``. Proportions are shifted by
576 the Haldane-Anscombe 0.5 so that 0 and 1 remain finite, and each point is
577 weighted by the inverse variance of its logit, ``n*p*(1-p)``.
579 This is a **descriptive** slope, not a fitness estimate. It absorbs any
580 change in who is sequencing, where, and how fast they report, and it
581 assumes the composition of the denominator is stable across the window.
583 Returns None unless there is something to fit. The success thresholds are
584 not decoration: with the continuity correction alone, a lineage observed
585 **zero** times in every week still yields p = 0.5/(n+1), which drifts purely
586 with the denominator. A shrinking denominator then manufactures a tight,
587 confident-looking positive slope for a lineage nobody has seen. Requiring
588 real observations is what stops that number from being printed.
589 """
590 usable = [(t, k, n) for t, k, n in points if n > 0]
591 if len(usable) < 3:
592 return None
593 if sum(k for _, k, _ in usable) < min_successes:
594 return None
595 if sum(1 for _, k, _ in usable if k > 0) < min_nonzero_weeks:
596 return None
598 rows = []
599 for t, k, n in usable:
600 p = (k + 0.5) / (n + 1.0)
601 weight = n * p * (1.0 - p)
602 if weight <= 0:
603 continue
604 rows.append((t, math.log(p / (1.0 - p)), weight))
605 if len(rows) < 3:
606 return None
608 sw = sum(w for _, _, w in rows)
609 mean_t = sum(w * t for t, _, w in rows) / sw
610 mean_y = sum(w * y for _, y, w in rows) / sw
611 sxx = sum(w * (t - mean_t) ** 2 for t, _, w in rows)
612 if sxx <= 0:
613 return None
614 sxy = sum(w * (t - mean_t) * (y - mean_y) for t, y, w in rows)
615 slope = sxy / sxx
616 intercept = mean_y - slope * mean_t
618 # Quasi-binomial dispersion, clamped at 1. With inverse-variance weights the
619 # model's own scale is 1, so an estimate below 1 is underdispersion -- almost
620 # always a short series that happens to sit near the line, not evidence the
621 # slope is better determined than binomial sampling allows. Letting it
622 # through would report a CI narrower than the data can support. Above 1 it
623 # is kept, so real overdispersion widens the interval as it should.
624 residual = sum(w * (y - intercept - slope * t) ** 2 for t, y, w in rows)
625 dof = len(rows) - 2
626 dispersion = max(residual / dof, 1.0) if dof > 0 else 1.0
627 stderr = math.sqrt(dispersion / sxx)
628 return {
629 "slope_per_week": slope,
630 "stderr": stderr,
631 "dispersion": dispersion,
632 "ci_low": slope - Z_95 * stderr,
633 "ci_high": slope + Z_95 * stderr,
634 "n_weeks": float(len(rows)),
635 }
638def unalias(name: str, aliases: dict[str, Any]) -> str:
639 """Expand a Pango alias to its full dotted path under A or B.
641 ``PQ.17`` unaliases to ``XDV.1.5.1.1.8.1.17``; applied repeatedly it walks
642 back to ``B.1.1.529...``. Recombinants stop the walk: their alias entry is a
643 list of parents, not a path, so ``XFG.1.1`` expands no further.
644 """
645 head, _, tail = name.partition(".")
646 target = aliases.get(head)
647 if not isinstance(target, str) or not target:
648 return name
649 return f"{target}.{tail}" if tail else target
652def unalias_full(name: str, aliases: dict[str, Any], limit: int = 20) -> str:
653 """Apply :func:`unalias` until it reaches a fixed point."""
654 current = name
655 for _ in range(limit):
656 expanded = unalias(current, aliases)
657 if expanded == current:
658 return current
659 current = expanded
660 return current
663def recombinant_parents(name: str, aliases: dict[str, Any]) -> list[str]:
664 """Parents of a recombinant lineage, or [] for an ordinary one."""
665 head = name.partition(".")[0]
666 target = aliases.get(head)
667 if isinstance(target, list):
668 seen: list[str] = []
669 for parent in target:
670 if parent not in seen:
671 seen.append(parent)
672 return seen
673 return []
676def parent_chain(name: str, definition: dict, limit: int = 60) -> list[str]:
677 """Walk a LAPIS lineage definition from a lineage up to its root."""
678 chain: list[str] = []
679 current = name
680 for _ in range(limit):
681 parents = (definition.get(current) or {}).get("parents") or []
682 if not parents:
683 break
684 current = parents[0]
685 if current in chain:
686 break
687 chain.append(current)
688 return chain
691def children_map(definition: dict) -> dict[str, list[str]]:
692 """Invert a lineage definition into parent -> children.
694 Worth building once and passing around: the SARS-CoV-2 definition holds
695 thousands of entries, so rebuilding it per lineage turns a list of names into
696 quadratic work.
697 """
698 children: dict[str, list[str]] = {}
699 for lineage, entry in definition.items():
700 for parent in (entry or {}).get("parents") or []:
701 children.setdefault(parent, []).append(lineage)
702 return children
705def descendants(
706 name: str, definition: dict, children: dict[str, list[str]] | None = None
707) -> list[str]:
708 """Every lineage whose parent chain passes through ``name``.
710 Pass ``children`` from :func:`children_map` when resolving several names.
711 """
712 if children is None:
713 children = children_map(definition)
714 found: set[str] = set()
715 stack = list(children.get(name, []))
716 while stack:
717 node = stack.pop()
718 if node in found:
719 continue
720 found.add(node)
721 stack.extend(children.get(node, []))
722 return sorted(found)
725def top_values(rows: Sequence[dict], field: str, limit: int) -> list[str]:
726 """The ``limit`` most frequent non-null values of a grouped-by field.
728 Used to discover what is actually circulating before naming anything, so a
729 situation report starts from the data rather than from a remembered list.
730 """
731 ranked = sorted(
732 (r for r in rows if r.get(field) not in (None, "")),
733 key=lambda r: -int(r.get("count") or 0),
734 )
735 return [str(r[field]) for r in ranked[:limit]]
738# Every C0 and C1 control character, tabs and newlines included -- a single
739# newline in a label is enough to forge a row in TSV output.
740_CONTROL_CHARS = re.compile(r"[\x00-\x1f\x7f-\x9f]")
743def sanitize(value: object, limit: int = 400) -> str:
744 """Flatten a value to one printable line for table and TSV rendering.
746 Cell values reach here from a remote instance: lineage labels, field names,
747 and error ``detail`` strings. Two reasons to flatten them.
749 Correctness: a tab or newline inside a label silently corrupts TSV columns.
751 Safety: this output is read by an agent, so remote text is untrusted input.
752 Stripping control characters and collapsing line breaks stops a hostile or
753 malformed instance from forging rows, headers, or anything that reads as
754 structure rather than data. It is not a substitute for pointing
755 ``--base-url`` only at deployments you trust.
757 JSON output deliberately skips this — the encoder escapes control
758 characters already, so values stay faithful for machine consumers.
759 """
760 text = "" if value is None else str(value)
761 text = re.sub(r"\s{2,}", " ", _CONTROL_CHARS.sub(" ", text)).strip()
762 return text if len(text) <= limit else text[: limit - 3] + "..."
765def format_table(rows: Sequence[dict], columns: Sequence[str]) -> str:
766 """Render rows as an aligned plain-text table."""
767 header = list(columns)
768 body = [[sanitize(r.get(c)) for c in header] for r in rows]
769 widths = [
770 max(len(header[i]), *(len(r[i]) for r in body)) if body else len(header[i])
771 for i in range(len(header))
772 ]
773 lines = [" ".join(h.ljust(widths[i]) for i, h in enumerate(header)).rstrip()]
774 for row in body:
775 lines.append(" ".join(cell.ljust(widths[i]) for i, cell in enumerate(row)).rstrip())
776 return "\n".join(lines)
779def emit(rows: Sequence[dict], columns: Sequence[str], fmt: str) -> str:
780 """Render rows as ``table``, ``tsv``, or ``json``."""
781 if fmt == "json":
782 return json.dumps([{c: r.get(c) for c in columns} for r in rows], indent=2)
783 if fmt == "tsv":
784 lines = ["\t".join(columns)]
785 lines += ["\t".join(sanitize(r.get(c)) for c in columns) for r in rows]
786 return "\n".join(lines)
787 return format_table(rows, columns)