GitHub

Saving and Loading Pulses

Optimized pulses are valuable — they can take minutes or hours to compute. This guide shows how to save pulses to disk and reload them for warm-starting, analysis, or deployment to hardware.

Why Save Pulses?

  • Warm-starting: Load a solved pulse as the initial guess for a harder problem (e.g., minimum-time optimization, higher fidelity targets, or robustness sweeps).
  • Reproducibility: Archive optimized pulses alongside your scripts so results can be inspected or extended later.
  • Hardware deployment: Export the final pulse for upload to an AWG or other control electronics.
  • Sharing: Send a .jld2 file to a collaborator — they can rebuild the trajectory and continue optimizing without re-running the original solve.

Saving a Pulse

After solving, extract the optimized pulse with get_pulse and persist it with save. The save(filename, pulse) method writes the pulse under the JLD2 key "pulse", and load_pulse(filename) reads it back as the original pulse type.

using Piccolo

# Set up and solve a simple problem
H_drift = PAULIS[:Z]
H_drives = [PAULIS[:X], PAULIS[:Y]]
sys = QuantumSystem(H_drift, H_drives, [1.0, 1.0])

T, N = 10.0, 100
times = collect(range(0, T, length = N))
pulse = ZeroOrderPulse(0.1 * randn(2, N), times)
qtraj = UnitaryTrajectory(sys, pulse, GATES[:X])

qcp = SmoothPulseProblem(qtraj, N; Q = 100.0, R = 1e-2, ddu_bound = 1.0)
solve!(qcp; max_iter = 50)

# Extract the optimized pulse and save it
optimized_pulse = get_pulse(qcp.qtraj)
save("my_pulse.jld2", optimized_pulse)
    constructing SmoothPulseProblem for UnitaryTrajectory...
┌ Warning: Trajectory has timestep variable :Δt but no bounds on it.
│ Adding default lower bound of 0 to prevent negative timesteps.
│ 
│ Recommended: Add explicit bounds when creating the trajectory:
│   NamedTrajectory(...; Δt_bounds=(min, max))
│ Example:
│   NamedTrajectory(qtraj, N; Δt_bounds=(1e-3, 0.5))
│ 
│ Or use timesteps_all_equal=true in problem options to fix timesteps.
└ @ DirectTrajOpt.Problems ~/.julia/packages/DirectTrajOpt/bnQ8f/src/problems.jl:66

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:    38156
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    38584

Total number of variables............................:     1587
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     1188
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1287
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  9.5926273e+01 8.32e+00 2.37e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.4988583e+01 6.31e+00 9.99e+01  -0.3 4.14e+00   2.0 9.72e-01 2.42e-01f  1
   2  8.9396197e+01 4.38e+00 2.08e+02   0.2 3.15e+00   1.5 9.08e-01 3.06e-01f  1
   3  6.2554324e+01 2.72e+00 1.41e+02  -0.2 2.24e+00   1.0 3.48e-01 3.78e-01f  1
   4  6.2478376e+01 7.42e-01 4.92e+01  -0.1 1.37e+00   1.5 8.80e-01 7.27e-01f  1
⋮
  46  1.8530316e-02 6.47e-06 1.83e-02  -4.0 4.42e-03  -0.5 1.00e+00 1.00e+00h  1
  47  1.8198220e-02 5.00e-05 1.26e-01  -4.0 1.20e-02  -1.0 1.00e+00 1.00e+00h  1
  48  1.7457683e-02 2.93e-04 1.01e-01  -4.0 3.27e-02  -1.5 1.00e+00 1.00e+00h  1
  49  1.5549877e-02 1.74e-03 2.30e-01  -4.0 8.32e-02  -1.9 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  1.4891538e-02 2.74e-04 2.60e-02  -4.0 2.91e-02  -1.5 1.00e+00 1.00e+00h  1

Number of Iterations....: 50

                                   (scaled)                 (unscaled)
Objective...............:   1.4891537531292980e-02    1.4891537531292980e-02
Dual infeasibility......:   2.5988992622093744e-02    2.5988992622093744e-02
Constraint violation....:   2.7424580025146167e-04    2.7424580025146167e-04
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.0186575618335294e-04    1.0186575618335294e-04
Overall NLP error.......:   2.5988992622093744e-02    2.5988992622093744e-02


