✨ add + method for Dense and SparseHamiltonian (#111)

Project.toml

@@ -1,7 +1,7 @@
 name = "OpenQuantumBase"
 uuid = "5dd19120-8766-11e9-1ef9-27094038a7db"
 authors = ["neversakura <[email protected]>"]
-version = "0.7.7"
+version = "0.7.8"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/hamiltonian/dense_hamiltonian.jl

@@ -40,6 +40,18 @@ function DenseHamiltonian(funcs, mats; unit=:h, dimensionless_time=true)
     DenseHamiltonian{eltype(mats[1]),dimensionless_time}(funcs, mats, cache, hsize)
 end
 
+function Base.:+(h1::DenseHamiltonian, h2::DenseHamiltonian)
+    @assert size(h1) == size(h2) "The two Hamiltonians need to have the same size."
+    @assert isdimensionlesstime(h1) == isdimensionlesstime(h2) "The two Hamiltonians need to have the time arguments."
+
+    (m1, m2) = promote(h1.m, h2.m)
+    cache = similar(m1[1])
+    mats = [m1; m2]
+    funcs = [h1.f; h2.f]
+    hsize = size(h1)
+    DenseHamiltonian{eltype(m1[1]),isdimensionlesstime(h1)}(funcs, mats, cache, hsize)
+end
+
 """
     function (h::DenseHamiltonian)(s::Real)
 

src/hamiltonian/sparse_hamiltonian.jl

@@ -36,6 +36,18 @@ function SparseHamiltonian(funcs, mats; unit=:h, dimensionless_time=true)
     SparseHamiltonian{eltype(mats[1]),dimensionless_time}(funcs, mats, cache, size(mats[1]))
 end
 
+function Base.:+(h1::SparseHamiltonian, h2::SparseHamiltonian)
+    @assert size(h1) == size(h2) "The two Hamiltonians need to have the same size."
+    @assert isdimensionlesstime(h1) == isdimensionlesstime(h2) "The two Hamiltonians need to have the time arguments."
+
+    (m1, m2) = promote(h1.m, h2.m)
+    mats = [m1; m2]
+    cache = similar(sum(mats))
+    funcs = [h1.f; h2.f]
+    hsize = size(h1)
+    SparseHamiltonian{eltype(m1[1]),isdimensionlesstime(h1)}(funcs, mats, cache, hsize)
+end
+
 isdimensionlesstime(::SparseHamiltonian{T,B}) where {T,B} = B
 issparse(::SparseHamiltonian) = true
 

test/hamiltonian/dense_hamiltonian.jl

@@ -13,35 +13,41 @@ Hc = H |> get_cache
 @test !isconstant(H)
 @test isdimensionlesstime(H)
 
+H1 = DenseHamiltonian([A], [σx])
+H2 = DenseHamiltonian([B], [σz])
+H3 = H1 + H2
+@test H3.m == H.m
+@test H3.f == H.f
+
 # update_cache method
 C = similar(σz)
 update_cache!(C, H, 10, 0.5)
 @test C == -1im * π * (σx + σz)
 
 # update_vectorized_cache method
 C = get_cache(H)
-C = C⊗C
+C = C ⊗ C
 update_vectorized_cache!(C, H, 10, 0.5)
 temp = -1im * π * (σx + σz)
 @test C == σi ⊗ temp - transpose(temp) ⊗ σi
 
 # in-place update for matrices
-du = [1.0 + 0.0im 0; 0 0]
-ρ  = PauliVec[1][1] * PauliVec[1][1]'
+du = [1.0+0.0im 0; 0 0]
+ρ = PauliVec[1][1] * PauliVec[1][1]'
 H(du, ρ, 2, 0.5)
 @test du ≈ -1.0im * π * ((σx + σz) * ρ - ρ * (σx + σz))
 
 # eigen-decomposition
 w, v = eigen_decomp(H, 0.5)
 @test w ≈ [-1, 1] / √2
-@test abs(v[:, 1]'*[1-sqrt(2), 1] / sqrt(4-2*sqrt(2))) ≈ 1
-@test abs(v[:, 2]'*[1+sqrt(2), 1] / sqrt(4+2*sqrt(2))) ≈ 1
+@test abs(v[:, 1]' * [1 - sqrt(2), 1] / sqrt(4 - 2 * sqrt(2))) ≈ 1
+@test abs(v[:, 2]' * [1 + sqrt(2), 1] / sqrt(4 + 2 * sqrt(2))) ≈ 1
 
-Hrot= rotate(H, v)
+Hrot = rotate(H, v)
 @test evaluate(Hrot, 0.5) ≈ [-1 0; 0 1] / sqrt(2)
 
 # error message test
-@test_throws ArgumentError DenseHamiltonian([(s)->1-s, (s)->s], [σx, σz], unit=:hh)
+@test_throws ArgumentError DenseHamiltonian([(s) -> 1 - s, (s) -> s], [σx, σz], unit=:hh)
 
 Hnd = DenseHamiltonian([A, B], [σx, σz], dimensionless_time=false)
 @test !isdimensionlesstime(Hnd)
@@ -66,22 +72,22 @@ update_cache!(Hstc, Hst, 10, 0.5)
 @test Hstc == -0.5im * (σx + σz)
 
 # update_vectorized_cache method
-C = Hstc⊗Hstc
+C = Hstc ⊗ Hstc
 update_vectorized_cache!(C, Hst, 10, 0.5)
 temp = -0.5im * (σx + σz)
 @test C == σi ⊗ temp - transpose(temp) ⊗ σi
 
 # in-place update for matrices
-du = [1.0 + 0.0im 0; 0 0]
-ρ  = PauliVec[1][1] * PauliVec[1][1]'
+du = [1.0+0.0im 0; 0 0]
+ρ = PauliVec[1][1] * PauliVec[1][1]'
 Hst(du, ρ, 2, 0.5)
 @test du ≈ -0.5im * ((σx + σz) * ρ - ρ * (σx + σz))
 
 # eigen-decomposition
 w, v = eigen_decomp(Hst, 0.5)
-@test 2π*w ≈ [-1, 1] / √2
-@test abs(v[:, 1]'*[1-sqrt(2), 1] / sqrt(4-2*sqrt(2))) ≈ 1
-@test abs(v[:, 2]'*[1+sqrt(2), 1] / sqrt(4+2*sqrt(2))) ≈ 1
+@test 2π * w ≈ [-1, 1] / √2
+@test abs(v[:, 1]' * [1 - sqrt(2), 1] / sqrt(4 - 2 * sqrt(2))) ≈ 1
+@test abs(v[:, 2]' * [1 + sqrt(2), 1] / sqrt(4 + 2 * sqrt(2))) ≈ 1
 
-Hrot= rotate(Hst, v)
-@test 2π*evaluate(Hrot, 0.5) ≈ [-1 0; 0 1] / sqrt(2)
+Hrot = rotate(Hst, v)
+@test 2π * evaluate(Hrot, 0.5) ≈ [-1 0; 0 1] / sqrt(2)

test/hamiltonian/sparse_hamiltonian.jl

@@ -10,6 +10,12 @@ u = [1.0 + 0.0im, 1] / sqrt(2)
 H_sparse = SparseHamiltonian([A, B], [spσx, spσz])
 @test isdimensionlesstime(H_sparse)
 
+Hs1 = SparseHamiltonian([A], [spσx])
+Hs2 = SparseHamiltonian([B], [spσz])
+Hs3 = Hs1 + Hs2
+@test Hs3.m == H_sparse.m
+@test Hs3.f == H_sparse.f
+
 H_real = convert(Real, H_sparse)
 @test eltype(H_real) <: Real
 @test H_sparse(0.0) ≈ H_real(0.0)
@@ -47,7 +53,7 @@ vf = [0, 0, 0, 1.0]
 @test abs(v[1, 2]) ≈ 1
 
 # ## Test suite for eigen decomposition of `SparseHamiltonian`
-np = 5 
+np = 5
 Hd = standard_driver(np, sp=true)
 Hp = alt_sec_chain(1, 0.5, 1, np, sp=true)
 H_sparse = SparseHamiltonian([A, B], [Hd, Hp], unit=:ħ)
@@ -56,6 +62,6 @@ wl, vl = haml_eigs(H_sparse, 0.5, 3)
 @test w1[1:3] ≈ wl
 @test abs.(v1[:, 1:3]' * vl) |> Diagonal ≈ I
 
-w2, v2 = haml_eigs(H_sparse, 0.5, 3, lobpcg = false)
+w2, v2 = haml_eigs(H_sparse, 0.5, 3, lobpcg=false)
 @test w1 ≈ w2
 @test v1 ≈ v2