Feature/coherent state fid obj

Project.toml

@@ -19,12 +19,12 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 TrajectoryIndexingUtils = "6dad8b7f-dd9a-4c28-9b70-85b9a079bfc8"
 
 [compat]
-DirectTrajOpt = "0.5"
+DirectTrajOpt = "0.6"
 Distributions = "0.25"
 ExponentialAction = "0.2"
 LinearAlgebra = "1.10, 1.11, 1.12"
-NamedTrajectories = "0.6"
-PiccoloQuantumObjects = "0.7"
+NamedTrajectories = "0.7"
+PiccoloQuantumObjects = "0.8"
 Random = "1.10, 1.11, 1.12"
 Reexport = "1.2"
 SparseArrays = "1.10, 1.11, 1.12"

docs/literate/examples/minimum_time_problem.jl

@@ -0,0 +1,157 @@
+# # Minimum Time Optimization
+
+# In this example, we'll show how to optimize for the **fastest possible gate** while maintaining high fidelity. This is useful when you want to minimize gate duration to reduce decoherence effects.
+
+# ## Problem Setup
+
+# We'll start with a smooth pulse solution for an X gate, then optimize for minimum time.
+
+using QuantumCollocation
+using PiccoloQuantumObjects
+using NamedTrajectories
+
+using PiccoloPlots
+using CairoMakie
+
+# Define a simple qubit system:
+
+H_drift = 0.05 * PAULIS.Z
+H_drives = [PAULIS.X, PAULIS.Y]
+T_max = 50.0  # Start with generous time budget
+drive_bounds = [(-0.5, 0.5), (-0.5, 0.5)]
+
+sys = QuantumSystem(H_drift, H_drives, T_max, drive_bounds)
+
+# Target gate: X gate (π rotation about X axis)
+
+U_goal = GATES.X
+N = 101  # Use more timesteps for better time resolution
+
+# ## Step 1: Find a Smooth Pulse Solution
+
+# First, we solve for a smooth pulse with a fixed time to get a good starting point:
+
+prob_smooth = UnitarySmoothPulseProblem(sys, U_goal, N)
+solve!(prob_smooth; max_iter=100)
+
+fid_smooth = unitary_rollout_fidelity(prob_smooth.trajectory, sys)
+duration_smooth = get_duration(prob_smooth.trajectory)
+
+println("Smooth pulse fidelity: ", fid_smooth)
+println("Smooth pulse duration: ", duration_smooth)
+
+# Let's visualize the smooth pulse solution:
+
+plot_unitary_populations(prob_smooth.trajectory)
+
+# ## Step 2: Optimize for Minimum Time
+
+# Now we'll use the smooth pulse as a starting point and optimize for minimum time while constraining the fidelity:
+
+prob_mintime = UnitaryMinimumTimeProblem(
+    prob_smooth,  # Use smooth pulse as initial guess
+    U_goal;
+    final_fidelity=0.9999  # Require high fidelity
+)
+
+# Solve the minimum time problem:
+
+solve!(prob_mintime; max_iter=300)
+
+# ## Results
+
+# Let's compare the results:
+
+fid_mintime = unitary_rollout_fidelity(prob_mintime.trajectory, sys)
+duration_mintime = get_duration(prob_mintime.trajectory)
+
+println("\n=== Comparison ===")
+println("Smooth pulse:")
+println("  Duration: ", duration_smooth)
+println("  Fidelity: ", fid_smooth)
+println("\nMinimum time:")
+println("  Duration: ", duration_mintime)
+println("  Fidelity: ", fid_mintime)
+println("\nSpeedup: ", duration_smooth / duration_mintime, "×")
+
+@assert fid_mintime > 0.9999
+@assert duration_mintime < duration_smooth
+
+# Visualize the minimum time solution:
+
+plot_unitary_populations(prob_mintime.trajectory)
+
+# ## Understanding the Tradeoff
+
+# Let's look at how the control pulses changed:
+
+fig = Figure(size=(1000, 600))
+
+## Smooth pulse controls
+ax1 = Axis(fig[1, 1], xlabel="Time", ylabel="Control amplitude", title="Smooth Pulse")
+for i in 1:size(prob_smooth.trajectory.u, 1)
+    lines!(ax1, prob_smooth.trajectory.u[i, :], label="u$i")
+end
+axislegend(ax1)
+
+## Minimum time controls  
+ax2 = Axis(fig[1, 2], xlabel="Time", ylabel="Control amplitude", title="Minimum Time")
+for i in 1:size(prob_mintime.trajectory.u, 1)
+    lines!(ax2, prob_mintime.trajectory.u[i, :], label="u$i")
+end
+axislegend(ax2)
+
+fig
+
+# Notice how the minimum time solution:
+# - Uses fewer timesteps (compressed duration)
+# - Has more aggressive control amplitudes
+# - May have sharper transitions
+
+# ## Varying the Fidelity Constraint
+
+# Let's see how gate duration depends on the required fidelity:
+
+fidelities = [0.99, 0.995, 0.999, 0.9999]
+durations = Float64[]
+
+for target_fid in fidelities
+    prob = UnitaryMinimumTimeProblem(
+        prob_smooth,
+        U_goal;
+        final_fidelity=target_fid
+    )
+    solve!(prob; max_iter=200, verbose=false)
+    push!(durations, get_duration(prob.trajectory))
+    println("Fidelity ", target_fid, " → Duration ", durations[end])
+end
+
+# Plot the tradeoff:
+
+fig_tradeoff = Figure()
+ax = Axis(
+    fig_tradeoff[1, 1],
+    xlabel="Required Fidelity",
+    ylabel="Gate Duration",
+    title="Fidelity vs Duration Tradeoff"
+)
+scatter!(ax, fidelities, durations, markersize=15)
+lines!(ax, fidelities, durations)
+fig_tradeoff
+
+# ## Key Takeaways
+
+# 1. **Start with smooth pulse** - Good initialization is crucial
+# 2. **Higher fidelity = longer duration** - There's always a tradeoff
+# 3. **More timesteps = better resolution** - Use enough N for time optimization
+# 4. **Control bounds matter** - Tighter bounds → longer minimum time
+#
+# ## When to Use Minimum Time Problems
+#
+# Use `UnitaryMinimumTimeProblem` when:
+# - Decoherence is limiting your gate fidelity
+# - You want to maximize throughput in quantum circuits
+# - You need to compare against theoretical speed limits
+# - You're exploring the fundamental limits of your control system
+#
+# For most applications, `UnitarySmoothPulseProblem` with reasonable duration is sufficient and more robust.

