From 12159f322b8981e5def1b7d4eae7b9837467ba5d Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Thu, 25 Oct 2018 22:02:52 +0200 Subject: [PATCH] Implement parametric types (#238) * Start parametric typing for Ket/Bra * Fix tensor vararg for StateVector * Fix subspacebasis field parameters * Fix semiclassical state for kets * Update testing scripts * Fix typo in appveyor * Fix semiclassical_stochastic ket typing * Fix travis 1.0 testing * Parametric typing for CompositeBasis * Rename Operator to AbstractOperator * Proper parametric type for CompositeBasis * Parametric typing for operators * Parametrize basis dimensions * Revert "Parametrize basis dimensions" This reverts commit b451987c22717e859ef7a1d0dbe26b66437c95ea. * Update basis checks for states * Update basis checks for dense operators * Update basis checks for sparse operators * Update LazyProduct implementation * Update LazySum implementation * Update LazyTensor implementation * Add non-type parameters where needed to ensure correct basis dispatch * Update basis checks for schroedinger * Update basis checks for metrics * Update basis checks for phasespace * Update basis checks for master * Update basis checks for mcwf * Update basis checks for semiclassical * Update superoperators and steadystate * Update basis checks for stochastic solvers * Update timecorrelations basis checks --- src/QuantumOptics.jl | 2 +- src/bases.jl | 15 +-- src/manybody.jl | 16 +-- src/master.jl | 96 +++++++++--------- src/mcwf.jl | 158 ++++++++++++++--------------- src/metrics.jl | 50 +++++---- src/operators.jl | 138 ++++++++++++------------- src/operators_dense.jl | 97 ++++++++++-------- src/operators_lazyproduct.jl | 45 ++++---- src/operators_lazysum.jl | 51 ++++++---- src/operators_lazytensor.jl | 64 ++++++------ src/operators_sparse.jl | 84 ++++++++------- src/particle.jl | 65 ++++++++---- src/phasespace.jl | 38 +++---- src/printing.jl | 4 +- src/schroedinger.jl | 19 ++-- src/semiclassical.jl | 45 ++++---- src/spectralanalysis.jl | 8 +- src/spin.jl | 8 +- src/state_definitions.jl | 4 +- src/states.jl | 67 ++++++++---- src/steadystate.jl | 12 +-- src/stochastic_definitions.jl | 6 +- src/stochastic_master.jl | 45 ++++---- src/stochastic_schroedinger.jl | 36 +++---- src/stochastic_semiclassical.jl | 59 +++++------ src/subspace.jl | 15 ++- src/superoperators.jl | 69 +++++++------ src/timecorrelations.jl | 30 +++--- test/test_embed.jl | 4 +- test/test_fock.jl | 2 +- test/test_manybody.jl | 2 +- test/test_metrics.jl | 8 +- test/test_operators.jl | 8 +- test/test_operators_dense.jl | 8 +- test/test_operators_lazyproduct.jl | 14 +-- test/test_operators_lazysum.jl | 24 ++--- test/test_operators_lazytensor.jl | 32 +++--- test/test_operators_sparse.jl | 12 +-- test/test_particle.jl | 8 +- test/test_phasespace.jl | 2 +- test/test_semiclassical.jl | 14 +-- test/test_spectralanalysis.jl | 6 +- test/test_spin.jl | 2 +- test/test_states.jl | 16 +-- test/test_superoperators.jl | 3 - test/test_transformations.jl | 2 +- 47 files changed, 801 insertions(+), 712 deletions(-) diff --git a/src/QuantumOptics.jl b/src/QuantumOptics.jl index a60160b4..77757dbb 100644 --- a/src/QuantumOptics.jl +++ b/src/QuantumOptics.jl @@ -6,7 +6,7 @@ export bases, Basis, GenericBasis, CompositeBasis, basis, tensor, ⊗, permutesystems, @samebases, states, StateVector, Bra, Ket, basisstate, norm, dagger, normalize, normalize!, - operators, Operator, expect, variance, identityoperator, ptrace, embed, dense, tr, + operators, AbstractOperator, expect, variance, identityoperator, ptrace, embed, dense, tr, sparse, operators_dense, DenseOperator, projector, dm, operators_sparse, SparseOperator, diagonaloperator, diff --git a/src/bases.jl b/src/bases.jl index d3aaa164..8a65b275 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -71,14 +71,17 @@ Stores the subbases in a vector and creates the shape vector directly from the shape vectors of these subbases. Instead of creating a CompositeBasis directly `tensor(b1, b2...)` or `b1 ⊗ b2 ⊗ …` can be used. """ -mutable struct CompositeBasis <: Basis +mutable struct CompositeBasis{B<:Tuple{Vararg{Basis}}} <: Basis shape::Vector{Int} - bases::Vector{Basis} + bases::B end -CompositeBasis(bases::Vector{Basis}) = CompositeBasis(Int[prod(b.shape) for b in bases], bases) -CompositeBasis(bases::Basis...) = CompositeBasis(Basis[bases...]) +CompositeBasis(bases::B) where B<:Tuple{Vararg{Basis}} = CompositeBasis{B}(Int[prod(b.shape) for b in bases], bases) +CompositeBasis(shape::Vector{Int}, bases::Vector{B}) where B<:Basis = (tmp = (bases...,); CompositeBasis{typeof(tmp)}(shape, tmp)) +CompositeBasis(bases::Vector{B}) where B<:Basis = CompositeBasis((bases...,)) +CompositeBasis(bases::Basis...) = CompositeBasis((bases...,)) -==(b1::CompositeBasis, b2::CompositeBasis) = equal_shape(b1.shape, b2.shape) && equal_bases(b1.bases, b2.bases) +==(b1::T, b2::T) where T<:CompositeBasis = equal_shape(b1.shape, b2.shape) +==(b1::CompositeBasis, b2::CompositeBasis) = false """ tensor(x, y, z...) @@ -98,7 +101,7 @@ Any given CompositeBasis is expanded so that the resulting CompositeBasis never contains another CompositeBasis. """ tensor(b1::Basis, b2::Basis) = CompositeBasis(Int[prod(b1.shape); prod(b2.shape)], Basis[b1, b2]) -tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis(Int[b1.shape; b2.shape], Basis[b1.bases; b2.bases]) +tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis(Int[b1.shape; b2.shape], (b1.bases..., b2.bases...)) function tensor(b1::CompositeBasis, b2::Basis) N = length(b1.bases) shape = Vector{Int}(undef, N+1) diff --git a/src/manybody.jl b/src/manybody.jl index 779decdf..dcc2d617 100644 --- a/src/manybody.jl +++ b/src/manybody.jl @@ -20,16 +20,18 @@ The basis has to know the associated one-body basis `b` and which occupation sta should be included. The occupations_hash is used to speed up checking if two many-body bases are equal. """ -mutable struct ManyBodyBasis <: Basis +mutable struct ManyBodyBasis{B<:Basis,H} <: Basis shape::Vector{Int} - onebodybasis::Basis + onebodybasis::B occupations::Vector{Vector{Int}} occupations_hash::UInt - function ManyBodyBasis(onebodybasis::Basis, occupations::Vector{Vector{Int}}) + function ManyBodyBasis{B,H}(onebodybasis::Basis, occupations::Vector{Vector{Int}}) where {B<:Basis,H} + @assert isa(H, UInt) new([length(occupations)], onebodybasis, occupations, hash(hash.(occupations))) end end +ManyBodyBasis(onebodybasis::B, occupations::Vector{Vector{Int}}) where B<:Basis = ManyBodyBasis{B,hash(hash.(occupations))}(onebodybasis,occupations) """ fermionstates(Nmodes, Nparticles) @@ -228,7 +230,7 @@ where ``X`` is the N-particle operator, ``x`` is the one-body operator and ``|u⟩`` are the one-body states associated to the different modes of the N-particle basis. """ -function manybodyoperator(basis::ManyBodyBasis, op::T)::T where T<:Operator +function manybodyoperator(basis::ManyBodyBasis, op::T) where T<:AbstractOperator @assert op.basis_l == op.basis_r if op.basis_l == basis.onebodybasis result = manybodyoperator_1(basis, op) @@ -321,7 +323,7 @@ end Expectation value of the one-body operator `op` in respect to the many-body `state`. """ -function onebodyexpect(op::Operator, state::Ket) +function onebodyexpect(op::AbstractOperator, state::Ket) @assert isa(state.basis, ManyBodyBasis) @assert op.basis_l == op.basis_r if state.basis.onebodybasis == op.basis_l @@ -335,7 +337,7 @@ function onebodyexpect(op::Operator, state::Ket) result end -function onebodyexpect(op::Operator, state::Operator) +function onebodyexpect(op::AbstractOperator, state::AbstractOperator) @assert op.basis_l == op.basis_r @assert state.basis_l == state.basis_r @assert isa(state.basis_l, ManyBodyBasis) @@ -349,7 +351,7 @@ function onebodyexpect(op::Operator, state::Operator) end result end -onebodyexpect(op::Operator, states::Vector) = [onebodyexpect(op, state) for state=states] +onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state=states] function onebodyexpect_1(op::DenseOperator, state::Ket) N = length(state.basis) diff --git a/src/master.jl b/src/master.jl index 9336c39b..fb106050 100644 --- a/src/master.jl +++ b/src/master.jl @@ -17,14 +17,14 @@ Integrate the master equation with dmaster_h as derivative function. Further information can be found at [`master`](@ref). """ -function master_h(tspan, rho0::DenseOperator, H::Operator, J::Vector; +function master_h(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector; rates::DecayRates=nothing, Jdagger::Vector=dagger.(J), fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} check_master(rho0, H, J, Jdagger, rates) tmp = copy(rho0) - dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) + dmaster_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) integrate_master(tspan, dmaster_, rho0, fout; kwargs...) end @@ -39,15 +39,15 @@ H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k ``` Further information can be found at [`master`](@ref). """ -function master_nh(tspan, rho0::DenseOperator, Hnh::Operator, J::Vector; +function master_nh(tspan, rho0::T, Hnh::AbstractOperator{B,B}, J::Vector; rates::DecayRates=nothing, - Hnhdagger::Operator=dagger(Hnh), + Hnhdagger::AbstractOperator=dagger(Hnh), Jdagger::Vector=dagger.(J), fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} check_master(rho0, Hnh, J, Jdagger, rates) tmp = copy(rho0) - dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) + dmaster_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) integrate_master(tspan, dmaster_, rho0, fout; kwargs...) end @@ -83,15 +83,15 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster. be changed. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master(tspan, rho0::DenseOperator, H::Operator, J::Vector; +function master(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector; rates::DecayRates=nothing, Jdagger::Vector=dagger.(J), fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} isreducible = check_master(rho0, H, J, Jdagger, rates) if !isreducible tmp = copy(rho0) - dmaster_h_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) + dmaster_h_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) return integrate_master(tspan, dmaster_h_, rho0, fout; kwargs...) else Hnh = copy(H) @@ -110,7 +110,7 @@ function master(tspan, rho0::DenseOperator, H::Operator, J::Vector; end Hnhdagger = dagger(Hnh) tmp = copy(rho0) - dmaster_nh_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) + dmaster_nh_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) return integrate_master(tspan, dmaster_nh_, rho0, fout; kwargs...) end end @@ -128,12 +128,12 @@ The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, J or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look at [`master_dynamic`](@ref). """ -function master_nh_dynamic(tspan, rho0::DenseOperator, f::Function; +function master_nh_dynamic(tspan, rho0::T, f::Function; rates::DecayRates=nothing, fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tmp = copy(rho0) - dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp) + dmaster_(t, rho::T, drho::T) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp) integrate_master(tspan, dmaster_, rho0, fout; kwargs...) end @@ -162,36 +162,36 @@ operators: be changed. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master_dynamic(tspan, rho0::DenseOperator, f::Function; +function master_dynamic(tspan, rho0::T, f::Function; rates::DecayRates=nothing, fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tmp = copy(rho0) - dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp) + dmaster_(t, rho::T, drho::T) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp) integrate_master(tspan, dmaster_, rho0, fout; kwargs...) end # Automatically convert Ket states to density operators -master(tspan, psi0::Ket, H::Operator, J::Vector; kwargs...) = master(tspan, dm(psi0), H, J; kwargs...) -master_h(tspan, psi0::Ket, H::Operator, J::Vector; kwargs...) = master_h(tspan, dm(psi0), H, J; kwargs...) -master_nh(tspan, psi0::Ket, Hnh::Operator, J::Vector; kwargs...) = master_nh(tspan, dm(psi0), Hnh, J; kwargs...) -master_dynamic(tspan, psi0::Ket, f::Function; kwargs...) = master_dynamic(tspan, dm(psi0), f; kwargs...) -master_nh_dynamic(tspan, psi0::Ket, f::Function; kwargs...) = master_nh_dynamic(tspan, dm(psi0), f; kwargs...) +master(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master(tspan, dm(psi0), H, J; kwargs...) +master_h(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_h(tspan, dm(psi0), H, J; kwargs...) +master_nh(tspan, psi0::Ket{B}, Hnh::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_nh(tspan, dm(psi0), Hnh, J; kwargs...) +master_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_dynamic(tspan, dm(psi0), f; kwargs...) +master_nh_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_nh_dynamic(tspan, dm(psi0), f; kwargs...) # Recasting needed for the ODE solver is just providing the underlying data -function recast!(x::Array{ComplexF64, 2}, rho::DenseOperator) +function recast!(x::T, rho::DenseOperator{B,B,T}) where {B<:Basis,T<:Matrix{ComplexF64}} rho.data = x end -recast!(rho::DenseOperator, x::Array{ComplexF64, 2}) = nothing +recast!(rho::DenseOperator{B,B,T}, x::T) where {B<:Basis,T<:Matrix{ComplexF64}} = nothing -function integrate_master(tspan, df::Function, rho0::DenseOperator, - fout::Union{Nothing, Function}; kwargs...) +function integrate_master(tspan, df::Function, rho0::T, + fout::Union{Nothing, Function}; kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tspan_ = convert(Vector{Float64}, tspan) x0 = rho0.data - state = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data) - dstate = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data) + state = T(rho0.basis_l, rho0.basis_r, rho0.data) + dstate = T(rho0.basis_l, rho0.basis_r, rho0.data) integrate(tspan_, df, x0, state, dstate, fout; kwargs...) end @@ -205,9 +205,9 @@ end # the type of the given decay rate object which can either be nothing, a vector # or a matrix. -function dmaster_h(rho::DenseOperator, H::Operator, +function dmaster_h(rho::T, H::AbstractOperator{B,B}, rates::Nothing, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator) + drho::T, tmp::T) where {B<:Basis,T<:DenseOperator{B,B}} operators.gemm!(-1im, H, rho, 0, drho) operators.gemm!(1im, rho, H, 1, drho) for i=1:length(J) @@ -222,9 +222,9 @@ function dmaster_h(rho::DenseOperator, H::Operator, return drho end -function dmaster_h(rho::DenseOperator, H::Operator, +function dmaster_h(rho::T, H::AbstractOperator{B,B}, rates::Vector{Float64}, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator) + drho::T, tmp::T) where {B<:Basis,T<:DenseOperator{B,B}} operators.gemm!(-1im, H, rho, 0, drho) operators.gemm!(1im, rho, H, 1, drho) for i=1:length(J) @@ -239,9 +239,9 @@ function dmaster_h(rho::DenseOperator, H::Operator, return drho end -function dmaster_h(rho::DenseOperator, H::Operator, +function dmaster_h(rho::T, H::AbstractOperator{B,B}, rates::Matrix{Float64}, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator) + drho::T, tmp::T) where {B<:Basis,T<:DenseOperator{B,B}} operators.gemm!(-1im, H, rho, 0, drho) operators.gemm!(1im, rho, H, 1, drho) for j=1:length(J), i=1:length(J) @@ -256,9 +256,9 @@ function dmaster_h(rho::DenseOperator, H::Operator, return drho end -function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator, +function dmaster_nh(rho::T1, Hnh::T2, Hnh_dagger::T2, rates::Nothing, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator) + drho::T1, tmp::T1) where {B<:Basis,T1<:DenseOperator{B,B},T2<:AbstractOperator{B,B}} operators.gemm!(-1im, Hnh, rho, 0, drho) operators.gemm!(1im, rho, Hnh_dagger, 1, drho) for i=1:length(J) @@ -268,9 +268,9 @@ function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator, return drho end -function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator, +function dmaster_nh(rho::T1, Hnh::T2, Hnh_dagger::T2, rates::Vector{Float64}, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator) + drho::T1, tmp::T1) where {B<:Basis,T1<:DenseOperator{B,B},T2<:AbstractOperator{B,B}} operators.gemm!(-1im, Hnh, rho, 0, drho) operators.gemm!(1im, rho, Hnh_dagger, 1, drho) for i=1:length(J) @@ -280,9 +280,9 @@ function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator, return drho end -function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator, +function dmaster_nh(rho::T1, Hnh::T2, Hnh_dagger::T2, rates::Matrix{Float64}, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator) + drho::T1, tmp::T1) where {B<:Basis,T1<:DenseOperator{B,B},T2<:AbstractOperator{B,B}} operators.gemm!(-1im, Hnh, rho, 0, drho) operators.gemm!(1im, rho, Hnh_dagger, 1, drho) for j=1:length(J), i=1:length(J) @@ -292,9 +292,9 @@ function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator, return drho end -function dmaster_h_dynamic(t::Float64, rho::DenseOperator, f::Function, +function dmaster_h_dynamic(t::Float64, rho::T, f::Function, rates::DecayRates, - drho::DenseOperator, tmp::DenseOperator) + drho::T, tmp::T) where {B<:Basis,T<:DenseOperator{B,B}} result = f(t, rho) QO_CHECKS[] && @assert 3 <= length(result) <= 4 if length(result) == 3 @@ -307,9 +307,9 @@ function dmaster_h_dynamic(t::Float64, rho::DenseOperator, f::Function, dmaster_h(rho, H, rates_, J, Jdagger, drho, tmp) end -function dmaster_nh_dynamic(t::Float64, rho::DenseOperator, f::Function, +function dmaster_nh_dynamic(t::Float64, rho::T, f::Function, rates::DecayRates, - drho::DenseOperator, tmp::DenseOperator) + drho::T, tmp::T) where {B<:Basis,T<:DenseOperator{B,B}} result = f(t, rho) QO_CHECKS[] && @assert 4 <= length(result) <= 5 if length(result) == 4 @@ -323,21 +323,21 @@ function dmaster_nh_dynamic(t::Float64, rho::DenseOperator, f::Function, end -function check_master(rho0::DenseOperator, H::Operator, J::Vector, Jdagger::Vector, rates::DecayRates) +function check_master(rho0::DenseOperator{B,B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis + # TODO: clean up type checks by dispatch; make type of J known isreducible = true # test if all operators are sparse or dense - check_samebases(rho0, H) if !(isa(H, DenseOperator) || isa(H, SparseOperator)) isreducible = false end for j=J - @assert isa(j, Operator) + @assert isa(j, AbstractOperator{B,B}) if !(isa(j, DenseOperator) || isa(j, SparseOperator)) isreducible = false end check_samebases(rho0, j) end for j=Jdagger - @assert isa(j, Operator) + @assert isa(j, AbstractOperator{B,B}) if !(isa(j, DenseOperator) || isa(j, SparseOperator)) isreducible = false end diff --git a/src/mcwf.jl b/src/mcwf.jl index f4185291..9a1b9487 100644 --- a/src/mcwf.jl +++ b/src/mcwf.jl @@ -23,19 +23,19 @@ Calculate MCWF trajectory where the Hamiltonian is given in hermitian form. For more information see: [`mcwf`](@ref) """ -function mcwf_h(tspan, psi0::Ket, H::Operator, J::Vector; - seed=rand(UInt), rates::DecayRates=nothing, - fout=nothing, Jdagger::Vector=dagger.(J), - tmp::Ket=copy(psi0), - display_beforeevent=false, display_afterevent=false, - kwargs...) +function mcwf_h(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector; + seed=rand(UInt), rates::DecayRates=nothing, + fout=nothing, Jdagger::Vector=dagger.(J), + tmp::T=copy(psi0), + display_beforeevent=false, display_afterevent=false, + kwargs...) where {B<:Basis,T<:Ket{B}} check_mcwf(psi0, H, J, Jdagger, rates) - f(t, psi, dpsi) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates) - j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates) - return integrate_mcwf(f, j, tspan, psi0, seed, fout; - display_beforeevent=display_beforeevent, - display_afterevent=display_afterevent, - kwargs...) + f(t::Float64, psi::T, dpsi::T) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates) + j(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates) + integrate_mcwf(f, j, tspan, psi0, seed, fout; + display_beforeevent=display_beforeevent, + display_afterevent=display_afterevent, + kwargs...) end """ @@ -49,17 +49,17 @@ H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k For more information see: [`mcwf`](@ref) """ -function mcwf_nh(tspan, psi0::Ket, Hnh::Operator, J::Vector; - seed=rand(UInt), fout=nothing, - display_beforeevent=false, display_afterevent=false, - kwargs...) +function mcwf_nh(tspan, psi0::T, Hnh::AbstractOperator{B,B}, J::Vector; + seed=rand(UInt), fout=nothing, + display_beforeevent=false, display_afterevent=false, + kwargs...) where {B<:Basis,T<:Ket{B}} check_mcwf(psi0, Hnh, J, J, nothing) - f(t, psi, dpsi) = dmcwf_nh(psi, Hnh, dpsi) - j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, nothing) - return integrate_mcwf(f, j, tspan, psi0, seed, fout; - display_beforeevent=display_beforeevent, - display_afterevent=display_afterevent, - kwargs...) + f(t::Float64, psi::T, dpsi::T) = dmcwf_nh(psi, Hnh, dpsi) + j(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, nothing) + integrate_mcwf(f, j, tspan, psi0, seed, fout; + display_beforeevent=display_beforeevent, + display_afterevent=display_afterevent, + kwargs...) end """ @@ -97,21 +97,21 @@ operators. If they are not given they are calculated automatically. * `display_afterevent=false`: `fout` is called after every jump. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function mcwf(tspan, psi0::Ket, H::Operator, J::Vector; - seed=rand(UInt), rates::DecayRates=nothing, - fout=nothing, Jdagger::Vector=dagger.(J), - display_beforeevent=false, display_afterevent=false, - kwargs...) +function mcwf(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector; + seed=rand(UInt), rates::DecayRates=nothing, + fout=nothing, Jdagger::Vector=dagger.(J), + display_beforeevent=false, display_afterevent=false, + kwargs...) where {B<:Basis,T<:Ket{B}} isreducible = check_mcwf(psi0, H, J, Jdagger, rates) if !isreducible tmp = copy(psi0) - dmcwf_h_(t, psi, dpsi) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates) - j_h(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates) - return integrate_mcwf(dmcwf_h_, j_h, tspan, psi0, seed, - fout; - display_beforeevent=display_beforeevent, - display_afterevent=display_afterevent, - kwargs...) + dmcwf_h_(t::Float64, psi::T, dpsi::T) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates) + j_h(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates) + integrate_mcwf(dmcwf_h_, j_h, tspan, psi0, seed, + fout; + display_beforeevent=display_beforeevent, + display_afterevent=display_afterevent, + kwargs...) else Hnh = copy(H) if typeof(rates) == Nothing @@ -123,13 +123,13 @@ function mcwf(tspan, psi0::Ket, H::Operator, J::Vector; Hnh -= 0.5im*rates[i]*Jdagger[i]*J[i] end end - dmcwf_nh_(t, psi, dpsi) = dmcwf_nh(psi, Hnh, dpsi) - j_nh(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates) - return integrate_mcwf(dmcwf_nh_, j_nh, tspan, psi0, seed, - fout; - display_beforeevent=display_beforeevent, - display_afterevent=display_afterevent, - kwargs...) + dmcwf_nh_(t::Float64, psi::T, dpsi::T) = dmcwf_nh(psi, Hnh, dpsi) + j_nh(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates) + integrate_mcwf(dmcwf_nh_, j_nh, tspan, psi0, seed, + fout; + display_beforeevent=display_beforeevent, + display_afterevent=display_afterevent, + kwargs...) end end @@ -159,13 +159,13 @@ and therefore must not be changed. * `display_afterevent=false`: `fout` is called after every jump. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function mcwf_dynamic(tspan, psi0::Ket, f::Function; +function mcwf_dynamic(tspan, psi0::T, f::Function; seed=rand(UInt), rates::DecayRates=nothing, fout=nothing, display_beforeevent=false, display_afterevent=false, - kwargs...) + kwargs...) where {T<:Ket} tmp = copy(psi0) - dmcwf_(t, psi, dpsi) = dmcwf_h_dynamic(t, psi, f, rates, dpsi, tmp) - j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, f, psi_new, rates) + dmcwf_(t::Float64, psi::T, dpsi::T) = dmcwf_h_dynamic(t, psi, f, rates, dpsi, tmp) + j_(rng, t::Float64, psi::T, psi_new::T) = jump_dynamic(rng, t, psi, f, psi_new, rates) integrate_mcwf(dmcwf_, j_, tspan, psi0, seed, fout; display_beforeevent=display_beforeevent, @@ -180,12 +180,12 @@ Calculate MCWF trajectory where the dynamic Hamiltonian is given in non-hermitia For more information see: [`mcwf_dynamic`](@ref) """ -function mcwf_nh_dynamic(tspan, psi0::Ket, f::Function; +function mcwf_nh_dynamic(tspan, psi0::T, f::Function; seed=rand(UInt), rates::DecayRates=nothing, fout=nothing, display_beforeevent=false, display_afterevent=false, - kwargs...) - dmcwf_(t, psi, dpsi) = dmcwf_nh_dynamic(t, psi, f, dpsi) - j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, f, psi_new, rates) + kwargs...) where T<:Ket + dmcwf_(t::Float64, psi::T, dpsi::T) = dmcwf_nh_dynamic(t, psi, f, dpsi) + j_(rng, t::Float64, psi::T, psi_new::T) = jump_dynamic(rng, t, psi, f, psi_new, rates) integrate_mcwf(dmcwf_, j_, tspan, psi0, seed, fout; display_beforeevent=display_beforeevent, @@ -193,8 +193,8 @@ function mcwf_nh_dynamic(tspan, psi0::Ket, f::Function; kwargs...) end -function dmcwf_h_dynamic(t::Float64, psi::Ket, f::Function, rates::DecayRates, - dpsi::Ket, tmp::Ket) +function dmcwf_h_dynamic(t::Float64, psi::T, f::Function, rates::DecayRates, + dpsi::T, tmp::T) where T<:Ket result = f(t, psi) QO_CHECKS[] && @assert 3 <= length(result) <= 4 if length(result) == 3 @@ -207,7 +207,7 @@ function dmcwf_h_dynamic(t::Float64, psi::Ket, f::Function, rates::DecayRates, dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates) end -function dmcwf_nh_dynamic(t::Float64, psi::Ket, f::Function, dpsi::Ket) +function dmcwf_nh_dynamic(t::Float64, psi::T, f::Function, dpsi::T) where T<:Ket result = f(t, psi) QO_CHECKS[] && @assert 3 <= length(result) <= 4 H, J, Jdagger = result[1:3] @@ -215,7 +215,7 @@ function dmcwf_nh_dynamic(t::Float64, psi::Ket, f::Function, dpsi::Ket) dmcwf_nh(psi, H, dpsi) end -function jump_dynamic(rng, t::Float64, psi::Ket, f::Function, psi_new::Ket, rates::DecayRates) +function jump_dynamic(rng, t::Float64, psi::T, f::Function, psi_new::T, rates::DecayRates) where T<:Ket result = f(t, psi) QO_CHECKS[] && @assert 3 <= length(result) <= 4 J = result[2] @@ -224,7 +224,7 @@ function jump_dynamic(rng, t::Float64, psi::Ket, f::Function, psi_new::Ket, rate else rates_ = result[4] end - jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket, rates::DecayRates) + jump(rng, t, psi, J, psi_new, rates) end """ @@ -249,19 +249,19 @@ Integrate a single Monte Carlo wave function trajectory. * `kwargs`: Further arguments are passed on to the ode solver. """ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, - psi0::Ket, seed, fout::Function; + psi0::T, seed, fout::Function; display_beforeevent=false, display_afterevent=false, #TODO: Remove kwargs save_everystep=false, callback=nothing, alg=OrdinaryDiffEq.DP5(), - kwargs...) + kwargs...) where {B<:Basis,D<:Vector{ComplexF64},T<:Ket{B,D}} tmp = copy(psi0) - as_ket(x::Vector{ComplexF64}) = Ket(psi0.basis, x) - as_vector(psi::Ket) = psi.data + as_ket(x::D) = T(psi0.basis, x) + as_vector(psi::T) = psi.data rng = MersenneTwister(convert(UInt, seed)) jumpnorm = Ref(rand(rng)) - djumpnorm(x::Vector{ComplexF64}, t, integrator) = norm(as_ket(x))^2 - (1-jumpnorm[]) + djumpnorm(x::D, t::Float64, integrator) = norm(as_ket(x))^2 - (1-jumpnorm[]) if !display_beforeevent && !display_afterevent function dojump(integrator) @@ -275,14 +275,14 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, save_positions = (display_beforeevent,display_afterevent)) - return timeevolution.integrate(float(tspan), dmcwf, as_vector(psi0), + timeevolution.integrate(float(tspan), dmcwf, as_vector(psi0), copy(psi0), copy(psi0), fout; callback = cb, kwargs...) else # Temporary workaround until proper tooling for saving # TODO: Replace by proper call to timeevolution.integrate - function fout_(x::Vector{ComplexF64}, t::Float64, integrator) + function fout_(x::D, t::Float64, integrator) recast!(x, state) fout(t, state) end @@ -324,7 +324,7 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, save_positions = (false,false)) full_cb = OrdinaryDiffEq.CallbackSet(callback,cb,scb) - function df_(dx::Vector{ComplexF64}, x::Vector{ComplexF64}, p, t) + function df_(dx::D, x::D, p, t) recast!(x, state) recast!(dx, dstate) dmcwf(t, state, dstate) @@ -346,9 +346,9 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, end function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, - psi0::Ket, seed, fout::Nothing; - kwargs...) - function fout_(t, x) + psi0::T, seed, fout::Nothing; + kwargs...) where {T<:Ket} + function fout_(t::Float64, x::T) psi = copy(x) psi /= norm(psi) return psi @@ -368,7 +368,7 @@ Default jump function. * `J`: List of jump operators. * `psi_new`: Result of jump. """ -function jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket, rates::Nothing) +function jump(rng, t::Float64, psi::T, J::Vector, psi_new::T, rates::Nothing) where T<:Ket if length(J)==1 operators.gemv!(complex(1.), J[1], psi, complex(0.), psi_new) psi_new.data ./= norm(psi_new) @@ -386,7 +386,7 @@ function jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket, rates::Nothing return nothing end -function jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket, rates::Vector{Float64}) +function jump(rng, t::Float64, psi::T, J::Vector, psi_new::T, rates::Vector{Float64}) where T<:Ket if length(J)==1 operators.gemv!(complex(sqrt(rates[1])), J[1], psi, complex(0.), psi_new) psi_new.data ./= norm(psi_new) @@ -410,8 +410,8 @@ Evaluate non-hermitian Schroedinger equation. The non-hermitian Hamiltonian is given in two parts - the hermitian part H and the jump operators J. """ -function dmcwf_h(psi::Ket, H::Operator, - J::Vector, Jdagger::Vector, dpsi::Ket, tmp::Ket, rates::Nothing) +function dmcwf_h(psi::T, H::AbstractOperator{B,B}, + J::Vector, Jdagger::Vector, dpsi::T, tmp::T, rates::Nothing) where {B<:Basis,T<:Ket{B}} operators.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi) for i=1:length(J) operators.gemv!(complex(1.), J[i], psi, complex(0.), tmp) @@ -420,8 +420,8 @@ function dmcwf_h(psi::Ket, H::Operator, return dpsi end -function dmcwf_h(psi::Ket, H::Operator, - J::Vector, Jdagger::Vector, dpsi::Ket, tmp::Ket, rates::Vector{Float64}) +function dmcwf_h(psi::T, H::AbstractOperator{B,B}, + J::Vector, Jdagger::Vector, dpsi::T, tmp::T, rates::Vector{Float64}) where {B<:Basis,T<:Ket{B}} operators.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi) for i=1:length(J) operators.gemv!(complex(rates[i]), J[i], psi, complex(0.), tmp) @@ -436,7 +436,7 @@ Evaluate non-hermitian Schroedinger equation. The given Hamiltonian is already the non-hermitian version. """ -function dmcwf_nh(psi::Ket, Hnh::Operator, dpsi::Ket) +function dmcwf_nh(psi::T, Hnh::AbstractOperator{B,B}, dpsi::T) where {B<:Basis,T<:Ket{B}} operators.gemv!(complex(0,-1.), Hnh, psi, complex(0.), dpsi) return dpsi end @@ -446,25 +446,23 @@ end Check input of mcwf. """ -function check_mcwf(psi0::Ket, H::Operator, J::Vector, Jdagger::Vector, rates::DecayRates) +function check_mcwf(psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis + # TODO: replace type checks by dispatch; make types of J known isreducible = true - check_samebases(basis(psi0), basis(H)) if !(isa(H, DenseOperator) || isa(H, SparseOperator)) isreducible = false end for j=J - @assert isa(j, Operator) + @assert isa(j, AbstractOperator{B,B}) if !(isa(j, DenseOperator) || isa(j, SparseOperator)) isreducible = false end - check_samebases(H, j) end for j=Jdagger - @assert isa(j, Operator) + @assert isa(j, AbstractOperator{B,B}) if !(isa(j, DenseOperator) || isa(j, SparseOperator)) isreducible = false end - check_samebases(H, j) end @assert length(J) == length(Jdagger) if typeof(rates) == Matrix{Float64} @@ -488,13 +486,13 @@ corresponding set of jump operators is calculated. * `rates`: Matrix of decay rates. * `J`: Vector of jump operators. """ -function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where T <: Operator +function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where {B<:Basis,T<:AbstractOperator{B,B}} @assert length(J) == size(rates)[1] == size(rates)[2] d, v = eigen(rates) d, [sum([v[j, i]*J[j] for j=1:length(d)]) for i=1:length(d)] end -function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where T <: Union{LazySum, LazyTensor, LazyProduct} +function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where {B<:Basis,T<:Union{LazySum{B,B},LazyTensor{B,B},LazyProduct{B,B}}} @assert length(J) == size(rates)[1] == size(rates)[2] d, v = eigen(rates) d, [LazySum([v[j, i]*J[j] for j=1:length(d)]...) for i=1:length(d)] diff --git a/src/metrics.jl b/src/metrics.jl index e307ecd1..57b6023b 100644 --- a/src/metrics.jl +++ b/src/metrics.jl @@ -22,12 +22,11 @@ T(ρ) = Tr\\{\\sqrt{ρ^† ρ}\\}. Depending if `rho` is hermitian either [`tracenorm_h`](@ref) or [`tracenorm_nh`](@ref) is called. """ -function tracenorm(rho::DenseOperator) - check_samebases(rho) +function tracenorm(rho::DenseOperator{B,B}) where B<:Basis ishermitian(rho) ? tracenorm_h(rho) : tracenorm_nh(rho) end -function tracenorm(rho::T) where T<:Operator - throw(ArgumentError("tracenorm not implemented for $(T). Use dense operators instead.")) +function tracenorm(rho::T) where T<:AbstractOperator + throw(ArgumentError("tracenorm not implemented for $(typeof(rho)). Use dense operators instead.")) end """ @@ -43,13 +42,12 @@ T(ρ) = Tr\\{\\sqrt{ρ^† ρ}\\} = \\sum_i |λ_i| where ``λ_i`` are the eigenvalues of `rho`. """ -function tracenorm_h(rho::DenseOperator) - check_samebases(rho) +function tracenorm_h(rho::DenseOperator{B,B}) where B<:Basis s = eigvals(Hermitian(rho.data)) sum(abs.(s)) end -function tracenorm_h(rho::T) where T<:Operator - throw(ArgumentError("tracenorm_h not implemented for $(T). Use dense operators instead.")) +function tracenorm_h(rho::T) where T<:AbstractOperator + throw(ArgumentError("tracenorm_h not implemented for $(typeof(rho)). Use dense operators instead.")) end @@ -70,8 +68,8 @@ It uses the identity where ``σ_i`` are the singular values of `rho`. """ tracenorm_nh(rho::DenseOperator) = sum(svdvals(rho.data)) -function tracenorm_nh(rho::T) where T<:Operator - throw(ArgumentError("tracenorm_nh not implemented for $(T). Use dense operators instead.")) +function tracenorm_nh(rho::T) where T<:AbstractOperator + throw(ArgumentError("tracenorm_nh not implemented for $(typeof(rho)). Use dense operators instead.")) end @@ -89,9 +87,9 @@ T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\}. It calls [`tracenorm`](@ref) which in turn either uses [`tracenorm_h`](@ref) or [`tracenorm_nh`](@ref) depending if ``ρ-σ`` is hermitian or not. """ -tracedistance(rho::DenseOperator, sigma::DenseOperator) = 0.5*tracenorm(rho - sigma) -function tracedistance(rho::T, sigma::T) where T<:Operator - throw(ArgumentError("tracedistance not implemented for $(T). Use dense operators instead.")) +tracedistance(rho::T, sigma::T) where T<:DenseOperator = 0.5*tracenorm(rho - sigma) +function tracedistance(rho::AbstractOperator, sigma::AbstractOperator) + throw(ArgumentError("tracedistance not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead.")) end """ @@ -107,9 +105,9 @@ T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\} = \\frac{1}{2} \ where ``λ_i`` are the eigenvalues of `rho` - `sigma`. """ -tracedistance_h(rho::DenseOperator, sigma::DenseOperator) = 0.5*tracenorm_h(rho - sigma) -function tracedistance_h(rho::T, sigma::T) where T<:Operator - throw(ArgumentError("tracedistance_h not implemented for $(T). Use dense operators instead.")) +tracedistance_h(rho::T, sigma::T) where {B<:Basis,T<:DenseOperator{B,B}}= 0.5*tracenorm_h(rho - sigma) +function tracedistance_h(rho::AbstractOperator, sigma::AbstractOperator) + throw(ArgumentError("tracedistance_h not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead.")) end """ @@ -129,9 +127,9 @@ It uses the identity where ``σ_i`` are the singular values of `rho` - `sigma`. """ -tracedistance_nh(rho::DenseOperator, sigma::DenseOperator) = 0.5*tracenorm_nh(rho - sigma) -function tracedistance_nh(rho::T, sigma::T) where T<:Operator - throw(ArgumentError("tracedistance_nh not implemented for $(T). Use dense operators instead.")) +tracedistance_nh(rho::T, sigma::T) where {T<:DenseOperator} = 0.5*tracenorm_nh(rho - sigma) +function tracedistance_nh(rho::AbstractOperator, sigma::AbstractOperator) + throw(ArgumentError("tracedistance_nh not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead.")) end @@ -153,7 +151,7 @@ natural logarithm and ``0\\log(0) ≡ 0``. * `rho`: Density operator of which to calculate Von Neumann entropy. * `tol=1e-15`: Tolerance for rounding errors in the computed eigenvalues. """ -function entropy_vn(rho::DenseOperator; tol::Float64=1e-15) +function entropy_vn(rho::DenseOperator{B,B}; tol::Float64=1e-15) where B<:Basis evals = eigvals(rho.data) evals[abs.(evals) .< tol] .= 0.0 sum([d == 0 ? 0 : -d*log(d) for d=evals]) @@ -173,7 +171,7 @@ F(ρ, σ) = Tr\\left(\\sqrt{\\sqrt{ρ}σ\\sqrt{ρ}}\\right), where ``\\sqrt{ρ}=\\sum_n\\sqrt{λ_n}|ψ⟩⟨ψ|``. """ -fidelity(rho::DenseOperator, sigma::DenseOperator) = tr(sqrt(sqrt(rho.data)*sigma.data*sqrt(rho.data))) +fidelity(rho::T, sigma::T) where {B<:Basis,T<:DenseOperator{B,B}} = tr(sqrt(sqrt(rho.data)*sigma.data*sqrt(rho.data))) """ @@ -181,9 +179,7 @@ fidelity(rho::DenseOperator, sigma::DenseOperator) = tr(sqrt(sqrt(rho.data)*sigm Partial transpose of rho with respect to subspace specified by index. """ -function ptranspose(rho::DenseOperator, index::Int=1) - @assert rho.basis_l == rho.basis_r - @assert typeof(rho.basis_l) == CompositeBasis +function ptranspose(rho::DenseOperator{B,B}, index::Int=1) where B<:CompositeBasis # Define permutation N = length(rho.basis_l.bases) @@ -210,7 +206,7 @@ end Peres-Horodecki criterion of partial transpose. """ -PPT(rho::DenseOperator, index::Int) = all(real.(eigvals(ptranspose(rho, index).data)) .>= 0.0) +PPT(rho::DenseOperator{B,B}, index::Int) where B<:CompositeBasis = all(real.(eigvals(ptranspose(rho, index).data)) .>= 0.0) """ @@ -225,7 +221,7 @@ N(ρ) = \\|ρᵀ\\|, ``` where `ρᵀ` is the partial transpose. """ -negativity(rho::DenseOperator, index::Int) = 0.5*(tracenorm(ptranspose(rho, index)) - 1.0) +negativity(rho::DenseOperator{B,B}, index::Int) where B<:CompositeBasis = 0.5*(tracenorm(ptranspose(rho, index)) - 1.0) """ @@ -238,6 +234,6 @@ N(ρ) = \\log₂\\|ρᵀ\\|, ``` where `ρᵀ` is the partial transpose. """ -logarithmic_negativity(rho::DenseOperator, index::Int) = log(2, tracenorm(ptranspose(rho, index))) +logarithmic_negativity(rho::DenseOperator{B,B}, index::Int) where B<:CompositeBasis = log(2, tracenorm(ptranspose(rho, index))) end # module diff --git a/src/operators.jl b/src/operators.jl index c0c1a309..f9bcfaf1 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -1,6 +1,6 @@ module operators -export Operator, length, basis, dagger, ishermitian, tensor, embed, +export AbstractOperator, length, basis, dagger, ishermitian, tensor, embed, tr, ptrace, normalize, normalize!, expect, variance, exp, permutesystems, identityoperator, dense @@ -20,59 +20,59 @@ All deriving operator classes have to define the fields `basis_l` and `basis_r` defining the left and right side bases. For fast time evolution also at least the function -`gemv!(alpha, op::Operator, x::Ket, beta, result::Ket)` should be +`gemv!(alpha, op::AbstractOperator, x::Ket, beta, result::Ket)` should be implemented. Many other generic multiplication functions can be defined in terms of this function and are provided automatically. """ -abstract type Operator end +abstract type AbstractOperator{BL<:Basis,BR<:Basis} end # Common error messages -arithmetic_unary_error(funcname, x::Operator) = throw(ArgumentError("$funcname is not defined for this type of operator: $(typeof(x)).\nTry to convert to another operator type first with e.g. dense() or sparse().")) -arithmetic_binary_error(funcname, a::Operator, b::Operator) = throw(ArgumentError("$funcname is not defined for this combination of types of operators: $(typeof(a)), $(typeof(b)).\nTry to convert to a common operator type first with e.g. dense() or sparse().")) +arithmetic_unary_error(funcname, x::AbstractOperator) = throw(ArgumentError("$funcname is not defined for this type of operator: $(typeof(x)).\nTry to convert to another operator type first with e.g. dense() or sparse().")) +arithmetic_binary_error(funcname, a::AbstractOperator, b::AbstractOperator) = throw(ArgumentError("$funcname is not defined for this combination of types of operators: $(typeof(a)), $(typeof(b)).\nTry to convert to a common operator type first with e.g. dense() or sparse().")) addnumbererror() = throw(ArgumentError("Can't add or subtract a number and an operator. You probably want 'op + identityoperator(op)*x'.")) -length(a::Operator) = length(a.basis_l)::Int*length(a.basis_r)::Int -basis(a::Operator) = (check_samebases(a); a.basis_l) +length(a::AbstractOperator) = length(a.basis_l)::Int*length(a.basis_r)::Int +basis(a::AbstractOperator) = (check_samebases(a); a.basis_l) # Arithmetic operations -+(a::Operator, b::Operator) = arithmetic_binary_error("Addition", a, b) -+(a::Number, b::Operator) = addnumbererror() -+(a::Operator, b::Number) = addnumbererror() ++(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Addition", a, b) ++(a::Number, b::AbstractOperator) = addnumbererror() ++(a::AbstractOperator, b::Number) = addnumbererror() --(a::Operator) = arithmetic_unary_error("Negation", a) --(a::Operator, b::Operator) = arithmetic_binary_error("Subtraction", a, b) --(a::Number, b::Operator) = addnumbererror() --(a::Operator, b::Number) = addnumbererror() +-(a::AbstractOperator) = arithmetic_unary_error("Negation", a) +-(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Subtraction", a, b) +-(a::Number, b::AbstractOperator) = addnumbererror() +-(a::AbstractOperator, b::Number) = addnumbererror() -*(a::Operator, b::Operator) = arithmetic_binary_error("Multiplication", a, b) -^(a::Operator, b::Int) = Base.power_by_squaring(a, b) +*(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Multiplication", a, b) +^(a::AbstractOperator, b::Int) = Base.power_by_squaring(a, b) -dagger(a::Operator) = arithmetic_unary_error("Hermitian conjugate", a) -Base.adjoint(a::Operator) = dagger(a) +dagger(a::AbstractOperator) = arithmetic_unary_error("Hermitian conjugate", a) +Base.adjoint(a::AbstractOperator) = dagger(a) -conj(a::Operator) = arithmetic_unary_error("Complex conjugate", a) -conj!(a::Operator) = conj(a::Operator) +conj(a::AbstractOperator) = arithmetic_unary_error("Complex conjugate", a) +conj!(a::AbstractOperator) = conj(a::AbstractOperator) -dense(a::Operator) = arithmetic_unary_error("Conversion to dense", a) +dense(a::AbstractOperator) = arithmetic_unary_error("Conversion to dense", a) """ - ishermitian(op::Operator) + ishermitian(op::AbstractOperator) Check if an operator is Hermitian. """ -ishermitian(op::Operator) = arithmetic_unary_error(ishermitian, op) +ishermitian(op::AbstractOperator) = arithmetic_unary_error(ishermitian, op) """ - tensor(x::Operator, y::Operator, z::Operator...) + tensor(x::AbstractOperator, y::AbstractOperator, z::AbstractOperator...) Tensor product ``\\hat{x}⊗\\hat{y}⊗\\hat{z}⊗…`` of the given operators. """ -tensor(a::Operator, b::Operator) = arithmetic_binary_error("Tensor product", a, b) -tensor(op::Operator) = op -tensor(operators::Operator...) = reduce(tensor, operators) +tensor(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Tensor product", a, b) +tensor(op::AbstractOperator) = op +tensor(operators::AbstractOperator...) = reduce(tensor, operators) """ @@ -81,25 +81,25 @@ tensor(operators::Operator...) = reduce(tensor, operators) Tensor product of operators where missing indices are filled up with identity operators. """ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, - indices::Vector{Int}, operators::Vector{T}) where T<:Operator + indices::Vector{Int}, operators::Vector{T}) where T<:AbstractOperator N = length(basis_l.bases) @assert length(basis_r.bases) == N @assert length(indices) == length(operators) sortedindices.check_indices(N, indices) tensor([i ∈ indices ? operators[indexin(i, indices)[1]] : identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...) end -embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Int, op::Operator) = embed(basis_l, basis_r, Int[index], [op]) -embed(basis::CompositeBasis, index::Int, op::Operator) = embed(basis, basis, Int[index], [op]) -embed(basis::CompositeBasis, indices::Vector{Int}, operators::Vector{T}) where {T<:Operator} = embed(basis, basis, indices, operators) +embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Int, op::AbstractOperator) = embed(basis_l, basis_r, Int[index], [op]) +embed(basis::CompositeBasis, index::Int, op::AbstractOperator) = embed(basis, basis, Int[index], [op]) +embed(basis::CompositeBasis, indices::Vector{Int}, operators::Vector{T}) where {T<:AbstractOperator} = embed(basis, basis, indices, operators) """ embed(basis1[, basis2], operators::Dict) -`operators` is a dictionary `Dict{Vector{Int}, Operator}`. The integer vector +`operators` is a dictionary `Dict{Vector{Int}, AbstractOperator}`. The integer vector specifies in which subsystems the corresponding operator is defined. """ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, - operators::Dict{Vector{Int}, T}) where T<:Operator + operators::Dict{Vector{Int}, T}) where T<:AbstractOperator @assert length(basis_l.bases) == length(basis_r.bases) N = length(basis_l.bases) if length(operators) == 0 @@ -127,33 +127,33 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, return permutesystems(op, perm) end end -embed(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:Operator} = embed(basis_l, basis_r, Dict([i]=>op_i for (i, op_i) in operators); kwargs...) -embed(basis::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:Operator} = embed(basis, basis, operators; kwargs...) -embed(basis::CompositeBasis, operators::Dict{Vector{Int}, T}; kwargs...) where {T<:Operator} = embed(basis, basis, operators; kwargs...) +embed(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:AbstractOperator} = embed(basis_l, basis_r, Dict([i]=>op_i for (i, op_i) in operators); kwargs...) +embed(basis::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:AbstractOperator} = embed(basis, basis, operators; kwargs...) +embed(basis::CompositeBasis, operators::Dict{Vector{Int}, T}; kwargs...) where {T<:AbstractOperator} = embed(basis, basis, operators; kwargs...) """ - tr(x::Operator) + tr(x::AbstractOperator) Trace of the given operator. """ -tr(x::Operator) = arithmetic_unary_error("Trace", x) +tr(x::AbstractOperator) = arithmetic_unary_error("Trace", x) -ptrace(a::Operator, index::Vector{Int}) = arithmetic_unary_error("Partial trace", a) +ptrace(a::AbstractOperator, index::Vector{Int}) = arithmetic_unary_error("Partial trace", a) """ normalize(op) Return the normalized operator so that its `tr(op)` is one. """ -normalize(op::Operator) = op/tr(op) +normalize(op::AbstractOperator) = op/tr(op) """ normalize!(op) In-place normalization of the given operator so that its `tr(x)` is one. """ -normalize!(op::Operator) = throw(ArgumentError("normalize! is not defined for this type of operator: $(typeof(op)).\n You may have to fall back to the non-inplace version 'normalize()'.")) +normalize!(op::AbstractOperator) = throw(ArgumentError("normalize! is not defined for this type of operator: $(typeof(op)).\n You may have to fall back to the non-inplace version 'normalize()'.")) """ expect(op, state) @@ -162,27 +162,27 @@ Expectation value of the given operator `op` for the specified `state`. `state` can either be a (density) operator or a ket. """ -expect(op::Operator, state::Ket) = dagger(state) * op * state -expect(op::Operator, state::Operator) = tr(op*state) +expect(op::AbstractOperator, state::Ket) = dagger(state) * op * state +expect(op::AbstractOperator, state::AbstractOperator) = tr(op*state) """ expect(index, op, state) If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number. """ -function expect(indices::Vector{Int}, op::Operator, state::Operator) +function expect(indices::Vector{Int}, op::AbstractOperator, state::AbstractOperator) N = length(state.basis_l.shape) indices_ = sortedindices.complement(N, indices) expect(op, ptrace(state, indices_)) end -function expect(indices::Vector{Int}, op::Operator, state::Ket) +function expect(indices::Vector{Int}, op::AbstractOperator, state::Ket) N = length(state.basis.shape) indices_ = sortedindices.complement(N, indices) expect(op, ptrace(state, indices_)) end -expect(index::Int, op::Operator, state) = expect([index], op, state) -expect(op::Operator, states::Vector) = [expect(op, state) for state=states] -expect(indices::Vector{Int}, op::Operator, states::Vector) = [expect(indices, op, state) for state=states] +expect(index::Int, op::AbstractOperator, state) = expect([index], op, state) +expect(op::AbstractOperator, states::Vector) = [expect(op, state) for state=states] +expect(indices::Vector{Int}, op::AbstractOperator, states::Vector) = [expect(indices, op, state) for state=states] """ variance(op, state) @@ -191,12 +191,12 @@ Variance of the given operator `op` for the specified `state`. `state` can either be a (density) operator or a ket. """ -function variance(op::Operator, state::Ket) +function variance(op::AbstractOperator, state::Ket) x = op*state stateT = dagger(state) stateT*op*x - (stateT*x)^2 end -function variance(op::Operator, state::Operator) +function variance(op::AbstractOperator, state::AbstractOperator) expect(op*op, state) - expect(op, state)^2 end @@ -205,41 +205,41 @@ end If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number """ -function variance(indices::Vector{Int}, op::Operator, state::Operator) +function variance(indices::Vector{Int}, op::AbstractOperator, state::AbstractOperator) N = length(state.basis_l.shape) indices_ = sortedindices.complement(N, indices) variance(op, ptrace(state, indices_)) end -function variance(indices::Vector{Int}, op::Operator, state::Ket) +function variance(indices::Vector{Int}, op::AbstractOperator, state::Ket) N = length(state.basis.shape) indices_ = sortedindices.complement(N, indices) variance(op, ptrace(state, indices_)) end -variance(index::Int, op::Operator, state) = variance([index], op, state) -variance(op::Operator, states::Vector) = [variance(op, state) for state=states] -variance(indices::Vector{Int}, op::Operator, states::Vector) = [variance(indices, op, state) for state=states] +variance(index::Int, op::AbstractOperator, state) = variance([index], op, state) +variance(op::AbstractOperator, states::Vector) = [variance(op, state) for state=states] +variance(indices::Vector{Int}, op::AbstractOperator, states::Vector) = [variance(indices, op, state) for state=states] """ - exp(op::Operator) + exp(op::AbstractOperator) Operator exponential. """ -exp(op::Operator) = throw(ArgumentError("exp() is not defined for this type of operator: $(typeof(op)).\nTry to convert to dense operator first with dense().")) +exp(op::AbstractOperator) = throw(ArgumentError("exp() is not defined for this type of operator: $(typeof(op)).\nTry to convert to dense operator first with dense().")) -permutesystems(a::Operator, perm::Vector{Int}) = arithmetic_unary_error("Permutations of subsystems", a) +permutesystems(a::AbstractOperator, perm::Vector{Int}) = arithmetic_unary_error("Permutations of subsystems", a) """ identityoperator(a::Basis[, b::Basis]) Return an identityoperator in the given bases. """ -identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:Operator} = throw(ArgumentError("Identity operator not defined for operator type $T.")) -identityoperator(::Type{T}, b::Basis) where {T<:Operator} = identityoperator(T, b, b) -identityoperator(op::T) where {T<:Operator} = identityoperator(T, op.basis_l, op.basis_r) +identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:AbstractOperator} = throw(ArgumentError("Identity operator not defined for operator type $T.")) +identityoperator(::Type{T}, b::Basis) where {T<:AbstractOperator} = identityoperator(T, b, b) +identityoperator(op::T) where {T<:AbstractOperator} = identityoperator(T, op.basis_l, op.basis_r) one(b::Basis) = identityoperator(b) -one(op::Operator) = identityoperator(op) +one(op::AbstractOperator) = identityoperator(op) # Fast in-place multiplication @@ -265,7 +265,7 @@ gemm!() = error("Not Implemented.") # Helper functions to check validity of arguments -function check_ptrace_arguments(a::Operator, indices::Vector{Int}) +function check_ptrace_arguments(a::AbstractOperator, indices::Vector{Int}) if !isa(a.basis_l, CompositeBasis) || !isa(a.basis_r, CompositeBasis) throw(ArgumentError("Partial trace can only be applied onto operators with composite bases.")) end @@ -290,12 +290,12 @@ function check_ptrace_arguments(a::StateVector, indices::Vector{Int}) sortedindices.check_indices(length(basis(a).shape), indices) end -samebases(a::Operator) = samebases(a.basis_l, a.basis_r)::Bool -samebases(a::Operator, b::Operator) = samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool -check_samebases(a::Operator) = check_samebases(a.basis_l, a.basis_r) +samebases(a::AbstractOperator) = samebases(a.basis_l, a.basis_r)::Bool +samebases(a::AbstractOperator, b::AbstractOperator) = samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool +check_samebases(a::AbstractOperator) = check_samebases(a.basis_l, a.basis_r) -multiplicable(a::Operator, b::Ket) = multiplicable(a.basis_r, b.basis) -multiplicable(a::Bra, b::Operator) = multiplicable(a.basis, b.basis_l) -multiplicable(a::Operator, b::Operator) = multiplicable(a.basis_r, b.basis_l) +multiplicable(a::AbstractOperator, b::Ket) = multiplicable(a.basis_r, b.basis) +multiplicable(a::Bra, b::AbstractOperator) = multiplicable(a.basis, b.basis_l) +multiplicable(a::AbstractOperator, b::AbstractOperator) = multiplicable(a.basis_r, b.basis_l) end # module diff --git a/src/operators_dense.jl b/src/operators_dense.jl index 50714a5a..f4bed95d 100644 --- a/src/operators_dense.jl +++ b/src/operators_dense.jl @@ -16,62 +16,73 @@ Dense array implementation of Operator. The matrix consisting of complex floats is stored in the `data` field. """ -mutable struct DenseOperator <: Operator - basis_l::Basis - basis_r::Basis +mutable struct DenseOperator{BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR data::Matrix{ComplexF64} - DenseOperator(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new(b1, b2, data) : throw(DimensionMismatch()) + function DenseOperator{BL,BR,T}(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} + if !(length(b1) == size(data, 1) && length(b2) == size(data, 2)) + throw(DimensionMismatch()) + end + new(b1, b2, data) + end end +DenseOperator{BL,BR}(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{BL,BR,T}(b1, b2, data) +DenseOperator(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{BL,BR,T}(b1, b2, data) +DenseOperator(b1::Basis, b2::Basis, data) = DenseOperator(b1, b2, convert(Matrix{ComplexF64}, data)) DenseOperator(b::Basis, data) = DenseOperator(b, b, data) DenseOperator(b1::Basis, b2::Basis) = DenseOperator(b1, b2, zeros(ComplexF64, length(b1), length(b2))) +DenseOperator{B1,B2}(b1::B1, b2::B2) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(b1, b2, zeros(ComplexF64, length(b1), length(b2))) DenseOperator(b::Basis) = DenseOperator(b, b) -DenseOperator(op::Operator) = dense(op) +DenseOperator(op::AbstractOperator) = dense(op) Base.copy(x::DenseOperator) = DenseOperator(x.basis_l, x.basis_r, copy(x.data)) """ - dense(op::Operator) + dense(op::AbstractOperator) Convert an arbitrary Operator into a [`DenseOperator`](@ref). """ operators.dense(x::DenseOperator) = copy(x) -==(x::DenseOperator, y::DenseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && (x.data == y.data) +==(x::DenseOperator, y::DenseOperator) = false +==(x::T, y::T) where T<:DenseOperator = (x.data == y.data) # Arithmetic operations -+(a::DenseOperator, b::DenseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data+b.data)) - --(a::DenseOperator) = DenseOperator(a.basis_l, a.basis_r, -a.data) --(a::DenseOperator, b::DenseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data-b.data)) - -*(a::DenseOperator, b::Ket) = (check_multiplicable(a, b); Ket(a.basis_l, a.data*b.data)) -*(a::Bra, b::DenseOperator) = (check_multiplicable(a, b); Bra(b.basis_r, transpose(b.data)*a.data)) -*(a::DenseOperator, b::DenseOperator) = (check_multiplicable(a, b); DenseOperator(a.basis_l, b.basis_r, a.data*b.data)) ++(a::T, b::T) where T<:DenseOperator = T(a.basis_l, a.basis_r, a.data+b.data) ++(a::DenseOperator, b::DenseOperator) = throw(bases.IncompatibleBases()) + +-(a::T) where T<:DenseOperator = T(a.basis_l, a.basis_r, -a.data) +-(a::T, b::T) where T<:DenseOperator = T(a.basis_l, a.basis_r, a.data-b.data) +-(a::DenseOperator, b::DenseOperator) = throw(bases.IncompatibleBases()) + +*(a::DenseOperator{BL,BR}, b::Ket{BR}) where {BL<:Basis,BR<:Basis} = Ket{BL}(a.basis_l, a.data*b.data) +*(a::DenseOperator, b::Ket) = throw(bases.IncompatibleBases()) +*(a::Bra{BL}, b::DenseOperator{BL,BR}) where {BL<:Basis,BR<:Basis} = Bra{BR}(b.basis_r, transpose(b.data)*a.data) +*(a::Bra, b::DenseOperator) = throw(bases.IncompatibleBases()) +*(a::DenseOperator{B1,B2,T}, b::DenseOperator{B2,B3,T}) where {B1<:Basis,B2<:Basis,B3<:Basis,T<:Matrix{ComplexF64}} = DenseOperator{B1,B3,T}(a.basis_l, b.basis_r, a.data*b.data) +*(a::DenseOperator, b::DenseOperator) = throw(bases.IncompatibleBases()) *(a::DenseOperator, b::Number) = DenseOperator(a.basis_l, a.basis_r, complex(b)*a.data) *(a::Number, b::DenseOperator) = DenseOperator(b.basis_l, b.basis_r, complex(a)*b.data) -function *(op1::Operator, op2::DenseOperator) - check_multiplicable(op1, op2) - result = DenseOperator(op1.basis_l, op2.basis_r) +function *(op1::AbstractOperator{B1,B2}, op2::DenseOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} + result = DenseOperator{B1,B3}(op1.basis_l, op2.basis_r) operators.gemm!(Complex(1.), op1, op2, Complex(0.), result) return result end -function *(op1::DenseOperator, op2::Operator) - check_multiplicable(op1, op2) - result = DenseOperator(op1.basis_l, op2.basis_r) +function *(op1::DenseOperator{B1,B2}, op2::AbstractOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} + result = DenseOperator{B1,B3}(op1.basis_l, op2.basis_r) operators.gemm!(Complex(1.), op1, op2, Complex(0.), result) return result end -function *(op::Operator, psi::Ket) - check_multiplicable(op, psi) - result = Ket(op.basis_l) +function *(op::AbstractOperator{BL,BR}, psi::Ket{BR}) where {BL<:Basis,BR<:Basis} + result = Ket{BL}(op.basis_l) operators.gemv!(Complex(1.), op, psi, Complex(0.), result) return result end -function *(psi::Bra, op::Operator) - check_multiplicable(psi, op) - result = Bra(op.basis_r) +function *(psi::Bra{BL}, op::AbstractOperator{BL,BR}) where {BL<:Basis,BR<:Basis} + result = Bra{BR}(op.basis_r) operators.gemv!(Complex(1.), psi, op, Complex(0.), result) return result end @@ -81,7 +92,8 @@ end operators.dagger(x::DenseOperator) = DenseOperator(x.basis_r, x.basis_l, x.data') -operators.ishermitian(A::DenseOperator) = (A.basis_l == A.basis_r) && ishermitian(A.data) +operators.ishermitian(A::DenseOperator) = false +operators.ishermitian(A::DenseOperator{B,B}) where B<:Basis = ishermitian(A.data) operators.tensor(a::DenseOperator, b::DenseOperator) = DenseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data)) @@ -96,7 +108,7 @@ Outer product ``|x⟩⟨y|`` of the given states. operators.tensor(a::Ket, b::Bra) = DenseOperator(a.basis, b.basis, reshape(kron(b.data, a.data), prod(a.basis.shape), prod(b.basis.shape))) -operators.tr(op::DenseOperator) = (check_samebases(op); tr(op.data)) +operators.tr(op::DenseOperator{B,B}) where B<:Basis = tr(op.data) function operators.ptrace(a::DenseOperator, indices::Vector{Int}) operators.check_ptrace_arguments(a, indices) @@ -124,15 +136,11 @@ end operators.normalize!(op::DenseOperator) = (rmul!(op.data, 1.0/tr(op)); nothing) -function operators.expect(op::DenseOperator, state::Ket)# where T <: Union{Ket, Bra} - check_samebases(op.basis_r, state.basis) - check_samebases(op.basis_l, state.basis) +function operators.expect(op::DenseOperator{B,B}, state::Ket{B}) where B<:Basis state.data' * op.data * state.data end -function operators.expect(op::DenseOperator, state::Operator) - check_samebases(op.basis_r, state.basis_l) - check_samebases(op.basis_l, state.basis_r) +function operators.expect(op::DenseOperator{B1,B2}, state::AbstractOperator{B2,B2}) where {B1<:Basis,B2<:Basis} result = ComplexF64(0.) @inbounds for i=1:size(op.data, 1), j=1:size(op.data,2) result += op.data[i,j]*state.data[j,i] @@ -140,12 +148,11 @@ function operators.expect(op::DenseOperator, state::Operator) result end -function operators.exp(op::DenseOperator) - check_samebases(op) - return DenseOperator(op.basis_l, op.basis_r, exp(op.data)) +function operators.exp(op::T) where {B<:Basis,T<:DenseOperator{B,B}} + return T(op.basis_l, op.basis_r, exp(op.data)) end -function operators.permutesystems(a::DenseOperator, perm::Vector{Int}) +function operators.permutesystems(a::DenseOperator{B1,B2}, perm::Vector{Int}) where {B1<:CompositeBasis,B2<:CompositeBasis} @assert length(a.basis_l.bases) == length(a.basis_r.bases) == length(perm) @assert isperm(perm) data = reshape(a.data, [a.basis_l.shape; a.basis_r.shape]...) @@ -154,7 +161,7 @@ function operators.permutesystems(a::DenseOperator, perm::Vector{Int}) DenseOperator(permutesystems(a.basis_l, perm), permutesystems(a.basis_r, perm), data) end -operators.identityoperator(::Type{DenseOperator}, b1::Basis, b2::Basis) = DenseOperator(b1, b2, Matrix{ComplexF64}(I, length(b1), length(b2))) +operators.identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:DenseOperator} = DenseOperator(b1, b2, Matrix{ComplexF64}(I, length(b1), length(b2))) """ projector(a::Ket, b::Bra) @@ -267,13 +274,13 @@ operators.gemm!(alpha, a::Matrix{ComplexF64}, b::Matrix{ComplexF64}, beta, resul operators.gemv!(alpha, a::Matrix{ComplexF64}, b::Vector{ComplexF64}, beta, result::Vector{ComplexF64}) = BLAS.gemv!('N', convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result) operators.gemv!(alpha, a::Vector{ComplexF64}, b::Matrix{ComplexF64}, beta, result::Vector{ComplexF64}) = BLAS.gemv!('T', convert(ComplexF64, alpha), b, a, convert(ComplexF64, beta), result) -operators.gemm!(alpha, a::DenseOperator, b::DenseOperator, beta, result::DenseOperator) = operators.gemm!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data) -operators.gemv!(alpha, a::DenseOperator, b::Ket, beta, result::Ket) = operators.gemv!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data) -operators.gemv!(alpha, a::Bra, b::DenseOperator, beta, result::Bra) = operators.gemv!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data) +operators.gemm!(alpha, a::DenseOperator{B1,B2}, b::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = operators.gemm!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data) +operators.gemv!(alpha, a::DenseOperator{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis} = operators.gemv!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data) +operators.gemv!(alpha, a::Bra{B1}, b::DenseOperator{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis} = operators.gemv!(convert(ComplexF64, alpha), a.data, b.data, convert(ComplexF64, beta), result.data) # Multiplication for Operators in terms of their gemv! implementation -function operators.gemm!(alpha, M::Operator, b::DenseOperator, beta, result::DenseOperator) +function operators.gemm!(alpha, M::AbstractOperator{B1,B2}, b::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} for i=1:size(b.data, 2) bket = Ket(b.basis_l, b.data[:,i]) resultket = Ket(M.basis_l, result.data[:,i]) @@ -282,7 +289,7 @@ function operators.gemm!(alpha, M::Operator, b::DenseOperator, beta, result::Den end end -function operators.gemm!(alpha, b::DenseOperator, M::Operator, beta, result::DenseOperator) +function operators.gemm!(alpha, b::DenseOperator{B1,B2}, M::AbstractOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} for i=1:size(b.data, 1) bbra = Bra(b.basis_r, vec(b.data[i,:])) resultbra = Bra(M.basis_r, vec(result.data[i,:])) diff --git a/src/operators_lazyproduct.jl b/src/operators_lazyproduct.jl index 67dbd6dd..fe1bd27d 100644 --- a/src/operators_lazyproduct.jl +++ b/src/operators_lazyproduct.jl @@ -19,39 +19,44 @@ The factors of the product are stored in the `operators` field. Additionally a complex factor is stored in the `factor` field which allows for fast multiplication with numbers. """ -mutable struct LazyProduct <: Operator - basis_l::Basis - basis_r::Basis +mutable struct LazyProduct{BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator}}} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR factor::ComplexF64 - operators::Vector{Operator} + operators::T - function LazyProduct(operators::Vector{Operator}, factor::Number=1) - if length(operators) < 1 - throw(ArgumentError("LazyProduct needs at least one operator.")) - end + function LazyProduct{BL,BR,T}(operators::T, factor::Number=1) where {BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator}}} for i = 2:length(operators) check_multiplicable(operators[i-1], operators[i]) end new(operators[1].basis_l, operators[end].basis_r, factor, operators) end end -LazyProduct(operators::Vector, factor::Number=1) = LazyProduct(convert(Vector{Operator}, operators), factor) -LazyProduct(operators::Operator...) = LazyProduct(Operator[operators...]) +function LazyProduct(operators::T, factor::Number=1) where {T<:Tuple{Vararg{AbstractOperator}}} + BL = typeof(operators[1].basis_l) + BR = typeof(operators[end].basis_r) + LazyProduct{BL,BR,T}(operators, factor) +end +LazyProduct(operators::Vector{T}, factor::Number=1) where T<:AbstractOperator = LazyProduct((operators...,), factor) +LazyProduct(operators::AbstractOperator...) = LazyProduct((operators...,)) +LazyProduct() = throw(ArgumentError("LazyProduct needs at least one operator!")) -Base.copy(x::LazyProduct) = LazyProduct([copy(op) for op in x.operators], x.factor) +Base.copy(x::T) where T<:LazyProduct = T(([copy(op) for op in x.operators]...,), x.factor) operators.dense(op::LazyProduct) = op.factor*prod(dense.(op.operators)) +operators.dense(op::LazyProduct{B1,B2,T}) where {B1<:Basis,B2<:Basis,T<:Tuple{AbstractOperator}} = op.factor*dense(op.operators[1]) SparseArrays.sparse(op::LazyProduct) = op.factor*prod(sparse.(op.operators)) +SparseArrays.sparse(op::LazyProduct{B1,B2,T}) where {B1<:Basis,B2<:Basis,T<:Tuple{AbstractOperator}} = op.factor*sparse(op.operators[1]) -==(x::LazyProduct, y::LazyProduct) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.operators==y.operators && x.factor == y.factor - +==(x::LazyProduct, y::LazyProduct) = false +==(x::T, y::T) where T<:LazyProduct = (x.operators==y.operators && x.factor == y.factor) # Arithmetic operations --(a::LazyProduct) = LazyProduct(a.operators, -a.factor) +-(a::T) where T<:LazyProduct = T(a.operators, -a.factor) -*(a::LazyProduct, b::LazyProduct) = (check_multiplicable(a, b); LazyProduct([a.operators; b.operators], a.factor*b.factor)) -*(a::LazyProduct, b::Number) = LazyProduct(a.operators, a.factor*b) -*(a::Number, b::LazyProduct) = LazyProduct(b.operators, a*b.factor) +*(a::LazyProduct{B1,B2}, b::LazyProduct{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = LazyProduct((a.operators..., b.operators...), a.factor*b.factor) +*(a::T, b::Number) where T<:LazyProduct = T(a.operators, a.factor*b) +*(a::Number, b::T) where T<:LazyProduct = T(b.operators, a*b.factor) /(a::LazyProduct, b::Number) = LazyProduct(a.operators, a.factor/b) @@ -62,13 +67,13 @@ operators.tr(op::LazyProduct) = throw(ArgumentError("Trace of LazyProduct is not operators.ptrace(op::LazyProduct, indices::Vector{Int}) = throw(ArgumentError("Partial trace of LazyProduct is not defined. Try to convert to another operator type first with e.g. dense() or sparse().")) -operators.permutesystems(op::LazyProduct, perm::Vector{Int}) = LazyProduct(Operator[permutesystems(op_i, perm) for op_i in op.operators], op.factor) +operators.permutesystems(op::LazyProduct, perm::Vector{Int}) = LazyProduct(([permutesystems(op_i, perm) for op_i in op.operators]...,), op.factor) operators.identityoperator(::Type{LazyProduct}, b1::Basis, b2::Basis) = LazyProduct(identityoperator(b1, b2)) # Fast in-place multiplication -function operators.gemv!(alpha, a::LazyProduct, b::Ket, beta, result::Ket) +function operators.gemv!(alpha, a::LazyProduct{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis} tmp1 = Ket(a.operators[end].basis_l) operators.gemv!(a.factor, a.operators[end], b, 0, tmp1) for i=length(a.operators)-1:-1:2 @@ -79,7 +84,7 @@ function operators.gemv!(alpha, a::LazyProduct, b::Ket, beta, result::Ket) operators.gemv!(alpha, a.operators[1], tmp1, beta, result) end -function operators.gemv!(alpha, a::Bra, b::LazyProduct, beta, result::Bra) +function operators.gemv!(alpha, a::Bra{B1}, b::LazyProduct{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis} tmp1 = Bra(b.operators[1].basis_r) operators.gemv!(b.factor, a, b.operators[1], 0, tmp1) for i=2:length(b.operators)-1 diff --git a/src/operators_lazysum.jl b/src/operators_lazysum.jl index 94f4c263..3d059a4c 100644 --- a/src/operators_lazysum.jl +++ b/src/operators_lazysum.jl @@ -18,37 +18,44 @@ All operators have to be given in respect to the same bases. The field `factors` accounts for an additional multiplicative factor for each operator stored in the `operators` field. """ -mutable struct LazySum <: Operator - basis_l::Basis - basis_r::Basis +mutable struct LazySum{BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator{BL,BR}}}} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR factors::Vector{ComplexF64} - operators::Vector{Operator} - - function LazySum(factors::Vector{ComplexF64}, operators::Vector{Operator}) - @assert length(operators)>0 + operators::T + function LazySum{BL,BR,T}(factors::Vector{ComplexF64}, operators::T) where {BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator{BL,BR}}}} @assert length(operators)==length(factors) - for i = 2:length(operators) - @assert operators[1].basis_l == operators[i].basis_l - @assert operators[1].basis_r == operators[i].basis_r - end new(operators[1].basis_l, operators[1].basis_r, factors, operators) end end -LazySum(factors::Vector{T}, operators::Vector) where {T<:Number} = LazySum(complex(factors), Operator[op for op in operators]) -LazySum(operators::Operator...) = LazySum(ones(ComplexF64, length(operators)), Operator[operators...]) +function LazySum(factors::Vector{ComplexF64}, operators::T) where {BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator{BL,BR}}}} + for i = 2:length(operators) + operators[1].basis_l == operators[i].basis_l || throw(bases.IncompatibleBases()) + operators[1].basis_r == operators[i].basis_r || throw(bases.IncompatibleBases()) + end + LazySum{BL,BR,T}(factors, operators) +end +LazySum(operators::AbstractOperator...) = LazySum(ones(ComplexF64, length(operators)), (operators...,)) +LazySum(factors::Vector{T}, operators::Vector{T2}) where {T<:Number,B1<:Basis,B2<:Basis,T2<:AbstractOperator{B1,B2}} = LazySum(complex(factors), (operators...,)) +LazySum() = throw(ArgumentError("LazySum needs at least one operator!")) -Base.copy(x::LazySum) = LazySum(copy(x.factors), [copy(op) for op in x.operators]) +Base.copy(x::T) where T<:LazySum = T(copy(x.factors), ([copy(op) for op in x.operators]...,)) operators.dense(op::LazySum) = sum(op.factors .* dense.(op.operators)) +operators.dense(op::LazySum{B1,B2,T}) where {B1<:Basis,B2<:Basis,T<:Tuple{AbstractOperator{B1,B2}}} = op.factors[1] * dense(op.operators[1]) SparseArrays.sparse(op::LazySum) = sum(op.factors .* sparse.(op.operators)) +SparseArrays.sparse(op::LazySum{B1,B2,T}) where {B1<:Basis,B2<:Basis,T<:Tuple{AbstractOperator{B1,B2}}} = op.factors[1] * sparse(op.operators[1]) -==(x::LazySum, y::LazySum) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.operators==y.operators && x.factors==y.factors +==(x::T, y::T) where T<:LazySum = (x.operators==y.operators && x.factors==y.factors) +==(x::LazySum, y::LazySum) = false # Arithmetic operations -+(a::LazySum, b::LazySum) = (check_samebases(a,b); LazySum([a.factors; b.factors], [a.operators; b.operators])) ++(a::LazySum{B1,B2,T1}, b::LazySum{B1,B2,T2}) where {B1<:Basis,B2<:Basis,T1<:Tuple{Vararg{AbstractOperator{B1,B2}}},T2<:Tuple{Vararg{AbstractOperator{B1,B2}}}} = LazySum([a.factors; b.factors], (a.operators..., b.operators...)) ++(a::LazySum{B1,B2}, b::LazySum{B3,B4}) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(bases.IncompatibleBases()) --(a::LazySum) = LazySum(-a.factors, a.operators) --(a::LazySum, b::LazySum) = (check_samebases(a,b); LazySum([a.factors; -b.factors], [a.operators; b.operators])) +-(a::T) where T<:LazySum = T(-a.factors, a.operators) +-(a::LazySum{B1,B2,T1}, b::LazySum{B1,B2,T2}) where {B1<:Basis,B2<:Basis,T1<:Tuple{Vararg{AbstractOperator{B1,B2}}},T2<:Tuple{Vararg{AbstractOperator{B1,B2}}}} = LazySum([a.factors; -b.factors], (a.operators..., b.operators...)) +-(a::LazySum{B1,B2}, b::LazySum{B3,B4}) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(bases.IncompatibleBases()) *(a::LazySum, b::Number) = LazySum(b*a.factors, a.operators) *(a::Number, b::LazySum) = LazySum(a*b.factors, b.operators) @@ -62,26 +69,26 @@ operators.tr(op::LazySum) = sum(op.factors .* tr.(op.operators)) function operators.ptrace(op::LazySum, indices::Vector{Int}) operators.check_ptrace_arguments(op, indices) rank = length(op.basis_l.shape) - length(indices) - D = Operator[ptrace(op_i, indices) for op_i in op.operators] + D = ([ptrace(op_i, indices) for op_i in op.operators]...,) LazySum(op.factors, D) end operators.normalize!(op::LazySum) = (op.factors /= tr(op); nothing) -operators.permutesystems(op::LazySum, perm::Vector{Int}) = LazySum(op.factors, Operator[permutesystems(op_i, perm) for op_i in op.operators]) +operators.permutesystems(op::LazySum, perm::Vector{Int}) = LazySum(op.factors, ([permutesystems(op_i, perm) for op_i in op.operators]...,)) operators.identityoperator(::Type{LazySum}, b1::Basis, b2::Basis) = LazySum(identityoperator(b1, b2)) # Fast in-place multiplication -function operators.gemv!(alpha, a::LazySum, b::Ket, beta, result::Ket) +function operators.gemv!(alpha, a::LazySum{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis,B3<:Basis} operators.gemv!(alpha*a.factors[1], a.operators[1], b, beta, result) for i=2:length(a.operators) operators.gemv!(alpha*a.factors[i], a.operators[i], b, 1, result) end end -function operators.gemv!(alpha, a::Bra, b::LazySum, beta, result::Bra) +function operators.gemv!(alpha, a::Bra{B1}, b::LazySum{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis,B3<:Basis} operators.gemv!(alpha*b.factors[1], a, b.operators[1], beta, result) for i=2:length(b.operators) operators.gemv!(alpha*b.factors[i], a, b.operators[i], 1, result) diff --git a/src/operators_lazytensor.jl b/src/operators_lazytensor.jl index 081d4f4e..16716cc1 100644 --- a/src/operators_lazytensor.jl +++ b/src/operators_lazytensor.jl @@ -20,24 +20,24 @@ specifies in which subsystem the corresponding operator lives. Additionally, a complex factor is stored in the `factor` field which allows for fast multiplication with numbers. """ -mutable struct LazyTensor <: Operator - basis_l::CompositeBasis - basis_r::CompositeBasis +mutable struct LazyTensor{BL<:CompositeBasis,BR<:CompositeBasis,T<:Tuple{Vararg{AbstractOperator}}} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR factor::ComplexF64 indices::Vector{Int} - operators::Vector{Operator} + operators::T - function LazyTensor(op::LazyTensor, factor::Number) + function LazyTensor{BL,BR,T}(op::LazyTensor{BL,BR,T}, factor::Number) where {BL<:CompositeBasis,BR<:CompositeBasis,T<:Tuple{Vararg{AbstractOperator}}} new(op.basis_l, op.basis_r, factor, op.indices, op.operators) end function LazyTensor(basis_l::Basis, basis_r::Basis, indices::Vector{Int}, ops::Vector, factor::Number=1) - if typeof(basis_l) != CompositeBasis + if !isa(basis_l, CompositeBasis) basis_l = CompositeBasis(basis_l.shape, Basis[basis_l]) end - if typeof(basis_r) != CompositeBasis + if !isa(basis_r, CompositeBasis) basis_r = CompositeBasis(basis_r.shape, Basis[basis_r]) end N = length(basis_l.bases) @@ -45,7 +45,7 @@ mutable struct LazyTensor <: Operator sortedindices.check_indices(N, indices) @assert length(indices) == length(ops) for n=1:length(indices) - @assert isa(ops[n], Operator) + @assert isa(ops[n], AbstractOperator) @assert ops[n].basis_l == basis_l.bases[indices[n]] @assert ops[n].basis_r == basis_r.bases[indices[n]] end @@ -54,13 +54,14 @@ mutable struct LazyTensor <: Operator indices = indices[perm] ops = ops[perm] end - new(basis_l, basis_r, complex(factor), indices, ops) + ops_tup = (ops...,) + new{typeof(basis_l),typeof(basis_r),typeof(ops_tup)}(basis_l, basis_r, complex(factor), indices, ops_tup) end end - +LazyTensor(op::LazyTensor{BL,BR,T}, factor::Number) where {BL<:CompositeBasis,BR<:CompositeBasis,T<:Tuple{Vararg{AbstractOperator}}} = LazyTensor{BL,BR,T}(op::LazyTensor{BL,BR,T}, factor::Number) LazyTensor(basis::Basis, indices::Vector{Int}, ops::Vector, factor::Number=1) = LazyTensor(basis, basis, indices, ops, factor) -LazyTensor(basis_l::Basis, basis_r::Basis, index::Int, operator::Operator, factor::Number=1) = LazyTensor(basis_l, basis_r, [index], Operator[operator], factor) -LazyTensor(basis::Basis, index::Int, operators::Operator, factor::Number=1.) = LazyTensor(basis, basis, index, operators, factor) +LazyTensor(basis_l::Basis, basis_r::Basis, index::Int, operator::AbstractOperator, factor::Number=1) = LazyTensor(basis_l, basis_r, [index], AbstractOperator[operator], factor) +LazyTensor(basis::Basis, index::Int, operators::AbstractOperator, factor::Number=1.) = LazyTensor(basis, basis, index, operators, factor) Base.copy(x::LazyTensor) = LazyTensor(x.basis_l, x.basis_r, copy(x.indices), [copy(op) for op in x.operators], x.factor) @@ -77,7 +78,7 @@ suboperator(op::LazyTensor, index::Int) = op.operators[findfirst(isequal(index), Return the suboperators corresponding to the subsystems specified by `indices`. Fails if there is no corresponding operator (i.e. it would be an identity operater). """ -suboperators(op::LazyTensor, indices::Vector{Int}) = op.operators[[findfirst(isequal(i), op.indices) for i in indices]] +suboperators(op::LazyTensor, indices::Vector{Int}) = [op.operators[[findfirst(isequal(i), op.indices) for i in indices]]...] operators.dense(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op.indices, DenseOperator[dense(x) for x in op.operators]) SparseArrays.sparse(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op.indices, SparseOperator[sparse(x) for x in op.operators]) @@ -88,10 +89,9 @@ SparseArrays.sparse(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op # Arithmetic operations -(a::LazyTensor) = LazyTensor(a, -a.factor) -function *(a::LazyTensor, b::LazyTensor) - check_multiplicable(a, b) +function *(a::LazyTensor{B1,B2}, b::LazyTensor{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} indices = sortedindices.union(a.indices, b.indices) - ops = Vector{Operator}(undef, length(indices)) + ops = Vector{AbstractOperator}(undef, length(indices)) for n in 1:length(indices) i = indices[n] in_a = i in a.indices @@ -110,15 +110,13 @@ function *(a::LazyTensor, b::LazyTensor) end *(a::LazyTensor, b::Number) = LazyTensor(a, a.factor*b) *(a::Number, b::LazyTensor) = LazyTensor(b, a*b.factor) -function *(a::LazyTensor, b::DenseOperator) - check_multiplicable(a, b) - result = DenseOperator(a.basis_l, b.basis_r) +function *(a::LazyTensor{B1,B2}, b::DenseOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} + result = DenseOperator{B1,B3}(a.basis_l, b.basis_r) operators.gemm!(complex(1.), a, b, complex(1.), result) result end -function *(a::DenseOperator, b::LazyTensor) - check_multiplicable(a, b) - result = DenseOperator(a.basis_l, b.basis_r) +function *(a::DenseOperator{B1,B2}, b::LazyTensor{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} + result = DenseOperator{B1,B3}(a.basis_l, b.basis_r) operators.gemm!(complex(1.), a, b, complex(1.), result) result end @@ -126,9 +124,9 @@ end /(a::LazyTensor, b::Number) = LazyTensor(a, a.factor/b) -operators.dagger(op::LazyTensor) = LazyTensor(op.basis_r, op.basis_l, op.indices, Operator[dagger(x) for x in op.operators], conj(op.factor)) +operators.dagger(op::LazyTensor) = LazyTensor(op.basis_r, op.basis_l, op.indices, AbstractOperator[dagger(x) for x in op.operators], conj(op.factor)) -operators.tensor(a::LazyTensor, b::LazyTensor) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [a.indices; b.indices .+ length(a.basis_l.bases)], Operator[a.operators; b.operators], a.factor*b.factor) +operators.tensor(a::LazyTensor, b::LazyTensor) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [a.indices; b.indices .+ length(a.basis_l.bases)], [a.operators..., b.operators...], a.factor*b.factor) function operators.tr(op::LazyTensor) b = basis(op) @@ -164,7 +162,7 @@ function operators.ptrace(op::LazyTensor, indices::Vector{Int}) if rank==1 return factor * identityoperator(b_l, b_r) end - ops = Vector{Operator}(undef, length(remaining_indices)) + ops = Vector{AbstractOperator}(undef, length(remaining_indices)) for i in 1:length(ops) ops[i] = suboperator(op, remaining_indices[i]) end @@ -178,10 +176,10 @@ function operators.permutesystems(op::LazyTensor, perm::Vector{Int}) b_r = permutesystems(op.basis_r, perm) indices = [findfirst(isequal(i), perm) for i in op.indices] perm_ = sortperm(indices) - LazyTensor(b_l, b_r, indices[perm_], op.operators[perm_], op.factor) + LazyTensor(b_l, b_r, indices[perm_], [op.operators[perm_]...], op.factor) end -operators.identityoperator(::Type{LazyTensor}, b1::Basis, b2::Basis) = LazyTensor(b1, b2, Int[], Operator[]) +operators.identityoperator(::Type{LazyTensor}, b1::Basis, b2::Basis) = LazyTensor(b1, b2, Int[], AbstractOperator[]) # Recursively calculate result_{IK} = \\sum_J op_{IJ} h_{JK} @@ -307,22 +305,22 @@ function gemm(alpha::ComplexF64, h::LazyTensor, op::Matrix{ComplexF64}, beta::Co _gemm_recursive_lazy_dense(1, N_k, 1, 1, alpha*h.factor, shape, strides_k, strides_j, h.indices, h, op, result) end -operators.gemm!(alpha, h::LazyTensor, op::DenseOperator, beta, result::DenseOperator) = gemm(convert(ComplexF64, alpha), h, op.data, convert(ComplexF64, beta), result.data) -operators.gemm!(alpha, op::DenseOperator, h::LazyTensor, beta, result::DenseOperator) = gemm(convert(ComplexF64, alpha), op.data, h, convert(ComplexF64, beta), result.data) +operators.gemm!(alpha, h::LazyTensor{B1,B2}, op::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = gemm(convert(ComplexF64, alpha), h, op.data, convert(ComplexF64, beta), result.data) +operators.gemm!(alpha, op::DenseOperator{B1,B2}, h::LazyTensor{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = gemm(convert(ComplexF64, alpha), op.data, h, convert(ComplexF64, beta), result.data) -function operators.gemv!(alpha::ComplexF64, a::LazyTensor, b::Ket, beta::ComplexF64, result::Ket) +function operators.gemv!(alpha::ComplexF64, a::LazyTensor{B1,B2}, b::Ket{B2}, beta::ComplexF64, result::Ket{B1}) where {B1<:Basis,B2<:Basis} b_data = reshape(b.data, length(b.data), 1) result_data = reshape(result.data, length(result.data), 1) gemm(alpha, a, b_data, beta, result_data) end -function operators.gemv!(alpha::ComplexF64, a::Bra, b::LazyTensor, beta::ComplexF64, result::Bra) +function operators.gemv!(alpha::ComplexF64, a::Bra{B1}, b::LazyTensor{B1,B2}, beta::ComplexF64, result::Bra{B2}) where {B1<:Basis,B2<:Basis} a_data = reshape(a.data, 1, length(a.data)) result_data = reshape(result.data, 1, length(result.data)) gemm(alpha, a_data, b, beta, result_data) end -operators.gemv!(alpha, a::LazyTensor, b::Ket, beta, result::Ket) = operators.gemv!(convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result) -operators.gemv!(alpha, a::Bra, b::LazyTensor, beta, result::Bra) = operators.gemv!(convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result) +operators.gemv!(alpha, a::LazyTensor{B1,B2}, b::Ket{B2}, beta, result::Ket{B2}) where {B1<:Basis,B2<:Basis} = operators.gemv!(convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result) +operators.gemv!(alpha, a::Bra{B1}, b::LazyTensor{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis} = operators.gemv!(convert(ComplexF64, alpha), a, b, convert(ComplexF64, beta), result) end # module diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl index 7bd9fe36..4bc4aa2b 100644 --- a/src/operators_sparse.jl +++ b/src/operators_sparse.jl @@ -18,11 +18,11 @@ Sparse array implementation of Operator. The matrix is stored as the julia built-in type `SparseMatrixCSC` in the `data` field. """ -mutable struct SparseOperator <: Operator - basis_l::Basis - basis_r::Basis - data::SparseMatrixCSC{ComplexF64, Int} - function SparseOperator(b1::Basis, b2::Basis, data) +mutable struct SparseOperator{BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR + data::T + function SparseOperator{BL,BR,T}(b1::Basis, b2::Basis, data::T) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} if length(b1) != size(data, 1) || length(b2) != size(data, 2) throw(DimensionMismatch()) end @@ -30,54 +30,68 @@ mutable struct SparseOperator <: Operator end end +SparseOperator(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} = SparseOperator{BL,BR,T}(b1, b2, data) +SparseOperator{BL,BR}(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} = SparseOperator{BL,BR,T}(b1, b2, data) +SparseOperator(b1::Basis, b2::Basis, data) = SparseOperator(b1, b2, convert(SparseMatrixCSC{ComplexF64,Int}, data)) SparseOperator(b::Basis, data::SparseMatrixCSC{ComplexF64, Int}) = SparseOperator(b, b, data) SparseOperator(b::Basis, data::Matrix{ComplexF64}) = SparseOperator(b, sparse(data)) SparseOperator(op::DenseOperator) = SparseOperator(op.basis_l, op.basis_r, sparse(op.data)) SparseOperator(b1::Basis, b2::Basis) = SparseOperator(b1, b2, spzeros(ComplexF64, length(b1), length(b2))) +SparseOperator{BL,BR}(b1::BL, b2::BR) where {BL<:Basis,BR<:Basis} = SparseOperator{BL,BR}(b1, b2, spzeros(ComplexF64, length(b1), length(b2))) SparseOperator(b::Basis) = SparseOperator(b, b) Base.copy(x::SparseOperator) = SparseOperator(x.basis_l, x.basis_r, copy(x.data)) operators.dense(a::SparseOperator) = DenseOperator(a.basis_l, a.basis_r, Matrix(a.data)) """ - sparse(op::Operator) + sparse(op::AbstractOperator) Convert an arbitrary operator into a [`SparseOperator`](@ref). """ -sparse(a::Operator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead.")) +sparse(a::AbstractOperator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead.")) sparse(a::SparseOperator) = copy(a) sparse(a::DenseOperator) = SparseOperator(a.basis_l, a.basis_r, sparse(a.data)) -==(x::SparseOperator, y::SparseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && (x.data == y.data) +==(x::SparseOperator, y::SparseOperator) = false +==(x::T, y::T) where T<:SparseOperator = (x.data == y.data) # Arithmetic operations -+(a::SparseOperator, b::SparseOperator) = (check_samebases(a,b); SparseOperator(a.basis_l, a.basis_r, a.data+b.data)) -+(a::SparseOperator, b::DenseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data+b.data)) -+(a::DenseOperator, b::SparseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data+b.data)) - --(a::SparseOperator) = SparseOperator(a.basis_l, a.basis_r, -a.data) --(a::SparseOperator, b::SparseOperator) = (check_samebases(a,b); SparseOperator(a.basis_l, a.basis_r, a.data-b.data)) --(a::SparseOperator, b::DenseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data-b.data)) --(a::DenseOperator, b::SparseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data-b.data)) - -*(a::SparseOperator, b::SparseOperator) = (check_multiplicable(a, b); SparseOperator(a.basis_l, b.basis_r, a.data*b.data)) -*(a::SparseOperator, b::Number) = SparseOperator(a.basis_l, a.basis_r, complex(b)*a.data) -*(a::Number, b::SparseOperator) = SparseOperator(b.basis_l, b.basis_r, complex(a)*b.data) - -/(a::SparseOperator, b::Number) = SparseOperator(a.basis_l, a.basis_r, a.data/complex(b)) ++(a::SparseOperator, b::SparseOperator) = throw(bases.IncompatibleBases()) ++(a::T, b::T) where T<:SparseOperator = T(a.basis_l, a.basis_r, a.data+b.data) ++(a::SparseOperator, b::DenseOperator) = throw(bases.IncompatibleBases()) ++(a::SparseOperator{B1,B2}, b::DenseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data+b.data) ++(a::DenseOperator, b::SparseOperator) = throw(bases.IncompatibleBases()) ++(a::DenseOperator{B1,B2}, b::SparseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data+b.data) + +-(a::T) where T<:SparseOperator = T(a.basis_l, a.basis_r, -a.data) +-(a::SparseOperator, b::SparseOperator) = throw(bases.IncompatibleBases()) +-(a::T, b::T) where T<:SparseOperator = SparseOperator(a.basis_l, a.basis_r, a.data-b.data) +-(a::SparseOperator, b::DenseOperator) = throw(bases.IncompatibleBases()) +-(a::SparseOperator{B1,B2}, b::DenseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data-b.data) +-(a::DenseOperator, b::SparseOperator) = throw(bases.IncompatibleBases()) +-(a::DenseOperator{B1,B2}, b::SparseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data-b.data) + +*(a::SparseOperator{B1,B2}, b::SparseOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = SparseOperator{B1,B3}(a.basis_l, b.basis_r, a.data*b.data) +*(a::SparseOperator, b::SparseOperator) = throw(bases.IncompatibleBases()) +*(a::T, b::Number) where T<:SparseOperator = T(a.basis_l, a.basis_r, complex(b)*a.data) +*(a::Number, b::T) where T<:SparseOperator = T(b.basis_l, b.basis_r, complex(a)*b.data) + +/(a::T, b::Number) where T<:SparseOperator = T(a.basis_l, a.basis_r, a.data/complex(b)) operators.dagger(x::SparseOperator) = SparseOperator(x.basis_r, x.basis_l, x.data') -operators.ishermitian(A::SparseOperator) = (A.basis_l == A.basis_r) && ishermitian(A.data) +operators.ishermitian(A::SparseOperator) = false +operators.ishermitian(A::SparseOperator{B,B}) where B<:Basis = ishermitian(A.data) operators.tensor(a::SparseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data)) operators.tensor(a::DenseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data)) operators.tensor(a::SparseOperator, b::DenseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data)) -operators.tr(op::SparseOperator) = (check_samebases(op); tr(op.data)) +operators.tr(op::SparseOperator{B,B}) where B<:Basis = tr(op.data) +operators.tr(op::SparseOperator) = throw(bases.IncompatibleBases()) -operators.conj(op::SparseOperator) = SparseOperator(op.basis_l, op.basis_r, conj(op.data)) +operators.conj(op::T) where T<:SparseOperator = T(op.basis_l, op.basis_r, conj(op.data)) operators.conj!(op::SparseOperator) = conj!(op.data) function operators.ptrace(op::SparseOperator, indices::Vector{Int}) @@ -89,16 +103,12 @@ function operators.ptrace(op::SparseOperator, indices::Vector{Int}) SparseOperator(b_l, b_r, data) end -function operators.expect(op::SparseOperator, state::Ket)# where T <: Union{Ket, Bra} - check_samebases(op.basis_r, state.basis) - check_samebases(op.basis_l, state.basis) +function operators.expect(op::SparseOperator{B,B}, state::Ket{B}) where B<:Basis state.data' * op.data * state.data end -function operators.expect(op::SparseOperator, state::DenseOperator) - check_samebases(op.basis_r, state.basis_l) - check_samebases(op.basis_l, state.basis_r) +function operators.expect(op::SparseOperator{B1,B2}, state::DenseOperator{B2,B2}) where {B1<:Basis,B2<:Basis} result = ComplexF64(0.) @inbounds for colindex = 1:op.data.n for i=op.data.colptr[colindex]:op.data.colptr[colindex+1]-1 @@ -108,7 +118,7 @@ function operators.expect(op::SparseOperator, state::DenseOperator) result end -function operators.permutesystems(rho::SparseOperator, perm::Vector{Int}) +function operators.permutesystems(rho::SparseOperator{B1,B2}, perm::Vector{Int}) where {B1<:CompositeBasis,B2<:CompositeBasis} @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm) @assert isperm(perm) shape = [rho.basis_l.shape; rho.basis_r.shape] @@ -116,7 +126,7 @@ function operators.permutesystems(rho::SparseOperator, perm::Vector{Int}) SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data) end -operators.identityoperator(::Type{SparseOperator}, b1::Basis, b2::Basis) = SparseOperator(b1, b2, sparse(ComplexF64(1)*I, length(b1), length(b2))) +operators.identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:SparseOperator} = SparseOperator(b1, b2, sparse(ComplexF64(1)*I, length(b1), length(b2))) operators.identityoperator(b1::Basis, b2::Basis) = identityoperator(SparseOperator, b1, b2) operators.identityoperator(b::Basis) = identityoperator(b, b) @@ -132,9 +142,9 @@ end # Fast in-place multiplication implementations -operators.gemm!(alpha, M::SparseOperator, b::DenseOperator, beta, result::DenseOperator) = sparsematrix.gemm!(convert(ComplexF64, alpha), M.data, b.data, convert(ComplexF64, beta), result.data) -operators.gemm!(alpha, a::DenseOperator, M::SparseOperator, beta, result::DenseOperator) = sparsematrix.gemm!(convert(ComplexF64, alpha), a.data, M.data, convert(ComplexF64, beta), result.data) -operators.gemv!(alpha, M::SparseOperator, b::Ket, beta, result::Ket) = sparsematrix.gemv!(convert(ComplexF64, alpha), M.data, b.data, convert(ComplexF64, beta), result.data) -operators.gemv!(alpha, b::Bra, M::SparseOperator, beta, result::Bra) = sparsematrix.gemv!(convert(ComplexF64, alpha), b.data, M.data, convert(ComplexF64, beta), result.data) +operators.gemm!(alpha, M::SparseOperator{B1,B2}, b::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = sparsematrix.gemm!(convert(ComplexF64, alpha), M.data, b.data, convert(ComplexF64, beta), result.data) +operators.gemm!(alpha, a::DenseOperator{B1,B2}, M::SparseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = sparsematrix.gemm!(convert(ComplexF64, alpha), a.data, M.data, convert(ComplexF64, beta), result.data) +operators.gemv!(alpha, M::SparseOperator{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis} = sparsematrix.gemv!(convert(ComplexF64, alpha), M.data, b.data, convert(ComplexF64, beta), result.data) +operators.gemv!(alpha, b::Bra{B1}, M::SparseOperator{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis} = sparsematrix.gemv!(convert(ComplexF64, alpha), b.data, M.data, convert(ComplexF64, beta), result.data) end # module diff --git a/src/particle.jl b/src/particle.jl index ceb1c49c..eb731b3e 100644 --- a/src/particle.jl +++ b/src/particle.jl @@ -27,13 +27,17 @@ of ``x_{min}`` and ``x_{max}`` are due to the periodic boundary conditions more or less arbitrary and are chosen to be ``-\\pi/dp`` and ``\\pi/dp`` with ``dp=(p_{max}-p_{min})/N``. """ -mutable struct PositionBasis <: Basis +mutable struct PositionBasis{X1,X2} <: Basis shape::Vector{Int} xmin::Float64 xmax::Float64 N::Int - PositionBasis(xmin::Real, xmax::Real, N::Int) = new([N], xmin, xmax, N) + function PositionBasis{X1,X2}(xmin::Real, xmax::Real, N::Int) where {X1,X2} + @assert isa(X1, Real) && isa(X2, Real) + new([N], xmin, xmax, N) + end end +PositionBasis(xmin::Real, xmax::Real, N::Int) = PositionBasis{xmin,xmax}(xmin,xmax,N) """ MomentumBasis(pmin, pmax, Npoints) @@ -49,13 +53,17 @@ of ``p_{min}`` and ``p_{max}`` are due to the periodic boundary conditions more or less arbitrary and are chosen to be ``-\\pi/dx`` and ``\\pi/dx`` with ``dx=(x_{max}-x_{min})/N``. """ -mutable struct MomentumBasis <: Basis +mutable struct MomentumBasis{P1,P2} <: Basis shape::Vector{Int} pmin::Float64 pmax::Float64 N::Int - MomentumBasis(pmin::Real, pmax::Real, N::Int) = new([N], pmin, pmax, N) + function MomentumBasis{P1,P2}(pmin::Real, pmax::Real, N::Int) where {P1,P2} + @assert isa(P1, Real) && isa(P2, Real) + new([N], pmin, pmax, N) + end end +MomentumBasis(pmin::Real, pmax::Real, N::Int) = MomentumBasis{pmin,pmax}(pmin, pmax, N) PositionBasis(b::MomentumBasis) = (dp = (b.pmax - b.pmin)/b.N; PositionBasis(-pi/dp, pi/dp, b.N)) MomentumBasis(b::PositionBasis) = (dx = (b.xmax - b.xmin)/b.N; MomentumBasis(-pi/dx, pi/dx, b.N)) @@ -264,7 +272,7 @@ end Abstract type for all implementations of FFT operators. """ -abstract type FFTOperator <: Operator end +abstract type FFTOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} end const PlanFFT = FFTW.cFFTWPlan @@ -274,15 +282,24 @@ const PlanFFT = FFTW.cFFTWPlan Operator performing a fast fourier transformation when multiplied with a state that is a Ket or an Operator. """ -mutable struct FFTOperators <: FFTOperator - basis_l::Basis - basis_r::Basis +mutable struct FFTOperators{BL<:Basis,BR<:Basis} <: FFTOperator{BL,BR} + basis_l::BL + basis_r::BR fft_l!::PlanFFT fft_r!::PlanFFT fft_l2!::PlanFFT fft_r2!::PlanFFT mul_before::Array{ComplexF64} mul_after::Array{ComplexF64} + function FFTOperators(b1::BL, b2::BR, + fft_l!::PlanFFT, + fft_r!::PlanFFT, + fft_l2!::PlanFFT, + fft_r2!::PlanFFT, + mul_before::Array{ComplexF64}, + mul_after::Array{ComplexF64}) where {BL<:Basis,BR<:Basis} + new{BL,BR}(b1, b2, fft_l!, fft_r!, fft_l2!, fft_r2!, mul_before, mul_after) + end end """ @@ -291,13 +308,27 @@ end Operator that can only perform fast fourier transformations on Kets. This is much more memory efficient when only working with Kets. """ -mutable struct FFTKets <: FFTOperator - basis_l::Basis - basis_r::Basis +mutable struct FFTKets{BL<:Basis,BR<:Basis} <: FFTOperator{BL,BR} + basis_l::BL + basis_r::BR fft_l!::PlanFFT fft_r!::PlanFFT mul_before::Array{ComplexF64} mul_after::Array{ComplexF64} + function FFTKets{BL,BR}(b1::BL, b2::BR, + fft_l!::PlanFFT, + fft_r!::PlanFFT, + mul_before::Array{ComplexF64}, + mul_after::Array{ComplexF64}) where {BL<:Basis,BR<:Basis} + new(b1, b2, fft_l!, fft_r!, mul_before, mul_after) + end +end +function FFTKets(b1::BL, b2::BR, + fft_l!::PlanFFT, + fft_r!::PlanFFT, + mul_before::Array{ComplexF64}, + mul_after::Array{ComplexF64}) where {BL<:Basis,BR<:Basis} + FFTKets{BL,BR}(b1, b2, fft_l!, fft_r!, mul_before, mul_after) end """ @@ -352,8 +383,8 @@ end function transform(basis_l::CompositeBasis, basis_r::CompositeBasis; ket_only::Bool=false, index::Vector{Int}=Int[]) @assert length(basis_l.bases) == length(basis_r.bases) if length(index) == 0 - check_pos = typeof.(basis_l.bases) .== PositionBasis - check_mom = typeof.(basis_l.bases) .== MomentumBasis + check_pos = [isa.(basis_l.bases, PositionBasis)...] + check_mom = [isa.(basis_l.bases, MomentumBasis)...] if any(check_pos) && !any(check_mom) index = [1:length(basis_l.bases);][check_pos] elseif any(check_mom) && !any(check_pos) @@ -362,11 +393,11 @@ function transform(basis_l::CompositeBasis, basis_r::CompositeBasis; ket_only::B throw(IncompatibleBases()) end end - if all(typeof.(basis_l.bases[index]) .== PositionBasis) - @assert all(typeof.(basis_r.bases[index]) .== MomentumBasis) + if all(isa.(basis_l.bases[index], PositionBasis)) + @assert all(isa.(basis_r.bases[index], MomentumBasis)) transform_xp(basis_l, basis_r, index; ket_only=ket_only) - elseif all(typeof.(basis_l.bases[index]) .== MomentumBasis) - @assert all(typeof.(basis_r.bases[index]) .== PositionBasis) + elseif all(isa.(basis_l.bases[index], MomentumBasis)) + @assert all(isa.(basis_r.bases[index], PositionBasis)) transform_px(basis_l, basis_r, index; ket_only=ket_only) else throw(IncompatibleBases()) diff --git a/src/phasespace.jl b/src/phasespace.jl index 2cef8462..09d6ff98 100755 --- a/src/phasespace.jl +++ b/src/phasespace.jl @@ -17,15 +17,13 @@ function can either be evaluated on one point α or on a grid specified by the vectors `xvec` and `yvec`. Note that conversion from `x` and `y` to `α` is done via the relation ``α = \\frac{1}{\\sqrt{2}}(x + i y)``. """ -function qfunc(rho::Operator, alpha::Number) +function qfunc(rho::AbstractOperator{B,B}, alpha::Number) where B<:FockBasis b = basis(rho) - @assert isa(b, FockBasis) _qfunc_operator(rho, convert(ComplexF64, alpha), Ket(b), Ket(b)) end -function qfunc(rho::Operator, xvec::Vector{Float64}, yvec::Vector{Float64}) +function qfunc(rho::AbstractOperator{B,B}, xvec::Vector{Float64}, yvec::Vector{Float64}) where B<:FockBasis b = basis(rho) - @assert isa(b, FockBasis) Nx = length(xvec) Ny = length(yvec) tmp1 = Ket(b) @@ -37,9 +35,8 @@ function qfunc(rho::Operator, xvec::Vector{Float64}, yvec::Vector{Float64}) result end -function qfunc(psi::Ket, alpha::Number) +function qfunc(psi::Ket{B}, alpha::Number) where B<:FockBasis b = basis(psi) - @assert isa(b, FockBasis) _alpha = convert(ComplexF64, alpha) _conj_alpha = conj(_alpha) N = length(psi.basis) @@ -51,9 +48,8 @@ function qfunc(psi::Ket, alpha::Number) return abs2(s)*exp(-abs2(_alpha))/pi end -function qfunc(psi::Ket, xvec::Vector{Float64}, yvec::Vector{Float64}) +function qfunc(psi::Ket{B}, xvec::Vector{Float64}, yvec::Vector{Float64}) where B<:FockBasis b = basis(psi) - @assert isa(b, FockBasis) Nx = length(xvec) Ny = length(yvec) points = Nx*Ny @@ -75,11 +71,11 @@ function qfunc(psi::Ket, xvec::Vector{Float64}, yvec::Vector{Float64}) result end -function qfunc(state::Union{Ket, Operator}, x::Number, y::Number) +function qfunc(state::Union{Ket{B}, AbstractOperator{B,B}}, x::Number, y::Number) where B<:FockBasis qfunc(state, ComplexF64(x, y)/sqrt(2)) end -function _qfunc_operator(rho::Operator, alpha::ComplexF64, tmp1::Ket, tmp2::Ket) +function _qfunc_operator(rho::AbstractOperator{B,B}, alpha::ComplexF64, tmp1::Ket, tmp2::Ket) where B<:FockBasis coherentstate(basis(rho), alpha, tmp1) operators.gemv!(complex(1.), rho, tmp1, complex(0.), tmp2) a = dot(tmp1.data, tmp2.data) @@ -97,9 +93,8 @@ function can either be evaluated on one point α or on a grid specified by the vectors `xvec` and `yvec`. Note that conversion from `x` and `y` to `α` is done via the relation ``α = \\frac{1}{\\sqrt{2}}(x + i y)``. """ -function wigner(rho::DenseOperator, x::Number, y::Number) +function wigner(rho::DenseOperator{B,B}, x::Number, y::Number) where B<:FockBasis b = basis(rho) - @assert isa(b, FockBasis) N = b.N::Int _2α = complex(convert(Float64, x), convert(Float64, y))*sqrt(2) abs2_2α = abs2(_2α) @@ -114,9 +109,8 @@ function wigner(rho::DenseOperator, x::Number, y::Number) exp(-abs2_2α/2)/pi*real(w) end -function wigner(rho::DenseOperator, xvec::Vector{Float64}, yvec::Vector{Float64}) +function wigner(rho::DenseOperator{B,B}, xvec::Vector{Float64}, yvec::Vector{Float64}) where B<:FockBasis b = basis(rho) - @assert isa(b, FockBasis) N = b.N::Int _2α = [complex(x, y)*sqrt(2) for x=xvec, y=yvec] abs2_2α = abs2.(_2α) @@ -244,11 +238,10 @@ resolution `(Ntheta, Nphi)`. This version calculates the Husimi Q SU(2) function at a position given by θ and ϕ. """ -function qfuncsu2(psi::Ket, Ntheta::Int; Nphi::Int=2Ntheta) +function qfuncsu2(psi::Ket{B}, Ntheta::Int; Nphi::Int=2Ntheta) where B<:SpinBasis b = psi.basis psi_bra_data = psi.data' lb = float(b.spinnumber) - @assert isa(b, SpinBasis) result = Array{Float64}(undef, Ntheta,Nphi) @inbounds for i = 0:Ntheta-1, j = 0:Nphi-1 result[i+1,j+1] = (2*lb+1)/(4pi)*abs2(psi_bra_data*coherentspinstate(b,pi-i*pi/(Ntheta-1),j*2pi/(Nphi-1)-pi).data) @@ -256,10 +249,9 @@ function qfuncsu2(psi::Ket, Ntheta::Int; Nphi::Int=2Ntheta) return result end -function qfuncsu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta) +function qfuncsu2(rho::DenseOperator{B,B}, Ntheta::Int; Nphi::Int=2Ntheta) where B<:SpinBasis b = basis(rho) lb = float(b.spinnumber) - @assert isa(b, SpinBasis) result = Array{Float64}(undef, Ntheta,Nphi) @inbounds for i = 0:Ntheta-1, j = 0:Nphi-1 c = coherentspinstate(b,pi-i*1pi/(Ntheta-1),j*2pi/(Nphi-1)-pi) @@ -268,19 +260,17 @@ function qfuncsu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta) return result end -function qfuncsu2(psi::Ket, theta::Real, phi::Real) +function qfuncsu2(psi::Ket{B}, theta::Real, phi::Real) where B<:SpinBasis b = basis(psi) psi_bra_data = psi.data' lb = float(b.spinnumber) - @assert isa(b, SpinBasis) result = (2*lb+1)/(4pi)*abs2(psi_bra_data*coherentspinstate(b,theta,phi).data) return result end -function qfuncsu2(rho::DenseOperator, theta::Real, phi::Real) +function qfuncsu2(rho::DenseOperator{B,B}, theta::Real, phi::Real) where B<:SpinBasis b = basis(rho) lb = float(b.spinnumber) - @assert isa(b, SpinBasis) c = coherentspinstate(b,theta,phi) result = abs((2*lb+1)/(4pi)*c.data'*rho.data*c.data) return result @@ -301,7 +291,7 @@ decomposing the state into the basis of spherical harmonics. This version calculates the Wigner SU(2) function at a position given by θ and ϕ """ -function wignersu2(rho::DenseOperator, theta::Real, phi::Real) +function wignersu2(rho::DenseOperator{B,B}, theta::Real, phi::Real) where B<:SpinBasis N = length(basis(rho))-1 @@ -351,7 +341,7 @@ function wignersu2(rho::DenseOperator, theta::Real, phi::Real) return wignermap*sqrt((N+1)/(4pi)) end -function wignersu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta) +function wignersu2(rho::DenseOperator{B,B}, Ntheta::Int; Nphi::Int=2Ntheta) where B<:SpinBasis N = length(basis(rho))-1 diff --git a/src/printing.jl b/src/printing.jl index 85ad5d3a..1ee7aa1e 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -104,7 +104,7 @@ function show(stream::IO, x::Bra) end end -function summary(stream::IO, x::Operator) +function summary(stream::IO, x::AbstractOperator) print(stream, "$(typeof(x).name.name)(dim=$(length(x.basis_l))x$(length(x.basis_r)))\n") if bases.samebases(x) print(stream, " basis: ") @@ -117,7 +117,7 @@ function summary(stream::IO, x::Operator) end end -show(stream::IO, x::Operator) = summary(stream, x) +show(stream::IO, x::AbstractOperator) = summary(stream, x) function show(stream::IO, x::DenseOperator) summary(stream, x) diff --git a/src/schroedinger.jl b/src/schroedinger.jl index 45e25bec..14847ad8 100644 --- a/src/schroedinger.jl +++ b/src/schroedinger.jl @@ -21,11 +21,10 @@ Integrate Schroedinger equation. normalized nor permanent! It is still in use by the ode solver and therefore must not be changed. """ -function schroedinger(tspan, psi0::T, H::Operator; +function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}; fout::Union{Function,Nothing}=nothing, - kwargs...) where T<:StateVector + kwargs...) where {B<:Basis,T<:StateVector{B}} tspan_ = convert(Vector{Float64}, tspan) - check_schroedinger(psi0, H) dschroedinger_(t::Float64, psi::T, dpsi::T) = dschroedinger(psi, H, dpsi) x0 = psi0.data state = T(psi0.basis, psi0.data) @@ -60,16 +59,16 @@ function schroedinger_dynamic(tspan, psi0::T, f::Function; end -recast!(x::Vector{ComplexF64}, psi::StateVector) = (psi.data = x); -recast!(psi::StateVector, x::Vector{ComplexF64}) = nothing +recast!(x::D, psi::StateVector{B,D}) where {B<:Basis, D<:Vector{ComplexF64}} = (psi.data = x); +recast!(psi::StateVector{B,D}, x::D) where {B<:Basis, D<:Vector{ComplexF64}} = nothing -function dschroedinger(psi::Ket, H::Operator, dpsi::Ket) +function dschroedinger(psi::Ket{B}, H::AbstractOperator{B,B}, dpsi::Ket{B}) where B<:Basis operators.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi) return dpsi end -function dschroedinger(psi::Bra, H::Operator, dpsi::Bra) +function dschroedinger(psi::Bra{B}, H::AbstractOperator{B,B}, dpsi::Bra{B}) where B<:Basis operators.gemv!(complex(0,1.), psi, H, complex(0.), dpsi) return dpsi end @@ -77,17 +76,17 @@ end function dschroedinger_dynamic(t::Float64, psi0::T, f::Function, dpsi::T) where T<:StateVector H = f(t, psi0) - QO_CHECKS[] && check_schroedinger(psi0, H) + # QO_CHECKS[] && check_schroedinger(psi0, H) dschroedinger(psi0, H, dpsi) end -function check_schroedinger(psi::Ket, H::Operator) +function check_schroedinger(psi::Ket, H::AbstractOperator) check_multiplicable(H, psi) check_samebases(H) end -function check_schroedinger(psi::Bra, H::Operator) +function check_schroedinger(psi::Bra, H::AbstractOperator) check_multiplicable(psi, H) check_samebases(H) end diff --git a/src/semiclassical.jl b/src/semiclassical.jl index fcf86a95..613470c5 100644 --- a/src/semiclassical.jl +++ b/src/semiclassical.jl @@ -7,7 +7,7 @@ import ..timeevolution: integrate, recast! using ..bases, ..states, ..operators, ..operators_dense, ..timeevolution -const QuantumState = Union{Ket, DenseOperator} +const QuantumState{B} = Union{Ket{B}, DenseOperator{B,B}} const DecayRates = Union{Nothing, Vector{Float64}, Matrix{Float64}} """ @@ -16,15 +16,18 @@ Semi-classical state. It consists of a quantum part, which is either a `Ket` or a `DenseOperator` and a classical part that is specified as a complex vector of arbitrary length. """ -mutable struct State{T<:QuantumState} +mutable struct State{B<:Basis,T<:QuantumState{B},C<:Vector{ComplexF64}} quantum::T - classical::Vector{ComplexF64} + classical::C + function State(quantum::T, classical::C) where {B<:Basis,T<:QuantumState{B},C<:Vector{ComplexF64}} + new{B,T,C}(quantum, classical) + end end Base.length(state::State) = length(state.quantum) + length(state.classical) Base.copy(state::State) = State(copy(state.quantum), copy(state.classical)) -function ==(a::State{T}, b::State{T}) where T<:QuantumState +function ==(a::State, b::State) samebases(a.quantum, b.quantum) && length(a.classical)==length(b.classical) && (a.classical==b.classical) && @@ -33,9 +36,9 @@ end operators.expect(op, state::State) = expect(op, state.quantum) operators.variance(op, state::State) = variance(op, state.quantum) -operators.ptrace(state::State, indices::Vector{Int}) = State{DenseOperator}(ptrace(state.quantum, indices), state.classical) +operators.ptrace(state::State, indices::Vector{Int}) = State(ptrace(state.quantum, indices), state.classical) -operators_dense.dm(x::State{Ket}) = State{DenseOperator}(dm(x.quantum), x.classical) +operators_dense.dm(x::State{B,T}) where {B<:Basis,T<:Ket{B}} = State(dm(x.quantum), x.classical) """ @@ -57,11 +60,11 @@ Integrate time-dependent Schrödinger equation coupled to a classical system. normalized nor permanent! * `kwargs...`: Further arguments are passed on to the ode solver. """ -function schroedinger_dynamic(tspan, state0::State{Ket}, fquantum::Function, fclassical::Function; +function schroedinger_dynamic(tspan, state0::S, fquantum::Function, fclassical::Function; fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:Ket{B},S<:State{B,T}} tspan_ = convert(Vector{Float64}, tspan) - dschroedinger_(t, state::State{Ket}, dstate::State{Ket}) = dschroedinger_dynamic(t, state, fquantum, fclassical, dstate) + dschroedinger_(t::Float64, state::S, dstate::S) = dschroedinger_dynamic(t, state, fquantum, fclassical, dstate) x0 = Vector{ComplexF64}(undef, length(state0)) recast!(state0, x0) state = copy(state0) @@ -88,13 +91,13 @@ Integrate time-dependent master equation coupled to a classical system. permanent! * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master_dynamic(tspan, state0::State{DenseOperator}, fquantum, fclassical; - rates::Union{Vector{Float64}, Matrix{Float64}, Nothing}=nothing, +function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical; + rates::DecayRates=nothing, fout::Union{Function,Nothing}=nothing, - tmp::DenseOperator=copy(state0.quantum), - kwargs...) + tmp::T=copy(state0.quantum), + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tspan_ = convert(Vector{Float64}, tspan) - function dmaster_(t, state::State{DenseOperator}, dstate::State{DenseOperator}) + function dmaster_(t::Float64, state::S, dstate::S) where {B<:Basis,T<:DenseOperator{B,B},S<:State{B,T}} dmaster_h_dynamic(t, state, fquantum, fclassical, rates, dstate, tmp) end x0 = Vector{ComplexF64}(undef, length(state0)) @@ -104,33 +107,33 @@ function master_dynamic(tspan, state0::State{DenseOperator}, fquantum, fclassica integrate(tspan_, dmaster_, x0, state, dstate, fout; kwargs...) end -function master_dynamic(tspan, state0::State{Ket}, fquantum, fclassical; kwargs...) +function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical; kwargs...) where {B<:Basis,T<:Ket{B}} master_dynamic(tspan, dm(state0), fquantum, fclassical; kwargs...) end -function recast!(state::State, x::Vector{ComplexF64}) +function recast!(state::State{B,T,C}, x::C) where {B<:Basis,T<:QuantumState{B},C<:Vector{ComplexF64}} N = length(state.quantum) copyto!(x, 1, state.quantum.data, 1, N) copyto!(x, N+1, state.classical, 1, length(state.classical)) x end -function recast!(x::Vector{ComplexF64}, state::State) +function recast!(x::C, state::State{B,T,C}) where {B<:Basis,T<:QuantumState{B},C<:Vector{ComplexF64}} N = length(state.quantum) copyto!(state.quantum.data, 1, x, 1, N) copyto!(state.classical, 1, x, N+1, length(state.classical)) end -function dschroedinger_dynamic(t::Float64, state::State{Ket}, fquantum::Function, - fclassical::Function, dstate::State{Ket}) +function dschroedinger_dynamic(t::Float64, state::State{B,T}, fquantum::Function, + fclassical::Function, dstate::State{B,T}) where {B<:Basis,T<:Ket{B}} fquantum_(t, psi) = fquantum(t, state.quantum, state.classical) timeevolution.timeevolution_schroedinger.dschroedinger_dynamic(t, state.quantum, fquantum_, dstate.quantum) fclassical(t, state.quantum, state.classical, dstate.classical) end -function dmaster_h_dynamic(t::Float64, state::State{DenseOperator}, fquantum::Function, - fclassical::Function, rates::DecayRates, dstate::State{DenseOperator}, tmp::DenseOperator) +function dmaster_h_dynamic(t::Float64, state::State{B,T}, fquantum::Function, + fclassical::Function, rates::DecayRates, dstate::State{B,T}, tmp::T) where {B<:Basis,T<:DenseOperator{B,B}} fquantum_(t, rho) = fquantum(t, state.quantum, state.classical) timeevolution.timeevolution_master.dmaster_h_dynamic(t, state.quantum, fquantum_, rates, dstate.quantum, tmp) fclassical(t, state.quantum, state.classical, dstate.classical) diff --git a/src/spectralanalysis.jl b/src/spectralanalysis.jl index f097816c..8daaaf5c 100644 --- a/src/spectralanalysis.jl +++ b/src/spectralanalysis.jl @@ -9,7 +9,7 @@ using Arpack, LinearAlgebra const nonhermitian_warning = "The given operator is not hermitian. If this is due to a numerical error make the operator hermitian first by calculating (x+dagger(x))/2 first." """ - eigenstates(op::Operator[, n::Int; warning=true]) + eigenstates(op::AbstractOperator[, n::Int; warning=true]) Calculate the lowest n eigenvalues and their corresponding eigenstates. @@ -61,7 +61,7 @@ end """ - eigenenergies(op::Operator[, n::Int; warning=true]) + eigenenergies(op::AbstractOperator[, n::Int; warning=true]) Calculate the lowest n eigenvalues. @@ -93,8 +93,8 @@ eigenenergies(op::SparseOperator, n::Int=6; kwargs...) = eigenstates(op, n; kwar arithmetic_unary_error = operators.arithmetic_unary_error -eigenstates(op::Operator, n::Int=0) = arithmetic_unary_error("eigenstates", op) -eigenenergies(op::Operator, n::Int=0) = arithmetic_unary_error("eigenenergies", op) +eigenstates(op::AbstractOperator, n::Int=0) = arithmetic_unary_error("eigenstates", op) +eigenenergies(op::AbstractOperator, n::Int=0) = arithmetic_unary_error("eigenenergies", op) """ diff --git a/src/spin.jl b/src/spin.jl index e2df7419..f6fa1e31 100644 --- a/src/spin.jl +++ b/src/spin.jl @@ -16,18 +16,20 @@ The basis can be created for arbitrary spinnumbers by using a rational number, e.g. `SpinBasis(3//2)`. The Pauli operators are defined for all possible spin numbers. """ -mutable struct SpinBasis <: Basis +mutable struct SpinBasis{S} <: Basis shape::Vector{Int} spinnumber::Rational{Int} - function SpinBasis(spinnumber::Rational{Int}) + function SpinBasis{S}(spinnumber::Rational{Int}) where S + @assert isa(S, Rational{Int}) n = numerator(spinnumber) d = denominator(spinnumber) @assert d==2 || d==1 @assert n > 0 N = numerator(spinnumber*2 + 1) - new([N], spinnumber) + new{spinnumber}([N], spinnumber) end end +SpinBasis(spinnumber::Rational{Int}) = SpinBasis{spinnumber}(spinnumber) SpinBasis(spinnumber::Int) = SpinBasis(convert(Rational{Int}, spinnumber)) ==(b1::SpinBasis, b2::SpinBasis) = b1.spinnumber==b2.spinnumber diff --git a/src/state_definitions.jl b/src/state_definitions.jl index 6cf80897..a6fca772 100644 --- a/src/state_definitions.jl +++ b/src/state_definitions.jl @@ -29,7 +29,7 @@ randoperator(b::Basis) = randoperator(b, b) Thermal state ``exp(-H/T)/Tr[exp(-H/T)]``. """ -function thermalstate(H::Operator,T::Real) +function thermalstate(H::AbstractOperator,T::Real) return normalize(exp(-dense(H)/T)) end @@ -38,7 +38,7 @@ end Coherent thermal state ``D(α)exp(-H/T)/Tr[exp(-H/T)]D^†(α)``. """ -function coherentthermalstate(basis::FockBasis,H::Operator,T::Real,alpha::Number) +function coherentthermalstate(basis::B,H::AbstractOperator{B,B},T::Real,alpha::Number) where B<:FockBasis return displace(basis,alpha)*thermalstate(H,T)*dagger(displace(basis,alpha)) end diff --git a/src/states.jl b/src/states.jl index f30131f7..6c1883a6 100644 --- a/src/states.jl +++ b/src/states.jl @@ -19,17 +19,17 @@ in respect to a certain basis. These coefficients are stored in the `data` field and the basis is defined in the `basis` field. """ -abstract type StateVector end +abstract type StateVector{B<:Basis,T<:Vector{ComplexF64}} end """ Bra(b::Basis[, data]) Bra state defined by coefficients in respect to the basis. """ -mutable struct Bra <: StateVector - basis::Basis - data::Vector{ComplexF64} - function Bra(b::Basis, data) +mutable struct Bra{B<:Basis,T<:Vector{ComplexF64}} <: StateVector{B,T} + basis::B + data::T + function Bra{B,T}(b::B, data::T) where {B<:Basis,T<:Vector{ComplexF64}} if length(b) != length(data) throw(DimensionMismatch()) end @@ -42,10 +42,10 @@ end Ket state defined by coefficients in respect to the given basis. """ -mutable struct Ket <: StateVector - basis::Basis - data::Vector{ComplexF64} - function Ket(b::Basis, data) +mutable struct Ket{B<:Basis,T<:Vector{ComplexF64}} <: StateVector{B,T} + basis::B + data::T + function Ket{B,T}(b::B, data::T) where {B<:Basis,T<:Vector{ComplexF64}} if length(b) != length(data) throw(DimensionMismatch()) end @@ -53,22 +53,44 @@ mutable struct Ket <: StateVector end end +Bra{B}(b::B, data::T) where {B<:Basis,T<:Vector{ComplexF64}} = Bra{B,T}(b, data) +Ket{B}(b::B, data::T) where {B<:Basis,T<:Vector{ComplexF64}} = Ket{B,T}(b, data) + +Bra(b::B, data::T) where {B<:Basis,T<:Vector{ComplexF64}} = Bra{B,T}(b, data) +Ket(b::B, data::T) where {B<:Basis,T<:Vector{ComplexF64}} = Ket{B,T}(b, data) + +Bra{B}(b::B) where B<:Basis = Bra{B}(b, zeros(ComplexF64, length(b))) +Ket{B}(b::B) where B<:Basis = Ket{B}(b, zeros(ComplexF64, length(b))) Bra(b::Basis) = Bra(b, zeros(ComplexF64, length(b))) Ket(b::Basis) = Ket(b, zeros(ComplexF64, length(b))) +Ket(b::Basis, data) = Ket(b, convert(Vector{ComplexF64}, data)) +Bra(b::Basis, data) = Ket(b, convert(Vector{ComplexF64}, data)) + copy(a::T) where {T<:StateVector} = T(a.basis, copy(a.data)) length(a::StateVector) = length(a.basis)::Int basis(a::StateVector) = a.basis -==(x::T, y::T) where {T<:StateVector} = samebases(x, y) && x.data==y.data +==(x::T, y::T) where {T<:Ket} = samebases(x, y) && x.data==y.data +==(x::T, y::T) where {T<:Bra} = samebases(x, y) && x.data==y.data +==(x::Ket, y::Ket) = false +==(x::Bra, y::Bra) = false # Arithmetic operations -+(a::T, b::T) where {T<:StateVector} = (check_samebases(a, b); T(a.basis, a.data+b.data)) ++(a::Ket{B}, b::Ket{B}) where {B<:Basis} = Ket(a.basis, a.data+b.data) ++(a::Bra{B}, b::Bra{B}) where {B<:Basis} = Bra(a.basis, a.data+b.data) ++(a::Ket, b::Ket) = throw(bases.IncompatibleBases()) ++(a::Bra, b::Bra) = throw(bases.IncompatibleBases()) + +-(a::Ket{B}, b::Ket{B}) where {B<:Basis} = Ket(a.basis, a.data-b.data) +-(a::Bra{B}, b::Bra{B}) where {B<:Basis} = Bra(a.basis, a.data-b.data) +-(a::Ket, b::Ket) = throw(bases.IncompatibleBases()) +-(a::Bra, b::Bra) = throw(bases.IncompatibleBases()) -(a::T) where {T<:StateVector} = T(a.basis, -a.data) --(a::T, b::T) where {T<:StateVector} = (check_samebases(a, b); T(a.basis, a.data-b.data)) -*(a::Bra, b::Ket) = (check_multiplicable(a, b); sum(a.data.*b.data)) +*(a::Bra{B,D}, b::Ket{B,D}) where {B<:Basis,D<:Vector{ComplexF64}} = transpose(a.data)*b.data +*(a::Bra, b::Ket) = throw(bases.IncompatibleBases()) *(a::Number, b::T) where {T<:StateVector} = T(b.basis, a*b.data) *(a::T, b::Number) where {T<:StateVector} = T(a.basis, b*a.data) @@ -88,9 +110,12 @@ dagger(x::Ket) = Bra(x.basis, conj(x.data)) Tensor product ``|x⟩⊗|y⟩⊗|z⟩⊗…`` of the given states. """ -tensor(a::T, b::T) where {T<:StateVector} = T(tensor(a.basis, b.basis), kron(b.data, a.data)) +tensor(a::Ket, b::Ket) = Ket(tensor(a.basis, b.basis), kron(b.data, a.data)) +tensor(a::Bra, b::Bra) = Bra(tensor(a.basis, b.basis), kron(b.data, a.data)) tensor(state::StateVector) = state -tensor(states::T...) where {T<:StateVector} = reduce(tensor, states) +tensor(states::Ket...) = reduce(tensor, states) +tensor(states::Bra...) = reduce(tensor, states) +tensor(states::Vector{T}) where T<:StateVector = reduce(tensor, states) # Normalization functions """ @@ -112,13 +137,21 @@ In-place normalization of the given bra or ket so that `norm(x)` is one. """ normalize!(x::StateVector) = (rmul!(x.data, 1.0/norm(x)); nothing) -function permutesystems(state::T, perm::Vector{Int}) where T<:StateVector +function permutesystems(state::T, perm::Vector{Int}) where T<:Ket + @assert length(state.basis.bases) == length(perm) + @assert isperm(perm) + data = reshape(state.data, state.basis.shape...) + data = permutedims(data, perm) + data = reshape(data, length(data)) + Ket(permutesystems(state.basis, perm), data) +end +function permutesystems(state::T, perm::Vector{Int}) where T<:Bra @assert length(state.basis.bases) == length(perm) @assert isperm(perm) data = reshape(state.data, state.basis.shape...) data = permutedims(data, perm) data = reshape(data, length(data)) - T(permutesystems(state.basis, perm), data) + Bra(permutesystems(state.basis, perm), data) end # Creation of basis states. diff --git a/src/steadystate.jl b/src/steadystate.jl index 0af87e4c..e37785ef 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -1,6 +1,6 @@ module steadystate -using ..states, ..operators, ..operators_dense, ..superoperators +using ..states, ..operators, ..operators_dense, ..superoperators, ..bases using ..timeevolution using Arpack, LinearAlgebra @@ -26,13 +26,13 @@ Calculate steady state using long time master equation evolution. density operator `rho` is further used and therefore must not be changed. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master(H::Operator, J::Vector; - rho0::DenseOperator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), +function master(H::AbstractOperator{B,B}, J::Vector; + rho0::DenseOperator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), hmin=1e-7, tol=1e-3, rates::Union{Vector{Float64}, Matrix{Float64}, Nothing}=nothing, Jdagger::Vector=dagger.(J), fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where B<:Basis t,u = timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger, hmin=hmin, hmax=Inf, display_initialvalue=false, @@ -88,7 +88,7 @@ function liouvillianspectrum(L::SparseSuperOperator; nev::Int = min(10, length(L return d[indices], ops end -liouvillianspectrum(H::Operator, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) = liouvillianspectrum(liouvillian(H, J; rates=rates); kwargs...) +liouvillianspectrum(H::AbstractOperator{B,B}, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) where B<:Basis = liouvillianspectrum(liouvillian(H, J; rates=rates); kwargs...) """ steadystate.eigenvector(L) @@ -118,7 +118,7 @@ function eigenvector(L::SuperOperator; tol::Real = 1e-9, nev::Int = 2, which::Sy return ops[1]/tr(ops[1]) end -eigenvector(H::Operator, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) = eigenvector(liouvillian(H, J; rates=rates); kwargs...) +eigenvector(H::AbstractOperator, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) = eigenvector(liouvillian(H, J; rates=rates); kwargs...) end # module diff --git a/src/stochastic_definitions.jl b/src/stochastic_definitions.jl index ed4ad63d..56ef1174 100644 --- a/src/stochastic_definitions.jl +++ b/src/stochastic_definitions.jl @@ -44,8 +44,8 @@ and H_s = iCe^{-i\\theta}. ``` """ -function homodyne_carmichael(H0::Operator, C::Vector{T}, theta::Vector{R}; - normalize_expect::Bool=true) where {T <: Operator, R <: Real} +function homodyne_carmichael(H0::AbstractOperator, C::Vector{T}, theta::Vector{R}; + normalize_expect::Bool=true) where {T <: AbstractOperator, R <: Real} n = length(C) @assert n == length(theta) Hs = 1.0im*C .* exp.(-1.0im .* theta) @@ -67,7 +67,7 @@ function homodyne_carmichael(H0::Operator, C::Vector{T}, theta::Vector{R}; return fdeterm_un, fstoch end end -homodyne_carmichael(H0::Operator, C::Operator, theta::Real; kwargs...) = +homodyne_carmichael(H0::AbstractOperator, C::AbstractOperator, theta::Real; kwargs...) = homodyne_carmichael(H0, [C], [theta]; kwargs...) end # module diff --git a/src/stochastic_master.jl b/src/stochastic_master.jl index 459fd754..68d538e8 100644 --- a/src/stochastic_master.jl +++ b/src/stochastic_master.jl @@ -41,23 +41,23 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster. be changed. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master(tspan, rho0::DenseOperator, H::Operator, +function master(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector, C::Vector; rates::DecayRates=nothing, Jdagger::Vector=dagger.(J), Cdagger::Vector=dagger.(C), fout::Union{Function,Nothing}=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tmp = copy(rho0) n = length(C) - dmaster_stoch(dx::DiffArray, t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = + dmaster_stoch(dx::DiffArray, t::Float64, rho::T, drho::T, n::Int) = dmaster_stochastic(dx, rho, C, Cdagger, drho, n) isreducible = check_master(rho0, H, J, Jdagger, rates) && check_master_stoch(rho0, C, Cdagger) if !isreducible - dmaster_h_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = + dmaster_h_determ(t::Float64, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch, rho0, fout, n; kwargs...) else @@ -77,7 +77,7 @@ function master(tspan, rho0::DenseOperator, H::Operator, end Hnhdagger = dagger(Hnh) - dmaster_nh_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = + dmaster_nh_determ(t::Float64, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch, rho0, fout, n; kwargs...) end @@ -113,11 +113,11 @@ dynamic Hamiltonian and J. number if you want to avoid an initial calculation of function outputs! * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master_dynamic(tspan::Vector{Float64}, rho0::DenseOperator, fdeterm::Function, fstoch::Function; +function master_dynamic(tspan::Vector{Float64}, rho0::T, fdeterm::Function, fstoch::Function; rates::DecayRates=nothing, fout::Union{Function,Nothing}=nothing, noise_processes::Int=0, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tmp = copy(rho0) @@ -128,23 +128,23 @@ function master_dynamic(tspan::Vector{Float64}, rho0::DenseOperator, fdeterm::Fu n = noise_processes end - dmaster_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp) - dmaster_stoch(dx::DiffArray, t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = + dmaster_determ(t::Float64, rho::T, drho::T) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp) + dmaster_stoch(dx::DiffArray, t::Float64, rho::T, drho::T, n::Int) = dmaster_stoch_dynamic(dx, t, rho, fstoch, drho, n) integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; kwargs...) end master_dynamic(tspan::Vector{Float64}, psi0::Ket, args...; kwargs...) = master_dynamic(tspan, dm(psi0), args...; kwargs...) # Derivative functions -function dmaster_stochastic(dx::Vector{ComplexF64}, rho::DenseOperator, - C::Vector, Cdagger::Vector, drho::DenseOperator, ::Int) +function dmaster_stochastic(dx::Vector{ComplexF64}, rho::T, + C::Vector, Cdagger::Vector, drho::T, ::Int) where {B<:Basis,T<:DenseOperator{B,B}} recast!(dx, drho) operators.gemm!(1, C[1], rho, 0, drho) operators.gemm!(1, rho, Cdagger[1], 1, drho) drho.data .-= tr(drho)*rho.data end -function dmaster_stochastic(dx::Array{ComplexF64, 2}, rho::DenseOperator, - C::Vector, Cdagger::Vector, drho::DenseOperator, n::Int) +function dmaster_stochastic(dx::Array{ComplexF64, 2}, rho::T, + C::Vector, Cdagger::Vector, drho::T, n::Int) where {B<:Basis,T<:DenseOperator{B,B}} for i=1:n dx_i = @view dx[:, i] recast!(dx_i, drho) @@ -155,8 +155,8 @@ function dmaster_stochastic(dx::Array{ComplexF64, 2}, rho::DenseOperator, end end -function dmaster_stoch_dynamic(dx::DiffArray, t::Float64, rho::DenseOperator, - f::Function, drho::DenseOperator, n::Int) +function dmaster_stoch_dynamic(dx::DiffArray, t::Float64, rho::T, + f::Function, drho::T, n::Int) where {B<:Basis,T<:DenseOperator{B,B}} result = f(t, rho) QO_CHECKS[] && @assert 2 == length(result) C, Cdagger = result @@ -175,32 +175,31 @@ function integrate_master_stoch(tspan, df::Function, dg::Function, integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...) end -function check_master_stoch(rho0::DenseOperator, C::Vector, Cdagger::Vector) +function check_master_stoch(rho0::DenseOperator{B,B}, C::Vector, Cdagger::Vector) where B<:Basis + # TODO: replace type checks by dispatch; make types of C known @assert length(C) == length(Cdagger) isreducible = true for c=C - @assert isa(c, Operator) + @assert isa(c, AbstractOperator{B,B}) if !(isa(c, DenseOperator) || isa(c, SparseOperator)) isreducible = false end - check_samebases(rho0, c) end for c=Cdagger - @assert isa(c, Operator) + @assert isa(c, AbstractOperator{B,B}) if !(isa(c, DenseOperator) || isa(c, SparseOperator)) isreducible = false end - check_samebases(rho0, c) end isreducible end # TODO: Speed up by recasting to n-d arrays, remove vector methods -function recast!(x::Union{Vector{ComplexF64}, SubArray{ComplexF64, 1}}, rho::DenseOperator) +function recast!(x::Union{Vector{ComplexF64}, SubArray{ComplexF64, 1}}, rho::DenseOperator{B,B,T}) where {B<:Basis,T<:Matrix{ComplexF64}} rho.data = reshape(x, size(rho.data)) end -recast!(state::DenseOperator, x::SubArray{ComplexF64, 1}) = (x[:] = state.data) -recast!(state::DenseOperator, x::Vector{ComplexF64}) = nothing +recast!(state::DenseOperator{B,B}, x::SubArray{ComplexF64, 1}) where B<:Basis = (x[:] = state.data) +recast!(state::DenseOperator{B,B}, x::Vector{ComplexF64}) where B<:Basis = nothing end # module diff --git a/src/stochastic_schroedinger.jl b/src/stochastic_schroedinger.jl index f2723855..a3171a0c 100644 --- a/src/stochastic_schroedinger.jl +++ b/src/stochastic_schroedinger.jl @@ -32,11 +32,11 @@ Integrate stochastic Schrödinger equation. each time step taken by the solver. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function schroedinger(tspan, psi0::Ket, H::Operator, Hs::Vector; +function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector; fout::Union{Function,Nothing}=nothing, normalize_state::Bool=false, calback=nothing, - kwargs...) + kwargs...) where {B<:Basis,T<:Ket{B}} tspan_ = convert(Vector{Float64}, tspan) n = length(Hs) @@ -44,14 +44,14 @@ function schroedinger(tspan, psi0::Ket, H::Operator, Hs::Vector; x0 = psi0.data state = copy(psi0) - check_schroedinger(psi0, H) + # TODO: replace checks by dispatch for h=Hs - check_schroedinger(psi0, h) + @assert isa(h, AbstractOperator{B,B}) end - dschroedinger_determ(t::Float64, psi::Ket, dpsi::Ket) = dschroedinger(psi, H, dpsi) + dschroedinger_determ(t::Float64, psi::T, dpsi::T) = dschroedinger(psi, H, dpsi) dschroedinger_stoch(dx::DiffArray, - t::Float64, psi::Ket, dpsi::Ket, n::Int) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n) + t::Float64, psi::T, dpsi::T, n::Int) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n) if normalize_state norm_func(u::Vector{ComplexF64}, t::Float64, integrator) = normalize!(u) @@ -66,7 +66,7 @@ function schroedinger(tspan, psi0::Ket, H::Operator, Hs::Vector; ncb=ncb, kwargs...) end -schroedinger(tspan, psi0::Ket, H::Operator, Hs::Operator; kwargs...) = schroedinger(tspan, psi0, H, [Hs]; kwargs...) +schroedinger(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, Hs::AbstractOperator{B,B}; kwargs...) where B<:Basis = schroedinger(tspan, psi0, H, [Hs]; kwargs...) """ stochastic.schroedinger_dynamic(tspan, state0, fdeterm, fstoch[; fout, ...]) @@ -94,10 +94,10 @@ Integrate stochastic Schrödinger equation with dynamic Hamiltonian. each time step taken by the solver. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function schroedinger_dynamic(tspan, psi0::Ket, fdeterm::Function, fstoch::Function; +function schroedinger_dynamic(tspan, psi0::T, fdeterm::Function, fstoch::Function; fout::Union{Function,Nothing}=nothing, noise_processes::Int=0, normalize_state::Bool=false, - kwargs...) + kwargs...) where T<:Ket tspan_ = convert(Vector{Float64}, tspan) if noise_processes == 0 @@ -111,9 +111,9 @@ function schroedinger_dynamic(tspan, psi0::Ket, fdeterm::Function, fstoch::Funct x0 = psi0.data state = copy(psi0) - dschroedinger_determ(t::Float64, psi::Ket, dpsi::Ket) = dschroedinger_dynamic(t, psi, fdeterm, dpsi) + dschroedinger_determ(t::Float64, psi::T, dpsi::T) = dschroedinger_dynamic(t, psi, fdeterm, dpsi) dschroedinger_stoch(dx::DiffArray, - t::Float64, psi::Ket, dpsi::Ket, n::Int) = + t::Float64, psi::T, dpsi::T, n::Int) = dschroedinger_stochastic(dx, t, psi, fstoch, dpsi, n) if normalize_state @@ -132,13 +132,13 @@ function schroedinger_dynamic(tspan, psi0::Ket, fdeterm::Function, fstoch::Funct end -function dschroedinger_stochastic(dx::Vector{ComplexF64}, psi::Ket, Hs::Vector{T}, - dpsi::Ket, index::Int) where T <: Operator +function dschroedinger_stochastic(dx::D, psi::T1, Hs::Vector{T2}, + dpsi::T1, index::Int) where {D<:Vector{ComplexF64},B<:Basis,T1<:Ket{B},T2<:AbstractOperator{B,B}} recast!(dx, dpsi) dschroedinger(psi, Hs[index], dpsi) end -function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, psi::Ket, Hs::Vector{T}, - dpsi::Ket, n::Int) where T <: Operator +function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, psi::T1, Hs::Vector{T2}, + dpsi::T1, n::Int) where {B<:Basis,T1<:Ket{B},T2<:AbstractOperator{B,B}} for i=1:n dx_i = @view dx[:, i] recast!(dx_i, dpsi) @@ -147,7 +147,7 @@ function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, psi::Ket, Hs::Vector end end function dschroedinger_stochastic(dx::DiffArray, - t::Float64, psi::Ket, f::Function, dpsi::Ket, n::Int) + t::Float64, psi::T, f::Function, dpsi::T, n::Int) where T<:Ket ops = f(t, psi) if QO_CHECKS[] @inbounds for h=ops @@ -157,7 +157,7 @@ function dschroedinger_stochastic(dx::DiffArray, dschroedinger_stochastic(dx, psi, ops, dpsi, n) end -recast!(psi::StateVector, x::SubArray{ComplexF64, 1}) = (x .= psi.data) -recast!(x::SubArray{ComplexF64, 1}, psi::StateVector) = (psi.data = x) +recast!(psi::Ket, x::SubArray{ComplexF64, 1}) = (x .= psi.data) +recast!(x::SubArray{ComplexF64, 1}, psi::Ket) = (psi.data = x) end # module diff --git a/src/stochastic_semiclassical.jl b/src/stochastic_semiclassical.jl index bfd662b9..2a4e7c70 100644 --- a/src/stochastic_semiclassical.jl +++ b/src/stochastic_semiclassical.jl @@ -56,16 +56,17 @@ Integrate time-dependent Schrödinger equation coupled to a classical system. each time step taken by the solver. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function schroedinger_semiclassical(tspan, state0::State{Ket}, fquantum::Function, +function schroedinger_semiclassical(tspan, state0::S, fquantum::Function, fclassical::Function; fstoch_quantum::Union{Nothing, Function}=nothing, fstoch_classical::Union{Nothing, Function}=nothing, fout::Union{Function,Nothing}=nothing, noise_processes::Int=0, noise_prototype_classical=nothing, normalize_state::Bool=false, - kwargs...) + kwargs...) where {B<:Basis,T<:Ket{B},S<:State{B,T}} tspan_ = convert(Vector{Float64}, tspan) - dschroedinger_det(t::Float64, state::State{Ket}, dstate::State{Ket}) = semiclassical.dschroedinger_dynamic(t, state, fquantum, fclassical, dstate) + dschroedinger_det(t::Float64, state::S, dstate::S) = + semiclassical.dschroedinger_dynamic(t, state, fquantum, fclassical, dstate) if isa(fstoch_quantum, Nothing) && isa(fstoch_classical, Nothing) throw(ArgumentError("No stochastic functions provided!")) @@ -104,9 +105,9 @@ function schroedinger_semiclassical(tspan, state0::State{Ket}, fquantum::Functio ncb = nothing end - dschroedinger_stoch(dx::DiffArray, - t::Float64, state::State{Ket}, dstate::State{Ket}, n::Int) = + dschroedinger_stoch(dx::DiffArray, t::Float64, state::S, dstate::S, n::Int) = dschroedinger_stochastic(dx, t, state, fstoch_quantum, fstoch_classical, dstate, n) + integrate_stoch(tspan_, dschroedinger_det, dschroedinger_stoch, x0, state, dstate, fout, n; noise_prototype_classical = noise_prototype_classical, ncb=ncb, @@ -156,7 +157,7 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster. noise. See the documentation for details. * `kwargs...`: Further arguments are passed on to the ode solver. """ -function master_semiclassical(tspan::Vector{Float64}, rho0::State{DenseOperator}, +function master_semiclassical(tspan::Vector{Float64}, rho0::S, fquantum::Function, fclassical::Function; fstoch_quantum::Union{Function, Nothing}=nothing, fstoch_classical::Union{Function, Nothing}=nothing, @@ -165,7 +166,7 @@ function master_semiclassical(tspan::Vector{Float64}, rho0::State{DenseOperator} noise_processes::Int=0, noise_prototype_classical=nothing, nonlinear::Bool=true, - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B},S<:State{B,T}} tmp = copy(rho0.quantum) if isa(fstoch_quantum, Nothing) && isa(fstoch_classical, Nothing) @@ -188,24 +189,24 @@ function master_semiclassical(tspan::Vector{Float64}, rho0::State{DenseOperator} end end - dmaster_determ(t::Float64, rho::State{DenseOperator}, drho::State{DenseOperator}) = + dmaster_determ(t::Float64, rho::S, drho::S) = dmaster_h_dynamic(t, rho, fquantum, fclassical, rates, drho, tmp) - dmaster_stoch(dx::DiffArray, t::Float64, rho::State{DenseOperator}, - drho::State{DenseOperator}, n::Int) = + dmaster_stoch(dx::DiffArray, t::Float64, rho::S, + drho::S, n::Int) = dmaster_stoch_dynamic(dx, t, rho, fstoch_quantum, fstoch_classical, drho, n) integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; noise_prototype_classical=noise_prototype_classical, kwargs...) end -master_semiclassical(tspan::Vector{Float64}, psi0::State{Ket}, args...; kwargs...) = +master_semiclassical(tspan::Vector{Float64}, psi0::State{B,T}, args...; kwargs...) where {B<:Basis,T<:Ket{B}} = master_semiclassical(tspan, dm(psi0), args...; kwargs...) # Derivative functions function dschroedinger_stochastic(dx::Vector{ComplexF64}, t::Float64, - state::State{Ket}, fstoch_quantum::Function, fstoch_classical::Nothing, - dstate::State{Ket}, ::Int) + state::S, fstoch_quantum::Function, fstoch_classical::Nothing, + dstate::S, ::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}} H = fstoch_quantum(t, state.quantum, state.classical) recast!(dx, dstate) QO_CHECKS[] && check_schroedinger(state.quantum, H[1]) @@ -213,8 +214,8 @@ function dschroedinger_stochastic(dx::Vector{ComplexF64}, t::Float64, recast!(dstate, dx) end function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, - t::Float64, state::State{Ket}, fstoch_quantum::Function, - fstoch_classical::Nothing, dstate::State{Ket}, n::Int) + t::Float64, state::S, fstoch_quantum::Function, + fstoch_classical::Nothing, dstate::S, n::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}} H = fstoch_quantum(t, state.quantum, state.classical) for i=1:n dx_i = @view dx[:, i] @@ -225,13 +226,13 @@ function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, end end function dschroedinger_stochastic(dx::DiffArray, t::Float64, - state::State{Ket}, fstoch_quantum::Nothing, fstoch_classical::Function, - dstate::State{Ket}, ::Int) + state::S, fstoch_quantum::Nothing, fstoch_classical::Function, + dstate::S, ::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}} dclassical = @view dx[length(state.quantum)+1:end, :] fstoch_classical(t, state.quantum, state.classical, dclassical) end -function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, t::Float64, state::State{Ket}, fstoch_quantum::Function, - fstoch_classical::Function, dstate::State{Ket}, n::Int) +function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, t::Float64, state::S, fstoch_quantum::Function, + fstoch_classical::Function, dstate::S, n::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}} dschroedinger_stochastic(dx, t, state, fstoch_quantum, nothing, dstate, n) dx_i = @view dx[length(state.quantum)+1:end, n+1:end] @@ -239,8 +240,8 @@ function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, t::Float64, state::S end function dmaster_stoch_dynamic(dx::Vector{ComplexF64}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Nothing, dstate::State{DenseOperator}, ::Int) + state::S, fstoch_quantum::Function, + fstoch_classical::Nothing, dstate::S, ::Int) where {B<:Basis,T<:DenseOperator{B,B},S<:State{B,T}} result = fstoch_quantum(t, state.quantum, state.classical) QO_CHECKS[] && @assert length(result) == 2 C, Cdagger = result @@ -252,8 +253,8 @@ function dmaster_stoch_dynamic(dx::Vector{ComplexF64}, t::Float64, recast!(dstate, dx) end function dmaster_stoch_dynamic(dx::Array{ComplexF64, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Nothing, dstate::State{DenseOperator}, n::Int) + state::S, fstoch_quantum::Function, + fstoch_classical::Nothing, dstate::S, n::Int) where {B<:Basis,T<:DenseOperator{B,B},S<:State{B,T}} result = fstoch_quantum(t, state.quantum, state.classical) QO_CHECKS[] && @assert length(result) == 2 C, Cdagger = result @@ -268,14 +269,14 @@ function dmaster_stoch_dynamic(dx::Array{ComplexF64, 2}, t::Float64, end end function dmaster_stoch_dynamic(dx::DiffArray, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Nothing, - fstoch_classical::Function, dstate::State{DenseOperator}, ::Int) + state::S, fstoch_quantum::Nothing, + fstoch_classical::Function, dstate::S, ::Int) where {B<:Basis,T<:DenseOperator{B,B},S<:State{B,T}} dclassical = @view dx[length(state.quantum)+1:end, :] fstoch_classical(t, state.quantum, state.classical, dclassical) end function dmaster_stoch_dynamic(dx::Array{ComplexF64, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Function, dstate::State{DenseOperator}, n::Int) + state::S, fstoch_quantum::Function, + fstoch_classical::Function, dstate::S, n::Int) where {B<:Basis,T<:DenseOperator{B,B},S<:State{B,T}} dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, nothing, dstate, n) dx_i = @view dx[length(state.quantum)+1:end, n+1:end] @@ -283,9 +284,9 @@ function dmaster_stoch_dynamic(dx::Array{ComplexF64, 2}, t::Float64, end function integrate_master_stoch(tspan, df::Function, dg::Function, - rho0::State{DenseOperator}, fout::Union{Nothing, Function}, + rho0::State{B,T}, fout::Union{Nothing, Function}, n::Int; - kwargs...) + kwargs...) where {B<:Basis,T<:DenseOperator{B,B}} tspan_ = convert(Vector{Float64}, tspan) x0 = Vector{ComplexF64}(undef, length(rho0)) recast!(rho0, x0) diff --git a/src/subspace.jl b/src/subspace.jl index 0fb64a84..f21f53c0 100644 --- a/src/subspace.jl +++ b/src/subspace.jl @@ -13,24 +13,29 @@ using ..bases, ..states, ..operators, ..operators_dense A basis describing a subspace embedded a higher dimensional Hilbert space. """ -mutable struct SubspaceBasis <: Basis +mutable struct SubspaceBasis{B<:Basis,T<:Ket,H} <: Basis shape::Vector{Int} - superbasis::Basis - basisstates::Vector{Ket} + superbasis::B + basisstates::Vector{T} basisstates_hash::UInt - function SubspaceBasis(superbasis::Basis, basisstates::Vector{Ket}) + function SubspaceBasis{B,T,H}(superbasis::B, basisstates::Vector{T}) where {B<:Basis,T<:Ket,H} for state = basisstates if state.basis != superbasis throw(ArgumentError("The basis of the basisstates has to be the superbasis.")) end end + @assert isa(H, UInt) basisstates_hash = hash(hash.([hash.(x.data) for x=basisstates])) new(Int[length(basisstates)], superbasis, basisstates, basisstates_hash) end end -SubspaceBasis(basisstates::Vector{Ket}) = SubspaceBasis(basisstates[1].basis, basisstates) +function SubspaceBasis(superbasis::B, basisstates::Vector{T}) where {B<:Basis,T<:Ket} + basisstates_hash = hash(hash.([hash.(x.data) for x=basisstates])) + SubspaceBasis{B,T,basisstates_hash}(superbasis, basisstates) +end +SubspaceBasis(basisstates::Vector{T}) where T<:Ket = SubspaceBasis(basisstates[1].basis, basisstates) ==(b1::SubspaceBasis, b2::SubspaceBasis) = b1.superbasis==b2.superbasis && b1.basisstates_hash==b2.basisstates_hash diff --git a/src/superoperators.jl b/src/superoperators.jl index 86ed05ca..7f0b8bc4 100644 --- a/src/superoperators.jl +++ b/src/superoperators.jl @@ -25,18 +25,18 @@ A_{bl_1,bl_2} = S_{(bl_1,bl_2) ↔ (br_1,br_2)} B_{br_1,br_2} A_{br_1,br_2} = B_{bl_1,bl_2} S_{(bl_1,bl_2) ↔ (br_1,br_2)} ``` """ -abstract type SuperOperator end +abstract type SuperOperator{B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} end """ DenseSuperOperator(b1[, b2, data]) SuperOperator stored as dense matrix. """ -mutable struct DenseSuperOperator <: SuperOperator - basis_l::Tuple{Basis, Basis} - basis_r::Tuple{Basis, Basis} - data::Matrix{ComplexF64} - function DenseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}, data::Matrix{ComplexF64}) +mutable struct DenseSuperOperator{B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}} <: SuperOperator{B1,B2} + basis_l::B1 + basis_r::B2 + data::T + function DenseSuperOperator{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}} if length(basis_l[1])*length(basis_l[2]) != size(data, 1) || length(basis_r[1])*length(basis_r[2]) != size(data, 2) throw(DimensionMismatch()) @@ -44,6 +44,7 @@ mutable struct DenseSuperOperator <: SuperOperator new(basis_l, basis_r, data) end end +DenseSuperOperator(basis_l::BL, basis_r::BR, data::T) where {BL<:Tuple{Basis,Basis},BR<:Tuple{Basis,Basis},T<:Matrix{ComplexF64}} = DenseSuperOperator{BL,BR,T}(basis_l, basis_r, data) function DenseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}) Nl = length(basis_l[1])*length(basis_l[2]) @@ -58,11 +59,11 @@ end SuperOperator stored as sparse matrix. """ -mutable struct SparseSuperOperator <: SuperOperator - basis_l::Tuple{Basis, Basis} - basis_r::Tuple{Basis, Basis} - data::SparseMatrixCSC{ComplexF64, Int} - function SparseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}, data::SparseMatrixCSC{ComplexF64, Int}) +mutable struct SparseSuperOperator{B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:SparseMatrixCSC{ComplexF64,Int}} <: SuperOperator{B1,B2} + basis_l::B1 + basis_r::B2 + data::T + function SparseSuperOperator{B1,B2,T}(basis_l::B1, basis_r::B2, data::T) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:SparseMatrixCSC{ComplexF64,Int}} if length(basis_l[1])*length(basis_l[2]) != size(data, 1) || length(basis_r[1])*length(basis_r[2]) != size(data, 2) throw(DimensionMismatch()) @@ -70,6 +71,7 @@ mutable struct SparseSuperOperator <: SuperOperator new(basis_l, basis_r, data) end end +SparseSuperOperator(basis_l::B1, basis_r::B2, data::T) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T<:SparseMatrixCSC{ComplexF64,Int}} = SparseSuperOperator{B1,B2,T}(basis_l, basis_r, data) function SparseSuperOperator(basis_l::Tuple{Basis, Basis}, basis_r::Tuple{Basis, Basis}) Nl = length(basis_l[1])*length(basis_l[2]) @@ -80,6 +82,8 @@ end SuperOperator(basis_l, basis_r, data::SparseMatrixCSC{ComplexF64, Int}) = SparseSuperOperator(basis_l, basis_r, data) SuperOperator(basis_l, basis_r, data::Matrix{ComplexF64}) = DenseSuperOperator(basis_l, basis_r, data) +SuperOperator{B1,B2}(basis_l::B1, basis_r::B2, data::SparseMatrixCSC{ComplexF64, Int}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = SparseSuperOperator(basis_l, basis_r, data) +SuperOperator{B1,B2}(basis_l::B1, basis_r::B2, data::Matrix{ComplexF64}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = DenseSuperOperator(basis_l, basis_r, data) Base.copy(a::T) where {T<:SuperOperator} = T(a.basis_l, a.basis_r, copy(a.data)) @@ -90,41 +94,41 @@ operators.dense(a::DenseSuperOperator) = copy(a) sparse(a::DenseSuperOperator) = SparseSuperOperator(a.basis_l, a.basis_r, sparse(a.data)) sparse(a::SparseSuperOperator) = copy(a) -==(a::T, b::T) where {T<:SuperOperator} = samebases(a, b) && (a.data == b.data) +==(a::T, b::T) where {T<:SuperOperator} = a.data == b.data +==(a::SuperOperator, b::SuperOperator) = false Base.length(a::SuperOperator) = length(a.basis_l[1])*length(a.basis_l[2])*length(a.basis_r[1])*length(a.basis_r[2]) bases.samebases(a::SuperOperator, b::SuperOperator) = samebases(a.basis_l[1], b.basis_l[1]) && samebases(a.basis_l[2], b.basis_l[2]) && samebases(a.basis_r[1], b.basis_r[1]) && samebases(a.basis_r[2], b.basis_r[2]) bases.multiplicable(a::SuperOperator, b::SuperOperator) = multiplicable(a.basis_r[1], b.basis_l[1]) && multiplicable(a.basis_r[2], b.basis_l[2]) -bases.multiplicable(a::SuperOperator, b::Operator) = multiplicable(a.basis_r[1], b.basis_l) && multiplicable(a.basis_r[2], b.basis_r) +bases.multiplicable(a::SuperOperator, b::AbstractOperator) = multiplicable(a.basis_r[1], b.basis_l) && multiplicable(a.basis_r[2], b.basis_r) # Arithmetic operations -function *(a::SuperOperator, b::DenseOperator) - check_multiplicable(a, b) +function *(a::SuperOperator{B1,B2}, b::DenseOperator{BL,BR}) where {BL<:Basis,BR<:Basis,B1<:Tuple{Basis,Basis},B2<:Tuple{BL,BR}} data = a.data*reshape(b.data, length(b.data)) return DenseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2]))) end -function *(a::SuperOperator, b::SparseOperator) - check_multiplicable(a, b) +function *(a::SuperOperator{B1,B2}, b::SparseOperator{BL,BR}) where {BL<:Basis,BR<:Basis,B1<:Tuple{Basis,Basis},B2<:Tuple{BL,BR}} data = a.data*reshape(b.data, length(b.data)) return SparseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2]))) end -function *(a::SuperOperator, b::SuperOperator) - check_multiplicable(a, b) - return SuperOperator(a.basis_l, b.basis_r, a.data*b.data) +function *(a::SuperOperator{B1,B2}, b::SuperOperator{B2,B3}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},B3<:Tuple{Basis,Basis}} + return SuperOperator{B1,B3}(a.basis_l, b.basis_r, a.data*b.data) end -*(a::SuperOperator, b::Number) = SuperOperator(a.basis_l, a.basis_r, a.data*b) -*(a::Number, b::SuperOperator) = SuperOperator(b.basis_l, b.basis_r, a*b.data) +*(a::T, b::Number) where T<:SuperOperator = T(a.basis_l, a.basis_r, a.data*b) +*(a::Number, b::T) where T<:SuperOperator = T(b.basis_l, b.basis_r, a*b.data) -/(a::SuperOperator, b::Number) = SuperOperator(a.basis_l, a.basis_r, a.data/b) +/(a::T, b::Number) where T<:SuperOperator = T(a.basis_l, a.basis_r, a.data/b) -+(a::SuperOperator, b::SuperOperator) = (check_samebases(a, b); SuperOperator(a.basis_l, a.basis_r, a.data+b.data)) ++(a::SuperOperator{B1,B2}, b::SuperOperator{B1,B2}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = SuperOperator{B1,B2}(a.basis_l, a.basis_r, a.data+b.data) ++(a::SuperOperator, b::SuperOperator) = throw(bases.IncompatibleBases()) +-(a::SuperOperator{B1,B2}, b::SuperOperator{B1,B2}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} = SuperOperator{B1,B2}(a.basis_l, a.basis_r, a.data-b.data) -(a::SuperOperator) = SuperOperator(a.basis_l, a.basis_r, -a.data) --(a::SuperOperator, b::SuperOperator) = (check_samebases(a, b); SuperOperator(a.basis_l, a.basis_r, a.data-b.data)) +-(a::SuperOperator, b::SuperOperator) = throw(bases.IncompatibleBases()) """ spre(op) @@ -139,7 +143,7 @@ For operators ``A``, ``B`` the relation holds. `op` can be a dense or a sparse operator. """ -spre(op::Operator) = SuperOperator((op.basis_l, op.basis_l), (op.basis_r, op.basis_r), tensor(op, identityoperator(op)).data) +spre(op::AbstractOperator) = SuperOperator((op.basis_l, op.basis_l), (op.basis_r, op.basis_r), tensor(op, identityoperator(op)).data) """ Create a super-operator equivalent for left side operator multiplication. @@ -152,17 +156,16 @@ For operators ``A``, ``B`` the relation holds. `op` can be a dense or a sparse operator. """ -spost(op::Operator) = SuperOperator((op.basis_r, op.basis_r), (op.basis_l, op.basis_l), kron(permutedims(op.data), identityoperator(op).data)) +spost(op::AbstractOperator) = SuperOperator((op.basis_r, op.basis_r), (op.basis_l, op.basis_l), kron(permutedims(op.data), identityoperator(op).data)) -function _check_input(H::Operator, J::Vector, Jdagger::Vector, rates::Union{Vector{Float64}, Matrix{Float64}}) +function _check_input(H::AbstractOperator{B1,B2}, J::Vector, Jdagger::Vector, rates::Union{Vector{Float64}, Matrix{Float64}}) where {B1<:Basis,B2<:Basis} + # TODO: replace basis checks by dispatch; make type of J known for j=J - @assert typeof(j) <: Operator - check_samebases(H, j) + @assert isa(j, AbstractOperator{B1,B2}) end for j=Jdagger - @assert typeof(j) <: Operator - check_samebases(H, j) + @assert isa(j, AbstractOperator{B1,B2}) end @assert length(J)==length(Jdagger) if typeof(rates) == Matrix{Float64} @@ -193,7 +196,7 @@ S ρ = -\\frac{i}{ħ} [H, ρ] + \\sum_i J_i ρ J_i^† - \\frac{1}{2} J_i^† J_ """ function liouvillian(H::T, J::Vector{T}; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), - Jdagger::Vector{T}=dagger.(J)) where T<:Operator + Jdagger::Vector{T}=dagger.(J)) where T<:AbstractOperator _check_input(H, J, Jdagger, rates) L = spre(-1im*H) + spost(1im*H) if typeof(rates) == Matrix{Float64} diff --git a/src/timecorrelations.jl b/src/timecorrelations.jl index ca31897d..3723381d 100644 --- a/src/timecorrelations.jl +++ b/src/timecorrelations.jl @@ -2,7 +2,7 @@ module timecorrelations export correlation, spectrum, correlation2spectrum -using ..states, ..operators, ..operators_dense +using ..states, ..operators, ..operators_dense, ..bases using ..metrics, ..timeevolution, ..steadystate using FFTW @@ -33,11 +33,11 @@ criterion specified in [`steadystate.master`](@ref). * `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump * `kwargs...`: Further arguments are passed on to the ode solver. """ -function correlation(tspan::Vector{Float64}, rho0::DenseOperator, H::Operator, J::Vector, - op1::Operator, op2::Operator; +function correlation(tspan::Vector{Float64}, rho0::DenseOperator{B,B}, H::AbstractOperator{B,B}, J::Vector, + op1::AbstractOperator{B,B}, op2::AbstractOperator{B,B}; rates::Union{Vector{Float64}, Matrix{Float64}, Nothing}=nothing, Jdagger::Vector=dagger.(J), - kwargs...) + kwargs...) where B<:Basis function fout(t, rho) expect(op1, rho) end @@ -46,12 +46,12 @@ function correlation(tspan::Vector{Float64}, rho0::DenseOperator, H::Operator, J u end -function correlation(rho0::DenseOperator, H::Operator, J::Vector, - op1::Operator, op2::Operator; +function correlation(rho0::DenseOperator{B,B}, H::AbstractOperator{B,B}, J::Vector, + op1::AbstractOperator{B,B}, op2::AbstractOperator{B,B}; tol::Float64=1e-4, h0=10., rates::Union{Vector{Float64}, Matrix{Float64}, Nothing}=nothing, Jdagger::Vector=dagger.(J), - kwargs...) + kwargs...) where B<:Basis op2rho0 = op2*rho0 exp1 = expect(op1, op2rho0) function fout(t, rho) @@ -96,11 +96,11 @@ automatically. * `kwargs...`: Further arguments are passed on to the ode solver. """ function spectrum(omega_samplepoints::Vector{Float64}, - H::Operator, J::Vector, op::Operator; - rho0::DenseOperator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), + H::AbstractOperator{B,B}, J::Vector, op::AbstractOperator{B,B}; + rho0::DenseOperator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), tol::Float64=1e-4, - rho_ss::DenseOperator=steadystate.master(H, J; tol=tol, rho0=rho0)[end][end], - kwargs...) + rho_ss::DenseOperator{B,B}=steadystate.master(H, J; tol=tol, rho0=rho0)[end][end], + kwargs...) where B<:Basis domega = minimum(diff(omega_samplepoints)) dt = 2*pi/abs(omega_samplepoints[end] - omega_samplepoints[1]) T = 2*pi/domega @@ -110,11 +110,11 @@ function spectrum(omega_samplepoints::Vector{Float64}, return omega_samplepoints, S end -function spectrum(H::Operator, J::Vector, op::Operator; - rho0::DenseOperator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), +function spectrum(H::AbstractOperator{B,B}, J::Vector, op::AbstractOperator{B,B}; + rho0::DenseOperator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))), tol::Float64=1e-4, h0=10., - rho_ss::DenseOperator=steadystate.master(H, J; tol=tol)[end][end], - kwargs...) + rho_ss::DenseOperator{B,B}=steadystate.master(H, J; tol=tol)[end][end], + kwargs...) where B<:Basis tspan, exp_values = correlation(rho_ss, H, J, dagger(op), op, tol=tol, h0=h0, kwargs...) dtmin = minimum(diff(tspan)) T = tspan[end] - tspan[1] diff --git a/test/test_embed.jl b/test/test_embed.jl index d38b7ae8..2613aea6 100644 --- a/test/test_embed.jl +++ b/test/test_embed.jl @@ -24,7 +24,7 @@ op2 = DenseOperator(b2, b2, rand(ComplexF64, length(b2), length(b2))) op3 = DenseOperator(b3, b3, rand(ComplexF64, length(b3), length(b3))) -# Test Vector{Int}, Vector{Operator} +# Test Vector{Int}, Vector{AbstractOperator} x = embed(b, [1,2], [op1, op2]) y = op1 ⊗ op2 ⊗ I3 @test 0 ≈ abs(tracedistance_nh(x, y)) @@ -46,7 +46,7 @@ y = I1 ⊗ I2 ⊗ op3 @test 0 ≈ abs(tracedistance_nh(x, y)) -# Test Dict(Int=>Operator) +# Test Dict(Int=>AbstractOperator) x = embed(b, Dict(1 => sparse(op1), 2 => sparse(op2))) y = op1 ⊗ op2 ⊗ I3 @test 0 ≈ abs(tracedistance_nh(dense(x), y)) diff --git a/test/test_fock.jl b/test/test_fock.jl index af3c0311..977c7f2e 100644 --- a/test/test_fock.jl +++ b/test/test_fock.jl @@ -6,7 +6,7 @@ using Random, SparseArrays, LinearAlgebra Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) randstate(b) = normalize(Ket(b, rand(ComplexF64, length(b)))) randop(bl, br) = DenseOperator(bl, br, rand(ComplexF64, length(bl), length(br))) randop(b) = randop(b, b) diff --git a/test/test_manybody.jl b/test/test_manybody.jl index 7791adc2..27ad3363 100644 --- a/test/test_manybody.jl +++ b/test/test_manybody.jl @@ -6,7 +6,7 @@ using Random, SparseArrays, LinearAlgebra Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x1::StateVector, x2::StateVector) = norm(x2-x1) # Test state creation diff --git a/test/test_metrics.jl b/test/test_metrics.jl index 645bad1a..74766a5b 100644 --- a/test/test_metrics.jl +++ b/test/test_metrics.jl @@ -49,8 +49,8 @@ sigma = tensor(psi2, dagger(psi2)) @test tracedistance(sigma, sigma) ≈ 0. rho = spinup(b1) ⊗ dagger(coherentstate(b2, 0.1)) -@test_throws bases.IncompatibleBases tracedistance(rho, rho) -@test_throws bases.IncompatibleBases tracedistance_h(rho, rho) +@test_throws ArgumentError tracedistance(rho, rho) +@test_throws ArgumentError tracedistance_h(rho, rho) @test tracedistance_nh(rho, rho) ≈ 0. @@ -77,8 +77,8 @@ rho_pT2 = ptranspose(rho, 2) @test rho_pT1.data ≈ rho_pT1_an.data @test rho_pT2.data ≈ rho_pT1_an.data -@test_throws AssertionError ptranspose(e ⊗ dagger(psi1)) -@test_throws AssertionError ptranspose(dm(e)) +@test_throws MethodError ptranspose(e ⊗ dagger(psi1)) +@test_throws MethodError ptranspose(dm(e)) rho = rho ⊗ dm(g) @test PPT(rho, 1) == PPT(rho, 2) == false diff --git a/test/test_operators.jl b/test/test_operators.jl index 586f815b..e7afc32e 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -2,11 +2,11 @@ using Test using QuantumOptics using LinearAlgebra, SparseArrays, Random -mutable struct test_operators <: Operator - basis_l::Basis - basis_r::Basis +mutable struct test_operators{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR data::Matrix{ComplexF64} - test_operators(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new(b1, b2, data) : throw(DimensionMismatch()) + test_operators(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new{typeof(b1),typeof(b2)}(b1, b2, data) : throw(DimensionMismatch()) end @testset "operators" begin diff --git a/test/test_operators_dense.jl b/test/test_operators_dense.jl index 4cd5e3c0..bd8dd5b0 100644 --- a/test/test_operators_dense.jl +++ b/test/test_operators_dense.jl @@ -6,7 +6,7 @@ using Random, SparseArrays, LinearAlgebra Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x1::StateVector, x2::StateVector) = norm(x2-x1) b1a = GenericBasis(2) @@ -52,20 +52,20 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l))) xbra2 = Bra(b_l, rand(ComplexF64, length(b_l))) # Addition -@test_throws bases.IncompatibleBases op1 + dagger(op2) +@test_throws DimensionMismatch op1 + dagger(op2) @test 1e-14 > D(op1 + op_zero, op1) @test 1e-14 > D(op1 + op2, op2 + op1) @test 1e-14 > D(op1 + (op2 + op3), (op1 + op2) + op3) # Subtraction -@test_throws bases.IncompatibleBases op1 - dagger(op2) +@test_throws DimensionMismatch op1 - dagger(op2) @test 1e-14 > D(op1-op_zero, op1) @test 1e-14 > D(op1-op2, op1 + (-op2)) @test 1e-14 > D(op1-op2, op1 + (-1*op2)) @test 1e-14 > D(op1-op2-op3, op1-(op2+op3)) # Test multiplication -@test_throws bases.IncompatibleBases op1*op2 +@test_throws DimensionMismatch op1*op2 @test 1e-11 > D(op1*(x1 + 0.3*x2), op1*x1 + 0.3*op1*x2) @test 1e-11 > D((op1+op2)*(x1+0.3*x2), op1*x1 + 0.3*op1*x2 + op2*x1 + 0.3*op2*x2) diff --git a/test/test_operators_lazyproduct.jl b/test/test_operators_lazyproduct.jl index 4d43a94e..a80fcfbf 100644 --- a/test/test_operators_lazyproduct.jl +++ b/test/test_operators_lazyproduct.jl @@ -7,7 +7,7 @@ using LinearAlgebra, Random Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x1::StateVector, x2::StateVector) = norm(x2-x1) b1a = GenericBasis(2) @@ -66,7 +66,7 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l))) @test 1e-14 > D(-op1_, -op1) # Test multiplication -@test_throws bases.IncompatibleBases op1a*op1a +@test_throws DimensionMismatch op1a*op1a @test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2)) @test 1e-11 > D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_) @test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2) @@ -79,11 +79,11 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l))) # Test identityoperator Idense = identityoperator(DenseOperator, b_l) -I = identityoperator(LazyProduct, b_l) -@test isa(I, LazyProduct) -@test dense(I) == Idense -@test 1e-11 > D(I*x1, x1) -@test 1e-11 > D(xbra1*I, xbra1) +id = identityoperator(LazyProduct, b_l) +@test isa(id, LazyProduct) +@test dense(id) == Idense +@test 1e-11 > D(id*x1, x1) +@test 1e-11 > D(xbra1*id, xbra1) # Test tr and normalize op1 = randoperator(b_l) diff --git a/test/test_operators_lazysum.jl b/test/test_operators_lazysum.jl index 10a270b8..e4dc9117 100644 --- a/test/test_operators_lazysum.jl +++ b/test/test_operators_lazysum.jl @@ -7,7 +7,7 @@ using LinearAlgebra, Random Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x1::StateVector, x2::StateVector) = norm(x2-x1) b1a = GenericBasis(2) @@ -21,10 +21,10 @@ b_l = b1a⊗b2a⊗b3a b_r = b1b⊗b2b⊗b3b # Test creation -@test_throws AssertionError LazySum() +@test_throws ArgumentError LazySum() @test_throws AssertionError LazySum([1., 2.], [randoperator(b_l)]) -@test_throws AssertionError LazySum(randoperator(b_l, b_r), sparse(randoperator(b_l, b_l))) -@test_throws AssertionError LazySum(randoperator(b_l, b_r), sparse(randoperator(b_r, b_r))) +@test_throws bases.IncompatibleBases LazySum(randoperator(b_l, b_r), sparse(randoperator(b_l, b_l))) +@test_throws bases.IncompatibleBases LazySum(randoperator(b_l, b_r), sparse(randoperator(b_r, b_r))) # Test copy op1 = 2*LazySum(randoperator(b_l, b_r), sparse(randoperator(b_l, b_r))) @@ -85,16 +85,16 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l))) # Test identityoperator Idense = identityoperator(DenseOperator, b_r) -I = identityoperator(LazySum, b_r) -@test isa(I, LazySum) -@test dense(I) == Idense -@test 1e-11 > D(I*x1, x1) +id = identityoperator(LazySum, b_r) +@test isa(id, LazySum) +@test dense(id) == Idense +@test 1e-11 > D(id*x1, x1) Idense = identityoperator(DenseOperator, b_l) -I = identityoperator(LazySum, b_l) -@test isa(I, LazySum) -@test dense(I) == Idense -@test 1e-11 > D(xbra1*I, xbra1) +id = identityoperator(LazySum, b_l) +@test isa(id, LazySum) +@test dense(id) == Idense +@test 1e-11 > D(xbra1*id, xbra1) # Test tr and normalize op1 = randoperator(b_l) diff --git a/test/test_operators_lazytensor.jl b/test/test_operators_lazytensor.jl index 1d36abc3..6119218f 100644 --- a/test/test_operators_lazytensor.jl +++ b/test/test_operators_lazytensor.jl @@ -2,18 +2,18 @@ using Test using QuantumOptics using LinearAlgebra, SparseArrays, Random -mutable struct test_lazytensor <: Operator - basis_l::Basis - basis_r::Basis +mutable struct test_lazytensor{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} + basis_l::BL + basis_r::BR data::Matrix{ComplexF64} - test_lazytensor(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new(b1, b2, data) : throw(DimensionMismatch()) + test_lazytensor(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new{typeof(b1),typeof(b2)}(b1, b2, data) : throw(DimensionMismatch()) end @testset "operators-lazytensor" begin Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x1::StateVector, x2::StateVector) = norm(x2-x1) b1a = GenericBasis(2) @@ -91,7 +91,7 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l))) @test 1e-14 > D(-op1_, -op1) # Test multiplication -@test_throws bases.IncompatibleBases op1*op2 +@test_throws DimensionMismatch op1*op2 @test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2)) @test 1e-11 > D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_) @test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2) @@ -106,18 +106,18 @@ xbra2 = Bra(b_l, rand(ComplexF64, length(b_l))) # Test identityoperator Idense = identityoperator(DenseOperator, b_r) -I = identityoperator(LazyTensor, b_r) -@test isa(I, LazyTensor) -@test dense(I) == Idense -@test 1e-11 > D(I*x1, x1) -@test I == identityoperator(LazyTensor, b1b) ⊗ identityoperator(LazyTensor, b2b) ⊗ identityoperator(LazyTensor, b3b) +id = identityoperator(LazyTensor, b_r) +@test isa(id, LazyTensor) +@test dense(id) == Idense +@test 1e-11 > D(id*x1, x1) +@test id == identityoperator(LazyTensor, b1b) ⊗ identityoperator(LazyTensor, b2b) ⊗ identityoperator(LazyTensor, b3b) Idense = identityoperator(DenseOperator, b_l) -I = identityoperator(LazyTensor, b_l) -@test isa(I, LazyTensor) -@test dense(I) == Idense -@test 1e-11 > D(xbra1*I, xbra1) -@test I == identityoperator(LazyTensor, b1a) ⊗ identityoperator(LazyTensor, b2a) ⊗ identityoperator(LazyTensor, b3a) +id = identityoperator(LazyTensor, b_l) +@test isa(id, LazyTensor) +@test dense(id) == Idense +@test 1e-11 > D(xbra1*id, xbra1) +@test id == identityoperator(LazyTensor, b1a) ⊗ identityoperator(LazyTensor, b2a) ⊗ identityoperator(LazyTensor, b3a) # Test tr and normalize diff --git a/test/test_operators_sparse.jl b/test/test_operators_sparse.jl index a7beda85..782dcd55 100644 --- a/test/test_operators_sparse.jl +++ b/test/test_operators_sparse.jl @@ -4,14 +4,14 @@ using Random, SparseArrays, LinearAlgebra # Custom operator type for testing error msg -mutable struct TestOperator <: Operator; end +mutable struct TestOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR}; end @testset "operators-sparse" begin Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x1::StateVector, x2::StateVector) = norm(x2-x1) sprandop(b1, b2) = sparse(randoperator(b1, b2)) sprandop(b) = sprandop(b, b) @@ -57,13 +57,13 @@ xbra1 = dagger(randstate(b_l)) xbra2 = dagger(randstate(b_l)) # Addition -@test_throws bases.IncompatibleBases op1 + dagger(op2) +@test_throws DimensionMismatch op1 + dagger(op2) @test 1e-14 > D(op1+op2, op1_+op2_) @test 1e-14 > D(op1+op2, op1+op2_) @test 1e-14 > D(op1+op2, op1_+op2) # Subtraction -@test_throws bases.IncompatibleBases op1 - dagger(op2) +@test_throws DimensionMismatch op1 - dagger(op2) @test 1e-14 > D(op1-op2, op1_-op2_) @test 1e-14 > D(op1-op2, op1-op2_) @test 1e-14 > D(op1-op2, op1_-op2) @@ -71,7 +71,7 @@ xbra2 = dagger(randstate(b_l)) @test 1e-14 > D(op1+(-1*op2), op1_ - op2_) # Test multiplication -@test_throws bases.IncompatibleBases op1*op2 +@test_throws DimensionMismatch op1*op2 @test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2)) @test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2) @test 1e-11 > D((op1+op2)*(x1+0.3*x2), (op1_+op2_)*(x1+0.3*x2)) @@ -332,7 +332,7 @@ operators.gemm!(alpha, state, op, beta, result) dat = sprandop(b1, b1).data @test SparseOperator(b1, dat) == SparseOperator(b1, Matrix{ComplexF64}(dat)) -@test_throws ArgumentError sparse(TestOperator()) +@test_throws ArgumentError sparse(TestOperator{Basis,Basis}()) @test 2*SparseOperator(b1, dat) == SparseOperator(b1, dat)*2 @test copy(op1) == deepcopy(op1) diff --git a/test/test_particle.jl b/test/test_particle.jl index 1a464834..172a53e2 100644 --- a/test/test_particle.jl +++ b/test/test_particle.jl @@ -6,7 +6,7 @@ using FFTW, LinearAlgebra, Random Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) N = 200 xmin = -32.5 @@ -218,10 +218,10 @@ operators.gemv!(Complex(1.), LazyProduct(Txp, Tpx), psi0_bx, Complex(0.), psi_) @test 1e-12 > norm(Txp*(Tpx*psi0_bx) - psi0_bx) psi_ = deepcopy(psi0_bx) -I = dense(identityoperator(basis_momentum)) -operators.gemv!(Complex(1.), LazyProduct(Txp, I, Tpx), psi0_bx, Complex(0.), psi_) +id = dense(identityoperator(basis_momentum)) +operators.gemv!(Complex(1.), LazyProduct(Txp, id, Tpx), psi0_bx, Complex(0.), psi_) @test 1e-12 > norm(psi_ - psi0_bx) -@test 1e-12 > norm(Txp*I*(Tpx*psi0_bx) - psi0_bx) +@test 1e-12 > norm(Txp*id*(Tpx*psi0_bx) - psi0_bx) # Test dense FFT operator Txp_dense = DenseOperator(Txp) diff --git a/test/test_phasespace.jl b/test/test_phasespace.jl index b2cf7e88..5c7606d3 100644 --- a/test/test_phasespace.jl +++ b/test/test_phasespace.jl @@ -6,7 +6,7 @@ using Random, LinearAlgebra Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) # Test quasi-probability functions b = FockBasis(100) diff --git a/test/test_semiclassical.jl b/test/test_semiclassical.jl index 678f4d59..4a40175f 100644 --- a/test/test_semiclassical.jl +++ b/test/test_semiclassical.jl @@ -77,7 +77,7 @@ semiclassical.recast!(x, state_) # Test master -basis = SpinBasis(1//2) +spinbasis = SpinBasis(1//2) # Random 2 level Hamiltonian a1 = 0.5 @@ -86,18 +86,18 @@ c = 1.3 d = -4.7 data = [a1 c-1im*d; c+1im*d a2] -H = DenseOperator(basis, data) +H = DenseOperator(spinbasis, data) a = (a1 + a2)/2 b = (a1 - a2)/2 r = [c d b] -sigma_r = c*sigmax(basis) + d*sigmay(basis) + b*sigmaz(basis) +sigma_r = c*sigmax(spinbasis) + d*sigmay(spinbasis) + b*sigmaz(spinbasis) -U(t) = exp(-1im*a*t)*(cos(norm(r)*t)*one(basis) - 1im*sin(norm(r)*t)*sigma_r/norm(r)) +U(t) = exp(-1im*a*t)*(cos(norm(r)*t)*one(spinbasis) - 1im*sin(norm(r)*t)*sigma_r/norm(r)) # Random initial state -psi0 = randstate(basis) +psi0 = randstate(spinbasis) T = [0:0.5:1;] @@ -109,7 +109,7 @@ function fclassical(t, quantumstate, u, du) end state0 = semiclassical.State(psi0, ComplexF64[complex(2., 3.)]) -function f(t, state::semiclassical.State{Ket}) +function f(t, state::semiclassical.State{B,T}) where {B<:Basis,T<:Ket{B}} @test 1e-5 > norm(state.quantum - U(t)*psi0) @test 1e-5 > abs(state.classical[1] - state0.classical[1]*exp(-t)) end @@ -117,7 +117,7 @@ semiclassical.schroedinger_dynamic(T, state0, fquantum_schroedinger, fclassical; tout, state_t = semiclassical.schroedinger_dynamic(T, state0, fquantum_schroedinger, fclassical) f(T[end], state_t[end]) -function f(t, state::semiclassical.State{DenseOperator}) +function f(t, state::semiclassical.State{B,T}) where {B<:Basis,T<:DenseOperator{B,B}} @test 1e-5 > tracedistance(state.quantum, dm(U(t)*psi0)) @test 1e-5 > abs(state.classical[1] - state0.classical[1]*exp(-t)) end diff --git a/test/test_spectralanalysis.jl b/test/test_spectralanalysis.jl index e40b7dea..3f287585 100644 --- a/test/test_spectralanalysis.jl +++ b/test/test_spectralanalysis.jl @@ -2,7 +2,7 @@ using Test using QuantumOptics using LinearAlgebra, SparseArrays, Random -mutable struct SpectralanalysisTestOperator <: Operator end +mutable struct SpectralanalysisTestOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} end @testset "spectralanalysis" begin @@ -11,9 +11,9 @@ Random.seed!(0) sprandop(b) = sparse(DenseOperator(b, rand(ComplexF64, length(b), length(b)))) # Test diagonalization -@test_throws ArgumentError eigenstates(SpectralanalysisTestOperator()) +@test_throws ArgumentError eigenstates(SpectralanalysisTestOperator{Basis,Basis}()) @test_throws bases.IncompatibleBases eigenstates(DenseOperator(GenericBasis(3), GenericBasis(4))) -@test_throws ArgumentError eigenenergies(SpectralanalysisTestOperator()) +@test_throws ArgumentError eigenenergies(SpectralanalysisTestOperator{Basis,Basis}()) @test_throws bases.IncompatibleBases eigenenergies(DenseOperator(GenericBasis(3), GenericBasis(4))) diff --git a/test/test_spin.jl b/test/test_spin.jl index 53becbaa..2408f8a6 100644 --- a/test/test_spin.jl +++ b/test/test_spin.jl @@ -4,7 +4,7 @@ using LinearAlgebra @testset "spin" begin -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) # Test creation @test_throws AssertionError SpinBasis(1//3) diff --git a/test/test_states.jl b/test/test_states.jl index 31c00f2c..ab37ebf5 100644 --- a/test/test_states.jl +++ b/test/test_states.jl @@ -49,15 +49,15 @@ ket_b2 = randstate(b2) ket_b3 = randstate(b3) # Addition -@test_throws bases.IncompatibleBases bra_b1 + bra_b2 -@test_throws bases.IncompatibleBases ket_b1 + ket_b2 +@test_throws DimensionMismatch bra_b1 + bra_b2 +@test_throws DimensionMismatch ket_b1 + ket_b2 @test 1e-14 > D(bra_b1 + Bra(b1), bra_b1) @test 1e-14 > D(ket_b1 + Ket(b1), ket_b1) @test 1e-14 > D(bra_b1 + dagger(ket_b1), dagger(ket_b1) + bra_b1) # Subtraction -@test_throws bases.IncompatibleBases bra_b1 - bra_b2 -@test_throws bases.IncompatibleBases ket_b1 - ket_b2 +@test_throws DimensionMismatch bra_b1 - bra_b2 +@test_throws DimensionMismatch ket_b1 - ket_b2 @test 1e-14 > D(bra_b1 - Bra(b1), bra_b1) @test 1e-14 > D(ket_b1 - Ket(b1), ket_b1) @test 1e-14 > D(bra_b1 - dagger(ket_b1), -dagger(ket_b1) + bra_b1) @@ -86,9 +86,9 @@ idx = LinearIndices(shape)[1, 4, 3] # Norm -basis = FockBasis(1) -bra = Bra(basis, [3im, -4]) -ket = Ket(basis, [-4im, 3]) +bf = FockBasis(1) +bra = Bra(bf, [3im, -4]) +ket = Ket(bf, [-4im, 3]) @test 5 ≈ norm(bra) @test 5 ≈ norm(ket) @@ -122,7 +122,7 @@ x2 = basisstate(b2, 1) # Test permutating systems b1 = GenericBasis(2) b2 = GenericBasis(5) -b3 = GenericBasis(3) +b3 = FockBasis(3) psi1 = randstate(b1) psi2 = randstate(b2) diff --git a/test/test_superoperators.jl b/test/test_superoperators.jl index 58d756d0..e4530678 100644 --- a/test/test_superoperators.jl +++ b/test/test_superoperators.jl @@ -172,9 +172,6 @@ tout, ρt = timeevolution.master([0.,1.], ρ₀, H, J; reltol=1e-7) @test dense(spre(op1)) == spre(op1) -@test_throws bases.IncompatibleBases L*op1 -@test_throws bases.IncompatibleBases L*spre(sm) - @test L/2.0 == 0.5*L == L*0.5 @test -L == SparseSuperOperator(L.basis_l, L.basis_r, -L.data) diff --git a/test/test_transformations.jl b/test/test_transformations.jl index c79066ef..71f7f7ee 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -6,7 +6,7 @@ using Random, LinearAlgebra Random.seed!(0) -D(op1::Operator, op2::Operator) = abs(tracedistance_nh(dense(op1), dense(op2))) +D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) D(x::Ket, y::Ket) = norm(x-y) # Test transformations