Skills › Scientific Agent Skills › qutip › steady_state_spectrum_planner.py
steady_state_spectrum_planner.py
245 lines · skills/qutip/scripts/steady_state_spectrum_planner.py at 5f2bc03
1#!/usr/bin/env python3
2"""Plan bounded steady-state and spectrum calculations for QuTiP 5."""
4from __future__ import annotations
6import argparse
7import math
8from typing import Any
10from _common import (
11 MAX_ABS_FREQUENCY,
12 MAX_FREQUENCY_POINTS,
13 MAX_HILBERT_DIMENSION,
14 MAX_TIME,
15 QUTIP_VERSION,
16 CliError,
17 add_output_arguments,
18 bounded_int,
19 emit_json,
20 finite_float,
21 run_cli,
22)
25STEADY_METHODS = ("direct", "eigen", "svd", "power", "propagator")
26LINEAR_SOLVERS = (
27 "auto",
28 "solve",
29 "lstsq",
30 "spsolve",
31 "gmres",
32 "lgmres",
33 "bicgstab",
34 "splu",
35)
38def create_plan(args: argparse.Namespace) -> dict[str, Any]:
39 """Return a current API plan plus numerical acceptance checks."""
41 dimension = bounded_int(
42 args.dimension,
43 name="dimension",
44 minimum=2,
45 maximum=MAX_HILBERT_DIMENSION,
46 )
47 frequency_points = bounded_int(
48 args.frequency_points,
49 name="frequency_points",
50 minimum=2,
51 maximum=MAX_FREQUENCY_POINTS,
52 )
53 tau_points = bounded_int(
54 args.tau_points,
55 name="tau_points",
56 minimum=2,
57 maximum=MAX_FREQUENCY_POINTS,
58 )
59 frequency_min = finite_float(
60 args.frequency_min,
61 name="frequency_min",
62 minimum=-MAX_ABS_FREQUENCY,
63 maximum=MAX_ABS_FREQUENCY,
64 )
65 frequency_max = finite_float(
66 args.frequency_max,
67 name="frequency_max",
68 minimum=-MAX_ABS_FREQUENCY,
69 maximum=MAX_ABS_FREQUENCY,
70 )
71 if frequency_max <= frequency_min:
72 raise CliError("frequency_max must be greater than frequency_min")
73 tau_max = finite_float(
74 args.tau_max,
75 name="tau_max",
76 minimum=0.0,
77 maximum=MAX_TIME,
78 minimum_inclusive=False,
79 )
80 residual_tolerance = finite_float(
81 args.residual_tolerance,
82 name="residual_tolerance",
83 minimum=1.0e-14,
84 maximum=1.0e-3,
85 )
86 positivity_tolerance = finite_float(
87 args.positivity_tolerance,
88 name="positivity_tolerance",
89 minimum=1.0e-14,
90 maximum=1.0e-3,
91 )
93 warnings: list[str] = []
94 blockers: list[str] = []
95 steady_kwargs: dict[str, Any] = {"method": args.steady_method}
96 if args.linear_solver != "auto":
97 steady_kwargs["solver"] = args.linear_solver
98 if args.steady_method == "svd" and dimension > 16:
99 warnings.append("svd is dense and may be expensive beyond small systems")
100 if args.degenerate_possible:
101 warnings.append(
102 "a small residual does not choose a unique state from a degenerate nullspace"
103 )
104 if args.time_dependent:
105 blockers.append(
106 "static steadystate/spectrum is not generally valid for a time-dependent "
107 "Hamiltonian; use a periodic/Floquet or explicit long-time analysis"
108 )
109 if not args.stationary_confirmed:
110 blockers.append("stationarity has not been confirmed for spectral analysis")
112 dt = tau_max / (tau_points - 1)
113 nyquist = math.pi / dt
114 angular_resolution = 2.0 * math.pi / (tau_points * dt)
115 requested_abs_frequency = max(abs(frequency_min), abs(frequency_max))
116 if requested_abs_frequency > nyquist:
117 warnings.append(
118 "requested frequency range exceeds the FFT Nyquist angular frequency"
119 )
121 direct_plan = {
122 "api": "qutip.spectrum",
123 "signature": "spectrum(H, wlist, c_ops, a_op, b_op, solver='es')",
124 "solver_strategy": args.direct_solver,
125 "frequency_grid": {
126 "minimum": frequency_min,
127 "maximum": frequency_max,
128 "points": frequency_points,
129 "units": "angular frequency",
130 },
131 "checks": [
132 "steady-state uniqueness or nullspace structure",
133 "operator order and adjoints",
134 "positive/negative-frequency convention",
135 "compare es, pi, or solve strategy near singular features",
136 "expand frequency range and refine spacing",
137 ],
138 }
139 fft_plan = {
140 "apis": [
141 "qutip.correlation_2op_1t",
142 "qutip.spectrum_correlation_fft",
143 ],
144 "tau_grid": {
145 "start": 0.0,
146 "stop": tau_max,
147 "points": tau_points,
148 "step": dt,
149 },
150 "derived_angular_frequency_limits": {
151 "nyquist": nyquist,
152 "approximate_bin_spacing": angular_resolution,
153 },
154 "checks": [
155 "uniform strictly increasing tau grid",
156 "correlation tail decayed before tau_max",
157 "double tau_max to test resolution",
158 "halve dt to test aliasing",
159 "compare window choices and disclose them",
160 "verify transform sign/normalization on an analytic signal",
161 ],
162 }
164 selected: list[dict[str, Any]] = []
165 if args.spectrum_mode in {"direct", "both"}:
166 selected.append(direct_plan)
167 if args.spectrum_mode in {"fft", "both"}:
168 selected.append(fft_plan)
170 return {
171 "report_type": "qutip.steady_state_spectrum_plan",
172 "schema_version": 1,
173 "qutip_version": QUTIP_VERSION,
174 "model_size": {
175 "hilbert_dimension": dimension,
176 "liouvillian_shape": [dimension * dimension, dimension * dimension],
177 },
178 "steady_state": {
179 "api": "qutip.steadystate",
180 "kwargs": steady_kwargs,
181 "checks": {
182 "liouvillian_residual_norm_at_most": residual_tolerance,
183 "trace_error_at_most": residual_tolerance,
184 "minimum_eigenvalue_at_least": -positivity_tolerance,
185 "hermitian": True,
186 "compare_multiple_initial_states": bool(args.degenerate_possible),
187 },
188 },
189 "spectrum_mode": args.spectrum_mode,
190 "spectrum_plans": selected,
191 "warnings": warnings,
192 "blockers": blockers,
193 "status": "blocked_pending_model_evidence" if blockers else "ready_to_implement",
194 "safety": {
195 "planner_only": True,
196 "network": False,
197 "model_code_loaded": False,
198 "bounded_grids": True,
199 },
200 }
203def build_parser() -> argparse.ArgumentParser:
204 parser = argparse.ArgumentParser(
205 description=(
206 "Plan bounded QuTiP steady-state and direct/FFT spectrum checks "
207 "without importing QuTiP or running a model."
208 )
209 )
210 parser.add_argument("--dimension", type=int, default=2)
211 parser.add_argument("--steady-method", choices=STEADY_METHODS, default="direct")
212 parser.add_argument("--linear-solver", choices=LINEAR_SOLVERS, default="auto")
213 parser.add_argument(
214 "--spectrum-mode",
215 choices=("direct", "fft", "both"),
216 default="both",
217 )
218 parser.add_argument(
219 "--direct-solver",
220 choices=("es", "pi", "solve"),
221 default="es",
222 )
223 parser.add_argument("--frequency-min", type=float, default=-5.0)
224 parser.add_argument("--frequency-max", type=float, default=5.0)
225 parser.add_argument("--frequency-points", type=int, default=1001)
226 parser.add_argument("--tau-max", type=float, default=50.0)
227 parser.add_argument("--tau-points", type=int, default=2001)
228 parser.add_argument("--residual-tolerance", type=float, default=1.0e-9)
229 parser.add_argument("--positivity-tolerance", type=float, default=1.0e-9)
230 parser.add_argument("--stationary-confirmed", action="store_true")
231 parser.add_argument("--time-dependent", action="store_true")
232 parser.add_argument("--degenerate-possible", action="store_true")
233 add_output_arguments(parser)
234 return parser
237def main(argv: list[str] | None = None) -> int:
238 args = build_parser().parse_args(argv)
239 report = create_plan(args)
240 emit_json(report, output=args.output, force=args.force)
241 return 0
244if __name__ == "__main__":
245 raise SystemExit(run_cli(main))