Refactor QuantumSystem and OpenQuantumSystem constructors

Project.toml

@@ -1,7 +1,7 @@
 name = "PiccoloQuantumObjects"
 uuid = "5a402ddf-f93c-42eb-975e-5582dcda653d"
 authors = ["Aaron Trowbridge <[email protected]> and contributors"]
-version = "0.6.1"
+version = "0.7.0"
 
 [deps]
 ExponentialAction = "e24c0720-ea99-47e8-929e-571b494574d3"

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
@@ -8,4 +9,10 @@ PiccoloQuantumObjects = "5a402ddf-f93c-42eb-975e-5582dcda653d"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 
 [compat]
-Literate = "2.7.0"
+CairoMakie = "0.13"
+Documenter = "1.15"
+Literate = "2.20"
+LiveServer = "1.5"
+NamedTrajectories = "0.5"
+PiccoloDocsTemplate = "0.2"
+Revise = "3.11"

docs/literate/quantum_systems.jl

@@ -22,24 +22,30 @@ The [`QuantumSystem`](@ref) type is used to represent a quantum system with a dr
 Hamiltonian and a set of drive Hamiltonians,
 
 ```math
-H = H_{\text{drift}} + \sum_i a_i H_{\text{drives}}^{(i)}
+H(u, t) = H_{\text{drift}} + \sum_i u_i H_{\text{drives}}^{(i)}
 ```
 
+where ``u`` is the control vector and ``t`` is time.
+
 ```@docs; canonical = false
 QuantumSystem
 ```
 
 `QuantumSystem`'s are containers for quantum dynamics. Internally, they compute the
-necessary isomorphisms to perform the dynamics in a real vector space.
+necessary isomorphisms to perform the dynamics in a real vector space. All systems
+require explicit specification of `T_max` (maximum time) and `drive_bounds` (control bounds).
 
 =#
 
 H_drift = PAULIS[:Z]
 H_drives = [PAULIS[:X], PAULIS[:Y]]
-system = QuantumSystem(H_drift, H_drives)
+T_max = 10.0
+drive_bounds = [(-1.0, 1.0), (-1.0, 1.0)]
+system = QuantumSystem(H_drift, H_drives, T_max, drive_bounds)
 
-a_drives = [1, 0]
-system.H(a_drives)
+u_controls = [1.0, 0.0]
+t = 0.0
+system.H(u_controls, t)
 
 #=
 To extract the drift and drive Hamiltonians from a `QuantumSystem`, use the 
@@ -58,27 +64,13 @@ drives[2] |> sparse
 
 #=
 !!! note
-    We can also construct a `QuantumSystem` directly from a Hamiltonian function. Internally,
-    `ForwardDiff.jl` is used to compute the drives.
+    We can also construct a `QuantumSystem` directly from a Hamiltonian function.
+    The function must accept `(u, t)` arguments where `u` is the control vector and `t` is time.
 =#
 
-H(a) = PAULIS[:Z] + a[1] * PAULIS[:X] + a[2] * PAULIS[:Y]
-system = QuantumSystem(H, 2)
-get_drives(system)[1] |> sparse
-
-# _Create a noise model with a confusion matrix._
-function H(a; C::Matrix{Float64}=[1.0 0.0; 0.0 1.0])
-    b = C * a
-    return b[1] * PAULIS.X + b[2] * PAULIS.Y
-end
-
-C_matrix = [0.99 0.01; -0.01 1.01]
-system = QuantumSystem(a -> H(a, C=C_matrix), 2; params=Dict(:C => C_matrix))
-confused_drives = get_drives(system)
-confused_drives[1] |> sparse
-
-# 
-confused_drives[2] |> sparse
+H(u, t) = PAULIS[:Z] + u[1] * PAULIS[:X] + u[2] * PAULIS[:Y]
+system = QuantumSystem(H, 10.0, [(-1.0, 1.0), (-1.0, 1.0)])
+system.H([1.0, 0.0], 0.0) |> sparse
 
 #=
 ## Open quantum systems
