Merge pull request #55 from albertomercurio/main

Merge from main branch

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/CI.yml

@@ -4,9 +4,13 @@ on:
   push:
     branches:
       - 'main'
+    paths-ignore:
+      - 'docs/**'
   pull_request:
     branches:
       - 'main'
+    paths-ignore:
+      - 'docs/**'
     types:
       - opened
       - reopened
@@ -38,7 +42,7 @@ jobs:
           - Core
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
@@ -52,6 +56,6 @@ jobs:
           directories: src
       - uses: codecov/codecov-action@v4
         with:
-          verbose: true
-        env:
-          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: true

.github/workflows/documentation.yml

@@ -23,7 +23,7 @@ jobs:
     if: ${{ !github.event.pull_request.draft }}
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - uses: julia-actions/julia-buildpkg@v1

Project.toml

@@ -19,21 +19,19 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-DiffEqCallbacks = "2.24, <3.2.0"
+DiffEqCallbacks = "2.24, <3.2.0, >=3.5.0"
 FFTW = "1.5"
 Graphs = "1.7.4"
 IncompleteLU = "0.2"
 LinearMaps = "3"
 LinearSolve = "2"
 OrdinaryDiffEq = "6.30"
 Reexport = "1"
-SafeTestsets = "0.1"
 SpecialFunctions = "2.2"
 julia = "1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 
 [targets]
-test = ["Test", "SafeTestsets"]
+test = ["Test"]

src/QuantumToolbox.jl

@@ -36,6 +36,7 @@ include("spin_lattice.jl")
 include("arnoldi.jl")
 include("eigsolve.jl")
 include("negativity.jl")
+include("progress_bar.jl")
 
 export QuantumObject, Qobj, BraQuantumObject, KetQuantumObject, OperatorQuantumObject, OperatorBraQuantumObject, OperatorKetQuantumObject, SuperOperatorQuantumObject, TimeEvolutionSol
 export isket, isbra, isoper, isoperbra, isoperket, issuper, ket2dm

src/progress_bar.jl

@@ -0,0 +1,52 @@
+export ProgressBar
+
+struct ProgressBar{CT,T1<:Integer, T2<:Real}
+    counter::CT
+    max_counts::T1
+    enable::Bool
+    bar_width::T1
+    start_time::T2
+    lock::ReentrantLock
+end
+
+function ProgressBar(max_counts; enable=true, bar_width=30)
+    return ProgressBar(Ref{Int64}(0), max_counts, enable, bar_width, time(), ReentrantLock())
+end
+
+function next!(p::ProgressBar)
+
+    lock(p.lock)
+
+    p.counter[] += 1
+
+    !p.enable && return
+
+    counter = p.counter[]
+    max_counts = p.max_counts
+    bar_width = p.bar_width
+    start_time = p.start_time
+
+    percentage = counter / max_counts
+    percentage_100 = lpad(round(100 * percentage, digits=1), 5, " ")
+    progress = floor(Int, bar_width * percentage)
+
+    # Calculate the elapsed time in seconds
+    elapsed_time = floor(Int, time() - start_time)
+    # Convert the elapsed time into a string in hours, minutes and seconds
+    elapsed_time_str = string(elapsed_time ÷ 3600, "h ", lpad((elapsed_time % 3600) ÷ 60, 2, "0"), "m ", lpad(elapsed_time % 60, 2, "0"), "s")
+
+    # Calculate the estimated time of arrival
+    eta = floor(Int, elapsed_time ÷ counter * (max_counts - counter))
+    # convert eta into a string in hours, minutes and seconds
+    eta_str = string(eta ÷ 3600, "h ", lpad((eta % 3600) ÷ 60, 2, "0"), "m ", lpad(eta % 60, 2, "0"), "s")
+
+    # Construct the progress bar string
+    bar = "[" * repeat("=", progress) * repeat(" ", bar_width - progress) * "]"
+
+    print("\rProgress: $bar $percentage_100% --- Elapsed Time: $elapsed_time_str (ETA: $eta_str)")
+    flush(stdout)
+
+    unlock(p.lock)
+
+    p.counter[] >= p.max_counts ? print("\n") : nothing
+end

src/quantum_object.jl

@@ -78,31 +78,62 @@ mutable struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType} <: Ab
     dims::Vector{Int}
 end
 
