Skip to content

Commit

Permalink
feat: fixed time problems
Browse files Browse the repository at this point in the history
  • Loading branch information
aarontrowbridge committed Jun 28, 2023
1 parent 01931c6 commit 2b491ba
Show file tree
Hide file tree
Showing 9 changed files with 608 additions and 406 deletions.
343 changes: 180 additions & 163 deletions Manifest.toml

Large diffs are not rendered by default.

62 changes: 25 additions & 37 deletions examples/qram/pi_gate.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -18,83 +18,65 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loading data dict from data/limited_drives_T_100_dt_1.0_dda_0.0001_a_0.12566370614359174_max_iter_100_00000.jld2:\n",
"Loading data dict from data/limited_drives_T_100_dt_1.0_dda_0.002_a_0.02_max_iter_500_00002.jld2:\n",
"\n",
" duration = 72.61324981519805\n",
" final_fidelity = 0.9922185163356951\n"
" duration = 26.80948914618727\n",
" final_fidelity = 0.9550975523905912\n"
]
},
{
"data": {
"text/plain": [
"NamedTrajectory{Float64}([1.0 1.0 … 0.007575094597536527 0.0066429378333538545; 0.0 0.0 … 0.0 0.0; … ; 1.8981443306253186e-6 3.7667442303727987e-6 … 1.107924470744334e-10 -5.112400723147125e-21; 0.7334671698504867 0.7334671698504867 … 0.7334671698504867 0.7334671698504867], [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.05541568848816433, 0.0, 0.0, -0.0034719494553613253, -0.000302106228627415, 6.171191907215001e-21, -5.112400723147125e-21, 0.7334671698504867], 100, :Δt, 519, (Ũ⃗ = 512, a = 2, da = 2, dda = 2, Δt = 1, states = 516, controls = 3), (a = ([-0.12566370614359174, -0.12566370614359174], [0.12566370614359174, 0.12566370614359174]), dda = ([-0.0001, -0.0001], [0.0001, 0.0001]), Δt = ([0.5], [1.2])), (Ũ⃗ = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], a = [0.0, 0.0]), (a = [0.0, 0.0],), (Ũ⃗ = [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],), (Ũ⃗ = 1:512, a = 513:514, da = 515:516, dda = 517:518, Δt = 519:519, states = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 507, 508, 509, 510, 511, 512, 513, 514, 515, 516], controls = [517, 518, 519]), (:Ũ⃗, :a, :da, :dda, :Δt), (:Ũ⃗, :a, :da), (:dda, :Δt))"
"NamedTrajectory{Float64}([1.0 1.0 … 0.08078165468060403 0.07882587664153458; 0.0 0.0 … 0.0 0.0; … ; 0.00047588630074294944 0.00076657588428622 … 1.7132035748174996e-7 1.1444956860524965e-19; 0.27080292066855877 0.27080292066855877 … 0.27080292066855877 0.27080292066855877], [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.09915545055176467, 0.0, 0.0, -0.008633865997297213, -0.0031638252315717543, 1.1444956860524965e-19, 1.1444956860524965e-19, 0.27080292066855877], 100, :Δt, 519, (Ũ⃗ = 512, a = 2, da = 2, dda = 2, Δt = 1, states = 516, controls = 3), (a = ([-0.02, -0.02], [0.02, 0.02]), dda = ([-0.002, -0.002], [0.002, 0.002]), Δt = ([0.1], [1.2])), (Ũ⃗ = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], a = [0.0, 0.0]), (a = [0.0, 0.0],), (Ũ⃗ = [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],), (Ũ⃗ = 1:512, a = 513:514, da = 515:516, dda = 517:518, Δt = 519:519, states = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 507, 508, 509, 510, 511, 512, 513, 514, 515, 516], controls = [517, 518, 519]), (:Ũ⃗, :a, :da, :dda, :Δt), (:Ũ⃗, :a, :da), (:dda, :Δt))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"data_path = \"data/T_100_dt_1.0_dda_0.0001_a_0.12566370614359174_max_iter_100_00000.jld2\"\n",
"data_path = \"data/limited_drives_T_100_dt_1.0_dda_0.002_a_0.02_max_iter_500_00002.jld2\"\n",
"data = load_problem(data_path; return_data=true)\n",
"traj = data[\"trajectory\"]"
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"100-element Vector{Float64}:\n",
" 0.0\n",
" 0.7334671698504867\n",
" 1.4669343397009733\n",
" 2.20040150955146\n",
" 2.9338686794019466\n",
" 3.6673358492524333\n",
" 4.40080301910292\n",
" 5.134270188953407\n",
" 5.867737358803893\n",
" 6.60120452865438\n",
"\n",
" 66.74551245639417\n",
" 67.47897962624465\n",
" 68.21244679609514\n",
" 68.94591396594562\n",
" 69.6793811357961\n",
" 70.41284830564659\n",
" 71.14631547549708\n",
" 71.87978264534756\n",
" 72.61324981519805"
"2×100 Matrix{Float64}:\n",
" 0.0 0.00233826 0.0045334 0.00658378 … 0.00453315 0.00233811 0.0\n",
" 0.0 -0.000790894 -0.00154689 -0.00224667 0.00167523 0.000856786 0.0"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pulse = traj.a\n",
"ts = times(traj)"
"ts = times(traj)\n",
"pulse = traj.a"
]
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"# save pulse and times\n",
"pulse_path = \"data/pulse.hdf5\"\n",
"pulse_path = \"data/pulse_2_drive_shorter_time.hdf5\"\n",
"pulse_file = h5open(pulse_path, \"w\")\n",
"pulse_file[\"pulse\"] = pulse\n",
"pulse_file[\"ts\"] = ts\n",
Expand All @@ -103,15 +85,21 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×100 Matrix{Float64}:\n",
" 0.0 0.00254727 0.00505691 … 0.00505566 0.00254656 0.0\n",
" 0.0 -0.000218839 -0.000436656 0.000442186 0.000221585 0.0"
"8×100 Matrix{Float64}:\n",
" 0.0 3.94116e-5 3.67995e-5 … 0.00379268 0.00191292 0.0\n",
" 0.0 -0.00172882 -0.00339493 -0.00183975 -0.000952879 0.0\n",
" 0.0 0.000255824 0.000559932 0.000507485 0.000229223 0.0\n",
" 0.0 5.76738e-5 0.000119058 -0.000104015 -5.0021e-5 0.0\n",
" 0.0 -0.000697201 -0.001356 0.00152926 0.000782295 0.0\n",
" 0.0 0.000126126 0.000269827 … -2.13574e-5 -2.33794e-5 0.0\n",
" 0.0 0.00110387 0.00218492 0.00216998 0.00110303 0.0\n",
" 0.0 0.000466185 0.000899852 -0.000223202 -0.000123923 0.0"
]
},
"metadata": {},
Expand All @@ -120,7 +108,7 @@
],
"source": [
"# load pulse and times\n",
"pulse_path = \"data/pulse.hdf5\"\n",
"pulse_path = \"data/pulse_full.hdf5\"\n",
"pulse_file = h5open(pulse_path, \"r\")\n",
"pulse_file[\"pulse\"][:,:]"
]
Expand Down
Binary file not shown.
Binary file not shown.
11 changes: 8 additions & 3 deletions examples/transmon/pi_gate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,8 @@ system = QuantumSystem(
H_drives
)

