Skills › Scientific Agent Skills › uncertainty-and-units › check_plausibility.py
check_plausibility.py
911 lines · skills/uncertainty-and-units/scripts/check_plausibility.py at 5f2bc03
1#!/usr/bin/env python3
2"""Test a set of quantities against dimensionless groups and known physical scales.
4Unit bookkeeping proves a calculation is dimensionally consistent. It cannot
5say whether the answer is physically possible. A cell 2 m across, a Reynolds
6number of 4e7 in a capillary, and a diffusion time of 300 years across a
7membrane are all dimensionally impeccable.
9This CLI closes that gap three ways. It evaluates named dimensionless groups
10and reports the regime they place the system in; it computes characteristic
11scales such as a diffusion time or a settling velocity; and it compares a
12quantity against a curated band of values that quantity is actually observed to
13take. Every expression is checked for dimensional consistency first, so passing
14a kinematic viscosity where a dynamic one belongs is caught before any number
15is reported.
17Constants come from scipy.constants at run time rather than from literals in
18this file, so they track the CODATA release SciPy ships.
19"""
21from __future__ import annotations
23import argparse
24import math
25import sys
26from typing import Any
28import _common
29from _common import CliError
30from convert_units import build_registry, checked_unit
33MAX_ENTRIES = 64
35# Available to every expression below without being supplied on the command
36# line. The unit strings state what scipy.constants documents each value in.
37CONSTANT_UNITS: tuple[tuple[str, str, str], ...] = (
38 ("k_B", "k", "J/K"),
39 ("N_A", "N_A", "1/mol"),
40 ("R_gas", "R", "J/(mol*K)"),
41 ("c_light", "c", "m/s"),
42 ("h_planck", "h", "J*s"),
43 ("e_charge", "e", "C"),
44 ("epsilon_0", "epsilon_0", "F/m"),
45 ("g_earth", "g", "m/s**2"),
46 ("sigma_sb", "sigma", "W/(m**2*K**4)"),
47 ("m_u", "m_u", "kg"),
48)
50# Each regime is [lower, upper, label]; null is an open end.
51GROUPS: tuple[dict[str, Any], ...] = (
52 {
53 "name": "reynolds",
54 "symbol": "Re",
55 "expression": "density * velocity * length / viscosity",
56 "inputs": {
57 "density": "[mass] / [length] ** 3",
58 "velocity": "[length] / [time]",
59 "length": "[length]",
60 "viscosity": "[mass] / ([length] * [time])",
61 },
62 "regimes": [
63 [None, 2300.0, "laminar (circular pipe, length = diameter)"],
64 [2300.0, 4000.0, "transitional"],
65 [4000.0, None, "turbulent"],
66 ],
67 "note": "viscosity is dynamic; the thresholds are pipe-flow values and "
68 "differ for external and open-channel flow",
69 "source": "White, Fluid Mechanics, 8th ed., ch. 6",
70 },
71 {
72 "name": "peclet",
73 "symbol": "Pe",
74 "expression": "velocity * length / diffusivity",
75 "inputs": {
76 "velocity": "[length] / [time]",
77 "length": "[length]",
78 "diffusivity": "[length] ** 2 / [time]",
79 },
80 "regimes": [
81 [None, 1.0, "diffusion dominates transport"],
82 [1.0, None, "advection dominates transport"],
83 ],
84 "note": "mass-transfer Peclet number; the thermal form replaces the "
85 "diffusivity with the thermal diffusivity",
86 "source": "Deen, Analysis of Transport Phenomena, 2nd ed., ch. 9",
87 },
88 {
89 "name": "damkohler",
90 "symbol": "Da",
91 "expression": "rate_constant * length / velocity",
92 "inputs": {
93 "rate_constant": "1 / [time]",
94 "length": "[length]",
95 "velocity": "[length] / [time]",
96 },
97 "regimes": [
98 [None, 0.1, "slow reaction relative to residence time; low conversion"],
99 [0.1, 10.0, "reaction and transport comparable"],
100 [10.0, None, "fast reaction relative to residence time; high ideal conversion"],
101 ],
102 "note": "first-order Da_I = residence time / reaction time; conversion also "
103 "depends on reactor geometry and mixing, not an automatic mass-transfer verdict",
104 "source": "Fogler, Elements of Chemical Reaction Engineering, 6th ed.",
105 },
106 {
107 "name": "knudsen",
108 "symbol": "Kn",
109 "expression": "mean_free_path / length",
110 "inputs": {"mean_free_path": "[length]", "length": "[length]"},
111 "regimes": [
112 [None, 0.01, "continuum; Navier-Stokes with no-slip applies"],
113 [0.01, 0.1, "slip flow; no-slip boundary condition fails"],
114 [0.1, 10.0, "transition; continuum treatment invalid"],
115 [10.0, None, "free molecular"],
116 ],
117 "note": "compute the mean free path with the mean_free_path_gas scale",
118 "source": "Karniadakis et al., Microflows and Nanoflows, ch. 1",
119 },
120 {
121 "name": "mach",
122 "symbol": "Ma",
123 "expression": "velocity / sound_speed",
124 "inputs": {
125 "velocity": "[length] / [time]",
126 "sound_speed": "[length] / [time]",
127 },
128 "regimes": [
129 [None, 0.3, "incompressible treatment valid to about 5% in density"],
130 [0.3, 0.8, "subsonic compressible"],
131 [0.8, 1.2, "transonic"],
132 [1.2, None, "supersonic"],
133 ],
134 "note": "the 0.3 threshold is a 5% density-change criterion, not a hard limit",
135 "source": "Anderson, Modern Compressible Flow, 4th ed., ch. 1",
136 },
137 {
138 "name": "womersley",
139 "symbol": "Wo",
140 "expression": "radius * sqrt(angular_frequency * density / viscosity)",
141 "inputs": {
142 "radius": "[length]",
143 "angular_frequency": "1 / [time]",
144 "density": "[mass] / [length] ** 3",
145 "viscosity": "[mass] / ([length] * [time])",
146 },
147 "regimes": [
148 [None, 1.0, "quasi-steady; velocity profile stays parabolic"],
149 [1.0, 10.0, "transitional; profile flattens and lags pressure"],
150 [10.0, None, "inertia-dominated plug flow"],
151 ],
152 "note": "angular frequency, not frequency: use 2*pi*f",
153 "source": "Womersley, J. Physiol. 127:553 (1955)",
154 },
155 {
156 "name": "capillary",
157 "symbol": "Ca",
158 "expression": "viscosity * velocity / surface_tension",
159 "inputs": {
160 "viscosity": "[mass] / ([length] * [time])",
161 "velocity": "[length] / [time]",
162 "surface_tension": "[mass] / [time] ** 2",
163 },
164 "regimes": [
165 [None, 0.001, "small viscous-to-capillary stress ratio"],
166 [0.001, None, "viscous effects may deform the interface"],
167 ],
168 "note": "illustrative 1e-3 scale, not a universal deformation or breakup threshold; "
169 "geometry, wetting, viscosity ratio and other stresses matter",
170 "source": "Bruus, Theoretical Microfluidics, ch. 7",
171 },
172 {
173 "name": "weber",
174 "symbol": "We",
175 "expression": "density * velocity ** 2 * length / surface_tension",
176 "inputs": {
177 "density": "[mass] / [length] ** 3",
178 "velocity": "[length] / [time]",
179 "length": "[length]",
180 "surface_tension": "[mass] / [time] ** 2",
181 },
182 "regimes": [
183 [None, 1.0, "small inertial-to-capillary stress ratio"],
184 [1.0, 12.0, "deformation possible; breakup depends on other parameters"],
185 [12.0, None, "bag breakup possible in low-Ohnesorge aerodynamic conditions"],
186 ],
187 "note": "the approximate We=12 bag-breakup threshold assumes ambient-fluid density, "
188 "relative speed, droplet diameter and low Ohnesorge number",
189 "source": "Pilch and Erdman, Int. J. Multiphase Flow 13:741 (1987)",
190 },
191 {
192 "name": "bond",
193 "symbol": "Bo",
194 "expression": "density_difference * g_earth * length ** 2 / surface_tension",
195 "inputs": {
196 "density_difference": "[mass] / [length] ** 3",
197 "length": "[length]",
198 "surface_tension": "[mass] / [time] ** 2",
199 },
200 "regimes": [
201 [None, 1.0, "surface tension dominates gravity"],
202 [1.0, None, "gravity dominates; the interface flattens"],
203 ],
204 "note": "also called the Eotvos number",
205 "source": "de Gennes et al., Capillarity and Wetting Phenomena, ch. 2",
206 },
207 {
208 "name": "stokes_number",
209 "symbol": "Stk",
210 "expression": (
211 "particle_density * particle_diameter ** 2 * velocity "
212 "/ (18 * viscosity * length)"
213 ),
214 "inputs": {
215 "particle_density": "[mass] / [length] ** 3",
216 "particle_diameter": "[length]",
217 "velocity": "[length] / [time]",
218 "viscosity": "[mass] / ([length] * [time])",
219 "length": "[length]",
220 },
221 "regimes": [
222 [None, 0.1, "small inertial lag in the Stokes-drag model"],
223 [0.1, 1.0, "partial slip; sampling bias likely"],
224 [1.0, None, "large inertial lag; particles may leave streamlines"],
225 ],
226 "note": "the tracer assumption behind particle image velocimetry needs "
227 "Stk well below 0.1; also check particle Reynolds number, buoyancy and size",
228 "source": "Raffel et al., Particle Image Velocimetry, 3rd ed., ch. 2",
229 },
230 {
231 "name": "biot",
232 "symbol": "Bi",
233 "expression": "heat_transfer_coefficient * length / conductivity",
234 "inputs": {
235 "heat_transfer_coefficient": "[mass] / ([time] ** 3 * [temperature])",
236 "length": "[length]",
237 "conductivity": "[length] * [mass] / ([time] ** 3 * [temperature])",
238 },
239 "regimes": [
240 [None, 0.1, "lumped capacitance valid; body is nearly isothermal"],
241 [0.1, None, "internal gradients matter; solve the conduction equation"],
242 ],
243 "note": "length is volume divided by surface area for an irregular body",
244 "source": "Incropera et al., Fundamentals of Heat and Mass Transfer, ch. 5",
245 },
246 {
247 "name": "fourier",
248 "symbol": "Fo",
249 "expression": "thermal_diffusivity * time / length ** 2",
250 "inputs": {
251 "thermal_diffusivity": "[length] ** 2 / [time]",
252 "time": "[time]",
253 "length": "[length]",
254 },
255 "regimes": [
256 [None, 0.05, "early relative to the internal diffusion timescale"],
257 [0.05, 1.0, "transient"],
258 [1.0, None, "beyond internal diffusion timescale; equilibrium not established"],
259 ],
260 "note": "dimensionless time for conduction; the same form with a mass "
261 "diffusivity governs diffusion; equilibrium also depends on Biot number, "
262 "geometry and boundary conditions",
263 "source": "Incropera et al., Fundamentals of Heat and Mass Transfer, ch. 5",
264 },
265 {
266 "name": "schmidt",
267 "symbol": "Sc",
268 "expression": "viscosity / (density * diffusivity)",
269 "inputs": {
270 "viscosity": "[mass] / ([length] * [time])",
271 "density": "[mass] / [length] ** 3",
272 "diffusivity": "[length] ** 2 / [time]",
273 },
274 "regimes": [
275 [None, 10.0, "momentum and mass diffusion timescales comparable within an order"],
276 [10.0, None, "momentum diffuses much faster than solute"],
277 ],
278 "note": "about 1 for gases and 1e3 with small molecules in water",
279 "source": "Cussler, Diffusion, 3rd ed., ch. 9",
280 },
281 {
282 "name": "deborah",
283 "symbol": "De",
284 "expression": "relaxation_time / observation_time",
285 "inputs": {"relaxation_time": "[time]", "observation_time": "[time]"},
286 "regimes": [
287 [None, 1.0, "material responds as a viscous liquid"],
288 [1.0, None, "material responds elastically on this timescale"],
289 ],
290 "note": "the same material is a liquid or a solid depending on how long "
291 "you watch it",
292 "source": "Reiner, Physics Today 17:62 (1964)",
293 },
294)
296SCALES: tuple[dict[str, Any], ...] = (
297 {
298 "name": "diffusion_time",
299 "expression": "length ** 2 / diffusivity",
300 "inputs": {"length": "[length]", "diffusivity": "[length] ** 2 / [time]"},
301 "unit": "s",
302 "note": "the L-squared scaling is why diffusion crosses a membrane in "
303 "microseconds and a millimetre of tissue in minutes",
304 "source": "Berg, Random Walks in Biology, ch. 2",
305 },
306 {
307 "name": "thermal_diffusion_time",
308 "expression": "length ** 2 / thermal_diffusivity",
309 "inputs": {
310 "length": "[length]",
311 "thermal_diffusivity": "[length] ** 2 / [time]",
312 },
313 "unit": "s",
314 "note": "same scaling with the thermal diffusivity k/(rho*c_p)",
315 "source": "Incropera et al., Fundamentals of Heat and Mass Transfer, ch. 5",
316 },
317 {
318 "name": "thermal_energy",
319 "expression": "k_B * temperature",
320 "inputs": {"temperature": "[temperature]"},
321 "unit": "J",
322 "note": "the energy scale every biological interaction is measured "
323 "against; 4.14e-21 J at 300 K, or 2.5 kJ/mol",
324 "source": "Phillips et al., Physical Biology of the Cell, 2nd ed., ch. 5",
325 },
326 {
327 "name": "molar_thermal_energy",
328 "expression": "R_gas * temperature",
329 "inputs": {"temperature": "[temperature]"},
330 "unit": "kJ/mol",
331 "note": "compare a specified standard-state free-energy difference with RT; "
332 "smaller energies still change equilibrium probabilities",
333 "source": "Phillips et al., Physical Biology of the Cell, 2nd ed., ch. 6",
334 },
335 {
336 "name": "stokes_settling_velocity",
337 "expression": (
338 "density_difference * g_earth * particle_diameter ** 2 "
339 "/ (18 * viscosity)"
340 ),
341 "inputs": {
342 "density_difference": "[mass] / [length] ** 3",
343 "particle_diameter": "[length]",
344 "viscosity": "[mass] / ([length] * [time])",
345 },
346 "unit": "m/s",
347 "note": "valid only while the particle Reynolds number stays below about "
348 "0.1; check it with the reynolds group",
349 "source": "Batchelor, An Introduction to Fluid Dynamics, ch. 4",
350 },
351 {
352 "name": "mean_free_path_gas",
353 "expression": (
354 "k_B * temperature / (sqrt(2) * pi * molecular_diameter ** 2 * pressure)"
355 ),
356 "inputs": {
357 "temperature": "[temperature]",
358 "molecular_diameter": "[length]",
359 "pressure": "[mass] / ([length] * [time] ** 2)",
360 },
361 "unit": "m",
362 "note": "hard-sphere estimate; the collision diameter of N2 is about 0.37 nm",
363 "source": "Atkins and de Paula, Physical Chemistry, 12th ed., ch. 1",
364 },
365 {
366 "name": "debye_length",
367 "expression": (
368 "sqrt(epsilon_0 * relative_permittivity * k_B * temperature "
369 "/ (2 * N_A * e_charge ** 2 * ionic_strength))"
370 ),
371 "inputs": {
372 "relative_permittivity": "",
373 "temperature": "[temperature]",
374 "ionic_strength": "[substance] / [length] ** 3",
375 },
376 "unit": "nm",
377 "note": "Debye-Huckel continuum estimate using I = 0.5*sum(c_i*z_i**2) in mol/m^3; "
378 "about 0.96 nm at 100 mM and 0.8 nm at 150 mM in water at 25 C",
379 "source": "Israelachvili, Intermolecular and Surface Forces, 3rd ed., ch. 14",
380 },
381 {
382 "name": "capillary_length",
383 "expression": "sqrt(surface_tension / (density * g_earth))",
384 "inputs": {
385 "surface_tension": "[mass] / [time] ** 2",
386 "density": "[mass] / [length] ** 3",
387 },
388 "unit": "mm",
389 "note": "density is the density contrast; liquid density approximates it for "
390 "liquid-air interfaces. About 2.7 mm for water-air",
391 "source": "de Gennes et al., Capillarity and Wetting Phenomena, ch. 2",
392 },
393)
395# Observed ranges, deliberately generous: a value outside one of these is worth
396# a second look, not automatically wrong.
397BANDS: tuple[dict[str, Any], ...] = (
398 {"name": "bacterial_cell_diameter", "low": 0.2, "high": 10.0, "unit": "um",
399 "source": "Milo and Phillips, Cell Biology by the Numbers, ch. 1"},
400 {"name": "eukaryotic_cell_diameter", "low": 5.0, "high": 100.0, "unit": "um",
401 "source": "Milo and Phillips, Cell Biology by the Numbers, ch. 1"},
402 {"name": "cell_membrane_thickness", "low": 3.0, "high": 5.0, "unit": "nm",
403 "source": "Alberts et al., Molecular Biology of the Cell, 7th ed., ch. 10"},
404 {"name": "dna_base_pair_rise", "low": 0.32, "high": 0.36, "unit": "nm",
405 "source": "Bloomfield et al., Nucleic Acids: Structures and Properties"},
406 {"name": "ribosome_diameter", "low": 20.0, "high": 30.0, "unit": "nm",
407 "source": "Milo and Phillips, Cell Biology by the Numbers, ch. 1"},
408 {"name": "protein_molar_mass", "low": 5.0, "high": 1000.0, "unit": "kg/mol",
409 "source": "Illustrative protein molar-mass band; Milo and Phillips, Cell Biology by the Numbers"},
410 {"name": "human_capillary_diameter", "low": 5.0, "high": 10.0, "unit": "um",
411 "source": "Guyton and Hall, Textbook of Medical Physiology, 14th ed., ch. 16"},
412 {"name": "mammalian_body_temperature", "low": 306.0, "high": 315.0, "unit": "K",
413 "source": "Guyton and Hall, Textbook of Medical Physiology, 14th ed., ch. 74"},
414 {"name": "resting_heart_rate", "low": 0.7, "high": 3.0, "unit": "Hz",
415 "source": "Guyton and Hall, Textbook of Medical Physiology, 14th ed., ch. 9"},
416 {"name": "blood_plasma_osmolarity", "low": 275.0, "high": 300.0, "unit": "mol/m**3",
417 "source": "Guyton and Hall, Textbook of Medical Physiology, 14th ed., ch. 25"},
418 {"name": "diffusivity_small_molecule_water", "low": 3e-10, "high": 3e-9,
419 "unit": "m**2/s", "source": "Cussler, Diffusion, 3rd ed., app. A"},
420 {"name": "diffusivity_protein_water", "low": 1e-11, "high": 1.5e-10,
421 "unit": "m**2/s", "source": "Cussler, Diffusion, 3rd ed., app. A"},
422 {"name": "dynamic_viscosity_water", "low": 0.5, "high": 1.5, "unit": "mPa*s",
423 "source": "IAPWS R12-08, viscosity of ordinary water substance"},
424 {"name": "surface_tension_water", "low": 0.06, "high": 0.08, "unit": "N/m",
425 "source": "IAPWS R1-76(2014); illustrative near-ambient liquid-water band"},
426 {"name": "sound_speed_water", "low": 1400.0, "high": 1560.0, "unit": "m/s",
427 "source": "Del Grosso and Mader, J. Acoust. Soc. Am. 52:1442 (1972)"},
428 {"name": "sound_speed_air", "low": 320.0, "high": 350.0, "unit": "m/s",
429 "source": "Cramer, J. Acoust. Soc. Am. 93:2510 (1993)"},
430 {"name": "atmospheric_pressure_sea_level", "low": 95.0, "high": 105.0,
431 "unit": "kPa", "source": "ISO 2533 standard atmosphere"},
432 {"name": "earth_surface_gravity", "low": 9.76, "high": 9.84, "unit": "m/s**2",
433 "source": "WGS 84 normal gravity at the ellipsoid"},
434 {"name": "visible_wavelength", "low": 380.0, "high": 750.0, "unit": "nm",
435 "source": "Illustrative conventional window; CIE e-ILV states visibility has no precise limits"},
436 {"name": "noncovalent_bond_energy", "low": 1.0, "high": 40.0, "unit": "kJ/mol",
437 "source": "Israelachvili, Intermolecular and Surface Forces, 3rd ed., ch. 2"},
438 {"name": "covalent_bond_energy", "low": 150.0, "high": 1000.0, "unit": "kJ/mol",
439 "source": "Atkins and de Paula, Physical Chemistry, 12th ed., data section"},
440 {"name": "atp_hydrolysis_energy", "low": 40.0, "high": 60.0, "unit": "kJ/mol",
441 "source": "Milo and Phillips, Cell Biology by the Numbers; magnitude of cellular hydrolysis free energy, not signed standard-state delta G"},
442)
444CATALOGUE = {
445 "groups": {entry["name"]: entry for entry in GROUPS},
446 "scales": {entry["name"]: entry for entry in SCALES},
447 "bands": {entry["name"]: entry for entry in BANDS},
448}
451def parse_quantity(text: str, registry: Any) -> tuple[str, Any]:
452 """Parse `name=value[ unit]` into a named pint quantity."""
454 name, separator, payload = text.partition("=")
455 name = name.strip()
456 if not separator or not name.isidentifier() or name.startswith("_"):
457 raise CliError(
458 f"quantity {text!r} must look like velocity=1.2 m/s, with a "
459 "Python-identifier name"
460 )
461 fields = payload.strip().split(None, 1)
462 if not fields:
463 raise CliError(f"quantity {name!r} needs a magnitude")
464 magnitude = _common.finite_float(fields[0])
465 unit = checked_unit(fields[1], label=name) if len(fields) == 2 else "dimensionless"
466 try:
467 return name, registry.Quantity(magnitude, unit)
468 except Exception as exc: # pint raises several unrelated types here
469 raise CliError(f"cannot read {name!r} as a quantity: {exc}") from exc
472def constant_quantities(registry: Any) -> dict[str, Any]:
473 """Return the physical constants every expression may use."""
475 try:
476 from scipy import constants
477 except ImportError as exc:
478 raise CliError(
479 f"SciPy is unavailable; install with `{_common.PINNED_INSTALL}`"
480 ) from exc
481 resolved: dict[str, Any] = {}
482 for name, attribute, unit in CONSTANT_UNITS:
483 resolved[name] = registry.Quantity(float(getattr(constants, attribute)), unit)
484 return resolved
487def check_dimensionality(name: str, quantity: Any, expected: str) -> None:
488 """Reject an input whose dimensionality is not the one the formula needs."""
490 if not expected:
491 if not quantity.dimensionless:
492 raise CliError(
493 f"{name} must be dimensionless, but carries {quantity.units:~P}"
494 )
495 return
496 if not quantity.check(expected):
497 raise CliError(
498 f"{name} must have dimensionality {expected}, but {quantity.units:~P} "
499 f"is {quantity.dimensionality}"
500 )
503def evaluate(
504 entry: dict[str, Any], supplied: dict[str, Any], constants: dict[str, Any]
505) -> Any:
506 """Reduce one catalogue expression against the supplied quantities."""
508 tree = _common.parse_expression(entry["expression"])
509 needed = set(_common.expression_variables(tree))
510 required = sorted(needed - set(constants))
511 missing = [name for name in required if name not in supplied]
512 if missing:
513 wanted = ", ".join(
514 f"{name} [{entry['inputs'].get(name) or 'dimensionless'}]"
515 for name in missing
516 )
517 raise CliError(f"{entry['name']} needs: {wanted}")
518 for name in required:
519 check_dimensionality(name, supplied[name], entry["inputs"].get(name, ""))
520 quantity = supplied[name]
521 # Temperature must be absolute, even when the caller supplied degC.
522 physical = quantity.m_as("kelvin") if name == "temperature" else quantity.to_base_units().magnitude # audit-units: ignore UNIT003
523 allow_zero = name in {"velocity", "time", "angular_frequency", "rate_constant", "relaxation_time", "density_difference"}
524 signed = name == "density_difference" and entry["name"] == "stokes_settling_velocity"
525 if not math.isfinite(float(physical)) or (not signed and (physical < 0 or (physical == 0 and not allow_zero))):
526 raise CliError(f"{name} must be physically positive (non-negative for zero-flow/time limits)")
527 values = {name: (supplied[name].to("kelvin") if name == "temperature" else supplied[name]) for name in required}
528 values.update({name: constants[name] for name in needed & set(constants)})
529 return _common.reduce_expression(tree, values, {"sqrt": lambda q: q**0.5})
532def classify(value: float, regimes: list[list[Any]]) -> str:
533 """Return the label of the regime a dimensionless value falls in."""
535 for lower, upper, label in regimes:
536 if (lower is None or value >= lower) and (upper is None or value < upper):
537 return label
538 return "outside every tabulated regime"
541def orders_of_magnitude(value: float, reference: float) -> float | None:
542 """Return log10 of the ratio between a value and a reference.
544 None when the ratio has no logarithm, which keeps the report strict JSON
545 rather than emitting an infinity.
546 """
548 if value <= 0 or reference <= 0:
549 return None
550 return math.log10(value / reference)
553def _describe(quantity: Any) -> dict[str, Any]:
554 """Report a quantity as a magnitude paired with its explicit unit string."""
556 # audit-units: ignore UNIT003 -- the unit travels with the magnitude below
557 value = float(quantity.magnitude)
558 return {"value": value, "unit": f"{quantity.units:~P}"}
561def evaluate_group(
562 entry: dict[str, Any], supplied: dict[str, Any], constants: dict[str, Any]
563) -> dict[str, Any]:
564 """Compute one dimensionless group and place it in a regime."""
566 result = evaluate(entry, supplied, constants)
567 if not getattr(result, "dimensionless", False):
568 raise CliError(
569 f"{entry['name']} evaluated to {result.units:~P}, which is not "
570 "dimensionless; an input carries the wrong quantity"
571 )
572 value = float(result.m_as("dimensionless"))
573 if not math.isfinite(value):
574 raise CliError("dimensionless result is not finite")
575 return {
576 "name": entry["name"],
577 "symbol": entry["symbol"],
578 "expression": entry["expression"],
579 "value": value,
580 "regime": classify(value, entry["regimes"]),
581 "note": entry["note"],
582 "source": entry["source"],
583 }
586def evaluate_scale(
587 entry: dict[str, Any], supplied: dict[str, Any], constants: dict[str, Any]
588) -> dict[str, Any]:
589 """Compute one characteristic scale in its tabulated unit."""
591 result = evaluate(entry, supplied, constants)
592 try:
593 magnitude = float(result.m_as(entry["unit"]))
594 except Exception as exc: # pint raises DimensionalityError and friends
595 raise CliError(
596 f"{entry['name']} evaluated to {result.units:~P}, which cannot be "
597 f"expressed in {entry['unit']}: {exc}"
598 ) from exc
599 if not math.isfinite(magnitude):
600 raise CliError("characteristic scale is not finite")
601 return {
602 "name": entry["name"],
603 "expression": entry["expression"],
604 "value": magnitude,
605 "unit": entry["unit"],
606 "note": entry["note"],
607 "source": entry["source"],
608 }
611def evaluate_band(
612 entry: dict[str, Any], name: str, quantity: Any
613) -> dict[str, Any]:
614 """Compare a supplied quantity against a curated plausible range."""
616 try:
617 magnitude = float(quantity.m_as(entry["unit"]))
618 except Exception as exc: # pint raises DimensionalityError and friends
619 raise CliError(
620 f"{name} is {quantity.units:~P}, which cannot be compared with the "
621 f"{entry['name']} band in {entry['unit']}: {exc}"
622 ) from exc
623 if entry["low"] <= magnitude <= entry["high"]:
624 verdict, decades = "plausible", 0.0
625 else:
626 edge = entry["low"] if magnitude < entry["low"] else entry["high"]
627 decades = orders_of_magnitude(magnitude, edge)
628 # A non-positive magnitude has no ratio to the band, and no physical
629 # quantity in this table can take one.
630 verdict = (
631 "implausible" if decades is None or abs(decades) >= 1 else "questionable"
632 )
633 return {
634 "interpretation": "comparison with an illustrative typical-value band, not a physical limit",
635 "quantity": name,
636 "band": entry["name"],
637 "value": magnitude,
638 "unit": entry["unit"],
639 "range": [entry["low"], entry["high"]],
640 "verdict": verdict,
641 "orders_of_magnitude_outside": decades,
642 "source": entry["source"],
643 }
646def overall_verdict(bands: list[dict[str, Any]]) -> str:
647 """Reduce the band results to a single word."""
649 verdicts = {item["verdict"] for item in bands}
650 if "implausible" in verdicts:
651 return "implausible"
652 if "questionable" in verdicts:
653 return "questionable"
654 return "plausible"
657def exit_code(document: dict[str, Any], threshold: str) -> int:
658 """Return 1 when the verdict meets the configured failure threshold."""
660 order = {"plausible": 0, "questionable": 1, "implausible": 2}
661 if threshold == "none":
662 return 0
663 return 1 if order[document["verdict"]] >= order[threshold] else 0
666def build_parser() -> argparse.ArgumentParser:
667 parser = argparse.ArgumentParser(
668 description=(
669 "Evaluate dimensionless groups and characteristic scales from "
670 "unit-bearing quantities, and compare quantities against curated "
671 "physical ranges."
672 )
673 )
674 parser.add_argument(
675 "--quantity",
676 action="append",
677 default=[],
678 metavar="NAME=VALUE UNIT",
679 help="input quantity, e.g. 'velocity=1.2 m/s'; repeatable",
680 )
681 parser.add_argument(
682 "--group",
683 action="append",
684 default=[],
685 help="dimensionless group to evaluate; repeatable",
686 )
687 parser.add_argument(
688 "--scale",
689 action="append",
690 default=[],
691 help="characteristic scale to compute; repeatable",
692 )
693 parser.add_argument(
694 "--band",
695 action="append",
696 default=[],
697 metavar="BAND=NAME",
698 help="compare a supplied quantity against a curated range, e.g. "
699 "'eukaryotic_cell_diameter=diameter'; repeatable",
700 )
701 parser.add_argument(
702 "--list",
703 action="store_true",
704 dest="list_catalogue",
705 help="list the available groups, scales, and bands, then exit",
706 )
707 parser.add_argument(
708 "--fail-on",
709 choices=("none", "questionable", "implausible"),
710 default="implausible",
711 help="verdict at which the exit status becomes 1 (default implausible)",
712 )
713 parser.add_argument("--format", choices=("json", "markdown"), default="json")
714 parser.add_argument("--output", help="write the report to this file")
715 parser.add_argument(
716 "--force", action="store_true", help="overwrite an existing output file"
717 )
718 return parser
721def _catalogue_document() -> dict[str, Any]:
722 """Describe everything this tool knows how to compute."""
724 return {
725 "groups": [
726 {
727 "name": entry["name"],
728 "symbol": entry["symbol"],
729 "expression": entry["expression"],
730 "inputs": entry["inputs"],
731 }
732 for entry in GROUPS
733 ],
734 "scales": [
735 {
736 "name": entry["name"],
737 "expression": entry["expression"],
738 "inputs": entry["inputs"],
739 "unit": entry["unit"],
740 }
741 for entry in SCALES
742 ],
743 "bands": [
744 {
745 "name": entry["name"],
746 "range": [entry["low"], entry["high"]],
747 "unit": entry["unit"],
748 }
749 for entry in BANDS
750 ],
751 }
754def _lookup(kind: str, name: str) -> dict[str, Any]:
755 """Fetch a catalogue entry or explain what is available."""
757 entry = CATALOGUE[kind].get(name)
758 if entry is None:
759 available = ", ".join(sorted(CATALOGUE[kind]))
760 raise CliError(f"unknown {kind[:-1]} {name!r}; available: {available}")
761 return entry
764def run(arguments: argparse.Namespace) -> dict[str, Any]:
765 """Evaluate every requested group, scale, and band."""
767 if arguments.list_catalogue:
768 return _catalogue_document()
770 requested = len(arguments.group) + len(arguments.scale) + len(arguments.band)
771 if requested == 0:
772 raise CliError("supply at least one of --group, --scale, or --band")
773 if max(requested, len(arguments.quantity)) > MAX_ENTRIES:
774 raise CliError(f"at most {MAX_ENTRIES} entries are supported per run")
776 registry = build_registry()
777 supplied: dict[str, Any] = {}
778 for item in arguments.quantity:
779 name, quantity = parse_quantity(item, registry)
780 if name in supplied:
781 raise CliError(f"quantity {name!r} is supplied more than once")
782 supplied[name] = quantity
783 constants = constant_quantities(registry)
784 shadowed = sorted(set(supplied) & set(constants))
785 if shadowed:
786 raise CliError(
787 f"these names are reserved for physical constants: {', '.join(shadowed)}"
788 )
790 groups = [
791 evaluate_group(_lookup("groups", name), supplied, constants)
792 for name in arguments.group
793 ]
794 scales = [
795 evaluate_scale(_lookup("scales", name), supplied, constants)
796 for name in arguments.scale
797 ]
799 bands: list[dict[str, Any]] = []
800 for item in arguments.band:
801 band_name, separator, quantity_name = item.partition("=")
802 quantity_name = quantity_name.strip()
803 if not separator or not quantity_name:
804 raise CliError(
805 f"band {item!r} must look like eukaryotic_cell_diameter=diameter"
806 )
807 if quantity_name not in supplied:
808 raise CliError(f"band {item!r} references unsupplied quantity "
809 f"{quantity_name!r}")
810 entry = _lookup("bands", band_name.strip())
811 bands.append(evaluate_band(entry, quantity_name, supplied[quantity_name]))
813 warnings: list[str] = []
814 for item in bands:
815 if item["verdict"] == "plausible":
816 continue
817 decades = item["orders_of_magnitude_outside"]
818 distance = (
819 "outside" if decades is None else f"{abs(decades):.1f} decades outside"
820 )
821 warnings.append(
822 f"{item['quantity']} = {item['value']:.4g} {item['unit']} sits "
823 f"{distance} the {item['band']} range "
824 f"{item['range'][0]:.4g}-{item['range'][1]:.4g} {item['unit']}"
825 )
826 for item in groups:
827 if item["value"] < 0:
828 warnings.append(
829 f"{item['symbol']} is not positive, which no physical "
830 "configuration produces; check the sign of an input"
831 )
833 return {
834 "quantities": {
835 name: _describe(quantity) for name, quantity in supplied.items()
836 },
837 "groups": groups,
838 "scales": scales,
839 "bands": bands,
840 "verdict": overall_verdict(bands),
841 "scope": "verdict covers only requested bands; groups and scales are model diagnostics, not validity certificates",
842 "warnings": warnings,
843 }
846def render_markdown(document: dict[str, Any]) -> str:
847 """Render a human-readable plausibility report."""
849 lines = ["# Plausibility check", ""]
850 if document.get("groups"):
851 lines += [
852 "## Dimensionless groups",
853 "",
854 "| Group | Symbol | Value | Regime |",
855 "| --- | --- | --- | --- |",
856 ]
857 for item in document["groups"]:
858 lines.append(
859 f"| {item['name']} | {item['symbol']} | {item['value']:.4g} | "
860 f"{item['regime']} |"
861 )
862 lines.append("")
863 if document.get("scales"):
864 lines += ["## Characteristic scales", "", "| Scale | Value |", "| --- | --- |"]
865 for item in document["scales"]:
866 lines.append(f"| {item['name']} | {item['value']:.4g} {item['unit']} |")
867 lines.append("")
868 if document.get("bands"):
869 lines += [
870 "## Magnitude bands",
871 "",
872 "| Quantity | Band | Value | Range | Verdict |",
873 "| --- | --- | --- | --- | --- |",
874 ]
875 for item in document["bands"]:
876 lines.append(
877 f"| {item['quantity']} | {item['band']} | "
878 f"{item['value']:.4g} {item['unit']} | "
879 f"{item['range'][0]:.4g}-{item['range'][1]:.4g} {item['unit']} | "
880 f"{item['verdict']} |"
881 )
882 lines.append("")
883 lines.append(f"**Verdict: {document['verdict']}**")
884 if document.get("warnings"):
885 lines += ["", "## Warnings", ""]
886 lines += [f"- {message}" for message in document["warnings"]]
887 return "\n".join(lines) + "\n"
890def main(argv: list[str] | None = None) -> int:
891 arguments = build_parser().parse_args(argv)
892 try:
893 document = run(arguments)
894 if arguments.format == "markdown" and not arguments.list_catalogue:
895 _common.emit_text(
896 render_markdown(document),
897 output=arguments.output,
898 force=arguments.force,
899 )
900 else:
901 _common.emit_json(document, output=arguments.output, force=arguments.force)
902 except CliError as error:
903 print(f"error: {error}", file=sys.stderr)
904 return 2
905 if arguments.list_catalogue:
906 return 0
907 return exit_code(document, arguments.fail_on)
910if __name__ == "__main__":
911 raise SystemExit(main())