docs/literate/examples/multilevel_transmon.jl

@@ -35,22 +35,18 @@ using CairoMakie
 
 ## define the time parameters
 
-T₀ = 10.0   # total time in ns
 N = 50      # number of time steps
-Δt = T₀ / N # time step
+T₀ = 10.0     # total time in ns
 
 ## define the system parameters
 levels = 5
 δ = 0.2
 
 ## add a bound to the controls
-u_bounds = [0.2, 0.2]
+u_bound = 0.2
 
 ## create the system
-sys = TransmonSystem(levels=levels, δ=δ, T_max=T₀, u_bounds=u_bounds)
-
-## let's look at a drive
-get_drives(sys)[1] |> sparse
+sys = TransmonSystem(levels=levels, δ=δ, T_max=T₀, drive_bounds=fill(u_bound, 2))
 
 
 # Since this is a multilevel transmon and we want to implement an, let's say, $X$ gate on the qubit subspace, i.e., the first two levels we can utilize the `EmbeddedOperator` type to define the target operator.
@@ -75,16 +71,16 @@ get_subspace_identity(op) |> sparse
 # We can then pass this embedded operator to the `UnitarySmoothPulseProblem` template to create the problem
 
 ## create the problem
-prob = UnitarySmoothPulseProblem(sys, op, N, Δt)
+prob = UnitarySmoothPulseProblem(sys, op, N)
 
 ## solve the problem
-load_path = joinpath(dirname(Base.active_project()), "data/multilevel_transmon_example_89ee72.jld2") # hide
-prob.trajectory = load_traj(load_path) # hide
-nothing # hide
+# load_path = joinpath(dirname(Base.active_project()), "data/multilevel_transmon_example_89ee72.jld2") # hide
+# prob.trajectory = load_traj(load_path) # hide
+# nothing # hide
+# solve!(prob; max_iter=50)
 
 #=
 ```julia
-solve!(prob; max_iter=50)
 ```
 
 ```@raw html
@@ -179,7 +175,7 @@ plot_unitary_populations(prob.trajectory; fig_size=(900, 700))
 
 ## create the a leakage suppression problem, initializing with the previous solution
 
-prob_leakage = UnitarySmoothPulseProblem(sys, op, N, Δt;
+prob_leakage = UnitarySmoothPulseProblem(sys, op, N;
     u_guess=prob.trajectory.u[:, :],
     piccolo_options=PiccoloOptions(
         leakage_constraint=true,
@@ -189,13 +185,12 @@ prob_leakage = UnitarySmoothPulseProblem(sys, op, N, Δt;
 )
 
 ## solve the problem
-load_path = joinpath(dirname(Base.active_project()), "data/multilevel_transmon_example_leakage_89ee72.jld2") # hide
-prob_leakage.trajectory = load_traj(load_path) # hide
-nothing # hide
-
+# load_path = joinpath(dirname(Base.active_project()), "data/multilevel_transmon_example_leakage_89ee72.jld2") # hide
+# prob_leakage.trajectory = load_traj(load_path) # hide
+# nothing # hide
+# solve!(prob_leakage; max_iter=250, options=IpoptOptions(eval_hessian=false))
 #=
 ```julia
-solve!(prob_leakage; max_iter=250)
 ```
 
 ```@raw html