part_iii solution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Piccolo\n",
    "\n",
    "using CairoMakie\n",
    "using ForwardDiff\n",
    "using LinearAlgebra\n",
    "using Random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example III.1\n",
    "-----\n",
    "**How to build a quantum control problem**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every quantum control problem requires a quantum system describing the dynamics, a goal, and time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = QuantumSystem(0.01 * PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])\n",
    "U_goal = GATES[:X]\n",
    "\n",
    "## Number of timesteps\n",
    "T = 50\n",
    "\n",
    "## Duration of timestep\n",
    "Δt = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NamedTrajectories..jl** stores the problem data. These are containers that have a lot of fields to help with this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NamedTrajectory |> fieldnames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are still coming up with nice utilities to make named trajectories easy to build. Most of the time, this process is hidden from the user. We will use this utility for now, just to see how it might be done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_drives = length(system.H_drives)\n",
    "a_bounds = fill(1.0, n_drives)\n",
    "da_bounds = fill(1.0, n_drives)\n",
    "dda_bounds = fill(1.0, n_drives)\n",
    "\n",
    "## This will help us initialize the trajectory\n",
    "traj = initialize_unitary_trajectory(\n",
    "    U_goal,\n",
    "    T,\n",
    "    Δt,\n",
    "    n_drives,\n",
    "    (a = a_bounds, da = da_bounds, dda = dda_bounds)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main thing a NamedTrajectory supports is indexing by symbols and plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Inspect the control\n",
    "traj.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(traj, [:a, :Ũ⃗])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trajectories also enable simple construction of objectives.\n",
    "\n",
    "Notice that we are using the iso_vec_operator symbol in the objective below here. This is because named trajectories store floating point numbers for each timestep. We need to turn our complex numbers to real numbers, and then vectorize. Type `U\\tilde\\vec` and autocomplete will help you to get the symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = 100.0\n",
    "## Notice that we are using the iso_vec_operator symbol\n",
    "J = UnitaryInfidelityObjective(:Ũ⃗, traj, Q)\n",
    "\n",
    "## Loss functions are evaluated on trajectories\n",
    "Z⃗ = vec(traj)\n",
    "J.L(Z⃗, traj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at that! The infidelity loss is already zero. Are we done? No! We initialized our trajectory with good values for the states, but the controls and states are not dynamically feasible (the dynamics constraint is primal infeasible).\n",
    "\n",
    "Let's add some regularization to our objective because we know that is important."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 1e-2\n",
    "J += QuadraticRegularizer(:a, traj, R)\n",
    "J += QuadraticRegularizer(:da, traj, R)\n",
    "J += QuadraticRegularizer(:dda, traj, R)\n",
    "\n",
    "J.L(Z⃗, traj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the dynamics constraints, which are implemented using integrators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Integrators\n",
    "integrators = [\n",
    "    UnitaryPadeIntegrator(system, :Ũ⃗, :a, traj),\n",
    "    DerivativeIntegrator(:a, :da, traj),\n",
    "    DerivativeIntegrator(:da, :dda, traj)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Miscellaneous settings\n",
    "Random.seed!(1234)\n",
    "ipopt_options = IpoptOptions(print_level=4, max_iter=50, recalc_y=\"yes\", recalc_y_feas_tol=1e-2)\n",
    "\n",
    "prob = QuantumControlProblem(\n",
    "    system,\n",
    "    traj,\n",
    "    J,\n",
    "    integrators,\n",
    "    ipopt_options=ipopt_options,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve!(prob)\n",
    "println(\"Unitary fidelity: \", unitary_fidelity(prob))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(prob.trajectory, [:a, :Ũ⃗])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, it is much easier to just call our problem templates! Everything we just did---and more---happens in the following 6 lines.\n",
    "\n",
    "```Julia\n",
    "system = QuantumSystem(0.01 * PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])\n",
    "U_goal = GATES[:X]\n",
    "T = 50\n",
    "Δt = 0.2\n",
    "problem = UnitarySmoothPulseProblem(system, U_goal, T, Δt)\n",
    "solve!(problem)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "-----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise III.1\n",
    "**Inspect a gradient for correctness**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is often the case that after we work very hard to solve for an analytic gradient, we need to make sure we did a good job coding it up. Not only do we need to profile the calculation to make sure our code is appropriately efficient, but we also need to be sure that our gradient was correct to begin with.\n",
    "\n",
    "The next cell defines the following objective:\n",
    "\\begin{equation}\n",
    "    J(\\vec{\\mathbf{Z}}) = ||\\mathbf{a}_\\text{1:T/2} - \\mathbf{a}_\\text{T:T/2+1}||_2^2\n",
    "\\end{equation}\n",
    "This is a cost that forces symmetric controls, which could be very nice to have for certain problems. Unfortunately, we made a mistake when writing our gradient.\n",
    "\n",
    "Create a test trajectory, then use **ForwardDiff.jl** to show the analytic gradient is wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function CustomObjective(;\n",
    "\tname::Union{Nothing, Symbol}=nothing,\n",
    "\ttimes::Union{Nothing, AbstractVector{Int}}=nothing,\n",
    "\tR::Union{Nothing, AbstractVector{<:Real}}=nothing,\n",
    ")\n",
    "    @assert !isnothing(name) \"name must be specified\"\n",
    "    @assert !isnothing(times) \"times must be specified\"\n",
    "    @assert !isnothing(R) \"R must be specified\"\n",
    "\n",
    "    params = Dict(\n",
    "        :type => :CustomObjective,\n",
    "        :name => name,\n",
    "        :times => times,\n",
    "        :R => R,\n",
    "    )\n",
    "\n",
    "    @views function L(Z⃗::AbstractVector{<:Real}, Z::NamedTrajectory)\n",
    "        J = 0.0\n",
    "        for (t₁, t₂) ∈ zip(times[1:end÷2], reverse(times[end÷2 + 1:end]))\n",
    "            rₜ₁ = Z⃗[slice(t₁, Z.components[name], Z.dim)]\n",
    "            rₜ₂ = Z⃗[slice(t₂, Z.components[name], Z.dim)]\n",
    "            rₜ = rₜ₁ .- rₜ₂\n",
    "            J += 0.5 * rₜ' * (R .* rₜ)\n",
    "        end\n",
    "        return J\n",
    "    end\n",
    "\n",
    "    @views function ∇L(Z⃗::AbstractVector{<:Real}, Z::NamedTrajectory)\n",
    "        ∇ = zeros(length(Z⃗))\n",
    "        for t ∈ times\n",
    "            rₜ_slice = slice(t, Z.components[name], Z.dim)\n",
    "            rₜ = Z⃗[rₜ_slice] \n",
    "            ∇[rₜ_slice] .= R .* rₜ\n",
    "        end\n",
    "        return ∇\n",
    "    end\n",
    "\n",
    "    ∂²L = nothing\n",
    "    ∂²L_structure = nothing\n",
    "\n",
    "    return Objective(L, ∇L, ∂²L, ∂²L_structure, Dict[params])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- # \n",
    "## Might help to lower the timesteps\n",
    "T_test = T\n",
    "## Useful a_guess to initialize the trajectory\n",
    "a_symmetric = ones(n_drives, T_test)\n",
    "a_asymmetric = stack(repeat([range(0.0, 1.0, length=T_test)], n_drives), dims=1)\n",
    "# ----- #\n",
    "\n",
    "## initialize a random trajectory\n",
    "n_drives = 2\n",
    "a_bounds = fill(1.0, n_drives)\n",
    "da_bounds = fill(1.0, n_drives)\n",
    "dda_bounds = fill(1.0, n_drives)\n",
    "\n",
    "Z = initialize_unitary_trajectory(\n",
    "    U_goal,\n",
    "    T_test,\n",
    "    Δt,\n",
    "    n_drives,\n",
    "    (a = a_bounds, da = da_bounds, dda = dda_bounds),\n",
    "    a_guess = a_asymmetric,\n",
    "    system=system\n",
    ")\n",
    "\n",
    "Z⃗ = vec(Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Construct the custom objective\n",
    "obj = CustomObjective(;\n",
    "\tname=:a,\n",
    "\ttimes=1:T_test,\n",
    "\tR=fill(1.0, n_drives)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = Figure()\n",
    "ax = Axis(fig[1, 1])\n",
    "plot!(ax, ForwardDiff.gradient(Z⃗ -> obj.L(Z⃗, Z), Z⃗), label=\"∇L (ForwardDiff)\")\n",
    "plot!(obj.∇L(Z⃗, Z), label=\"∇L (Custom)\")\n",
    "axislegend(ax, position=:rb)\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Write the correct gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Someone should open up a PR for this nice new objective.** Symmetric controls would be nice to have!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise III.2\n",
    "**The problem template grand tour**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constrained optimization gives **Piccolo.jl** the ability to solve many different kinds of problems. These next few examples will walk you through implementing many of these types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## This should be familiar by now.\n",
    "Δ = 0.1\n",
    "system = system = QuantumSystem(Δ * PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])\n",
    "U_goal = GATES[:X]\n",
    "\n",
    "T = 50\n",
    "Δt = 0.2\n",
    "\n",
    "## Make sure to pass any provided options dictionaries to each problem! \n",
    "## ipopt_options = ...\n",
    "## piccolo_options = ...\n",
    "\n",
    "## Anything enclosed by `# ----- #` should remain untouched. Use the problem template to pass the test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smooth pulse problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- #\n",
    "Random.seed!(42);\n",
    "ipopt_options = IpoptOptions(print_level=1, max_iter=50)\n",
    "piccolo_options = PiccoloOptions(verbose=false)\n",
    "test_fidelity() = unitary_fidelity(smooth_pulse_prob) > 0.99 ? \"Pass\" : \"Fail\"\n",
    "# ----- #\n",
    "\n",
    "\n",
    "smooth_pulse_prob = UnitarySmoothPulseProblem(\n",
    "    system, U_goal, T, Δt, \n",
    "    ipopt_options=ipopt_options,\n",
    "    piccolo_options=piccolo_options\n",
    ")\n",
    "\n",
    "\n",
    "# ----- #\n",
    "solve!(smooth_pulse_prob)\n",
    "test_fidelity()\n",
    "# ----- #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## View your result\n",
    "plot_unitary_populations(smooth_pulse_prob.trajectory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minimum time problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- #\n",
    "Random.seed!(42);\n",
    "ipopt_options = IpoptOptions(print_level=1, max_iter=50, recalc_y=\"yes\", recalc_y_feas_tol=1e-2)\n",
    "piccolo_options = PiccoloOptions(verbose=false)\n",
    "test_time() = get_duration(mintime_prob.trajectory) < get_duration(smooth_pulse_prob.trajectory) ? \"Pass\" : \"Fail\"\n",
    "# ----- #\n",
    "\n",
    "\n",
    "mintime_prob = UnitaryMinimumTimeProblem(\n",
    "    smooth_pulse_prob, final_fidelity=0.9999,\n",
    "    ipopt_options=ipopt_options,\n",
    "    piccolo_options=piccolo_options,\n",
    ")\n",
    "\n",
    "\n",
    "# ----- #\n",
    "solve!(mintime_prob)\n",
    "test_time()\n",
    "# ----- #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_unitary_populations(mintime_prob.trajectory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bang-bang control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- #\n",
    "Random.seed!(42)\n",
    "piccolo_options = piccolo_options=PiccoloOptions(verbose=false)\n",
    "ipopt_options = ipopt_options=IpoptOptions(print_level=1, max_iter=50)\n",
    "function test_sparsity()\n",
    "    r = 1e-3\n",
    "    if 2sum(bangbang_prob.trajectory.da .> r) < sum(smooth_pulse_prob.trajectory.da .> r)\n",
    "        \"Pass\"\n",
    "    else\n",
    "        \"Fail. Try increasing the bang bang regularization parameter.\"\n",
    "    end\n",
    "end\n",
    "# ----- #\n",
    "\n",
    "\n",
    "bangbang_prob = UnitaryBangBangProblem(\n",
    "    system, U_goal, T, Δt, R_bang_bang=10.0,\n",
    "    ipopt_options=ipopt_options,\n",
    "    piccolo_options=piccolo_options\n",
    ")\n",
    "\n",
    "\n",
    "# ----- #\n",
    "solve!(bangbang_prob)\n",
    "test_sparsity()\n",
    "# ----- #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_unitary_populations(bangbang_prob.trajectory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(bangbang_prob.trajectory, [:da])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling over quantum systems\n",
    "*Solving problems with shared controls.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleable_system(Δ) = QuantumSystem(Δ * PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## It doesn't work if the drift changes\n",
    "unitary_fidelity(smooth_pulse_prob.trajectory, sampleable_system(Δ)) |> println\n",
    "unitary_fidelity(smooth_pulse_prob.trajectory, sampleable_system(0.0)) |> println"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the problem for a bunch of different systems--but make everyone use the same controls!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- # \n",
    "Random.seed!(42)\n",
    "ipopt_options = IpoptOptions(print_level=1, max_iter=50)\n",
    "piccolo_options = PiccoloOptions(verbose=false)\n",
    "function test_robustness()\n",
    "    test_vals = range(-0.2, 0.2, length=10)\n",
    "    test_fn(Δ) = unitary_fidelity(\n",
    "        sample_prob.trajectory, \n",
    "        sampleable_system(Δ), \n",
    "        unitary_name=:Ũ⃗1\n",
    "    )\n",
    "    for Δ in test_vals\n",
    "        if test_fn(Δ) < 0.99\n",
    "            return \"Fail. The fidelity is below 0.99 for Δ = $Δ. Consider decreasing the curvature regularization.\"\n",
    "        end\n",
    "    end\n",
    "    return \"Pass\"\n",
    "end\n",
    "# ----- #\n",
    "\n",
    "\n",
    "sample_prob = UnitarySamplingProblem(\n",
    "    [sampleable_system(-0.2), \n",
    "     sampleable_system(0.1),\n",
    "     sampleable_system(0.2)],\n",
    "    U_goal,\n",
    "    T,\n",
    "    Δt,\n",
    "    R_dda=1e-6,\n",
    "    ipopt_options=ipopt_options,\n",
    "    piccolo_options=piccolo_options\n",
    ")\n",
    "\n",
    "\n",
    "# ----- #\n",
    "solve!(sample_prob)\n",
    "test_robustness()\n",
    "# ----- #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(sample_prob.trajectory, [:a, :Ũ⃗1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sweep_sample(Δ) = unitary_fidelity(sample_prob.trajectory, sampleable_system(Δ), unitary_name=:Ũ⃗1)\n",
    "sweep_default(Δ) = unitary_fidelity(smooth_pulse_prob.trajectory, sampleable_system(Δ))\n",
    "     \n",
    "f = Figure()\n",
    "ax = Axis(f[1,1], yscale=log10, limits=(nothing, (10^-6, 1)))\n",
    "\n",
    "Δs = range(-.2, .2, length=200)\n",
    "lines!(ax, Δs, 1 .- sweep_sample.(Δs), label=\"Robust infidelity\")\n",
    "lines!(ax, Δs, 1 .- sweep_default.(Δs), label=\"Default infidelity\")\n",
    "axislegend(ax, position=:lb)\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling over gates\n",
    "*Solving independent problems, together.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleable_U_goal(θ) = exp(-im * θ * PAULIS[:X])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- #\n",
    "Random.seed!(42)\n",
    "Q_symb = :a\n",
    "ipopt_options = IpoptOptions(print_level=1, max_iter=50)\n",
    "piccolo_options = PiccoloOptions(verbose=false, free_time=false)\n",
    "test_direct_sum() = nothing\n",
    "# ----- #\n",
    "\n",
    "\n",
    "θs = range(0, π, length=5)\n",
    "probs = QuantumControlProblem[]\n",
    "\n",
    "## Create a UnitarySmoothPulseProblem for each θ in θs and store them in probs\n",
    "for θ in θs\n",
    "    prob_θ = UnitarySmoothPulseProblem(\n",
    "        system, sampleable_U_goal(θ), T, Δt,\n",
    "        piccolo_options=piccolo_options,\n",
    "        ipopt_options=ipopt_options\n",
    "    )\n",
    "    solve!(prob_θ)\n",
    "    push!(probs, prob_θ)\n",
    "end\n",
    "\n",
    "## Create a UnitaryDirectSumProblem with the problems in probs \n",
    "direct_sum_prob = UnitaryDirectSumProblem(\n",
    "    probs, 0.9999, Q_symb=Q_symb, \n",
    "    ipopt_options=ipopt_options, \n",
    "    piccolo_options=piccolo_options\n",
    ")\n",
    "\n",
    "\n",
    "# ----- #\n",
    "solve!(direct_sum_prob)\n",
    "test_direct_sum()\n",
    "# ----- #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## After solve, the controls are close by\n",
    "control_symbols = [Symbol(\"a$i\") for i in eachindex(θs)]\n",
    "plot(direct_sum_prob.trajectory, control_symbols, merge_components=true, ignored_lables=control_symbols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EmbeddedOperators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gates are defined in a computational subspace, but the physical system might not care. We need to model this, and adapt our costs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- #\n",
    "a = annihilate(4)\n",
    "TRANSMON = Dict(\n",
    "    :I => Matrix{ComplexF64}(I, size(a)),\n",
    "    :X => (a + a') / 2,\n",
    "    :Y => (a - a') / 2im,\n",
    "    :N => a'a,\n",
    "    :A => a'a'a*a\n",
    ");\n",
    "transmon = QuantumSystem(0.1 * (TRANSMON[:A]), [TRANSMON[:X], TRANSMON[:Y]])\n",
    "test_embedding() = unembed(U_goal_embedded) == GATES[:X] ? \"Pass\" : \"Fail\"\n",
    "# ----- #\n",
    "\n",
    "\n",
    "## create an EmbeddedOperator for GATES[:X]\n",
    "U_goal_embedded = EmbeddedOperator(GATES[:X], transmon)\n",
    "\n",
    "\n",
    "# ----- #\n",
    "test_embedding()\n",
    "# ----- #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Solve the problem with the EmbeddedOperator."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.4",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}