Skills › Scientific Agent Skills › qutip › solver_config_planner.py
solver_config_planner.py
326 lines · skills/qutip/scripts/solver_config_planner.py at 5f2bc03
1#!/usr/bin/env python3
2"""Create a bounded QuTiP 5 solver and validation plan without simulation."""
4from __future__ import annotations
6import argparse
7from typing import Any
9from _common import (
10 DEFAULT_SEED,
11 MAX_TIME,
12 MAX_TIME_POINTS,
13 MAX_TRAJECTORIES,
14 QUTIP_VERSION,
15 CliError,
16 add_output_arguments,
17 bounded_int,
18 emit_json,
19 finite_float,
20 run_cli,
21)
24MODELS = (
25 "closed",
26 "lindblad",
27 "quantum-jump",
28 "bloch-redfield",
29 "diffusive",
30 "periodic-closed",
31 "periodic-open",
32 "heom",
33 "piqs",
34)
37def create_plan(args: argparse.Namespace) -> dict[str, Any]:
38 """Select a current solver and produce a model-specific checklist."""
40 t_final = finite_float(
41 args.t_final,
42 name="t_final",
43 minimum=0.0,
44 maximum=MAX_TIME,
45 minimum_inclusive=False,
46 )
47 time_points = bounded_int(
48 args.time_points,
49 name="time_points",
50 minimum=2,
51 maximum=MAX_TIME_POINTS,
52 )
53 trajectories = bounded_int(
54 args.trajectories,
55 name="trajectories",
56 minimum=1,
57 maximum=MAX_TRAJECTORIES,
58 )
59 seed = bounded_int(args.seed, name="seed", minimum=0, maximum=2**63 - 1)
60 collapse_channels = bounded_int(
61 args.collapse_channels,
62 name="collapse_channels",
63 minimum=0,
64 maximum=64,
65 )
66 period = None
67 if args.period is not None:
68 period = finite_float(
69 args.period,
70 name="period",
71 minimum=0.0,
72 maximum=MAX_TIME,
73 minimum_inclusive=False,
74 )
76 method = "bdf" if args.stiff else "adams"
77 options: dict[str, Any] = {
78 "method": method,
79 "atol": 1.0e-10,
80 "rtol": 1.0e-8,
81 "store_states": bool(args.store_states),
82 "store_final_state": True,
83 "progress_bar": "",
84 }
85 required_inputs = [
86 "unit convention and hbar convention",
87 "ordered subsystem dimensions",
88 "initial-state norm/trace/positivity audit",
89 "Hamiltonian Hermiticity audit",
90 "time-grid scale justification",
91 ]
92 convergence = [
93 "tighten atol and rtol",
94 "compare at least one alternative integrator",
95 "increase output-grid density",
96 "sweep every Hilbert-space truncation",
97 ]
98 warnings: list[str] = []
99 call_arguments: dict[str, Any] = {
100 "tlist": {"start": 0.0, "stop": t_final, "points": time_points},
101 "options": options,
102 }
104 if args.model == "closed":
105 solver = "sesolve"
106 required_inputs.append("evidence that no dissipative channel is modeled")
107 if args.initial_state != "ket":
108 warnings.append(
109 "sesolve requires a ket; use mesolve for a density-matrix initial state"
110 )
111 if collapse_channels:
112 warnings.append("closed model conflicts with nonzero collapse_channels")
113 elif args.model == "lindblad":
114 solver = "mesolve"
115 required_inputs.extend(
116 [
117 "Markovian Lindblad approximation",
118 "each collapse operator written as sqrt(rate) times its operator",
119 ]
120 )
121 if collapse_channels == 0:
122 warnings.append("no collapse channels supplied; mesolve may defer to sesolve")
123 elif args.model == "quantum-jump":
124 solver = "mcsolve"
125 required_inputs.extend(
126 [
127 "physical meaning of the jump unravelling",
128 "nonzero collapse-operator list",
129 ]
130 )
131 convergence.extend(["increase ntraj", "repeat or pair seeds deliberately"])
132 call_arguments.update({"ntraj": trajectories, "seeds": seed})
133 options["keep_runs_results"] = bool(args.store_states)
134 if collapse_channels == 0:
135 warnings.append("mcsolve requires a nonzero jump channel for this plan")
136 elif args.model == "bloch-redfield":
137 solver = "brmesolve"
138 required_inputs.extend(
139 [
140 "Hermitian bath-coupling operators",
141 "angular-frequency bath spectra including negative-frequency behavior",
142 "Born-Markov and bath-stationarity assumptions",
143 "justified sec_cutoff",
144 ]
145 )
146 convergence.extend(
147 ["vary sec_cutoff", "track minimum density-matrix eigenvalue over time"]
148 )
149 call_arguments["sec_cutoff"] = 0.1
150 if not args.weak_coupling_confirmed:
151 warnings.append("weak-coupling assumption has not been confirmed")
152 elif args.model == "diffusive":
153 # Unmonitored channels mix an initially pure state; SSE has no c_ops.
154 solver = (
155 "ssesolve"
156 if args.initial_state == "ket" and collapse_channels == 0
157 else "smesolve"
158 )
159 options.pop("atol")
160 options.pop("rtol")
161 options["method"] = "platen"
162 convergence = [
163 "decrease stochastic dt",
164 "compare an alternative stochastic integrator",
165 "increase output-grid density",
166 "sweep every Hilbert-space truncation",
167 "increase ntraj",
168 ]
169 required_inputs.extend(
170 [
171 "separate monitored sc_ops from unmonitored c_ops",
172 "homodyne versus heterodyne measurement model",
173 "measurement efficiency and record convention",
174 ]
175 )
176 options["dt"] = t_final / max(time_points - 1, 1) / 2.0
177 options["store_measurement"] = False
178 options["keep_runs_results"] = bool(args.store_states)
179 call_arguments.update(
180 {"ntraj": trajectories, "seeds": seed, "heterodyne": False}
181 )
182 elif args.model == "periodic-closed":
183 solver = "FloquetBasis + fsesolve"
184 call_arguments["floquet_basis_options"] = {
185 key: options.pop(key) for key in ("method", "atol", "rtol")
186 }
187 required_inputs.append(
188 "pass floquet_basis_options to FloquetBasis, not fsesolve result options"
189 )
190 required_inputs.extend(
191 [
192 "verified H(t + T) equals H(t)",
193 "quasi-energy branch convention",
194 ]
195 )
196 convergence.extend(
197 ["compare one-period propagator with direct evolution", "sweep precompute grid"]
198 )
199 call_arguments["T"] = period
200 if period is None:
201 warnings.append("a positive --period is required")
202 elif args.model == "periodic-open":
203 solver = "FloquetBasis + fmmesolve"
204 required_inputs.extend(
205 [
206 "verified H(t + T) equals H(t)",
207 "paired coupling operators and bath spectrum callbacks",
208 "weak-coupling Floquet-Markov assumptions",
209 "temperature/frequency unit convention",
210 ]
211 )
212 convergence.extend(
213 [
214 "sweep Floquet precompute grid",
215 "compare with direct mesolve in a valid limit",
216 ]
217 )
218 call_arguments["T"] = period
219 if period is None:
220 warnings.append("a positive --period is required")
221 if not args.weak_coupling_confirmed:
222 warnings.append("weak-coupling assumption has not been confirmed")
223 elif args.model == "heom":
224 solver = "qutip.solver.heom.HEOMSolver"
225 required_inputs.extend(
226 [
227 "bath correlation expansion and coupling operator",
228 "temperature, cutoff, and energy units",
229 "factorized or continued ADO initial condition",
230 ]
231 )
232 convergence.extend(
233 [
234 "increase hierarchy max_depth",
235 "increase bath expansion terms",
236 "compare Matsubara and Pade/environment approximations",
237 ]
238 )
239 call_arguments.update({"max_depth": "model-specific", "bath": "required"})
240 else:
241 solver = "qutip.piqs.Dicke.liouvillian + mesolve"
242 required_inputs.extend(
243 [
244 "permutation symmetry of constituents and dynamics",
245 "Dicke versus uncoupled basis",
246 "separate local and collective rates",
247 ]
248 )
249 convergence.extend(
250 ["compare a small-N full-space model", "verify symmetry-sector assumptions"]
251 )
253 if args.time_dependent:
254 required_inputs.extend(
255 [
256 "QobjEvo with trusted Pythonic callables or numeric arrays",
257 "coefficient sampling/envelope convergence",
258 ]
259 )
260 if args.stiff and args.model == "diffusive":
261 warnings.append(
262 "--stiff does not select bdf for stochastic dynamics; reduce dt and "
263 "review a supported implicit stochastic method"
264 )
265 elif args.stiff:
266 convergence.append("compare bdf with lsoda on representative observables")
268 return {
269 "report_type": "qutip.solver_config_plan",
270 "schema_version": 1,
271 "qutip_version": QUTIP_VERSION,
272 "model": args.model,
273 "initial_state": args.initial_state,
274 "recommended_solver": solver,
275 "call_configuration": call_arguments,
276 "required_inputs_and_assumptions": required_inputs,
277 "convergence_plan": convergence,
278 "warnings": warnings,
279 "status": "needs_review" if warnings else "ready_for_model_construction",
280 "safety": {
281 "network": False,
282 "model_code_loaded": False,
283 "planner_only": True,
284 },
285 }
288def build_parser() -> argparse.ArgumentParser:
289 parser = argparse.ArgumentParser(
290 description=(
291 "Plan a QuTiP 5 solver configuration and model-specific numerical "
292 "checks without importing QuTiP or running a simulation."
293 )
294 )
295 parser.add_argument("--model", choices=MODELS, required=True)
296 parser.add_argument(
297 "--initial-state",
298 choices=("ket", "density"),
299 default="ket",
300 )
301 parser.add_argument(
302 "--collapse-channels", type=int, default=0,
303 help="number of c_ops (unmonitored channels for diffusive models)",
304 )
305 parser.add_argument("--time-dependent", action="store_true")
306 parser.add_argument("--stiff", action="store_true")
307 parser.add_argument("--weak-coupling-confirmed", action="store_true")
308 parser.add_argument("--period", type=float)
309 parser.add_argument("--t-final", type=float, default=10.0)
310 parser.add_argument("--time-points", type=int, default=201)
311 parser.add_argument("--trajectories", type=int, default=400)
312 parser.add_argument("--seed", type=int, default=DEFAULT_SEED)
313 parser.add_argument("--store-states", action="store_true")
314 add_output_arguments(parser)
315 return parser
318def main(argv: list[str] | None = None) -> int:
319 args = build_parser().parse_args(argv)
320 report = create_plan(args)
321 emit_json(report, output=args.output, force=args.force)
322 return 0
325if __name__ == "__main__":
326 raise SystemExit(run_cli(main))