Skills › Scientific Agent Skills › qutip › qobj_model_validator.py
qobj_model_validator.py
327 lines · skills/qutip/scripts/qobj_model_validator.py at 5f2bc03
1#!/usr/bin/env python3
2"""Validate a bounded numeric QuTiP model described by strict JSON."""
4from __future__ import annotations
6import argparse
7import math
8import re
9from collections.abc import Mapping
10from typing import Any
12from _common import (
13 MAX_MODEL_OBJECTS,
14 QUTIP_VERSION,
15 CliError,
16 add_output_arguments,
17 bounded_dimensions,
18 emit_json,
19 finite_float,
20 load_json_object,
21 load_qutip,
22 run_cli,
23 validate_keys,
24)
27ROLES = {"hamiltonian", "initial_state", "collapse_operator", "observable"}
28UNIT_CONVENTIONS = {"hbar=1-angular-frequency", "dimensionless-scaled"}
29NAME_PATTERN = re.compile(r"[A-Za-z][A-Za-z0-9_.-]{0,63}\Z")
30MAX_ABS_ENTRY = 1.0e12
33def _name(value: Any, *, context: str) -> str:
34 if not isinstance(value, str) or not NAME_PATTERN.fullmatch(value):
35 raise CliError(
36 f"{context} must match {NAME_PATTERN.pattern!r} and be at most 64 characters"
37 )
38 return value
41def _complex_scalar(value: Any, *, name: str) -> complex:
42 if isinstance(value, Mapping):
43 validate_keys(
44 value,
45 allowed={"real", "imag"},
46 required={"real", "imag"},
47 context=name,
48 )
49 real = finite_float(value["real"], name=f"{name}.real")
50 imag = finite_float(value["imag"], name=f"{name}.imag")
51 result = complex(real, imag)
52 else:
53 result = complex(finite_float(value, name=name), 0.0)
54 if abs(result) > MAX_ABS_ENTRY:
55 raise CliError(f"{name} magnitude exceeds {MAX_ABS_ENTRY}")
56 return result
59def _numeric_data(
60 value: Any,
61 *,
62 role: str,
63 dimension: int,
64 context: str,
65) -> tuple[list[complex] | list[list[complex]], bool]:
66 if not isinstance(value, list) or not value:
67 raise CliError(f"{context}.data must be a non-empty JSON array")
69 is_matrix = all(isinstance(row, list) for row in value)
70 if any(isinstance(row, list) for row in value) and not is_matrix:
71 raise CliError(f"{context}.data cannot mix scalar entries and matrix rows")
73 if not is_matrix:
74 if role != "initial_state":
75 raise CliError(f"{context}.data must be a square matrix for role {role}")
76 if len(value) != dimension:
77 raise CliError(
78 f"{context}.data ket length is {len(value)}; expected {dimension}"
79 )
80 return [
81 _complex_scalar(item, name=f"{context}.data[{index}]")
82 for index, item in enumerate(value)
83 ], True
85 if len(value) != dimension:
86 raise CliError(
87 f"{context}.data has {len(value)} rows; expected {dimension}"
88 )
89 matrix: list[list[complex]] = []
90 for row_index, row in enumerate(value):
91 if len(row) != dimension:
92 raise CliError(
93 f"{context}.data[{row_index}] has {len(row)} entries; "
94 f"expected {dimension}"
95 )
96 matrix.append(
97 [
98 _complex_scalar(
99 item,
100 name=f"{context}.data[{row_index}][{column_index}]",
101 )
102 for column_index, item in enumerate(row)
103 ]
104 )
105 return matrix, False
108def _state_summary(qobj: Any, *, tolerance: float) -> tuple[dict[str, Any], bool]:
109 if qobj.isket:
110 norm = float(qobj.norm())
111 valid = abs(norm - 1.0) <= tolerance
112 return {
113 "representation": "ket",
114 "norm": norm,
115 "normalization_error": abs(norm - 1.0),
116 "normalized_within_tolerance": valid,
117 }, valid
119 trace = complex(qobj.tr())
120 hermitian = bool(qobj.isherm)
121 if hermitian:
122 eigenvalues = [float(value) for value in qobj.eigenenergies()]
123 minimum_eigenvalue = min(eigenvalues)
124 else:
125 minimum_eigenvalue = None
126 valid = (
127 hermitian
128 and abs(trace - 1.0) <= tolerance
129 and minimum_eigenvalue is not None
130 and minimum_eigenvalue >= -tolerance
131 )
132 return {
133 "representation": "density_matrix",
134 "is_hermitian": hermitian,
135 "trace": trace,
136 "trace_error": abs(trace - 1.0),
137 "minimum_eigenvalue": minimum_eigenvalue,
138 "positive_within_tolerance": (
139 minimum_eigenvalue is not None and minimum_eigenvalue >= -tolerance
140 ),
141 }, valid
144def validate_model(document: Mapping[str, Any], qutip_module: Any | None = None) -> dict[str, Any]:
145 """Validate one bounded model and return a portable report."""
147 validate_keys(
148 document,
149 allowed={"schema_version", "unit_convention", "tolerance", "objects"},
150 required={"schema_version", "unit_convention", "objects"},
151 context="model",
152 )
153 if type(document["schema_version"]) is not int or document["schema_version"] != 1:
154 raise CliError("schema_version must be the integer 1")
155 convention = document["unit_convention"]
156 if convention not in UNIT_CONVENTIONS:
157 raise CliError(
158 "unit_convention must be one of: "
159 + ", ".join(sorted(UNIT_CONVENTIONS))
160 )
161 tolerance = finite_float(
162 document.get("tolerance", 1.0e-9),
163 name="tolerance",
164 minimum=1.0e-14,
165 maximum=1.0e-3,
166 )
167 objects = document["objects"]
168 if not isinstance(objects, list) or not 1 <= len(objects) <= MAX_MODEL_OBJECTS:
169 raise CliError(
170 f"objects must contain from 1 through {MAX_MODEL_OBJECTS} entries"
171 )
173 qutip = qutip_module or load_qutip()
174 summaries: list[dict[str, Any]] = []
175 names: set[str] = set()
176 dimensions_seen: list[list[int]] = []
177 hamiltonians = 0
178 initial_states = 0
179 valid = True
181 for index, raw in enumerate(objects):
182 context = f"objects[{index}]"
183 if not isinstance(raw, Mapping):
184 raise CliError(f"{context} must be a JSON object")
185 validate_keys(
186 raw,
187 allowed={"name", "role", "dims", "data", "rate"},
188 required={"name", "role", "dims", "data"},
189 context=context,
190 )
191 name = _name(raw["name"], context=f"{context}.name")
192 if name in names:
193 raise CliError(f"object names must be unique; duplicate {name!r}")
194 names.add(name)
195 role = raw["role"]
196 if role not in ROLES:
197 raise CliError(f"{context}.role must be one of: {', '.join(sorted(ROLES))}")
198 dimensions = bounded_dimensions(raw["dims"], name=f"{context}.dims")
199 dimension = math.prod(dimensions)
200 dimensions_seen.append(dimensions)
201 data, is_ket = _numeric_data(
202 raw["data"],
203 role=role,
204 dimension=dimension,
205 context=context,
206 )
207 if is_ket:
208 qobj = qutip.Qobj(
209 data,
210 dims=[dimensions, [1] * len(dimensions)],
211 )
212 else:
213 qobj = qutip.Qobj(data, dims=[dimensions, dimensions])
215 summary: dict[str, Any] = {
216 "name": name,
217 "role": role,
218 "dims": qobj.dims,
219 "shape": list(qobj.shape),
220 "qobj_type": qobj.type,
221 }
223 object_valid = True
224 if role == "hamiltonian":
225 hamiltonians += 1
226 if "rate" in raw:
227 raise CliError(f"{context}.rate is valid only for collapse_operator")
228 object_valid = bool(qobj.isoper and qobj.isherm)
229 summary["is_hermitian"] = bool(qobj.isherm)
230 elif role == "initial_state":
231 initial_states += 1
232 if "rate" in raw:
233 raise CliError(f"{context}.rate is valid only for collapse_operator")
234 state_summary, object_valid = _state_summary(
235 qobj,
236 tolerance=tolerance,
237 )
238 summary.update(state_summary)
239 elif role == "observable":
240 if "rate" in raw:
241 raise CliError(f"{context}.rate is valid only for collapse_operator")
242 object_valid = bool(qobj.isoper and qobj.isherm)
243 summary["is_hermitian"] = bool(qobj.isherm)
244 else:
245 if is_ket:
246 raise CliError(f"{context}.data must be an operator matrix")
247 if "rate" not in raw:
248 raise CliError(f"{context}.rate is required for collapse_operator")
249 rate = finite_float(
250 raw["rate"],
251 name=f"{context}.rate",
252 minimum=0.0,
253 maximum=10_000.0,
254 )
255 scaled = math.sqrt(rate) * qobj
256 object_valid = bool(qobj.isoper)
257 summary.update(
258 {
259 "rate": rate,
260 "scaling": "sqrt(rate) * operator",
261 "scaled_operator_norm": float(scaled.norm()),
262 }
263 )
265 summary["valid"] = object_valid
266 valid = valid and object_valid
267 summaries.append(summary)
269 if hamiltonians != 1:
270 raise CliError(f"model must contain exactly one hamiltonian; found {hamiltonians}")
271 if initial_states != 1:
272 raise CliError(
273 f"model must contain exactly one initial_state; found {initial_states}"
274 )
276 reference_dims = dimensions_seen[0]
277 compatible_dims = all(item == reference_dims for item in dimensions_seen)
278 valid = valid and compatible_dims
280 return {
281 "report_type": "qutip.qobj_model_validation",
282 "schema_version": 1,
283 "qutip_version": QUTIP_VERSION,
284 "unit_convention": convention,
285 "tolerance": tolerance,
286 "limits": {
287 "hilbert_dimension": 64,
288 "model_objects": MAX_MODEL_OBJECTS,
289 },
290 "checks": {
291 "exactly_one_hamiltonian": True,
292 "exactly_one_initial_state": True,
293 "all_dimensions_compatible": compatible_dims,
294 "numeric_data_only": True,
295 "collapse_rates_nonnegative": True,
296 },
297 "objects": summaries,
298 "valid": valid,
299 "interpretation": (
300 "Structural/numerical preflight only; physical assumptions and "
301 "approximation validity require independent review."
302 ),
303 }
306def build_parser() -> argparse.ArgumentParser:
307 parser = argparse.ArgumentParser(
308 description=(
309 "Validate a bounded strict-JSON QuTiP model containing one "
310 "Hamiltonian, one initial state, and optional collapse/observable objects."
311 )
312 )
313 parser.add_argument("model", help="local strict-JSON model path")
314 add_output_arguments(parser)
315 return parser
318def main(argv: list[str] | None = None) -> int:
319 args = build_parser().parse_args(argv)
320 document = load_json_object(args.model)
321 report = validate_model(document)
322 emit_json(report, output=args.output, force=args.force)
323 return 0 if report["valid"] else 1
326if __name__ == "__main__":
327 raise SystemExit(run_cli(main))