Number of objective function evaluations             = 51
Number of objective gradient evaluations             = 51
Number of equality constraint evaluations            = 51
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 51
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 50
Total seconds in IPOPT                               = 20.785

EXIT: Maximum Number of Iterations Exceeded.
    initializing optimizer...
      building evaluator: 3 integrators, 0 nonlinear constraints
      dynamics constraints: 1188, nonlinear constraints: 0
        integrator 1 jacobian structure: 0.225s
        integrator 2 jacobian structure: 0.014s
        integrator 3 jacobian structure: 0.002s
      jacobian structure: 38016 nonzeros (0.261s)
        integrator 1 hessian structure: 0.069s
        integrator 2 hessian structure: 0.039s
        integrator 3 hessian structure: 0.006s
        computing objective hessian structure (CompositeObjective)...
          sub-objective 1 (KnotPointObjective): 0.232s
          sub-objective 2 (QuadraticRegularizer): 0.199s
          sub-objective 3 (QuadraticRegularizer): 0.001s
          sub-objective 4 (QuadraticRegularizer): 0.017s
          sub-objective 5 (NullObjective): 0.005s
        objective hessian structure: 0.491s
      hessian structure: 38944 nonzeros (0.606s)
      linear index maps built (0.026s)
      evaluator ready (total: 0.994s)
    evaluator created (1.282s)
    NL constraint bounds extracted (0.017s)
    NLP block data built (0.0s)
    Ipopt optimizer configured (0.006s)
    variables set (0.084s)
        applying constraint: initial value of Ũ⃗
        applying constraint: initial value of u
        applying constraint: final value of u
        applying constraint: bounds on Ũ⃗
        applying constraint: bounds on u
        applying constraint: bounds on du
        applying constraint: bounds on ddu
        applying constraint: bounds on Δt
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0
    linear constraints added: 10 (0.322s)
    optimizer initialization complete (total: 1.71s)