@@ -95,7 +87,9 @@ OpenQuantumSystem
 H_drives = [PAULIS[:X]]
 a = annihilate(2)
 dissipation_operators = [a'a, a]
-system = OpenQuantumSystem(H_drives, dissipation_operators=dissipation_operators)
+T_max = 10.0
+drive_bounds = [(-1.0, 1.0)]
+system = OpenQuantumSystem(H_drives, T_max, drive_bounds, dissipation_operators=dissipation_operators)
 system.dissipation_operators[1] |> sparse
 
 # 
@@ -110,41 +104,6 @@ system.dissipation_operators[2] |> sparse
 
 get_drift(system) |> sparse
 
-#=
-## Time Dependent Quantum Systems
-A [`TimeDependentQuantumSystem`](@ref) is a `QuantumSystem` with time-dependent Hamiltonians.
-```@docs; canonical = false
-TimeDependentQuantumSystem
-```
-
-A function `H(a, t)` or carrier and phase kwargs are used to specify time-dependent drives,
-```math
-    H(a, t) = H_{\text{drift}} + \sum_i a_i \cos(\omega_i t + \phi_i) H_{\text{drives}}^{(i)}
-```
-=#
-# _Create a time-dependent Hamiltonian with a time-dependent drive._
-H(a, t) = PAULIS.Z + a[1] * cos(t) * PAULIS.X
-system = TimeDependentQuantumSystem(H, 1)
-
-# _The drift Hamiltonian is the Z operator, but its now a function of time!_
-get_drift(system)(0.0) |> sparse
-
-# _The drive Hamiltonian is the X operator, but its now a function of time!_
-get_drives(system)[1](0.0) |> sparse
-
-# _Change the time to π._
-get_drives(system)[1](π) |> sparse
-
-# _Similar matrix constructors exist, but with carrier and phase kwargs._
-system = TimeDependentQuantumSystem(PAULIS.Z, [PAULIS.X], carriers=[1.0], phases=[0.0])
-
-# _This is the same as before, t=0.0:_
-get_drives(system)[1](0.0) |> sparse
-
-# _and at π:_
-get_drives(system)[1](π) |> sparse
-
-
 #=
 ## Composite quantum systems
 
@@ -155,10 +114,10 @@ lifted to the full Hilbert space.
 
 =#
 
-system_1 = QuantumSystem([PAULIS[:X]])
-system_2 = QuantumSystem([PAULIS[:Y]])
+system_1 = QuantumSystem([PAULIS[:X]], 1.0, [(-1.0, 1.0)])
+system_2 = QuantumSystem([PAULIS[:Y]], 1.0, [(-1.0, 1.0)])
 H_drift = PAULIS[:Z] ⊗ PAULIS[:Z]
-system = CompositeQuantumSystem(H_drift, [system_1, system_2]);
+system = CompositeQuantumSystem(H_drift, Matrix{ComplexF64}[], [system_1, system_2], 1.0, Float64[]);
 
 # _The drift Hamiltonian is the ZZ coupling._
 get_drift(system) |> sparse
@@ -214,11 +173,11 @@ is_reachable
 =#
 
 # _Y can be reached by commuting Z and X._
-system = QuantumSystem(PAULIS[:Z], [PAULIS[:X]])
+system = QuantumSystem(PAULIS[:Z], [PAULIS[:X]], 1.0, [(-1.0, 1.0)])
 is_reachable(PAULIS[:Y], system)
 
 # _Y cannot be reached by X alone._
-system = QuantumSystem([PAULIS[:X]])
+system = QuantumSystem([PAULIS[:X]], 1.0, [(-1.0, 1.0)])
 is_reachable(PAULIS[:Y], system)
 
 #=
@@ -231,7 +190,7 @@ direct_sum
 =#
 
 # _Create a pair of non-interacting qubits._
-system_1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
-system_2 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
+system_1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
+system_2 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
 system = direct_sum(system_1, system_2)
 get_drift(system) |> sparse

docs/literate/quickstart.jl

@@ -0,0 +1,143 @@
+# # Quickstart
+
+# This quickstart guide will help you get up and running with PiccoloQuantumObjects.jl.
+
+# ## Installation
+
+# ```julia
+# using Pkg
+# Pkg.add("PiccoloQuantumObjects")
+# ```
+
+# ## Basic Usage
+
+using PiccoloQuantumObjects
+using LinearAlgebra
+using SparseArrays
+using NamedTrajectories
+
+# ### Creating a Quantum System
+
+# Define Hamiltonian components
+H_drift = PAULIS[:Z]
+H_drives = [PAULIS[:X], PAULIS[:Y]]
+
+# Specify time and control bounds
+T_max = 10.0  # Maximum evolution time
+drive_bounds = [(-1.0, 1.0), (-1.0, 1.0)]  # Control bounds for each drive
+
+# Create the quantum system
+system = QuantumSystem(H_drift, H_drives, T_max, drive_bounds)
+
+# ### Working with Quantum States
+
+# Create quantum states
+ψ_ground = ket_from_string("g", [2])  # Ground state
+
+#
+ψ_excited = ket_from_string("e", [2])  # Excited state
+
+# Calculate fidelity
+fidelity(ψ_ground, ψ_excited)
+
+# ### Simulating Evolution (Rollouts)
+
+# Initial state
+ψ_init = ComplexF64[1.0, 0.0]
+
+# Define control sequence
+T = 10  # Number of time steps
+controls = rand(2, T)  # Random controls for two drives
+Δt = fill(0.1, T)  # Time step duration
+
+# Perform rollout
+ψ̃_rollout = rollout(ψ_init, controls, Δt, system)
+
+# Check final state fidelity
+ψ_goal = ComplexF64[0.0, 1.0]
+rollout_fidelity(ψ_init, ψ_goal, controls, Δt, system)
+
+# ### Open Quantum Systems
+
+# Add dissipation operators
+a = annihilate(2)
+dissipation_operators = [a'a, a]
+
+# Create open quantum system
+open_system = OpenQuantumSystem(
+    H_drives, 
+    T_max, 
+    drive_bounds, 
+    dissipation_operators=dissipation_operators
+)
+
+# ### Composite Systems
+
+# Create subsystems
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [(-1.0, 1.0)])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [(-1.0, 1.0)])
+
+# Define coupling
+H_coupling = 0.1 * kron(PAULIS[:Z], PAULIS[:Z])
+
+# Create composite system
+composite_sys = CompositeQuantumSystem(
+    H_coupling, 
+    Matrix{ComplexF64}[], 
+    [sys1, sys2], 
+    10.0, 
+    Float64[]
+)
+
+# ## Visualization
+
+# PiccoloQuantumObjects.jl integrates with NamedTrajectories.jl for plotting trajectories.
+
+using CairoMakie
+
+# ### Plotting Controls and States
+
+# Create a trajectory with controls and states
+T_plot = 50
+controls_plot = 0.5 * sin.(2π * (1:T_plot) / T_plot)
+Δt_plot = fill(T_max / T_plot, T_plot)
+
+# Perform a rollout to get the state evolution
+ψ̃_traj = rollout(ψ_init, hcat(controls_plot, -controls_plot)', Δt_plot, system)
+
+# Create a NamedTrajectory for plotting
+traj = NamedTrajectory(
+    (
+        ψ̃ = ψ̃_traj,
+        a = hcat(controls_plot, -controls_plot)',
+        Δt = Δt_plot
+    );
+    timestep=:Δt,
+    controls=:a,
+    initial=(ψ̃ = ket_to_iso(ψ_init),),
+    goal=(ψ̃ = ket_to_iso(ψ_goal),)
+)
+
+# Plot the trajectory 
+plot(traj)
+
+# ### Plotting State Populations
+
+# Use transformations to plot populations directly from isomorphic states
+plot(
+    traj, 
+    [:a],
+    transformations=[
+        :ψ̃ => (ψ̃ -> abs2.(iso_to_ket(ψ̃)))
+    ],
+    transformation_labels=["Populations"]
+)
+
+#=
+## Next Steps
+
+- Explore the [Quantum Systems](@ref) manual for detailed system construction
+- Learn about [Quantum Objects](@ref) for working with states and operators
+- See [Rollouts](@ref) for simulation and fidelity calculations
+- Understand [Isomorphisms](@ref) for the underlying mathematical transformations
+=#

docs/literate/rollouts.jl

@@ -81,7 +81,7 @@ T = 10
 ψ_init = ComplexF64[1.0, 0.0]
 controls = rand(2, T)
 Δt = fill(0.1, T)
-system = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
+system = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
 ψ̃_rollout = rollout(ψ_init, controls, Δt, system) 
 ψ̃_rollout |> size
 

docs/make.jl

@@ -3,11 +3,12 @@ using PiccoloDocsTemplate
 
 pages = [
     "Home" => "index.md",
+    "Quickstart" => "generated/quickstart.md",
     "Manual" => [
-        "Isomorphisms" => "generated/isomorphisms.md",
-        "Quantum Objects" => "generated/quantum_objects.md",
         "Quantum Systems" => "generated/quantum_systems.md",
+        "Quantum Objects" => "generated/quantum_objects.md",
         "Rollouts" => "generated/rollouts.md",
+        "Isomorphisms" => "generated/isomorphisms.md",
     ],
     "Library" => "lib.md",
 ]

ext/QuantumToolboxExt.jl

@@ -71,10 +71,10 @@ get_c_ops(sys::OpenQuantumSystem) = Qobj.(sys.dissipation_operators)
     # X gate
     traj = named_trajectory_type_2()
 
-    sys = QuantumSystem([PAULIS.X, PAULIS.Y])
+    sys = QuantumSystem([PAULIS.X, PAULIS.Y], 3.92, [(-1.0, 1.0), (-1.0, 1.0)])
 
     open_sys = OpenQuantumSystem(
-        [PAULIS.X, PAULIS.Y], dissipation_operators=[annihilate(2)]
+        [PAULIS.X, PAULIS.Y], 1.0, [1.0, 1.0], dissipation_operators=[annihilate(2)]
     )
 
     ψ0 = Qobj([1; 0])