Skills › Scientific Agent Skills › uncertainty-and-units › _common.py
_common.py
685 lines · skills/uncertainty-and-units/scripts/_common.py at 5f2bc03
1#!/usr/bin/env python3
2"""Shared, standard-library-first helpers for the uncertainty and unit CLIs."""
4from __future__ import annotations
6import ast
7import json
8import math
9import os
10import stat
11import tempfile
12from pathlib import Path
13from typing import Any, Callable, Iterable
16MAX_INPUT_BYTES = 4 * 1024 * 1024
17MAX_REPORT_BYTES = 4 * 1024 * 1024
18MAX_COMPONENTS = 256
19MAX_VARIABLES = 64
20MAX_TRIALS = 5_000_000
21# A useful initial trial count, not a convergence guarantee. JCGM 101:2008
22# 7.9 and 8.2 require numerical stabilization for a validation claim.
23RECOMMENDED_TRIALS = 1_000_000
24DEFAULT_TRIALS = RECOMMENDED_TRIALS
25DEFAULT_SEED = 20_260_726
26MAX_EXPRESSION_CHARS = 2_000
27MAX_EXPRESSION_NODES = 400
28MAX_EXPRESSION_DEPTH = 32
29MAX_LITERAL_EXPONENT = 64
31PINNED_INSTALL = (
32 'uv pip install "pint==0.26.1" "uncertainties==3.2.3" '
33 '"numpy==2.5.3" "scipy==1.18.1"'
34)
36# Standard uncertainty of a Type B component equals its half-width divided by
37# these factors (JCGM 100:2008 4.3.7-4.3.9).
38DISTRIBUTION_DIVISORS: dict[str, float] = {
39 "normal": 1.0,
40 "rectangular": math.sqrt(3.0),
41 "triangular": math.sqrt(6.0),
42 "arcsine": math.sqrt(2.0),
43 "exact": 1.0,
44}
46ALLOWED_FUNCTIONS = (
47 "abs",
48 "acos",
49 "acosh",
50 "asin",
51 "asinh",
52 "atan",
53 "atan2",
54 "atanh",
55 "cos",
56 "cosh",
57 "degrees",
58 "erf",
59 "exp",
60 "expm1",
61 "fabs",
62 "hypot",
63 "log",
64 "log10",
65 "log1p",
66 "radians",
67 "sin",
68 "sinh",
69 "sqrt",
70 "tan",
71 "tanh",
72)
74ALLOWED_CONSTANTS = {"pi": math.pi, "e": math.e, "tau": math.tau}
76_ALLOWED_NODES = (
77 ast.Expression,
78 ast.BinOp,
79 ast.UnaryOp,
80 ast.Call,
81 ast.Name,
82 ast.Constant,
83 ast.Load,
84 ast.Add,
85 ast.Sub,
86 ast.Mult,
87 ast.Div,
88 ast.Pow,
89 ast.USub,
90 ast.UAdd,
91)
94class CliError(ValueError):
95 """An expected command-line validation error."""
98def bounded_int(minimum: int, maximum: int) -> Callable[[str], int]:
99 """Return an argparse converter for a bounded integer."""
101 def convert(value: str) -> int:
102 try:
103 parsed = int(value)
104 except ValueError as exc:
105 raise CliError(f"expected an integer, got {value!r}") from exc
106 if not minimum <= parsed <= maximum:
107 raise CliError(
108 f"expected an integer from {minimum} through {maximum}, got {parsed}"
109 )
110 return parsed
112 return convert
115def finite_float(value: str) -> float:
116 """Parse a finite floating-point value."""
118 try:
119 parsed = float(value)
120 except (TypeError, ValueError) as exc:
121 raise CliError(f"expected a number, got {value!r}") from exc
122 if not math.isfinite(parsed):
123 raise CliError("value must be finite")
124 return parsed
127def non_negative_float(value: str) -> float:
128 """Parse a finite, non-negative floating-point value."""
130 parsed = finite_float(value)
131 if parsed < 0:
132 raise CliError("value must not be negative")
133 return parsed
136def positive_float(value: str) -> float:
137 """Parse a finite, strictly positive floating-point value."""
139 parsed = finite_float(value)
140 if parsed <= 0:
141 raise CliError("value must be greater than zero")
142 return parsed
145def probability(value: str) -> float:
146 """Parse a coverage probability strictly between zero and one."""
148 parsed = finite_float(value)
149 if not 0 < parsed < 1:
150 raise CliError("coverage probability must be strictly between 0 and 1")
151 return parsed
154def as_finite(value: Any, *, label: str) -> float:
155 """Coerce a JSON scalar to a finite float."""
157 if isinstance(value, bool) or not isinstance(value, (int, float)):
158 raise CliError(f"{label} must be a number")
159 number = float(value)
160 if not math.isfinite(number):
161 raise CliError(f"{label} must be finite")
162 return number
165def as_degrees_of_freedom(value: Any, *, label: str) -> float:
166 """Coerce a degrees-of-freedom entry, treating null as infinite."""
168 if value is None:
169 return math.inf
170 if isinstance(value, str) and value.strip().lower() in {"inf", "infinite"}:
171 return math.inf
172 number = as_finite(value, label=label)
173 if number <= 0:
174 raise CliError(f"{label} must be greater than zero")
175 return number
178def normalize_distribution(value: Any, *, label: str) -> str:
179 """Validate a supported probability-density shape."""
181 if value is None:
182 return "normal"
183 if not isinstance(value, str):
184 raise CliError(f"{label} must be a string")
185 name = value.strip().lower()
186 aliases = {"uniform": "rectangular", "u-shaped": "arcsine", "gaussian": "normal"}
187 name = aliases.get(name, name)
188 if name not in DISTRIBUTION_DIVISORS:
189 allowed = ", ".join(sorted(DISTRIBUTION_DIVISORS))
190 raise CliError(f"{label} must be one of: {allowed}")
191 return name
194def checked_input_file(
195 value: str | os.PathLike[str],
196 *,
197 suffixes: Iterable[str],
198 max_bytes: int = MAX_INPUT_BYTES,
199) -> Path:
200 """Return a bounded regular local file, rejecting URLs and symlinks."""
202 raw = os.fspath(value)
203 if "://" in raw:
204 raise CliError("network URLs are not accepted; provide a local file")
205 path = Path(raw).expanduser()
206 if path.is_symlink():
207 raise CliError(f"input must not be a symlink: {path}")
208 try:
209 info = path.stat()
210 except OSError as exc:
211 raise CliError(f"cannot access input file {path}: {exc}") from exc
212 if not stat.S_ISREG(info.st_mode):
213 raise CliError(f"input is not a regular file: {path}")
214 if info.st_size > max_bytes:
215 raise CliError(f"input is {info.st_size} bytes; limit is {max_bytes} bytes")
216 allowed = {suffix.lower() for suffix in suffixes}
217 if path.suffix.lower() not in allowed:
218 raise CliError(f"input suffix must be one of: {', '.join(sorted(allowed))}")
219 return path.resolve()
222def checked_output_file(
223 value: str | os.PathLike[str],
224 *,
225 suffixes: Iterable[str],
226 force: bool = False,
227) -> Path:
228 """Validate an explicit local output without following symlinks."""
230 raw = os.fspath(value)
231 if "://" in raw:
232 raise CliError("network URLs are not accepted as output paths")
233 path = Path(raw).expanduser()
234 if path.name in {"", ".", ".."}:
235 raise CliError("output must name a file")
236 if path.is_symlink():
237 raise CliError(f"output must not be a symlink: {path}")
238 allowed = {suffix.lower() for suffix in suffixes}
239 if path.suffix.lower() not in allowed:
240 raise CliError(f"output suffix must be one of: {', '.join(sorted(allowed))}")
241 parent = path.parent
242 if not parent.exists() or not parent.is_dir() or parent.is_symlink():
243 raise CliError(f"output parent must be an existing regular directory: {parent}")
244 if path.exists():
245 if not path.is_file():
246 raise CliError(f"output exists and is not a regular file: {path}")
247 if not force:
248 raise CliError(f"refusing to overwrite existing output: {path}")
249 return parent.resolve() / path.name
252def atomic_write_bytes(path: Path, payload: bytes, *, force: bool = False) -> None:
253 """Write bytes through a private same-directory temporary file."""
255 destination = checked_output_file(path, suffixes={path.suffix.lower()}, force=force)
256 descriptor, temporary_name = tempfile.mkstemp(
257 prefix=f".{destination.name}.", suffix=".tmp", dir=destination.parent
258 )
259 temporary = Path(temporary_name)
260 try:
261 with os.fdopen(descriptor, "wb") as handle:
262 handle.write(payload)
263 handle.flush()
264 os.fsync(handle.fileno())
265 os.chmod(temporary, 0o600)
266 if destination.exists() and not force:
267 raise CliError(f"refusing to overwrite existing output: {destination}")
268 os.replace(temporary, destination)
269 finally:
270 temporary.unlink(missing_ok=True)
273def _json_bytes(document: Any) -> bytes:
274 """Serialize deterministic strict JSON."""
276 try:
277 payload = (
278 json.dumps(
279 document,
280 indent=2,
281 sort_keys=True,
282 ensure_ascii=False,
283 allow_nan=False,
284 )
285 + "\n"
286 ).encode("utf-8")
287 except (TypeError, ValueError) as exc:
288 raise CliError(f"report is not strict JSON: {exc}") from exc
289 if len(payload) > MAX_REPORT_BYTES:
290 raise CliError(
291 f"report is {len(payload)} bytes; limit is {MAX_REPORT_BYTES} bytes"
292 )
293 return payload
296def emit_json(
297 document: Any,
298 *,
299 output: str | os.PathLike[str] | None = None,
300 force: bool = False,
301) -> None:
302 """Print deterministic JSON or write it atomically."""
304 payload = _json_bytes(document)
305 if output is None:
306 print(payload.decode("utf-8"), end="")
307 return
308 destination = checked_output_file(output, suffixes={".json"}, force=force)
309 atomic_write_bytes(destination, payload, force=force)
312def emit_text(
313 text: str,
314 *,
315 output: str | os.PathLike[str] | None = None,
316 force: bool = False,
317) -> None:
318 """Print text or write it atomically to Markdown."""
320 payload = text.encode("utf-8")
321 if len(payload) > MAX_REPORT_BYTES:
322 raise CliError(
323 f"report is {len(payload)} bytes; limit is {MAX_REPORT_BYTES} bytes"
324 )
325 if output is None:
326 print(text, end="" if text.endswith("\n") else "\n")
327 return
328 destination = checked_output_file(output, suffixes={".md"}, force=force)
329 atomic_write_bytes(destination, payload, force=force)
332def load_json(value: str | os.PathLike[str]) -> Any:
333 """Load bounded strict JSON from a local file."""
335 path = checked_input_file(value, suffixes={".json"})
337 def reject_constant(constant: str) -> None:
338 raise CliError(f"non-standard JSON constant is not allowed: {constant}")
340 def unique_object(pairs: list[tuple[str, Any]]) -> dict[str, Any]:
341 result: dict[str, Any] = {}
342 for key, value in pairs:
343 if key in result:
344 raise CliError(f"duplicate JSON key: {key!r}")
345 result[key] = value
346 return result
348 try:
349 with path.open("r", encoding="utf-8") as handle:
350 return json.load(handle, parse_constant=reject_constant,
351 object_pairs_hook=unique_object)
352 except (OSError, json.JSONDecodeError) as exc:
353 raise CliError(f"cannot read valid JSON from {path.name}: {exc}") from exc
356def read_text_file(
357 value: str | os.PathLike[str], *, suffixes: Iterable[str] = (".py",)
358) -> tuple[str, Path]:
359 """Read a bounded local UTF-8 text file."""
361 path = checked_input_file(value, suffixes=suffixes)
362 try:
363 return path.read_text(encoding="utf-8"), path
364 except (OSError, UnicodeDecodeError) as exc:
365 raise CliError(f"cannot read {path.name} as UTF-8 text: {exc}") from exc
368# --- bounded expression handling -------------------------------------------
369#
370# Measurement models arrive as text. Nothing here compiles or executes that
371# text: the string is parsed to an AST, every node is checked against a
372# whitelist, and the tree is reduced by an explicit walk over the six operators
373# and the named functions below.
376def parse_expression(text: str) -> ast.Expression:
377 """Parse a measurement model into a validated, bounded expression tree."""
379 if not isinstance(text, str) or not text.strip():
380 raise CliError("expression must be a non-empty string")
381 if len(text) > MAX_EXPRESSION_CHARS:
382 raise CliError(
383 f"expression is {len(text)} characters; limit is {MAX_EXPRESSION_CHARS}"
384 )
385 try:
386 tree = ast.parse(text, mode="eval")
387 except (SyntaxError, ValueError) as exc:
388 raise CliError(f"cannot parse expression: {exc}") from exc
390 nodes = list(ast.walk(tree))
391 if len(nodes) > MAX_EXPRESSION_NODES:
392 raise CliError(
393 f"expression has {len(nodes)} nodes; limit is {MAX_EXPRESSION_NODES}"
394 )
395 for node in nodes:
396 if not isinstance(node, _ALLOWED_NODES):
397 raise CliError(
398 f"expression may not contain {type(node).__name__}; allowed syntax is "
399 "names, numbers, + - * / **, and whitelisted functions"
400 )
401 if isinstance(node, ast.Constant):
402 if isinstance(node.value, bool) or not isinstance(node.value, (int, float)):
403 raise CliError("expression constants must be numbers")
404 if not math.isfinite(float(node.value)):
405 raise CliError("expression constants must be finite")
406 if isinstance(node, ast.Name):
407 if node.id.startswith("_"):
408 raise CliError("expression names must not start with an underscore")
409 if isinstance(node, ast.Call):
410 if not isinstance(node.func, ast.Name):
411 raise CliError("only direct calls to whitelisted functions are allowed")
412 if node.func.id not in ALLOWED_FUNCTIONS:
413 allowed = ", ".join(ALLOWED_FUNCTIONS)
414 raise CliError(
415 f"function {node.func.id!r} is not allowed; allowed: {allowed}"
416 )
417 if node.keywords:
418 raise CliError("function calls may not use keyword arguments")
419 arities = {"atan2": (2,), "hypot": (2,), "log": (1, 2)}
420 if len(node.args) not in arities.get(node.func.id, (1,)):
421 raise CliError(f"wrong number of arguments for {node.func.id}")
422 if isinstance(node, ast.BinOp) and isinstance(node.op, ast.Pow):
423 exponent = node.right
424 if isinstance(exponent, ast.Constant) and isinstance(
425 exponent.value, (int, float)
426 ):
427 if abs(float(exponent.value)) > MAX_LITERAL_EXPONENT:
428 raise CliError(
429 f"literal exponents are limited to +/-{MAX_LITERAL_EXPONENT}"
430 )
432 _check_depth(tree.body, 1)
433 return tree
436def _check_depth(node: ast.AST, depth: int) -> None:
437 """Reject deeply nested expressions."""
439 if depth > MAX_EXPRESSION_DEPTH:
440 raise CliError(f"expression nesting exceeds {MAX_EXPRESSION_DEPTH} levels")
441 for child in ast.iter_child_nodes(node):
442 _check_depth(child, depth + 1)
445def expression_variables(tree: ast.Expression) -> list[str]:
446 """Return the free variable names of a parsed expression, in sorted order."""
448 called = {
449 node.func.id
450 for node in ast.walk(tree)
451 if isinstance(node, ast.Call) and isinstance(node.func, ast.Name)
452 }
453 names = {
454 node.id
455 for node in ast.walk(tree)
456 if isinstance(node, ast.Name)
457 and node.id not in called
458 and node.id not in ALLOWED_CONSTANTS
459 }
460 if len(names) > MAX_VARIABLES:
461 raise CliError(f"expression uses more than {MAX_VARIABLES} variables")
462 return sorted(names)
465def reduce_expression(
466 tree: ast.Expression,
467 variables: dict[str, Any],
468 functions: dict[str, Callable[..., Any]],
469) -> Any:
470 """Reduce a validated expression tree against concrete values."""
472 def walk(node: ast.AST) -> Any:
473 if isinstance(node, ast.Expression):
474 return walk(node.body)
475 if isinstance(node, ast.Constant):
476 return node.value
477 if isinstance(node, ast.Name):
478 if node.id in variables:
479 return variables[node.id]
480 if node.id in ALLOWED_CONSTANTS:
481 return ALLOWED_CONSTANTS[node.id]
482 raise CliError(f"no value supplied for variable {node.id!r}")
483 if isinstance(node, ast.UnaryOp):
484 operand = walk(node.operand)
485 return -operand if isinstance(node.op, ast.USub) else +operand
486 if isinstance(node, ast.BinOp):
487 left = walk(node.left)
488 right = walk(node.right)
489 if isinstance(node.op, ast.Add):
490 return left + right
491 if isinstance(node.op, ast.Sub):
492 return left - right
493 if isinstance(node.op, ast.Mult):
494 return left * right
495 if isinstance(node.op, ast.Div):
496 return left / right
497 return left**right
498 if isinstance(node, ast.Call):
499 name = node.func.id # type: ignore[union-attr]
500 handler = functions.get(name)
501 if handler is None:
502 raise CliError(f"function {name!r} is unavailable in this mode")
503 return handler(*[walk(argument) for argument in node.args])
504 raise CliError(f"unsupported expression node {type(node).__name__}")
506 try:
507 return walk(tree)
508 except CliError:
509 raise
510 except ZeroDivisionError as exc:
511 raise CliError("expression divides by zero at the supplied values") from exc
512 except (ArithmeticError, ValueError, TypeError) as exc:
513 raise CliError(f"cannot compute the expression: {exc}") from exc
516def scalar_functions() -> dict[str, Callable[..., Any]]:
517 """Return uncertainty-aware scalar functions with analytic derivatives."""
519 try:
520 from uncertainties import umath, wrap
521 except ImportError as exc:
522 raise CliError(
523 f"the uncertainties package is unavailable; install with `{PINNED_INSTALL}`"
524 ) from exc
525 # abs(UFloat) and umath.fabs are deprecated in uncertainties 3.2.3.
526 # At zero the derivative is undefined; do not invent a zero uncertainty.
527 absolute = wrap(math.fabs, [lambda x: 1.0 if x > 0 else (-1.0 if x < 0 else math.nan)])
528 mapping: dict[str, Callable[..., Any]] = {"abs": absolute, "fabs": absolute}
529 for name in ALLOWED_FUNCTIONS:
530 if name in {"abs", "fabs"}:
531 continue
532 handler = getattr(umath, name, None)
533 if handler is not None:
534 mapping[name] = handler
535 return mapping
538def array_functions() -> dict[str, Callable[..., Any]]:
539 """Return vectorized functions for Monte Carlo sampling."""
541 try:
542 import numpy as np
543 except ImportError as exc:
544 raise CliError(
545 f"NumPy is unavailable; install with `{PINNED_INSTALL}`"
546 ) from exc
547 mapping: dict[str, Callable[..., Any]] = {
548 "abs": np.abs,
549 "acos": np.arccos,
550 "acosh": np.arccosh,
551 "asin": np.arcsin,
552 "asinh": np.arcsinh,
553 "atan": np.arctan,
554 "atan2": np.arctan2,
555 "atanh": np.arctanh,
556 "cos": np.cos,
557 "cosh": np.cosh,
558 "degrees": np.degrees,
559 "exp": np.exp,
560 "expm1": np.expm1,
561 "fabs": np.fabs,
562 "hypot": np.hypot,
563 # np.log's second positional argument is an output buffer, NOT a base.
564 "log": lambda value, base=None: (
565 np.log(value) if base is None else np.log(value) / np.log(base) # audit-units: ignore UNC003 -- MC arrays
566 ),
567 "log10": np.log10,
568 "log1p": np.log1p,
569 "radians": np.radians,
570 "sin": np.sin,
571 "sinh": np.sinh,
572 "sqrt": np.sqrt,
573 "tan": np.tan,
574 "tanh": np.tanh,
575 }
577 def erf(value: Any) -> Any:
578 try:
579 from scipy.special import erf as scipy_erf
580 except ImportError as exc:
581 raise CliError(
582 f"erf needs SciPy; install with `{PINNED_INSTALL}`"
583 ) from exc
584 return scipy_erf(value)
586 mapping["erf"] = erf
587 return mapping
590def json_dof(value: float | None) -> float | None:
591 """Render infinite degrees of freedom as JSON null."""
593 if value is None or not math.isfinite(value):
594 return None
595 return float(value)
598def format_dof(value: float | None) -> str:
599 """Render degrees of freedom for a report table."""
601 if value is None or not math.isfinite(value):
602 return "inf"
603 return f"{value:g}"
606def welch_satterthwaite(
607 combined_uncertainty: float, terms: Iterable[tuple[float, float]]
608) -> float:
609 """Return effective degrees of freedom from |c_i u_i| and nu_i pairs.
611 JCGM 100:2008 equation G.2b. Components with infinite degrees of freedom
612 contribute nothing to the denominator.
613 """
615 if combined_uncertainty <= 0:
616 return math.inf
617 denominator = 0.0
618 for contribution, dof in terms:
619 if not math.isfinite(dof) or dof <= 0:
620 continue
621 denominator += (contribution**4) / dof
622 if denominator <= 0:
623 return math.inf
624 return (combined_uncertainty**4) / denominator
627def coverage_factor(dof: float, coverage_probability: float) -> float:
628 """Return the Student-t coverage factor for the given effective dof."""
630 if not 0 < coverage_probability < 1:
631 raise CliError("coverage probability must be strictly between 0 and 1")
632 tail = 0.5 * (1.0 + coverage_probability)
633 if not math.isfinite(dof):
634 from statistics import NormalDist
636 return float(NormalDist().inv_cdf(tail))
637 try:
638 from scipy.stats import t as student_t
639 except ImportError as exc:
640 raise CliError(
641 "a Student-t coverage factor needs SciPy; install with "
642 f"`{PINNED_INSTALL}`, or supply infinite degrees of freedom"
643 ) from exc
644 return float(student_t.ppf(tail, dof))
647def numerical_tolerance(combined_uncertainty: float, significant_digits: int) -> float:
648 """Return the JCGM 101:2008 clause 8 numerical tolerance for u_c.
650 Writing u_c to `significant_digits` digits as c x 10**exponent, the
651 tolerance is half of the unit in that last retained digit.
652 """
654 if significant_digits not in (1, 2):
655 raise CliError("numerical tolerance is defined for 1 or 2 significant digits")
656 if combined_uncertainty <= 0 or not math.isfinite(combined_uncertainty):
657 raise CliError("combined standard uncertainty must be finite and positive")
658 rounded = float(f"{combined_uncertainty:.{significant_digits}g}")
659 decade = math.log10(rounded) # audit-units: ignore UNC003 -- plain float
660 exponent = math.floor(decade) - (significant_digits - 1)
661 return 0.5 * (10.0**exponent)
664def shortest_coverage_interval(
665 sorted_sample: Any, coverage_probability: float
666) -> tuple[float, float]:
667 """Return the shortest coverage interval of a sorted Monte Carlo sample."""
669 try:
670 import numpy as np
671 except ImportError as exc:
672 raise CliError(
673 f"NumPy is unavailable; install with `{PINNED_INSTALL}`"
674 ) from exc
675 total = int(sorted_sample.size)
676 if total < 2:
677 raise CliError("a coverage interval needs at least two Monte Carlo trials")
678 # JCGM 101:2008 7.7.2 rounds p*M to the nearest integer.
679 inside = int(math.floor(coverage_probability * total + 0.5))
680 inside = min(max(inside, 1), total - 1)
681 lower = sorted_sample[: total - inside]
682 upper = sorted_sample[inside:]
683 widths = upper - lower
684 index = int(np.argmin(widths))
685 return float(lower[index]), float(upper[index])