Add find_degenerate function and fix the Hamiltonian interface (#94)

* 🐛`Hamiltonian` interface does not try to promote the elements of `mats`

* 🎨new `find_degenerate` function
1. Round the energies with `digits` option instead of error tolerance
2. export `find_degenerate` function

* 🔖tag version v0.7.2

Project.toml

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

src/OpenQuantumBase.jl

@@ -114,7 +114,7 @@ export q_translate, standard_driver, local_field_term, two_local_term,
     single_clause, q_translate_state, collective_operator, hamming_weight_operator
 
 export check_positivity, check_unitary, check_density_matrix, partial_trace,
-    matrix_decompose, low_level_matrix, fidelity, inst_population, gibbs_state, purity, check_pure_state
+    matrix_decompose, low_level_matrix, fidelity, inst_population, gibbs_state, purity, check_pure_state, find_degenerate
 
 export construct_interpolations, gradient, log_uniform
 

src/hamiltonian/interface.jl

@@ -15,12 +15,6 @@ function ConstantHamiltonian(mat::SparseMatrixCSC; unit=:h, static=true)
     ConstantSparseHamiltonian(mat, size(mat))
 end
 
-function Hamiltonian(f, mats; unit=:h, dimensionless_time=true, static=true)
-    @warn "Elements in `mats` have different types. Will attempt to promote them to a common type."
-    mats = promote(mats)
-    Hamiltonian(f, mats, unit=unit, dimensionless_time=dimensionless_time, static=static)
-end
-
 function Hamiltonian(f, mats::AbstractVector{T}; unit=:h, dimensionless_time=true, static=true) where {T<:Matrix}
     hsize = size(mats[1])
     # use static array for size smaller than 100

src/math_util.jl

@@ -229,6 +229,8 @@ $(SIGNATURES)
 
 Calculate the indices of the degenerate values in list `w`. `w` must be sorted in ascending order.
 
+If the `digits` keyword argument is provided, it rounds the energies to the specified number of digits after the decimal place (or before if negative), in base 10. The default value is 8.
+
 # Examples
 ```julia-repl
 julia> w = [1,1,1,2,3,4,4,5]
@@ -238,19 +240,24 @@ julia> find_degenerate(w)
  [6, 7]
 ```
 """
-function find_degenerate(w; atol::Real=1e-6, rtol::Real=0)
-    dw = w[2:end] - w[1:end-1]
-    inds = findall((x)->isapprox(x, 0.0, atol=atol, rtol=rtol), dw)
-    res = Array{Array{Int, 1}, 1}()
-    temp = Int[]
-    for i in eachindex(inds)
-        push!(temp, inds[i])
-        if (i==length(inds)) || !(inds[i]+1 == inds[i+1])
-            push!(temp, inds[i]+1)
-            push!(res, temp)
-            temp = Int[]
+function find_degenerate(w; digits::Integer=8)
+    w = round.(w, digits = digits)
+    res = []
+    tmp = []
+    flag = true
+    for i in 2:length(w)
+        if w[i] == w[i-1] && flag
+            push!(tmp, i-1, i)
+            flag = false
+        elseif w[i] == w[i-1]
+            push!(tmp, i)
+        elseif !flag
+            push!(res, tmp)
+            tmp = []
+            flag = true
         end
     end
+    isempty(tmp) ? nothing : push!(res, tmp)
     res
 end
 

src/projection/projection.jl

@@ -43,8 +43,7 @@ function project_to_lowlevel(
     coupling,
     dH;
     lvl=2,
-    atol::Real=1e-6,
-    rtol::Real=0,
+    digits::Integer=6,
     direction=:forward,
     refs=zeros(0, 0),
 ) where {T <: Real,S <: Real}
@@ -63,7 +62,7 @@ function project_to_lowlevel(
         # this is needed for StaticArrays
         w = w[1:lvl]
         v = v[:, 1:lvl]
-        d_inds = find_degenerate(w, atol=atol, rtol=rtol)
+        d_inds = find_degenerate(w, digits=digits)
         if !isempty(d_inds)
             @warn "Degenerate energy levels detected at" _s_axis[1]
             @warn "With" d_inds
@@ -81,7 +80,7 @@ function project_to_lowlevel(
         # this is needed for StaticArrays
         w = w[1:lvl]
         v = v[:, 1:lvl]
-        d_inds = find_degenerate(w, atol=atol, rtol=rtol)
+        d_inds = find_degenerate(w, digits=digits)
         if !isempty(d_inds)
             @warn "Possible degenerate detected at" s
             @warn "With levels" d_inds
@@ -98,14 +97,13 @@ function project_to_lowlevel(
     coupling,
     dH;
     lvl=2,
-    atol::Real=1e-6,
-    rtol::Real=0,
+    digits::Integer=6,
     direction=:forward,
     refs=zeros(0, 0),
 ) where {T <: Complex,S <: Real}
     @warn "The projection method only works with real Hamitonians. Convert the complex Hamiltonian to real one."
     H_real = convert(Real, H)
-    project_to_lowlevel(H_real, s_axis, coupling, dH, lvl=lvl, atol=atol, rtol=rtol, direction=direction, refs=refs)
+    project_to_lowlevel(H_real, s_axis, coupling, dH, lvl=lvl, digits=digits, direction=direction, refs=refs)
 end
 
 function update_refs!(refs, v, lvl, d_inds)

test/utilities/math_util.jl

@@ -15,7 +15,7 @@ m =
 # == units conversion test ===
 @test temperature_2_freq(1e3) ≈ 20.8366176361328 atol = 1e-4 rtol = 1e-4
 @test freq_2_temperature(20) ≈ 959.8489324422699 atol = 1e-4 rtol = 1e-4
-@test temperature_2_freq(1e3) ≈ 1 /  temperature_2_β(1e3) / 2 / π
+@test temperature_2_freq(1e3) ≈ 1 / temperature_2_β(1e3) / 2 / π
 @test β_2_temperature(0.47) ≈ 16.251564065921915
 # === unitary test ===
 u_res = exp(-1.0im * 5 * 0.5 * σx)
@@ -57,37 +57,37 @@ w, v = eigen_decomp(H_check, [0.5])
 @test log_uniform(1, 10, 3) == [1, 10^0.5, 10]
 
 v = sqrt.([0.4, 0.6])
-ρ1 = v*v'
+ρ1 = v * v'
 ρ2 = [0.5 0; 0 0.5]
-ρ3 = ones(3, 3)/3
+ρ3 = ones(3, 3) / 3
 @test ρ1 == partial_trace(ρ1 ⊗ ρ2 ⊗ ρ2, [1])
 @test ρ2 == partial_trace(ρ1 ⊗ ρ2 ⊗ ρ2, [2])
-@test ρ3 ≈ partial_trace(ρ1⊗ρ2⊗ρ3, [2,2,3], [3])
-@test_throws ArgumentError partial_trace(ρ1⊗ρ2⊗ρ3, [3,2,3], [3])
-@test_throws ArgumentError partial_trace(rand(4,5), [2,2], [1])
+@test ρ3 ≈ partial_trace(ρ1 ⊗ ρ2 ⊗ ρ3, [2, 2, 3], [3])
+@test_throws ArgumentError partial_trace(ρ1 ⊗ ρ2 ⊗ ρ3, [3, 2, 3], [3])
+@test_throws ArgumentError partial_trace(rand(4, 5), [2, 2], [1])
 @test purity(ρ1) ≈ 1
 @test purity(ρ2) == 0.5
 @test check_pure_state(ρ1)
 @test !check_pure_state(ρ2)
 @test !check_pure_state([0.4 0.5; 0.5 0.6])
 
-ρ = PauliVec[1][1]*PauliVec[1][1]'
-σ = PauliVec[3][1]*PauliVec[3][1]'
+ρ = PauliVec[1][1] * PauliVec[1][1]'
+σ = PauliVec[3][1] * PauliVec[3][1]'
 @test fidelity(ρ, σ) ≈ 0.5
 @test fidelity(ρ, ρ) ≈ 1
 @test check_density_matrix(ρ)
 @test !check_density_matrix(σx)
 
-w = [1,1,1,2,3,4,4,5]
-@test OpenQuantumBase.find_degenerate(w) == [[1,2,3],[6,7]]
-@test isempty(OpenQuantumBase.find_degenerate([1,2,3]))
+w = [-0.000000003, 0, 1, 1.000000001, 1.000000004, 1.000000006, 2, 3, 4, 4, 5, 5]
+@test find_degenerate(w, digits=8) == [[1, 2], [3, 4, 5], [9, 10], [11, 12]]
+@test isempty(find_degenerate([1, 2, 3]))
 
 gibbs = gibbs_state(σz, 12)
-@test gibbs[2,2] ≈ 1/(1+exp(-temperature_2_β(12)*2))
-@test gibbs[1,1] ≈ 1 - 1/(1+exp(-temperature_2_β(12)*2))
+@test gibbs[2, 2] ≈ 1 / (1 + exp(-temperature_2_β(12) * 2))
+@test gibbs[1, 1] ≈ 1 - 1 / (1 + exp(-temperature_2_β(12) * 2))
 
-@test low_level_matrix(σz⊗σz+0.1σz⊗σi, 2) == [0.0+0.0im   0.0+0.0im   0.0+0.0im  0.0+0.0im
-0.0+0.0im  -0.9+0.0im   0.0+0.0im  0.0+0.0im
-0.0+0.0im   0.0+0.0im  -1.1+0.0im  0.0+0.0im
-0.0+0.0im   0.0+0.0im   0.0+0.0im  0.0+0.0im]
-@test_logs (:warn, "Subspace dimension bigger than total dimension.") low_level_matrix(σz⊗σz+0.1σz⊗σi, 5)
+@test low_level_matrix(σz ⊗ σz + 0.1σz ⊗ σi, 2) == [0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
+    0.0+0.0im -0.9+0.0im 0.0+0.0im 0.0+0.0im
+    0.0+0.0im 0.0+0.0im -1.1+0.0im 0.0+0.0im
+    0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
+@test_logs (:warn, "Subspace dimension bigger than total dimension.") low_level_matrix(σz ⊗ σz + 0.1σz ⊗ σi, 5)