From 6aac8b4d5fa3794e1fcc92f3d8fa3fe791125b6e Mon Sep 17 00:00:00 2001 From: Jack Champagne Date: Sat, 18 Apr 2026 12:41:35 -0400 Subject: [PATCH 1/3] bench: add Ipopt+MadNLP allocation profile script Adds benchmark/alloc_profile.jl which runs Profile.Allocs-traced solves on the shared bilinear toy problem for both in-tree NLP backends, saving AllocProfileResult JLD2s via HarmoniqsBenchmarks' benchmark_memory! / save_alloc_profile API. benchmark/Project.toml pins HarmoniqsBenchmarks to the feat/alloc-profile branch until that PR merges. --- .gitignore | 5 +- benchmark/Project.toml | 17 +++++ benchmark/alloc_profile.jl | 133 +++++++++++++++++++++++++++++++++++++ 3 files changed, 154 insertions(+), 1 deletion(-) create mode 100644 benchmark/Project.toml create mode 100644 benchmark/alloc_profile.jl diff --git a/.gitignore b/.gitignore index 2117a39..60d314a 100644 --- a/.gitignore +++ b/.gitignore @@ -33,10 +33,13 @@ Manifest.toml docs/Manifest.toml # Project specific ignores below -# generated example artifacts +# generated example artifacts /examples/**/plots/ /examples/**/trajectories/ +# benchmark output artifacts +/benchmark/results/ + # external pkgs and configs pardiso.lic /.CondaPkg/ diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000..f08b1cf --- /dev/null +++ b/benchmark/Project.toml @@ -0,0 +1,17 @@ +[deps] +DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" +ExponentialAction = "e24c0720-ea99-47e8-929e-571b494574d3" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +HarmoniqsBenchmarks = "f45d0b76-2d23-4568-9599-481e0da131db" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +NamedTrajectories = "538bc3a1-5ab9-4fc3-b776-35ca1e893e08" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[sources] +DirectTrajOpt = {path = ".."} +# TODO: drop rev pin once HarmoniqsBenchmarks.jl#1 (feat/alloc-profile) merges +HarmoniqsBenchmarks = {url = "https://github.com/harmoniqs/HarmoniqsBenchmarks.jl", rev = "feat/alloc-profile"} diff --git a/benchmark/alloc_profile.jl b/benchmark/alloc_profile.jl new file mode 100644 index 0000000..0d21197 --- /dev/null +++ b/benchmark/alloc_profile.jl @@ -0,0 +1,133 @@ +# ============================================================================= +# Ipopt + MadNLP allocation profile — bilinear toy problem +# +# Runs `solve!` once per solver under Profile.Allocs via benchmark_memory! +# from HarmoniqsBenchmarks.jl and saves the sampled trace to +# benchmark/results/allocs/ for hot-path triage. The Piccolissimo alloc- +# profile testitem covers the Altissimo side; this script is the sibling +# for the in-tree NLP solvers. +# +# Uses the same `bilinear_dynamics_and_trajectory` fixture the main test +# suite uses, so the profiled problem is deterministic and small (N=10, +# 4-state × 2-control) — we care about allocation *patterns*, not absolute +# counts on a production-size problem. +# +# Run: +# julia --project=benchmark benchmark/alloc_profile.jl +# ============================================================================= + +using Random +using NamedTrajectories +using SparseArrays +using LinearAlgebra +using DirectTrajOpt +using MathOptInterface +const MOI = MathOptInterface +using Ipopt +using MadNLP +using HarmoniqsBenchmarks + +# Resolve the MadNLPSolverExt extension module so MadNLPOptions is accessible +# (matches the pattern used in Piccolissimo.jl/benchmark/benchmarks.jl). +const MadNLPSolverExt = [ + mod for mod in reverse(Base.loaded_modules_order) + if Symbol(mod) == :MadNLPSolverExt +][1] + +# Pull in the bilinear fixture without duplicating it. +include(joinpath(@__DIR__, "..", "test", "test_utils.jl")) + +Random.seed!(42) + +const RESULTS_DIR = joinpath(@__DIR__, "results", "allocs") +mkpath(RESULTS_DIR) + +# ---------------------------------------------------------------------------- +# Problem builder — wraps the shared fixture with a QuadraticRegularizer-style +# objective so both Ipopt and MadNLP see the same NLP. +# ---------------------------------------------------------------------------- +function build_problem(; N = 10) + G, traj = bilinear_dynamics_and_trajectory(; N = N) + + integrators = [ + BilinearIntegrator(G, :x, :u, traj), + DerivativeIntegrator(:u, :du, traj), + DerivativeIntegrator(:du, :ddu, traj), + ] + + J = TerminalObjective(x -> norm(x - traj.goal.x)^2, :x, traj) + J += QuadraticRegularizer(:u, traj, 1.0) + + prob = DirectTrajOptProblem(traj, J, integrators) + return prob, traj +end + +# ---------------------------------------------------------------------------- +# Profile one solver. Warmup runs on a throwaway deepcopy so JIT/compile +# allocations stay out of the recorded trace. +# ---------------------------------------------------------------------------- +function profile_solver(; solver_name, options_ctor, N = 10, sample_rate = 1.0) + prob_warmup, traj = build_problem(; N = N) + prob_profiled, _ = build_problem(; N = N) + + state_dim = traj.dims[:x] + ctrl_dim = sum(traj.dims[cn] for cn in traj.control_names if cn != traj.timestep; init = 0) + + println("\n[$(solver_name)] JIT warmup on throwaway problem copy...") + DirectTrajOpt.solve!(prob_warmup; options = options_ctor()) + + println("[$(solver_name)] Profiling allocations (sample_rate=$(sample_rate))...") + profile = benchmark_memory!( + package = "DirectTrajOpt", + solver = solver_name, + benchmark_name = "bilinear_N$(N)_$(lowercase(solver_name))", + N = traj.N, + state_dim = state_dim, + control_dim = ctrl_dim, + sample_rate = sample_rate, + warmup = false, + runner = "local", + ) do + DirectTrajOpt.solve!(prob_profiled; options = options_ctor()) + end + + mb = profile.total_bytes / (1024 * 1024) + println("[$(solver_name)] captured $(profile.total_count) samples, $(round(mb; digits=2)) MB total") + + path = save_alloc_profile(RESULTS_DIR, profile.benchmark_name, profile) + println("[$(solver_name)] saved to $(path)") + return profile, path +end + +# ---------------------------------------------------------------------------- +# Entry points +# +# sample_rate default is 0.01 because Ipopt/MadNLP generate orders of magnitude +# more fine-grained allocations than the solve's wall-time budget accommodates +# at sample_rate=1.0 (an N=10 bilinear toy can hang for 15+ minutes at 1.0). +# 0.01 still gives statistically useful traces for hot-path triage. +# ---------------------------------------------------------------------------- +function main(; N = 10, sample_rate = 0.01) + ipopt_profile, ipopt_path = profile_solver(; + solver_name = "Ipopt", + options_ctor = () -> IpoptOptions(max_iter = 50, print_level = 0), + N = N, + sample_rate = sample_rate, + ) + + madnlp_profile, madnlp_path = profile_solver(; + solver_name = "MadNLP", + options_ctor = () -> MadNLPSolverExt.MadNLPOptions(max_iter = 50, print_level = MadNLP.ERROR), + N = N, + sample_rate = sample_rate, + ) + + println("\nDone.") + println(" Ipopt profile: $(ipopt_path) ($(ipopt_profile.total_count) samples)") + println(" MadNLP profile: $(madnlp_path) ($(madnlp_profile.total_count) samples)") + return (ipopt = ipopt_profile, madnlp = madnlp_profile) +end + +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end From 8b5006183ad59f37e8ed16f2039d3b292149904a Mon Sep 17 00:00:00 2001 From: Jack Champagne Date: Sun, 19 Apr 2026 02:09:45 -0400 Subject: [PATCH 2/3] bench: analyzer + coerce MadNLP.ERROR to Int MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit analyze_allocs.jl ranks alloc hotspots by scaled bytes across all profile JLD2s in benchmark/results/allocs/ with noise-frame filtering. Companion to the existing alloc_profile.jl driver. alloc_profile.jl: coerce MadNLP.ERROR to Int for MadNLPOptions — symmetric to Piccolissimo commit 472fe4d. --- benchmark/alloc_profile.jl | 2 +- benchmark/analyze_allocs.jl | 136 ++++++++++++++++++++++++++++++++++++ 2 files changed, 137 insertions(+), 1 deletion(-) create mode 100644 benchmark/analyze_allocs.jl diff --git a/benchmark/alloc_profile.jl b/benchmark/alloc_profile.jl index 0d21197..ae5bf8d 100644 --- a/benchmark/alloc_profile.jl +++ b/benchmark/alloc_profile.jl @@ -117,7 +117,7 @@ function main(; N = 10, sample_rate = 0.01) madnlp_profile, madnlp_path = profile_solver(; solver_name = "MadNLP", - options_ctor = () -> MadNLPSolverExt.MadNLPOptions(max_iter = 50, print_level = MadNLP.ERROR), + options_ctor = () -> MadNLPSolverExt.MadNLPOptions(max_iter = 50, print_level = Int(MadNLP.ERROR)), N = N, sample_rate = sample_rate, ) diff --git a/benchmark/analyze_allocs.jl b/benchmark/analyze_allocs.jl new file mode 100644 index 0000000..9d1ff53 --- /dev/null +++ b/benchmark/analyze_allocs.jl @@ -0,0 +1,136 @@ +using HarmoniqsBenchmarks +using Printf + +const DEFAULT_RESULTS_DIR = joinpath(@__DIR__, "results", "allocs") +results_dir() = isempty(ARGS) ? DEFAULT_RESULTS_DIR : ARGS[1] + +# Noise filters — frames / types from Profile.Allocs itself or the Julia +# toplevel/runtime that do not tell us anything about user-code hotpaths. +const NOISE_FRAME_PATTERNS = [ + "Profile.Allocs", + "gc-alloc-profiler", + "gc-stock.c", + "gc.c:", + "jl_apply", + "jl_toplevel_", + "ijl_toplevel_", + "jl_interpret_toplevel_thunk", + "jl_repl_entrypoint", + "interpreter.c", + "_include(", + "include_string(", + "loading.jl", + "client.jl", + "_start() at sys.so", + "ip:0x", + "_start at ", + " at Base.jl:", + "true_main at jlapi.c", + "__libc_start_main", + "loader_exe.c", + "jl_system_image_data", + "macro expansion at Allocs.jl", + "boot.jl:", + "jl_f__call_latest", +] + +const WRAPPER_FRAME_PATTERNS = [ + "alloc_profile.jl", + "benchmark_memory!", + "HarmoniqsBenchmarks", +] + +const NOISE_TYPE_PATTERNS = [ + "Profile.Allocs", +] + +_is_noise_frame(f) = any(p -> occursin(p, f), NOISE_FRAME_PATTERNS) +_is_noise_type(t) = any(p -> occursin(p, t), NOISE_TYPE_PATTERNS) + +function _first_user_frame(stack) + for f in stack + _is_noise_frame(f) && continue + any(p -> occursin(p, f), WRAPPER_FRAME_PATTERNS) && continue + return f + end + return isempty(stack) ? "" : stack[end] +end + +_is_wrapper_frame(f) = any(p -> occursin(p, f), WRAPPER_FRAME_PATTERNS) + +function top_frames(profile; k = 25, scale_to_total = true, drop_wrappers = true) + by_frame = Dict{String, Tuple{Int, Int}}() + for s in profile.samples + _is_noise_type(s.type_name) && continue + for frame in s.stacktrace + _is_noise_frame(frame) && continue + drop_wrappers && _is_wrapper_frame(frame) && continue + cnt, bytes = get(by_frame, frame, (0, 0)) + by_frame[frame] = (cnt + 1, bytes + s.size_bytes) + end + end + ranked = sort(collect(by_frame); by = x -> -x[2][2])[1:min(k, length(by_frame))] + scale = scale_to_total ? (1 / profile.sample_rate) : 1.0 + println("\nTop $(length(ranked)) user frames by allocated bytes (scaled ×$(Int(scale))):") + println(rpad(" bytes", 14), rpad("samples", 10), "frame") + for (frame, (cnt, bytes)) in ranked + @printf " %-12s %-8d %s\n" _fmt_bytes(bytes * scale) cnt _truncate(frame, 140) + end +end + +function top_leaf_callsites(profile; k = 25, scale_to_total = true) + by_leaf = Dict{String, Tuple{Int, Int}}() + for s in profile.samples + _is_noise_type(s.type_name) && continue + leaf = _first_user_frame(s.stacktrace) + cnt, bytes = get(by_leaf, leaf, (0, 0)) + by_leaf[leaf] = (cnt + 1, bytes + s.size_bytes) + end + ranked = sort(collect(by_leaf); by = x -> -x[2][2])[1:min(k, length(by_leaf))] + scale = scale_to_total ? (1 / profile.sample_rate) : 1.0 + println("\nTop $(length(ranked)) leaf call sites by allocated bytes (scaled ×$(Int(scale))):") + println(rpad(" bytes", 14), rpad("samples", 10), "leaf") + for (leaf, (cnt, bytes)) in ranked + @printf " %-12s %-8d %s\n" _fmt_bytes(bytes * scale) cnt _truncate(leaf, 140) + end +end + +function top_types(profile; k = 15, scale_to_total = true) + by_type = Dict{String, Tuple{Int, Int}}() + for s in profile.samples + _is_noise_type(s.type_name) && continue + cnt, bytes = get(by_type, s.type_name, (0, 0)) + by_type[s.type_name] = (cnt + 1, bytes + s.size_bytes) + end + ranked = sort(collect(by_type); by = x -> -x[2][2])[1:min(k, length(by_type))] + scale = scale_to_total ? (1 / profile.sample_rate) : 1.0 + println("\nTop $(length(ranked)) allocated types (scaled ×$(Int(scale))):") + println(rpad(" bytes", 14), rpad("samples", 10), "type") + for (t, (cnt, bytes)) in ranked + @printf " %-12s %-8d %s\n" _fmt_bytes(bytes * scale) cnt _truncate(t, 120) + end +end + +_fmt_bytes(b) = b >= 1 << 30 ? @sprintf("%.2f GB", b / (1 << 30)) : + b >= 1 << 20 ? @sprintf("%.2f MB", b / (1 << 20)) : + b >= 1 << 10 ? @sprintf("%.2f KB", b / (1 << 10)) : + @sprintf("%d B", Int(round(b))) + +_truncate(s, n) = length(s) <= n ? s : string(first(s, n - 1), "…") + +function main() + dir = results_dir() + files = sort(filter(f -> endswith(f, "_allocs.jld2"), readdir(dir; join = true))) + isempty(files) && (println("no *_allocs.jld2 files under $dir"); return) + for path in files + profile = load_alloc_profile(path) + println("=" ^ 100) + println(basename(path)) + @printf " solver=%s N=%d sample_rate=%g samples=%d total=%s\n" profile.solver profile.N profile.sample_rate profile.total_count _fmt_bytes(profile.total_bytes) + top_types(profile; k = 10) + top_leaf_callsites(profile; k = 20) + top_frames(profile; k = 20) + end +end + +main() From 9ba555767db864335ff35f174527dfcdec7a8796 Mon Sep 17 00:00:00 2001 From: Jack Champagne Date: Sun, 19 Apr 2026 02:09:48 -0400 Subject: [PATCH 3/3] =?UTF-8?q?debug:=20diagnostic=20for=20Ipopt=E2=86=94M?= =?UTF-8?q?adNLP=20global-phase=20basin=20divergence?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Identified by Gennadi. abs2(tr(U_goal'·U))/n² is phase-insensitive, so +im·U_goal and -im·U_goal are equally valid optima. With free Δt, solvers fall into different basins depending on init + step sequence; with fixed Δt, both consistently land in the same basin. --- scripts/diagnose_phase_bifurcation.jl | 201 ++++++++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 scripts/diagnose_phase_bifurcation.jl diff --git a/scripts/diagnose_phase_bifurcation.jl b/scripts/diagnose_phase_bifurcation.jl new file mode 100644 index 0000000..d2c39a2 --- /dev/null +++ b/scripts/diagnose_phase_bifurcation.jl @@ -0,0 +1,201 @@ +# Diagnose global phase bifurcation between Ipopt and MadNLP +# +# Run this AFTER the main comparison script has populated: +# qcp_ipopt, qcp_madnlp, integrator, qtraj +# +# Or run standalone — it sets everything up. + +import Ipopt +import MadNLP + +using DirectTrajOpt +using Piccolo +using Piccolissimo +using LinearAlgebra + +function _get_MadNLPSolverExt() + cur_mods = reverse(Base.loaded_modules_order) + candidate_mods = [n for n = cur_mods if Symbol(n) == :MadNLPSolverExt] + length(candidate_mods) == 1 || error("Expected 1 MadNLPSolverExt, found $(length(candidate_mods))") + return candidate_mods[1] +end +const MadNLPSolverExt = _get_MadNLPSolverExt() + +# ---------------------------------------------------------------------------- # +# Setup: same problem for both solvers +# ---------------------------------------------------------------------------- # + +H_drift = 0.5 * PAULIS.Z +H_drives = [PAULIS.X, PAULIS.Y] +drive_bounds = [1., 1.] +sys = QuantumSystem(H_drift, H_drives, drive_bounds) +T = 10.0 +N = 100 +U_goal = GATES.X + +qtraj = UnitaryTrajectory(sys, U_goal, T) +integrator = HermitianExponentialIntegrator(qtraj, N) + +# Build two identical problems +qcp = SmoothPulseProblem(qtraj, N; Q=1e2, R=1e-2, ddu_bound=1e0, Δt_bounds=(0.5, 1.5), integrator=integrator) +qcp_ipopt = deepcopy(qcp) +qcp_madnlp = deepcopy(qcp) + +# ---------------------------------------------------------------------------- # +# Solve +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("Solving with Ipopt...") +println("="^70) +solve!(qcp_ipopt; options=IpoptOptions(max_iter=100)) + +println("\n" * "="^70) +println("Solving with MadNLP...") +println("="^70) +solve!(qcp_madnlp; options=MadNLPSolverExt.MadNLPOptions(max_iter=100)) + +# ---------------------------------------------------------------------------- # +# Extract final unitaries from the solver trajectory data +# ---------------------------------------------------------------------------- # + +function get_final_unitary(qcp) + traj = get_trajectory(qcp) + # The state Ũ⃗ at the final knot point + Ũ⃗_final = traj[:Ũ⃗][:, end] + return iso_vec_to_operator(Ũ⃗_final) +end + +U_ipopt = get_final_unitary(qcp_ipopt) +U_madnlp = get_final_unitary(qcp_madnlp) + +# ---------------------------------------------------------------------------- # +# Diagnose: what phase did each solver converge to? +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("PHASE ANALYSIS") +println("="^70) + +# The trace inner product tr(U_goal' * U_final) +# For a perfect solution: tr(X' * (α·X)) = α·tr(I) = 2α +# where α = ±i (from Schrödinger evolution: U = exp(-iHt)) +trace_ipopt = tr(U_goal' * U_ipopt) +trace_madnlp = tr(U_goal' * U_madnlp) + +println("\ntr(U_goal' * U_final):") +println(" Ipopt: $trace_ipopt") +println(" MadNLP: $trace_madnlp") + +println("\nPhase angle (arg of trace):") +println(" Ipopt: $(angle(trace_ipopt)) rad ($(rad2deg(angle(trace_ipopt)))°)") +println(" MadNLP: $(angle(trace_madnlp)) rad ($(rad2deg(angle(trace_madnlp)))°)") + +println("\nPhase sign diagnostic (im * tr / |tr|):") +phase_ipopt = trace_ipopt / abs(trace_ipopt) +phase_madnlp = trace_madnlp / abs(trace_madnlp) +println(" Ipopt: $phase_ipopt") +println(" MadNLP: $phase_madnlp") + +println("\nFidelity (abs2-based, phase-insensitive):") +n = size(U_goal, 1) +fid_ipopt = abs2(trace_ipopt) / n^2 +fid_madnlp = abs2(trace_madnlp) / n^2 +println(" Ipopt: $fid_ipopt") +println(" MadNLP: $fid_madnlp") + +println("\nFidelity via Piccolo API:") +println(" Ipopt: $(fidelity(qcp_ipopt))") +println(" MadNLP: $(fidelity(qcp_madnlp))") + +# ---------------------------------------------------------------------------- # +# Check if solvers converged to different phase basins +# ---------------------------------------------------------------------------- # + +phase_diff = angle(trace_ipopt) - angle(trace_madnlp) +# Normalize to [-π, π] +phase_diff = mod(phase_diff + π, 2π) - π + +println("\n" * "="^70) +println("PHASE DIFFERENCE BETWEEN SOLVERS") +println("="^70) +println(" Δφ = $(phase_diff) rad ($(rad2deg(phase_diff))°)") +if abs(phase_diff) > π/2 + println(" ⚠ SOLVERS CONVERGED TO DIFFERENT PHASE BASINS") + println(" This is the bifurcation Gennadi identified.") + println(" The abs2 in the objective makes ±(im·U_goal) equivalent.") +else + println(" ✓ Solvers converged to the same phase basin this time.") + println(" Re-run to potentially observe bifurcation (depends on numerics).") +end + +# ---------------------------------------------------------------------------- # +# Show the actual final unitaries for inspection +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("FINAL UNITARIES") +println("="^70) + +println("\nU_goal = X gate:") +display(U_goal) + +println("\n\nU_final (Ipopt):") +display(U_ipopt) + +println("\n\nU_final (MadNLP):") +display(U_madnlp) + +println("\n\nU_final / (im * U_goal) — should be ≈ ±1 * I if converged:") +println("\n Ipopt: U / (i·X) =") +display(U_ipopt / (im * U_goal)) +println("\n MadNLP: U / (i·X) =") +display(U_madnlp / (im * U_goal)) + +# ---------------------------------------------------------------------------- # +# Constraint violations (for completeness) +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("CONSTRAINT VIOLATIONS") +println("="^70) + +deltas_ipopt = zeros(integrator.dim) +deltas_madnlp = zeros(integrator.dim) +DirectTrajOpt.evaluate!(deltas_ipopt, integrator, get_trajectory(qcp_ipopt)) +DirectTrajOpt.evaluate!(deltas_madnlp, integrator, get_trajectory(qcp_madnlp)) + +println(" max |constraint| Ipopt: $(maximum(abs.(deltas_ipopt)))") +println(" max |constraint| MadNLP: $(maximum(abs.(deltas_madnlp)))") + +# ---------------------------------------------------------------------------- # +# Datavec comparison +# ---------------------------------------------------------------------------- # + +println("\n" * "="^70) +println("TRAJECTORY DIVERGENCE") +println("="^70) + +data_ipopt = get_trajectory(qcp_ipopt).data +data_madnlp = get_trajectory(qcp_madnlp).data +max_diff = maximum(abs.(data_ipopt .- data_madnlp)) +println(" max |data_ipopt - data_madnlp| = $max_diff") + +println("\n" * "="^70) +println("DIAGNOSIS COMPLETE") +println("="^70) +println(""" +The objective `abs2(tr(U_goal' * U)) / n^2` at: + qc2/src/quantum_objectives.jl:49 + +is phase-insensitive: abs2(z) = abs2(-z) = abs2(iz) = ... +Both +im*U_goal and -im*U_goal are equally valid optima. + +When Δt is a free variable, the expanded search space creates +multiple basins of attraction. The solver's early numerical +trajectory (from initialization + step selection) determines +which basin it falls into. + +With fixed Δt, the landscape is more constrained and both +solvers consistently find the same basin. +""")