From 1b4dae63252be4b87cc7b2ac6b19d015d360c98b Mon Sep 17 00:00:00 2001
From: Aaron Trowbridge <aaron.j.trowbridge@gmail.com>
Date: Thu, 7 Nov 2024 17:17:09 -0500
Subject: [PATCH] density operators update

---
 src/isomorphisms.jl    |  1 +
 src/quantum_systems.jl | 51 +++++++++++++++++++++++++++++++++++++++---
 2 files changed, 49 insertions(+), 3 deletions(-)

diff --git a/src/isomorphisms.jl b/src/isomorphisms.jl
index 72c001d..0beebf9 100644
--- a/src/isomorphisms.jl
+++ b/src/isomorphisms.jl
@@ -15,6 +15,7 @@ export ad_vec
 export ⊗
 
 using LinearAlgebra
+using SparseArrays
 using TestItemRunner
 
 ⊗(xs::AbstractVecOrMat...) = kron(xs...)
diff --git a/src/quantum_systems.jl b/src/quantum_systems.jl
index 43c3335..01755bf 100644
--- a/src/quantum_systems.jl
+++ b/src/quantum_systems.jl
@@ -112,16 +112,42 @@ function QuantumSystem(H::Function, n_drives::Int)
 end
 
 
+function L_function(Ls::AbstractVector{<:AbstractMatrix})
+    return sum([conj(L) ⊗ L - 1 / 2 * ad_vec(L'L, anti=true) for L in Ls])
+end
+
+function QuantumSystem(
+    H_drift::AbstractMatrix,
+    H_drives::Vector{<:AbstractMatrix},
+    dissipation_operators::Vector{<:AbstractMatrix}
+)
+    H_drift = sparse(H_drift)
+    H_drives = sparse.(H_drives)
+
+    H = a -> H_drift + sum(a .* H_drives)
+
+    𝒟̃ = sparse(iso(L_function(dissipation_operators)))
+
+    G = a -> Isomorphisms.G(ad_vec(H(a))) + 𝒟̃
+
+    ∂Gs = Isomorphisms.G.(ad_vec.(H_drives))
+    ∂G = a -> ∂Gs
+
+    levels = size(H_drift, 1)
 
-function Base.copy(sys::QuantumSystem)
     return QuantumSystem(
-        copy(sys.H_drift),
-        copy.(sys.H_drives)
+        H,
+        G,
+        ∂G,
+        levels,
+        length(H_drives)
     )
+
 end
 
 
 
+
 # ============================================================================= #
 
 @testitem "System creation" begin
@@ -132,5 +158,24 @@ end
     system = QuantumSystem(H_drift, H_drives)
 end
 
+@testitem "System creation with dissipation" begin
+    H_drift = GATES[:Z]
+    H_drives = [GATES[:X], GATES[:Y]]
+    dissipation_operators = [GATES[:Z], GATES[:X]]
+
+    system = QuantumSystem(H_drift, H_drives, dissipation_operators)
+
+    # test jacobians
+    a = randn(system.n_drives)
+    ∂G = system.∂G(a)
+    @test length(∂G) == system.n_drives
+    @test all(∂G .≈ QuantumSystems.generator_jacobian(system.G)(a))
+end
+
+
+
+
+
+
 
 end