Skills › Scientific Agent Skills › uncertainty-and-units › propagate_uncertainty.py
propagate_uncertainty.py
701 lines · skills/uncertainty-and-units/scripts/propagate_uncertainty.py at 5f2bc03
1#!/usr/bin/env python3
2"""Propagate uncertainty through a measurement model two ways and compare them.
4The GUM uncertainty framework (JCGM 100:2008) linearizes the model about the
5best estimates. A Monte Carlo run (JCGM 101:2008) propagates the distributions
6themselves. This fixed-trial CLI reports an endpoint comparison inspired by
7JCGM 101 clause 8. It does not implement the adaptive stabilization procedure
8required for a full validation claim, or validate the measurement model.
9"""
11from __future__ import annotations
13import argparse
14import math
15import sys
16from typing import Any
18import _common
19from _common import CliError
22def parse_variable(text: str) -> dict[str, Any]:
23 """Parse `name=value,u[,distribution[,dof]]` into a variable record."""
25 if "=" not in text:
26 raise CliError(
27 f"variable {text!r} must look like name=value,standard_uncertainty"
28 )
29 name, _, payload = text.partition("=")
30 name = name.strip()
31 fields = [item.strip() for item in payload.split(",")]
32 if len(fields) < 2 or len(fields) > 4:
33 raise CliError(
34 f"variable {name!r} takes value,standard_uncertainty[,distribution[,dof]]"
35 )
36 record: dict[str, Any] = {
37 "name": name,
38 "value": _common.finite_float(fields[0]),
39 "standard_uncertainty": _common.non_negative_float(fields[1]),
40 }
41 if len(fields) >= 3 and fields[2]:
42 record["distribution"] = fields[2]
43 if len(fields) == 4 and fields[3]:
44 if fields[3].lower() in {"inf", "infinite"}:
45 record["dof"] = None
46 else:
47 record["dof"] = _common.finite_float(fields[3])
48 return record
51def parse_correlation(text: str) -> tuple[str, str, float]:
52 """Parse `a,b=r` into a correlation coefficient between two inputs."""
54 pair, _, value = text.partition("=")
55 names = [item.strip() for item in pair.split(",")]
56 if len(names) != 2 or not all(names) or not value.strip():
57 raise CliError(f"correlation {text!r} must look like name_a,name_b=0.4")
58 coefficient = _common.finite_float(value)
59 if not -1.0 <= coefficient <= 1.0:
60 raise CliError("correlation coefficients must lie in [-1, 1]")
61 if names[0] == names[1]:
62 raise CliError("a variable cannot be correlated with itself")
63 return names[0], names[1], coefficient
66def normalize_variables(records: list[Any]) -> list[dict[str, Any]]:
67 """Validate variable records from either the CLI or a JSON spec."""
69 if not records:
70 raise CliError("at least one input variable is required")
71 if len(records) > _common.MAX_VARIABLES:
72 raise CliError(f"at most {_common.MAX_VARIABLES} variables are supported")
73 seen: set[str] = set()
74 normalized: list[dict[str, Any]] = []
75 for record in records:
76 if not isinstance(record, dict):
77 raise CliError("each variable must be a JSON object")
78 name = record.get("name")
79 if not isinstance(name, str) or not name.isidentifier():
80 raise CliError(f"variable name {name!r} must be a Python identifier")
81 if name.startswith("_"):
82 raise CliError("variable names must not start with an underscore")
83 if name in _common.ALLOWED_CONSTANTS or name in _common.ALLOWED_FUNCTIONS:
84 raise CliError(f"variable name {name!r} shadows a built-in symbol")
85 if name in seen:
86 raise CliError(f"variable {name!r} is defined more than once")
87 seen.add(name)
88 uncertainty = record.get("standard_uncertainty", record.get("u"))
89 entry = {
90 "name": name,
91 "value": _common.as_finite(record.get("value"), label=f"{name}.value"),
92 "standard_uncertainty": _common.as_finite(
93 uncertainty, label=f"{name}.standard_uncertainty"
94 ),
95 "distribution": _common.normalize_distribution(
96 record.get("distribution"), label=f"{name}.distribution"
97 ),
98 "dof": _common.as_degrees_of_freedom(
99 record.get("dof"), label=f"{name}.dof"
100 ),
101 "unit": record.get("unit") if isinstance(record.get("unit"), str) else None,
102 }
103 if entry["standard_uncertainty"] < 0:
104 raise CliError(f"{name}: standard uncertainty must not be negative")
105 if entry["distribution"] == "exact" and entry["standard_uncertainty"] != 0:
106 raise CliError(f"{name}: an exact input must have zero uncertainty")
107 normalized.append(entry)
108 return normalized
111def gum_framework(
112 tree: Any,
113 variables: list[dict[str, Any]],
114 correlations: dict[tuple[str, str], float],
115 coverage_probability: float,
116) -> dict[str, Any]:
117 """Evaluate the model and propagate uncertainty by first-order expansion."""
119 try:
120 from uncertainties import ufloat
121 except ImportError as exc:
122 raise CliError(
123 "the uncertainties package is unavailable; install with "
124 f"`{_common.PINNED_INSTALL}`"
125 ) from exc
127 if correlations:
128 correlation_factor(variables, correlations)
129 if any(math.isfinite(entry["dof"]) for entry in variables):
130 raise CliError(
131 "correlated inputs with finite degrees of freedom need a justified "
132 "joint uncertainty model; this CLI cannot apply independent-input "
133 "Welch-Satterthwaite to them"
134 )
135 handles = {
136 entry["name"]: ufloat(entry["value"], entry["standard_uncertainty"])
137 for entry in variables
138 }
139 result = _common.reduce_expression(tree, handles, _common.scalar_functions())
140 value = float(getattr(result, "nominal_value", result))
141 derivatives = getattr(result, "derivatives", {})
143 inputs: list[dict[str, Any]] = []
144 for entry in variables:
145 sensitivity = float(derivatives.get(handles[entry["name"]], 0.0))
146 contribution = sensitivity * entry["standard_uncertainty"]
147 inputs.append({**entry, "sensitivity": sensitivity, "contribution": contribution})
149 independent = sum(item["contribution"] ** 2 for item in inputs)
150 by_name = {item["name"]: item for item in inputs}
151 covariance_term = 0.0
152 for (first, second), coefficient in correlations.items():
153 left = by_name[first]
154 right = by_name[second]
155 covariance_term += (
156 2.0
157 * left["sensitivity"]
158 * right["sensitivity"]
159 * coefficient
160 * left["standard_uncertainty"]
161 * right["standard_uncertainty"]
162 )
163 variance = independent + covariance_term
164 if not all(math.isfinite(number) for number in (value, independent, covariance_term)):
165 raise CliError("model value or linearized uncertainty is not finite")
166 if variance < -1e-12 * independent:
167 raise CliError(
168 "the supplied correlations give a negative combined variance; "
169 "check the correlation matrix for consistency"
170 )
171 combined = math.sqrt(max(variance, 0.0)) # audit-units: ignore UNC003 -- plain float
173 # Budget percentages are taken against the sum of squared contributions,
174 # not against u_c**2. With correlated inputs the covariance term can make
175 # u_c**2 smaller than that sum, and percentages of it would exceed 100.
176 for item in inputs:
177 item["variance_fraction"] = (
178 (item["contribution"] ** 2) / independent if independent > 0 else 0.0
179 )
181 effective_dof = _common.welch_satterthwaite(
182 combined, [(abs(item["contribution"]), item["dof"]) for item in inputs]
183 )
184 factor = _common.coverage_factor(effective_dof, coverage_probability)
185 expanded = factor * combined
186 for item in inputs:
187 item["dof"] = _common.json_dof(item["dof"])
188 return {
189 "value": value,
190 "combined_standard_uncertainty": combined,
191 "independent_variance": independent,
192 "covariance_term": covariance_term,
193 "effective_degrees_of_freedom": (
194 None if not math.isfinite(effective_dof) else effective_dof
195 ),
196 "coverage_probability": coverage_probability,
197 "coverage_factor": factor,
198 "expanded_uncertainty": expanded,
199 "coverage_interval": [value - expanded, value + expanded],
200 "inputs": inputs,
201 }
204def _draw(rng: Any, entry: dict[str, Any], trials: int) -> Any:
205 """Sample one input from its assigned probability density."""
207 import numpy as np
209 value = entry["value"]
210 uncertainty = entry["standard_uncertainty"]
211 if uncertainty == 0:
212 return np.full(trials, value)
213 shape = entry["distribution"]
214 if shape == "normal":
215 return rng.normal(value, uncertainty, trials)
216 half_width = uncertainty * _common.DISTRIBUTION_DIVISORS[shape]
217 if shape == "rectangular":
218 return rng.uniform(value - half_width, value + half_width, trials)
219 if shape == "triangular":
220 return rng.triangular(value - half_width, value, value + half_width, trials)
221 if shape == "arcsine":
222 # audit-units: ignore UNC003 -- ordinary array of draws, not a ufloat
223 return value + half_width * np.cos(rng.uniform(0.0, math.pi, trials))
224 raise CliError(f"cannot sample distribution {shape!r}")
227def correlation_factor(
228 variables: list[dict[str, Any]],
229 correlations: dict[tuple[str, str], float],
230) -> tuple[list[str], Any]:
231 """Validate a correlation matrix and factor it, including singular PSD cases."""
233 import numpy as np
235 involved = sorted({name for pair in correlations for name in pair})
236 by_name = {entry["name"]: entry for entry in variables}
237 for name in involved:
238 entry = by_name[name]
239 if entry["standard_uncertainty"] <= 0:
240 raise CliError(f"{name}: a correlated input needs a positive uncertainty")
242 index = {name: position for position, name in enumerate(involved)}
243 size = len(involved)
244 matrix = np.eye(size)
245 for (first, second), coefficient in correlations.items():
246 if first == second or not math.isfinite(coefficient) or not -1 <= coefficient <= 1:
247 raise CliError("correlations need distinct inputs and finite r in [-1, 1]")
248 matrix[index[first], index[second]] = coefficient
249 matrix[index[second], index[first]] = coefficient
250 eigenvalues, eigenvectors = np.linalg.eigh(matrix)
251 tolerance = 1e-12 * max(size, 1)
252 if eigenvalues.min() < -tolerance:
253 raise CliError(
254 "the correlation matrix is not positive semidefinite; the supplied "
255 "coefficients cannot come from a single joint distribution"
256 )
257 # audit-units: ignore UNC003 -- numeric eigenvalues, not uncertain variables
258 factor = eigenvectors * np.sqrt(np.maximum(eigenvalues, 0.0))
259 return involved, factor
262def _draw_correlated(
263 rng: Any,
264 variables: list[dict[str, Any]],
265 correlations: dict[tuple[str, str], float],
266 trials: int,
267) -> dict[str, Any]:
268 """Sample a jointly normal distribution, allowing perfect correlations."""
270 involved, factor = correlation_factor(variables, correlations)
271 by_name = {entry["name"]: entry for entry in variables}
272 for name in involved:
273 if by_name[name]["distribution"] != "normal":
274 raise CliError(
275 f"{name}: correlated Monte Carlo sampling requires normal marginals "
276 "and a jointly normal model; other joint distributions need a "
277 "separate sampler, not a change to the declared input distribution"
278 )
279 index = {name: position for position, name in enumerate(involved)}
280 size = len(involved)
282 standard = rng.standard_normal((size, trials))
283 correlated = factor @ standard
284 samples: dict[str, Any] = {}
285 for name in involved:
286 entry = by_name[name]
287 samples[name] = (
288 entry["value"] + entry["standard_uncertainty"] * correlated[index[name]]
289 )
290 for entry in variables:
291 if entry["name"] not in samples:
292 samples[entry["name"]] = _draw(rng, entry, trials)
293 return samples
296def monte_carlo(
297 tree: Any,
298 variables: list[dict[str, Any]],
299 correlations: dict[tuple[str, str], float],
300 coverage_probability: float,
301 trials: int,
302 seed: int,
303) -> dict[str, Any]:
304 """Propagate the input distributions by Monte Carlo sampling."""
306 try:
307 import numpy as np
308 except ImportError as exc:
309 raise CliError(
310 f"NumPy is unavailable; install with `{_common.PINNED_INSTALL}`"
311 ) from exc
313 if trials * max(len(variables), 1) > 40_000_000:
314 raise CliError(
315 f"{trials} trials across {len(variables)} inputs exceeds the sampling "
316 "budget; lower --trials"
317 )
319 rng = np.random.default_rng(seed)
320 if correlations:
321 samples = _draw_correlated(rng, variables, correlations, trials)
322 else:
323 samples = {
324 entry["name"]: _draw(rng, entry, trials) for entry in variables
325 }
327 drawn = _common.reduce_expression(tree, samples, _common.array_functions())
328 drawn = np.asarray(drawn, dtype=float)
329 if drawn.shape != (trials,):
330 drawn = np.broadcast_to(drawn, (trials,)).astype(float)
331 finite = int(np.count_nonzero(np.isfinite(drawn)))
332 if finite < trials:
333 raise CliError(
334 f"{trials - finite} of {trials} Monte Carlo trials produced a "
335 "non-finite result; the model is undefined over part of the input "
336 "distribution"
337 )
339 tail = 0.5 * (1.0 - coverage_probability)
340 ordered = np.sort(drawn)
341 symmetric = (
342 float(np.quantile(ordered, tail)),
343 float(np.quantile(ordered, 1.0 - tail)),
344 )
345 shortest = _common.shortest_coverage_interval(ordered, coverage_probability)
346 return {
347 "trials": trials,
348 "seed": seed,
349 "mean": float(np.mean(drawn)),
350 "standard_uncertainty": float(np.std(drawn, ddof=1)),
351 "median": float(np.median(drawn)),
352 "coverage_probability": coverage_probability,
353 "probabilistically_symmetric_interval": list(symmetric),
354 "shortest_coverage_interval": list(shortest),
355 "numerical_stability_verified": False,
356 }
359def validate_linearization(
360 framework: dict[str, Any], sampling: dict[str, Any], significant_digits: int
361) -> dict[str, Any]:
362 """Compare fixed-trial endpoints; do not imply adaptive JCGM validation."""
364 combined = sampling["standard_uncertainty"]
365 if combined <= 0:
366 return {
367 "significant_digits": significant_digits,
368 "numerical_tolerance": None,
369 "gum_framework_validated": None,
370 "endpoint_agreement": None,
371 "note": "Monte Carlo standard uncertainty is zero; no comparison is defined",
372 }
373 tolerance = _common.numerical_tolerance(combined, significant_digits)
374 guf_low, guf_high = framework["coverage_interval"]
375 mc_low, mc_high = sampling["probabilistically_symmetric_interval"]
376 low_gap = abs(guf_low - mc_low)
377 high_gap = abs(guf_high - mc_high)
378 validated = low_gap <= tolerance and high_gap <= tolerance
379 return {
380 "significant_digits": significant_digits,
381 "numerical_tolerance": tolerance,
382 "endpoint_difference_low": low_gap,
383 "endpoint_difference_high": high_gap,
384 "endpoint_agreement": validated,
385 "gum_framework_validated": None,
386 "numerical_stability_verified": False,
387 "note": (
388 ("endpoints agree" if validated else "endpoints disagree")
389 + " at this trial count; this is a fixed-trial diagnostic. Check Monte "
390 "Carlo stabilization (JCGM 101 7.9 and 8.2, tolerance delta/5), input "
391 "distributions, and model validity before selecting a reported interval"
392 ),
393 }
396def collect_warnings(
397 framework: dict[str, Any],
398 sampling: dict[str, Any] | None,
399 correlations: dict[tuple[str, str], float],
400) -> list[str]:
401 """Flag conditions that make the headline numbers misleading."""
403 warnings: list[str] = []
404 dof = framework["effective_degrees_of_freedom"]
405 if dof is not None and dof < 6:
406 warnings.append(
407 f"effective degrees of freedom is {dof:.1f}; the coverage factor is "
408 "dominated by a small Type A sample and is unstable"
409 )
410 if correlations:
411 warnings.append(
412 "correlated sampling assumes a jointly normal distribution; covariance "
413 "alone does not specify a general joint distribution"
414 )
415 dominant = max(framework["inputs"], key=lambda item: item["variance_fraction"])
416 if dominant["variance_fraction"] > 0.9 and len(framework["inputs"]) > 1:
417 warnings.append(
418 f"{dominant['name']} contributes "
419 f"{dominant['variance_fraction'] * 100:.1f}% of the summed squared "
420 "contributions; the other inputs barely affect the result"
421 )
422 if framework["covariance_term"] != 0:
423 share = framework["covariance_term"] / framework["independent_variance"]
424 warnings.append(
425 f"the covariance term changes the combined variance by {share * 100:+.0f}% "
426 "of the summed squared contributions, so the budget percentages below "
427 "do not add up to u_c"
428 )
429 for item in framework["inputs"]:
430 if item["value"] != 0:
431 relative = item["standard_uncertainty"] / abs(item["value"])
432 if relative > 0.3:
433 warnings.append(
434 f"{item['name']} has a relative standard uncertainty of "
435 f"{relative * 100:.0f}%; check nonlinearity and physical support "
436 "of the input distribution (a linear model can still be exact)"
437 )
438 if sampling is not None:
439 if any(item["dof"] is not None for item in framework["inputs"]):
440 warnings.append(
441 "finite input degrees of freedom affect the GUM coverage factor "
442 "but are not sampled: this CLI uses the declared fixed PDFs, not "
443 "Student-t input distributions or uncertainty in estimated variances"
444 )
445 if sampling["trials"] < _common.RECOMMENDED_TRIALS:
446 warnings.append(
447 f"{sampling['trials']} Monte Carlo trials is below the "
448 f"{_common.RECOMMENDED_TRIALS} initial trial-count guideline; "
449 "endpoint agreement can be affected by sampling noise"
450 )
451 combined = framework["combined_standard_uncertainty"]
452 if combined > 0:
453 shift = abs(sampling["mean"] - framework["value"]) / combined
454 if shift > 0.1:
455 warnings.append(
456 f"the Monte Carlo mean differs from the model value by "
457 f"{shift:.2f} u_c, which indicates a nonlinear model"
458 )
459 return warnings
462def render_markdown(document: dict[str, Any]) -> str:
463 """Render a human-readable propagation report."""
465 framework = document["gum_framework"]
466 unit = f" {document['unit']}" if document.get("unit") else ""
467 lines = [
468 f"# Uncertainty propagation: {document.get('measurand') or 'measurand'}",
469 "",
470 f"Model: `{document['expression']}`",
471 "",
472 "## GUM uncertainty framework (JCGM 100:2008)",
473 "",
474 f"- Value: {framework['value']:.6g}{unit}",
475 f"- Combined standard uncertainty u_c: "
476 f"{framework['combined_standard_uncertainty']:.6g}{unit}",
477 f"- Effective degrees of freedom: "
478 f"{_common.format_dof(framework['effective_degrees_of_freedom'])}",
479 f"- Coverage factor k: {framework['coverage_factor']:.4g} "
480 f"at p = {framework['coverage_probability']:.3g}",
481 f"- Expanded uncertainty U: {framework['expanded_uncertainty']:.6g}{unit}",
482 f"- Coverage interval: [{framework['coverage_interval'][0]:.6g}, "
483 f"{framework['coverage_interval'][1]:.6g}]{unit}",
484 "",
485 "## Uncertainty budget",
486 "",
487 "| Input | Value | u(x) | Distribution | c = dy/dx | c*u(x) | % variance |",
488 "| --- | --- | --- | --- | --- | --- | --- |",
489 ]
490 for item in framework["inputs"]:
491 lines.append(
492 f"| {item['name']} | {item['value']:.6g} | "
493 f"{item['standard_uncertainty']:.4g} | {item['distribution']} | "
494 f"{item['sensitivity']:.4g} | {item['contribution']:.4g} | "
495 f"{item['variance_fraction'] * 100:.1f} |"
496 )
497 sampling = document.get("monte_carlo")
498 if sampling:
499 lines += [
500 "",
501 "## Monte Carlo propagation (JCGM 101:2008)",
502 "",
503 f"- Trials: {sampling['trials']} (seed {sampling['seed']})",
504 f"- Mean: {sampling['mean']:.6g}{unit}",
505 f"- Standard uncertainty: {sampling['standard_uncertainty']:.6g}{unit}",
506 f"- Probabilistically symmetric interval: "
507 f"[{sampling['probabilistically_symmetric_interval'][0]:.6g}, "
508 f"{sampling['probabilistically_symmetric_interval'][1]:.6g}]{unit}",
509 f"- Shortest coverage interval: "
510 f"[{sampling['shortest_coverage_interval'][0]:.6g}, "
511 f"{sampling['shortest_coverage_interval'][1]:.6g}]{unit}",
512 ]
513 validation = document.get("validation")
514 if validation and validation.get("endpoint_agreement") is not None:
515 verdict = "AGREE" if validation["endpoint_agreement"] else "DISAGREE"
516 lines += [
517 "",
518 "## Fixed-trial interval comparison (not full JCGM validation)",
519 "",
520 f"- Numerical tolerance delta: {validation['numerical_tolerance']:.4g}",
521 f"- Endpoint differences: "
522 f"{validation['endpoint_difference_low']:.4g} (low), "
523 f"{validation['endpoint_difference_high']:.4g} (high)",
524 f"- Verdict: **{verdict}** - {validation['note']}",
525 ]
526 if document.get("warnings"):
527 lines += ["", "## Warnings", ""]
528 lines += [f"- {message}" for message in document["warnings"]]
529 return "\n".join(lines) + "\n"
532def build_parser() -> argparse.ArgumentParser:
533 parser = argparse.ArgumentParser(
534 description=(
535 "Propagate uncertainty through a measurement model with both the GUM "
536 "framework and Monte Carlo sampling, then check the linearization."
537 )
538 )
539 parser.add_argument(
540 "--expression",
541 help="measurement model, e.g. 'm / (pi * (d / 2) ** 2 * h)'",
542 )
543 parser.add_argument(
544 "--variable",
545 action="append",
546 default=[],
547 metavar="NAME=VALUE,U[,DIST[,DOF]]",
548 help="input estimate and its standard uncertainty; repeatable",
549 )
550 parser.add_argument(
551 "--correlation",
552 action="append",
553 default=[],
554 metavar="A,B=R",
555 help="correlation coefficient between two inputs; repeatable",
556 )
557 parser.add_argument("--spec", help="JSON file holding the whole model")
558 parser.add_argument("--measurand", help="name of the output quantity")
559 parser.add_argument("--unit", help="unit label for the report only")
560 parser.add_argument(
561 "--coverage",
562 type=_common.probability,
563 default=0.95,
564 help="coverage probability for the expanded uncertainty (default 0.95)",
565 )
566 parser.add_argument(
567 "--trials",
568 type=_common.bounded_int(1_000, _common.MAX_TRIALS),
569 default=_common.DEFAULT_TRIALS,
570 help=f"Monte Carlo trials (default {_common.DEFAULT_TRIALS})",
571 )
572 parser.add_argument(
573 "--seed",
574 type=_common.bounded_int(0, 2**32 - 1),
575 default=_common.DEFAULT_SEED,
576 help="Monte Carlo seed",
577 )
578 parser.add_argument(
579 "--significant-digits",
580 type=_common.bounded_int(1, 2),
581 default=2,
582 help="significant digits retained in u_c for the clause 8 tolerance",
583 )
584 parser.add_argument(
585 "--no-monte-carlo",
586 action="store_true",
587 help="skip sampling and report only the linearized result",
588 )
589 parser.add_argument("--format", choices=("json", "markdown"), default="json")
590 parser.add_argument("--output", help="write the report to this file")
591 parser.add_argument(
592 "--force", action="store_true", help="overwrite an existing output file"
593 )
594 return parser
597def run(arguments: argparse.Namespace) -> dict[str, Any]:
598 """Build the model from CLI flags or a spec file and propagate it."""
600 spec: dict[str, Any] = {}
601 if arguments.spec:
602 loaded = _common.load_json(arguments.spec)
603 if not isinstance(loaded, dict):
604 raise CliError("the spec file must hold a JSON object")
605 spec = loaded
607 expression = arguments.expression or spec.get("expression")
608 if not expression:
609 raise CliError("supply --expression or a spec file containing 'expression'")
611 records: list[Any] = list(spec.get("variables", []))
612 records += [parse_variable(item) for item in arguments.variable]
613 variables = normalize_variables(records)
615 tree = _common.parse_expression(expression)
616 required = set(_common.expression_variables(tree))
617 supplied = {entry["name"] for entry in variables}
618 missing = sorted(required - supplied)
619 if missing:
620 raise CliError(f"no value supplied for: {', '.join(missing)}")
621 unused = sorted(supplied - required)
622 if unused:
623 raise CliError(
624 f"these variables do not appear in the model: {', '.join(unused)}"
625 )
627 pairs: list[tuple[str, str, float]] = []
628 for item in spec.get("correlations", []):
629 if not isinstance(item, (list, tuple)) or len(item) != 3:
630 raise CliError("each spec correlation must be [name_a, name_b, r]")
631 pairs.append(
632 (str(item[0]), str(item[1]), _common.as_finite(item[2], label="correlation"))
633 )
634 pairs += [parse_correlation(item) for item in arguments.correlation]
636 correlations: dict[tuple[str, str], float] = {}
637 for first, second, coefficient in pairs:
638 if first == second:
639 raise CliError("a variable cannot be correlated with itself")
640 for name in (first, second):
641 if name not in supplied:
642 raise CliError(f"correlation references unknown variable {name!r}")
643 if not -1.0 <= coefficient <= 1.0:
644 raise CliError("correlation coefficients must lie in [-1, 1]")
645 key = (first, second) if first < second else (second, first)
646 if key in correlations and correlations[key] != coefficient:
647 raise CliError(f"conflicting correlations given for {key[0]} and {key[1]}")
648 if coefficient != 0.0:
649 correlations[key] = coefficient
651 coverage = spec.get("coverage_probability", arguments.coverage)
652 coverage = _common.as_finite(coverage, label="coverage_probability")
653 if not 0 < coverage < 1:
654 raise CliError("coverage probability must be strictly between 0 and 1")
656 framework = gum_framework(tree, variables, correlations, coverage)
657 document: dict[str, Any] = {
658 "measurand": arguments.measurand or spec.get("measurand"),
659 "unit": arguments.unit or spec.get("unit"),
660 "expression": expression,
661 "gum_framework": framework,
662 "correlations": [
663 {"inputs": list(key), "coefficient": value}
664 for key, value in sorted(correlations.items())
665 ],
666 }
667 sampling = None
668 if not arguments.no_monte_carlo:
669 sampling = monte_carlo(
670 tree, variables, correlations, coverage, arguments.trials, arguments.seed
671 )
672 document["monte_carlo"] = sampling
673 document["validation"] = validate_linearization(
674 framework, sampling, arguments.significant_digits
675 )
676 document["warnings"] = collect_warnings(framework, sampling, correlations)
677 return document
680def main(argv: list[str] | None = None) -> int:
681 arguments = build_parser().parse_args(argv)
682 try:
683 document = run(arguments)
684 if arguments.format == "markdown":
685 _common.emit_text(
686 render_markdown(document),
687 output=arguments.output,
688 force=arguments.force,
689 )
690 else:
691 _common.emit_json(
692 document, output=arguments.output, force=arguments.force
693 )
694 except CliError as error:
695 print(f"error: {error}", file=sys.stderr)
696 return 2
697 return 0
700if __name__ == "__main__":
701 raise SystemExit(main())