[ Info: Loaded cached trajectory from saving_loading_base_b6b244c.jld2

Bundling a pulse with metadata — fidelity, gate name, system parameters, anything you want to remember — uses jldsave from JLD2 directly:

using JLD2

jldsave(
    "my_pulse_with_metadata.jld2";
    pulse = optimized_pulse,
    fidelity = fidelity(qcp),
    duration = sum(get_timesteps(get_trajectory(qcp))),
    gate = "X",
    system_drift = H_drift,
)

Loading a Pulse

load_pulse reads back a pulse written with save:

loaded_pulse = load_pulse("my_pulse.jld2")
ZeroOrderPulse
  drives: 2
  duration: 10.534068838745865

It is also the right loader for a metadata bundle saved with jldsave(...; pulse=...):

loaded_with_metadata = load_pulse("my_pulse_with_metadata.jld2")
ZeroOrderPulse
  drives: 2
  duration: 10.534068838745865

To inspect every key in a metadata bundle, use JLD2.load:

bundle = load("my_pulse_with_metadata.jld2")
keys(bundle)
KeySet for a Dict{String, Any} with 5 entries. Keys:
  "duration"
  "system_drift"
  "gate"
  "fidelity"
  "pulse"

The loaded pulse is a regular pulse object — you can evaluate it, check its duration, and pass it to any trajectory constructor:

duration(loaded_pulse)
10.534068838745865
n_drives(loaded_pulse)
2
loaded_pulse(5.0)
2-element Vector{Float64}:
  0.08212236543160609
 -0.08425623181764982

Warm-Starting from a Saved Pulse

The most common use case: load a previously optimized pulse and use it as the starting point for a new optimization. This works because the pulse already carries good control values — the optimizer just needs to refine them.

# Load the saved pulse
loaded = load_pulse("my_pulse.jld2")

# Build a new trajectory using the loaded pulse
qtraj_warm = UnitaryTrajectory(sys, loaded, GATES[:X])

# Optimize — starting from a good initial guess converges faster
qcp_warm = SmoothPulseProblem(qtraj_warm, N; Q = 1000.0, R = 1e-3, ddu_bound = 1.0)
solve!(qcp_warm; max_iter = 20)
fidelity(qcp_warm)
0.9999924087336114

Warm-Starting for Minimum-Time Optimization

A particularly powerful pattern: solve a fixed-duration problem first, save the pulse, then load it for time-optimal compression.

# Step 1: Solve and save (in one script)
qcp = SmoothPulseProblem(qtraj, N; Q=100.0, Δt_bounds=(0.01, 0.5))
solve!(qcp; max_iter=200)
save("gate_pulse.jld2", get_pulse(qcp.qtraj))

# Step 2: Load and minimize time (in another script)
saved_pulse = load_pulse("gate_pulse.jld2")
qtraj = UnitaryTrajectory(sys, saved_pulse, U_goal)
qcp = SmoothPulseProblem(qtraj, N; Q=100.0, Δt_bounds=(0.01, 0.5))
qcp_mintime = MinimumTimeProblem(qcp; final_fidelity=0.999)
solve!(qcp_mintime; max_iter=100)

Saving and Loading Trajectories

If you need to save the full optimization state (controls, quantum states at every timestep, and timesteps), save the NamedTrajectory instead. The trajectory-aware save writes under the JLD2 key "traj" and pairs with load_traj:

traj = get_trajectory(qcp)
save("my_trajectory.jld2", traj)

Load it back:

loaded_traj = load_traj("my_trajectory.jld2")
typeof(loaded_traj)
NamedTrajectory{Float64, (:Ũ⃗, :Δt, :t, :u, :du, :ddu), NTuple{6, Int64}, (:Ũ⃗, :u, :du, :ddu, :Δt), NTuple{5, Tuple{Vector{Float64}, Vector{Float64}}}, (:Ũ⃗, :u), Tuple{Vector{Float64}, Vector{Float64}}, (:u,), Tuple{Vector{Float64}}, (:Ũ⃗,), Tuple{Vector{Float64}}, (:Ũ⃗, :Δt, :t, :u, :du, :ddu), NTuple{6, UnitRange{Int64}}, NTuple{6, Symbol}, Tuple{Symbol, Symbol}, NTuple{4, Symbol}, (), Tuple{}, (), Tuple{}, Tuple{}}

The trajectory contains all NLP decision variables. You can inspect controls, states, and timesteps:

loaded_traj[:u][:, 1:3]  # First 3 timesteps of controls
2×3 Matrix{Float64}:
 0.0  -0.0628788  -0.0664955
 0.0  -0.0350897  -0.0495202

Organizing Saved Pulses

For projects with many gates, use a consistent directory structure:

my_project/
├── scripts/
│   ├── optimize_X.jl
│   ├── optimize_CZ.jl
│   └── mintime_CZ.jl      # loads from data/CZ.jld2
└── data/
    ├── X.jld2
    ├── CZ.jld2
    └── CZ_mintime.jld2

Each optimization script saves its result to data/, and downstream scripts (like minimum-time solves) load from there. This creates a clear dependency chain and makes results easy to find.

Tips

Always Save After a Successful Solve

Optimization can be expensive. Get in the habit of saving immediately after solve! returns a good result:

solve!(qcp; max_iter=300)
if fidelity(qcp) > 0.999
    save("data/my_gate.jld2", get_pulse(qcp.qtraj))
end

Include Metadata

Future-you will thank present-you for saving context:

using Dates
jldsave("data/CZ.jld2";
    pulse = get_pulse(qcp.qtraj),
    fidelity = fidelity(qcp),
    duration = sum(get_timesteps(get_trajectory(qcp))),
    gate_name = "CZ",
    system_config = "2-atom Rydberg, 8.7 μm spacing",
    date = string(Dates.now()),
)

Converting Pulse Types Before Saving

If you optimized with ZeroOrderPulse but need a smooth pulse for hardware, convert before saving:

zop = get_pulse(qcp.qtraj)  # ZeroOrderPulse from SmoothPulseProblem

# Convert to cubic spline
ctrl = hcat([zop(t) for t in times]...)
tgts = zeros(size(ctrl))
for k in 1:(N-1)
    tgts[:, k] = (ctrl[:, k+1] - ctrl[:, k]) / (times[k+1] - times[k])
end
tgts[:, N] = tgts[:, N-1]
smooth_pulse = CubicSplinePulse(ctrl, tgts, times)

save("data/CZ_smooth.jld2", smooth_pulse)
# Clean up temporary files
rm("my_pulse.jld2"; force = true)
rm("my_pulse_with_metadata.jld2"; force = true)
rm("my_trajectory.jld2"; force = true)

See Also

This page was generated using Literate.jl.