-function QuantumObject(A::AbstractVector{T}; type::Type{ObjType}=KetQuantumObject, dims=nothing) where
-{T,ObjType<:Union{BraQuantumObject,KetQuantumObject}}
-
-    dims === nothing ? dims = [length(A)] : nothing
-    prod(dims) != length(A) && throw(DimensionMismatch("The dims parameter does not fit the dimension of the Vector."))
-    # if !(norm(A) ≈ 1)
-    #     @warn "The norm of the input data is not one."
-    # end
-    ObjType <: KetQuantumObject ? QuantumObject(A, type, dims) : QuantumObject(transpose(A), type, dims)
-end
-
-function QuantumObject(A::AbstractMatrix{T}; type::Type{ObjType}=OperatorQuantumObject, dims=nothing) where
-{T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+function QuantumObject(A::AbstractArray{T, N}; type::Type{ObjType}=Nothing, dims=nothing) where {T,N,ObjType<:Union{Nothing,QuantumObjectType}}
 
-    n = checksquare(A)
-    if type <: SuperOperatorQuantumObject
-        sqrt(n) - isqrt(n) != 0 ? throw(DomainError(n, "The dimension of the matrix is not compatible with the SuperOperator type")) : nothing
-        dims === nothing ? dims = [isqrt(n)] : nothing
-        prod(dims) != isqrt(n) && throw(DimensionMismatch("The dims parameter does not fit the dimension of the Matrix."))
+    # only accept 1D- and 2D-array
+    if N == 1
+        Size = (length(A), 1)
     else
-        dims === nothing ? dims = [n] : nothing
-        prod(dims) != n && throw(DimensionMismatch("The dims parameter does not fit the dimension of the Matrix."))
+        Size = size(A)
+        N > 2 ? throw(DomainError(Size, "The dimension of the array is not compatible with Quantum Object")) : nothing
+    end
+
+    # decide QuantumObjectType from the size of A
+    if type == Nothing
+        if Size[1] == Size[2]
+            type = OperatorQuantumObject
+        elseif Size[2] == 1
+            type = KetQuantumObject
+        elseif Size[1] == 1
+            type = BraQuantumObject
+        else
+            throw(DomainError(Size, "The dimension of the array is not compatible with Quantum Object"))
+        end
+    end
+
+    # decide dims from the size of A and the given type
+    if dims === nothing 
+        if (type <: KetQuantumObject) || (type <: OperatorQuantumObject)
+            dims = [Size[1]]
+        elseif type <: SuperOperatorQuantumObject
+            dims = [isqrt(Size[1])]
+        elseif type <: BraQuantumObject
+            dims = [Size[2]]
+        end
     end
 
-    QuantumObject(A, type, dims)
+    _check_QuantumObject(type, prod(dims), Size[1], Size[2])
+    return QuantumObject(A, type, dims)
+end
+
+function _check_QuantumObject(type::Type{KetQuantumObject}, prod_dims::Int, m::Int, n::Int)
+    (n != 1) ? throw(DomainError((m, n), "The dimension of the array is not compatible with Ket type")) : nothing
+    prod_dims != m ? throw(DimensionMismatch("The dims parameter does not fit the dimension of the Array.")) : nothing
+end
+
+function _check_QuantumObject(type::Type{BraQuantumObject}, prod_dims::Int, m::Int, n::Int)
+    (m != 1) ? throw(DomainError((m, n), "The dimension of the array is not compatible with Bra type")) : nothing
+    prod_dims != n ? throw(DimensionMismatch("The dims parameter does not fit the dimension of the Array.")) : nothing
+end
+
+function _check_QuantumObject(type::Type{OperatorQuantumObject}, prod_dims::Int, m::Int, n::Int)
+    (m != n) ? throw(DomainError((m, n), "The dimension of the array is not compatible with Operator type")) : nothing
+    prod_dims != m ? throw(DimensionMismatch("The dims parameter does not fit the dimension of the Array.")) : nothing
+end
+
+function _check_QuantumObject(type::Type{SuperOperatorQuantumObject}, prod_dims::Int, m::Int, n::Int)
+    (m != n) ? throw(DomainError((m, n), "The dimension of the array is not compatible with SuperOperator type")) : nothing
+    prod_dims != sqrt(m) ? throw(DimensionMismatch("The dims parameter does not fit the dimension of the Array.")) : nothing
 end
 
 function QuantumObject(A::QuantumObject{<:AbstractArray}; type::Type{ObjType}=A.type, dims=A.dims) where

src/time_evolution/mcsolve.jl

@@ -63,7 +63,7 @@ function _mcsolve_prob_func(prob, i, repeat)
 
     prm = merge(internal_params, (expvals = similar(internal_params.expvals),
                 cache_mc = similar(internal_params.cache_mc), weights_mc = similar(internal_params.weights_mc), 
-                cumsum_weights_mc = similar(internal_params.weights_mc), random_n = Ref(rand()), progr_mc = ODEProgress(0), jump_times_which_idx = Ref(1),
+                cumsum_weights_mc = similar(internal_params.weights_mc), random_n = Ref(rand()), progr_mc = ProgressBar(size(internal_params.expvals, 2), enable=false), jump_times_which_idx = Ref(1),
                 jump_times = similar(internal_params.jump_times), jump_which = similar(internal_params.jump_which)))
 
     remake(prob, p=prm)
@@ -149,7 +149,7 @@ function mcsolveProblem(H::QuantumObject{MT1,OperatorQuantumObject},
 
     params2 = (expvals = expvals, e_ops_mc = e_ops2,
                 is_empty_e_ops_mc = isempty(e_ops),
-                progr_mc = ODEProgress(0), seeds = seeds,
+                progr_mc = ProgressBar(length(t_l), enable=false), seeds = seeds,
                 random_n = Ref(rand()), c_ops = get_data.(c_ops), cache_mc = cache_mc, 
                 weights_mc = weights_mc, cumsum_weights_mc = cumsum_weights_mc,
                 jump_times = jump_times, jump_which = jump_which,

src/time_evolution/mesolve.jl

@@ -34,6 +34,7 @@ end
         e_ops::AbstractVector=[],
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
+        progress_bar::Bool=true,
         kwargs...)
 
 Generates the ODEProblem for the master equation time evolution of an open quantum system.
@@ -47,6 +48,7 @@ Generates the ODEProblem for the master equation time evolution of an open quant
 - `e_ops::AbstractVector=[]`: The list of the operators for which the expectation values are calculated.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing`: The time-dependent Hamiltonian or Liouvillian.
 - `params::NamedTuple=NamedTuple()`: The parameters of the time evolution.
+- `progress_bar::Bool=true`: Whether to show the progress bar.
 - `kwargs...`: The keyword arguments for the ODEProblem.
 
 # Returns
@@ -60,6 +62,7 @@ function mesolveProblem(H::QuantumObject{MT1,HOpType},
     e_ops::Vector{QuantumObject{Te, OperatorQuantumObject}}=QuantumObject{MT1, OperatorQuantumObject}[],
     H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
     params::NamedTuple=NamedTuple(),
+    progress_bar::Bool=true,
     kwargs...) where {MT1<:AbstractMatrix,T2,Tc<:AbstractMatrix,Te<:AbstractMatrix,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
@@ -72,7 +75,7 @@ function mesolveProblem(H::QuantumObject{MT1,HOpType},
     ρ0 = mat2vec(ket2dm(ψ0).data)
     L = liouvillian(H, c_ops).data
 
-    progr = ODEProgress(0)
+    progr = ProgressBar(length(t_l), enable=progress_bar)
     expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
     e_ops2 = @. mat2vec(adjoint(get_data(e_ops)))
 
@@ -105,6 +108,7 @@ end
         e_ops::AbstractVector=[],
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
+        progress_bar::Bool=true,
         kwargs...)
 
 Time evolution of an open quantum system using master equation.
@@ -118,6 +122,7 @@ Time evolution of an open quantum system using master equation.
 - `e_ops::AbstractVector`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Named Tuple of parameters to pass to the solver.
+- `progress_bar::Bool`: Whether to show the progress bar.
 - `kwargs...`: Additional keyword arguments to pass to the solver.
 
 # Returns
@@ -131,13 +136,14 @@ function mesolve(H::QuantumObject{MT1,HOpType},
     e_ops::Vector{QuantumObject{Te, OperatorQuantumObject}}=QuantumObject{MT1, OperatorQuantumObject}[],
     H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
     params::NamedTuple=NamedTuple(),
+    progress_bar::Bool=true,
     kwargs...) where {MT1<:AbstractMatrix,T2,Tc<:AbstractMatrix,Te<:AbstractMatrix,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
     COpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
 
     prob = mesolveProblem(H, ψ0, t_l, c_ops; alg=alg, e_ops=e_ops,
-            H_t=H_t, params=params, kwargs...)
+            H_t=H_t, params=params, progress_bar=progress_bar, kwargs...)
 
     return mesolve(prob, alg; kwargs...)
 end

src/time_evolution/sesolve.jl

@@ -32,6 +32,7 @@ end
         e_ops::AbstractVector=[],
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
+        progress_bar::Bool=true,
         kwargs...)
 
 Generates the ODEProblem for the Schrödinger time evolution of a quantum system.
@@ -44,6 +45,7 @@ Generates the ODEProblem for the Schrödinger time evolution of a quantum system
 - `e_ops::AbstractVector`: The list of operators to be evaluated during the evolution.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
 - `params::NamedTuple`: The parameters of the system.
+- `progress_bar::Bool`: Whether to show the progress bar.
 - `kwargs...`: The keyword arguments passed to the `ODEProblem` constructor.
 
 # Returns
@@ -56,6 +58,7 @@ function sesolveProblem(H::QuantumObject{MT1,OperatorQuantumObject},
     e_ops::Vector{QuantumObject{MT2, OperatorQuantumObject}}=QuantumObject{MT1, OperatorQuantumObject}[],
     H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
     params::NamedTuple=NamedTuple(),
+    progress_bar::Bool=true,
     kwargs...) where {MT1<:AbstractMatrix,T2,MT2<:AbstractMatrix}
 
     H.dims != ψ0.dims && throw(ErrorException("The two operators don't have the same Hilbert dimension."))
@@ -65,7 +68,7 @@ function sesolveProblem(H::QuantumObject{MT1,OperatorQuantumObject},
     ϕ0 = get_data(ψ0)
     U = -1im * get_data(H)
 
-    progr = ODEProgress(0)
+    progr = ProgressBar(length(t_l), enable=progress_bar)
     expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
     e_ops2 = get_data.(e_ops)
 
@@ -98,6 +101,7 @@ end
         e_ops::AbstractVector=[],
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
+        progress_bar::Bool=true,
         kwargs...)
 
 Time evolution of a closed quantum system using the Schrödinger equation.
@@ -110,6 +114,7 @@ Time evolution of a closed quantum system using the Schrödinger equation.
 - `e_ops::AbstractVector`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Dictionary of parameters to pass to the solver.
+- `progress_bar::Bool`: Whether to show the progress bar.
 - `kwargs...`: Additional keyword arguments to pass to the solver.
 
 - Returns
@@ -122,10 +127,11 @@ function sesolve(H::QuantumObject{MT1,OperatorQuantumObject},
     e_ops::Vector{QuantumObject{MT2, OperatorQuantumObject}}=QuantumObject{MT1, OperatorQuantumObject}[],
     H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
     params::NamedTuple=NamedTuple(),
+    progress_bar::Bool=true,
     kwargs...) where {MT1<:AbstractMatrix,T2,MT2<:AbstractMatrix}
 
     prob = sesolveProblem(H, ψ0, t_l; alg=alg, e_ops=e_ops,
-            H_t=H_t, params=params, kwargs...)
+            H_t=H_t, params=params, progress_bar=progress_bar, kwargs...)
 
     return sesolve(prob, alg; kwargs...)
 end

src/time_evolution/time_evolution.jl

@@ -26,22 +26,6 @@ end
 
 LiouvillianDirectSolver(;tol=1e-16) = LiouvillianDirectSolver(tol)
 
-# It is needed to keep track of the index for saving the expectation values
-mutable struct ODEProgress{T<:Integer}
-    counter::T
-    # max_counts::T
-    # start_time::T2
-end
-
-function next!(p::ODEProgress)
-    p.counter += 1
-    # prog = round(100 * p.counter / p.max_counts, digits=2)
-    # eta = round((time() - p.start_time) / p.counter * (p.max_counts - p.counter), digits=2)
-    # # convert eta into a string in hours, minutes and seconds
-    # eta_str = string(round(Int, eta / 3600), "h ", round(Int, (eta % 3600) / 60), "m ", round(Int, eta % 60), "s")
-    # print("\rProgress: $prog%  ---  ETA: $eta_str")
-end
-
 abstract type LindbladJumpCallbackType end
 
 struct ContinuousLindbladJumpCallback <: LindbladJumpCallbackType

src/time_evolution/time_evolution_dynamical.jl

@@ -467,7 +467,7 @@ function _dsf_mcsolve_prob_func(prob, i, repeat)
     prm = merge(internal_params, (U = copy(internal_params.U), e_ops_mc = copy(internal_params.e_ops_mc),
                 c_ops = copy(internal_params.c_ops), expvals = similar(internal_params.expvals), 
                 cache_mc = similar(internal_params.cache_mc), weights_mc = similar(internal_params.weights_mc), 
-                cumsum_weights_mc = similar(internal_params.weights_mc), random_n = Ref(rand()), progr_mc = ODEProgress(0), jump_times_which_idx = Ref(1),
+                cumsum_weights_mc = similar(internal_params.weights_mc), random_n = Ref(rand()), progr_mc = ProgressBar(size(internal_params.expvals, 2), enable=false), jump_times_which_idx = Ref(1),
                 jump_times = similar(internal_params.jump_times), jump_which = similar(internal_params.jump_which),
                 αt_list = copy(internal_params.αt_list), dsf_cache1 = similar(internal_params.dsf_cache1),
                 dsf_cache2 = similar(internal_params.dsf_cache2), expv_cache = copy(internal_params.expv_cache),