T = 400
Δt = 0.1
T = 100
Δt = 0.4
Q = 200.
R = 0.01
R_L1 = 1.0
Expand All @@ -36,7 +36,10 @@ max_iter = 500

time = T * Δt

free_time = false

prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt;
free_time=free_time,
a_bounds=a_bounds,
dda_bound=dda_bound,
max_iter=max_iter,
Expand All @@ -47,7 +50,9 @@ prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt;

solve!(prob)

println("fidelity = ", unitary_fidelity(prob.trajectory[end].Ũ⃗, prob.trajectory.goal.Ũ⃗))
Ũ⃗_final = unitary_rollout(prob.trajectory.initial.Ũ⃗, prob.trajectory.a, prob.trajectory.timestep, system)[:, end]

println("fidelity = ", unitary_fidelity(Ũ⃗_final, prob.trajectory.goal.Ũ⃗))

plot_dir = joinpath(@__DIR__, "plots/pi_gate/")

Expand Down
95 changes: 64 additions & 31 deletions src/dynamics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,22 @@ function QuantumDynamics(
println(" constructing knot point dynamics functions...")
end

@assert all([
!isnothing(state(integrator)) &&
!isnothing(controls(integrator)) &&
!isnothing(timestep(integrator))
for integrator integrators
])
free_time = traj.timestep isa Symbol

if free_time
@assert all([
!isnothing(state(integrator)) &&
!isnothing(controls(integrator)) &&
!isnothing(timestep(integrator))
for integrator integrators
])
else
@assert all([
!isnothing(state(integrator)) &&
!isnothing(controls(integrator))
for integrator integrators
])
end

for integrator integrators
if integrator isa QuantumIntegrator && controls(integrator) isa Tuple
Expand Down Expand Up @@ -91,10 +101,16 @@ function QuantumDynamics(

∂[integrator_comps, 1:2traj.dim] = ∂P(zₜ, zₜ₊₁)
else
x_comps, u_comps, Δt_comps = comps(integrator, traj)

∂xₜf, ∂xₜ₊₁f, ∂uₜf, ∂Δtₜf =
Integrators.jacobian(integrator, zₜ, zₜ₊₁, traj)
if free_time
x_comps, u_comps, Δt_comps = comps(integrator, traj)
∂xₜf, ∂xₜ₊₁f, ∂uₜf, ∂Δtₜf =
Integrators.jacobian(integrator, zₜ, zₜ₊₁, traj)
∂[integrator_comps, Δt_comps] = ∂Δtₜf
else
x_comps, u_comps = comps(integrator, traj)
∂xₜf, ∂xₜ₊₁f, ∂uₜf =
Integrators.jacobian(integrator, zₜ, zₜ₊₁, traj)
end

∂[integrator_comps, x_comps] = ∂xₜf
∂[integrator_comps, x_comps .+ traj.dim] = ∂xₜ₊₁f
Expand All @@ -107,20 +123,26 @@ function QuantumDynamics(
∂[integrator_comps, u_comps] = ∂uₜf
end

∂[integrator_comps, Δt_comps] = ∂Δtₜf
end

elseif integrator isa DerivativeIntegrator

x_comps, dx_comps, Δt_comps = comps(integrator, traj)

∂xₜf, ∂xₜ₊₁f, ∂dxₜf, ∂Δtₜf =
Integrators.jacobian(integrator, zₜ, zₜ₊₁, traj)
if free_time
x_comps, dx_comps, Δt_comps = comps(integrator, traj)
∂xₜf, ∂xₜ₊₁f, ∂dxₜf, ∂Δtₜf =
Integrators.jacobian(integrator, zₜ, zₜ₊₁, traj)
else
x_comps, dx_comps = comps(integrator, traj)
∂xₜf, ∂xₜ₊₁f, ∂dxₜf =
Integrators.jacobian(integrator, zₜ, zₜ₊₁, traj)
end

∂[integrator_comps, x_comps] = ∂xₜf
∂[integrator_comps, x_comps .+ traj.dim] = ∂xₜ₊₁f
∂[integrator_comps, dx_comps] = ∂dxₜf
∂[integrator_comps, Δt_comps] = ∂Δtₜf
if free_time
∂[integrator_comps, Δt_comps] = ∂Δtₜf
end
else
error("integrator type not supported: $(typeof(integrator))")
end
Expand All @@ -147,10 +169,16 @@ function QuantumDynamics(
μ∂²[1:2traj.dim, 1:2traj.dim] = sparse(μ∂²P(zₜ, zₜ₊₁, μₜ[integrator_comps]))

else
x_comps, u_comps, Δt_comps = comps(integrator, traj)

μ∂uₜ∂xₜf, μ∂²uₜf, μ∂Δtₜ∂xₜf, μ∂Δtₜ∂uₜf, μ∂²Δtₜf, μ∂xₜ₊₁∂uₜf, μ∂xₜ₊₁∂Δtₜf =
hessian_of_the_lagrangian(integrator, zₜ, zₜ₊₁, μₜ[integrator_comps], traj)
if free_time
x_comps, u_comps, Δt_comps = comps(integrator, traj)
μ∂uₜ∂xₜf, μ∂²uₜf, μ∂Δtₜ∂xₜf, μ∂Δtₜ∂uₜf, μ∂²Δtₜf, μ∂xₜ₊₁∂uₜf, μ∂xₜ₊₁∂Δtₜf =
hessian_of_the_lagrangian(integrator, zₜ, zₜ₊₁, μₜ[integrator_comps], traj)
else
x_comps, u_comps = comps(integrator, traj)
μ∂uₜ∂xₜf, μ∂²uₜf, μ∂xₜ₊₁∂uₜf =
hessian_of_the_lagrangian(integrator, zₜ, zₜ₊₁, μₜ[integrator_comps], traj)
end

if u_comps isa Tuple
for (uᵢ_comps, μ∂uₜᵢ∂xₜf) zip(u_comps, μ∂uₜ∂xₜf)
Expand All @@ -159,32 +187,37 @@ function QuantumDynamics(
for (uᵢ_comps, μ∂²uₜᵢf) zip(u_comps, μ∂²uₜf)
μ∂²[uᵢ_comps, uᵢ_comps] += μ∂²uₜᵢf
end
for (uᵢ_comps, μ∂Δtₜ∂uₜᵢf) zip(u_comps, μ∂Δtₜ∂uₜf)
μ∂²[uᵢ_comps, Δt_comps] += μ∂Δtₜ∂uₜᵢf
if free_time
for (uᵢ_comps, μ∂Δtₜ∂uₜᵢf) zip(u_comps, μ∂Δtₜ∂uₜf)
μ∂²[uᵢ_comps, Δt_comps] += μ∂Δtₜ∂uₜᵢf
end
end
for (uᵢ_comps, μ∂xₜ₊₁∂uₜᵢf) zip(u_comps, μ∂xₜ₊₁∂uₜf)
μ∂²[uᵢ_comps, x_comps .+ traj.dim] += μ∂xₜ₊₁∂uₜᵢf
end
else
μ∂²[x_comps, u_comps] += μ∂uₜ∂xₜf
μ∂²[u_comps, u_comps] += μ∂²uₜf
μ∂²[u_comps, Δt_comps] += μ∂Δtₜ∂uₜf
if free_time
μ∂²[u_comps, Δt_comps] += μ∂Δtₜ∂uₜf
end
μ∂²[u_comps, x_comps .+ traj.dim] += μ∂xₜ₊₁∂uₜf
end

μ∂²[x_comps, Δt_comps] += μ∂Δtₜ∂xₜf
μ∂²[Δt_comps, x_comps .+ traj.dim] += μ∂xₜ₊₁∂Δtₜf
μ∂²[Δt_comps, Δt_comps] .+= μ∂²Δtₜf
if free_time
μ∂²[x_comps, Δt_comps] += μ∂Δtₜ∂xₜf
μ∂²[Δt_comps, x_comps .+ traj.dim] += μ∂xₜ₊₁∂Δtₜf
μ∂²[Δt_comps, Δt_comps] .+= μ∂²Δtₜf
end
end

elseif integrator isa DerivativeIntegrator
if free_time
x_comps, dx_comps, Δt_comps = comps(integrator, traj)

x_comps, dx_comps, Δt_comps = comps(integrator, traj)

μ∂dxₜ∂Δtₜf = -μₜ[integrator_comps]

μ∂²[dx_comps, Δt_comps] += μ∂dxₜ∂Δtₜf
μ∂dxₜ∂Δtₜf = -μₜ[integrator_comps]

μ∂²[dx_comps, Δt_comps] += μ∂dxₜ∂Δtₜf
end
end
end

Expand Down
Loading

0 comments on commit 2b491ba

Please sign in to comment.