diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2660f07..709ed21 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -10,9 +10,9 @@ jobs: fail-fast: false matrix: version: - - '1.3' - '1.4' - '1.5' + - '1.6' os: - ubuntu-latest - macOS-latest diff --git a/.github/workflows/deploy-nightly.yml b/.github/workflows/deploy-nightly.yml index 51db94c..e4f27c0 100644 --- a/.github/workflows/deploy-nightly.yml +++ b/.github/workflows/deploy-nightly.yml @@ -14,6 +14,7 @@ jobs: - '1.3' - '1.4' - '1.5' + - '1.6' - 'nightly' os: - ubuntu-latest diff --git a/README.md b/README.md index 41cbdb6..c53681f 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,7 @@ # OpenQuantumSystems.jl -| Service | Master | Develop | -| :------------ | :------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | :--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| Stable | [![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI/badge.svg?branch=master)](https://github.com/detrin/OpenQuantumSystems.jl/actions) | [![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI/badge.svg?branch=devel)](https://github.com/detrin/OpenQuantumSystems.jl/actions) | -| Nightly Julia | [![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI-nightly-julia/badge.svg?branch=master)](https://github.com/detrin/OpenQuantumSystems.jl/actions) | [![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI-nightly-julia/badge.svg?branch=devel)](https://github.com/detrin/OpenQuantumSystems.jl/actions) | -| Nigtly deploy | [![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/Deploy%20Nightly/badge.svg?branch=master)](https://github.com/detrin/OpenQuantumSystems.jl/actions) | [![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/Deploy%20Nightly/badge.svg?branch=devel)](https://github.com/detrin/OpenQuantumSystems.jl/actions) | -| coveralls | [![Coverage Status](https://coveralls.io/repos/detrin/OpenQuantumSystems.jl/badge.svg?branch=master)](https://coveralls.io/r/detrin/OpenQuantumSystems.jl?branch=master) | [![Coverage Status](https://coveralls.io/repos/github/detrin/OpenQuantumSystems.jl/badge.svg?branch=devel)](https://coveralls.io/github/detrin/OpenQuantumSystems.jl?branch=devel) | -| codecov | [![codecov](https://codecov.io/gh/detrin/OpenQuantumSystems.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/detrin/OpenQuantumSystems.jl) | [![codecov](https://codecov.io/gh/detrin/OpenQuantumSystems.jl/branch/devel/graph/badge.svg)](https://codecov.io/gh/detrin/OpenQuantumSystems.jl) | +[![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI/badge.svg?branch=master)](https://github.com/detrin/OpenQuantumSystems.jl/actions) +[![codecov](https://codecov.io/gh/detrin/OpenQuantumSystems.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/detrin/OpenQuantumSystems.jl) **OpenQuantumSystems.jl** is a numerical framework written in [Julia] that makes it easy to simulate various kinds of quantum systems. It is inspired by the @@ -40,3 +35,16 @@ or create a pull request. 4. - GPU support for Schrodinger equation and possibly state decomposition. - Anharmonic oscillators. - Double excited states. + +### Branches status + +Master + +[![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI-nightly-julia/badge.svg?branch=master)](https://github.com/detrin/OpenQuantumSystems.jl/actions) +[![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/Deploy%20Nightly/badge.svg?branch=master)](https://github.com/detrin/OpenQuantumSystems.jl/actions) + +Develop + +[![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI/badge.svg?branch=devel)](https://github.com/detrin/OpenQuantumSystems.jl/actions) +[![Actions Status](https://github.com/detrin/OpenQuantumSystems.jl/workflows/CI-nightly-julia/badge.svg?branch=devel)](https://github.com/detrin/OpenQuantumSystems.jl/actions) +[![codecov](https://codecov.io/gh/detrin/OpenQuantumSystems.jl/branch/devel/graph/badge.svg)](https://codecov.io/gh/detrin/OpenQuantumSystems.jl) diff --git a/src/OpenQuantumSystems.jl b/src/OpenQuantumSystems.jl index 3512a6d..df39b84 100644 --- a/src/OpenQuantumSystems.jl +++ b/src/OpenQuantumSystems.jl @@ -102,6 +102,7 @@ export bases, avg_gate_fidelity +include("core.jl") include("sortedindices.jl") include("bases.jl") include("states.jl") diff --git a/src/core.jl b/src/core.jl new file mode 100644 index 0000000..0e9b094 --- /dev/null +++ b/src/core.jl @@ -0,0 +1,2 @@ + +const ComputableType = Union{AbstractFloat, Complex} diff --git a/src/operators.jl b/src/operators.jl index 406617b..e805329 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -25,14 +25,11 @@ abstract type DataOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} end # Common error messages -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::AbstractOperator) = length(a.basis_l)::Int * length(a.basis_r)::Int +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::AbstractOperator) = length(a.basis_l)::Int*length(a.basis_r)::Int basis(a::AbstractOperator) = (check_samebases(a); a.basis_l) # Ensure scalar broadcasting @@ -48,8 +45,7 @@ Base.broadcastable(x::AbstractOperator) = Ref(x) -(a::Number, b::AbstractOperator) = addnumbererror() -(a::AbstractOperator, b::Number) = addnumbererror() -*(a::AbstractOperator, b::AbstractOperator) = - arithmetic_binary_error("Multiplication", a, b) +*(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Multiplication", a, b) ^(a::AbstractOperator, b::Int) = Base.power_by_squaring(a, b) @@ -76,8 +72,7 @@ ishermitian(op::AbstractOperator) = arithmetic_unary_error(ishermitian, op) Tensor product ``\\hat{x}⊗\\hat{y}⊗\\hat{z}⊗…`` of the given operators. """ -tensor(a::AbstractOperator, b::AbstractOperator) = - arithmetic_binary_error("Tensor product", a, b) +tensor(a::AbstractOperator, b::AbstractOperator) = arithmetic_binary_error("Tensor product", a, b) tensor(op::AbstractOperator) = op tensor(operators::AbstractOperator...) = reduce(tensor, operators) @@ -87,12 +82,8 @@ tensor(operators::AbstractOperator...) = 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, - operators::Vector{T}, -) where {T<:AbstractOperator} +function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, + indices::Vector, operators::Vector{T}) where T<:AbstractOperator @assert check_embed_indices(indices) @@ -110,10 +101,7 @@ function embed( (opsb.basis_r == basis_r.bases[idxsb]) || throw(IncompatibleBases()) end - embed_op = tensor([ - i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : - identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i = 1:N - ]...) + embed_op = tensor([i ∈ indices_sb ? ops_sb[indexin(i, indices_sb)[1]] : identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...) # Embed all joint-subspace operators. idxop_comp = [x for x in zip(indices, operators) if typeof(x[1]) <: Array] @@ -130,12 +118,8 @@ end Embed operator acting on a joint Hilbert space where missing indices are filled up with identity operators. """ -function embed( - basis_l::CompositeBasis, - basis_r::CompositeBasis, - indices::Vector{Int}, - op::T, -) where {T<:AbstractOperator} +function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, + indices::Vector{Int}, op::T) where T<:AbstractOperator N = length(basis_l.bases) @assert length(basis_r.bases) == N @@ -143,9 +127,7 @@ function embed( reduce(tensor, basis_r.bases[indices]) == op.basis_r || throw(IncompatibleBases()) index_order = [idx for idx in 1:length(basis_l.bases) if idx ∉ indices] - all_operators = AbstractOperator[ - identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in index_order - ] + all_operators = AbstractOperator[identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in index_order] for idx in indices pushfirst!(index_order, idx) @@ -177,32 +159,22 @@ function embed( expand_dims_r = dims_r[expand_order] # Prepare the permutation to the correctly ordered basis. - perm_order_l = [indexin(idx, expand_order)[1] for idx = 1:length(dims_l)] - perm_order_r = [indexin(idx, expand_order)[1] for idx = 1:length(dims_r)] + perm_order_l = [indexin(idx, expand_order)[1] for idx in 1:length(dims_l)] + perm_order_r = [indexin(idx, expand_order)[1] for idx in 1:length(dims_r)] # Perform the index expansion, the permutation, and the index collapse. - M = ( - reshape(permuted_op.data, tuple([expand_dims_l; expand_dims_r]...)) |> - x -> - permutedims(x, [perm_order_l; perm_order_r .+ length(dims_l)]) |> - x -> sparse(reshape(x, (prod(dims_l), prod(dims_r)))) - ) + M = (reshape(permuted_op.data, tuple([expand_dims_l; expand_dims_r]...)) |> + x -> permutedims(x, [perm_order_l; perm_order_r .+ length(dims_l)]) |> + x -> sparse(reshape(x, (prod(dims_l), prod(dims_r))))) unpermuted_op.data = M return unpermuted_op end -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, - operators::Vector{T}, -) where {T<:AbstractOperator} = embed(basis, basis, indices, operators) -embed(basis::CompositeBasis, indices::Vector{Int}, op::AbstractOperator) = - embed(basis, basis, indices, op) +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, operators::Vector{T}) where {T<:AbstractOperator} = embed(basis, basis, indices, operators) +embed(basis::CompositeBasis, indices::Vector{Int}, op::AbstractOperator) = embed(basis, basis, indices, op) """ embed(basis1[, basis2], operators::Dict) @@ -210,11 +182,8 @@ embed(basis::CompositeBasis, indices::Vector{Int}, op::AbstractOperator) = `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<:AbstractOperator} +function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, + 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 @@ -224,47 +193,27 @@ function embed( operator_list = [operator_list...;] indices_flat = [indices...;] start_indices_flat = [i[1] for i in indices] - complement_indices_flat = Int[i for i = 1:N if i ∉ indices_flat] + complement_indices_flat = Int[i for i=1:N if i ∉ indices_flat] operators_flat = T[] - if all([minimum(I):maximum(I);] == I for I in indices) - for i = 1:N + if all([minimum(I):maximum(I);]==I for I in indices) + for i in 1:N if i in complement_indices_flat - push!( - operators_flat, - identityoperator(T, basis_l.bases[i], basis_r.bases[i]), - ) + push!(operators_flat, identityoperator(T, basis_l.bases[i], basis_r.bases[i])) elseif i in start_indices_flat push!(operators_flat, operator_list[indexin(i, start_indices_flat)[1]]) end end return tensor(operators_flat...) else - complement_operators = [ - identityoperator(T, basis_l.bases[i], basis_r.bases[i]) - for i in complement_indices_flat - ] + complement_operators = [identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in complement_indices_flat] op = tensor([operator_list; complement_operators]...) perm = sortperm([indices_flat; complement_indices_flat]) return permutesystems(op, perm) end end -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...) +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...) """ @@ -281,15 +230,14 @@ ptrace(a::AbstractOperator, index::Vector{Int}) = arithmetic_unary_error("Partia Return the normalized operator so that its `tr(op)` is one. """ -normalize(op::AbstractOperator) = 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::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()'.")) +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) @@ -298,40 +246,27 @@ Expectation value of the given operator `op` for the specified `state`. `state` can either be a (density) operator or a ket. """ -expect(op::AbstractOperator{B,B}, state::Ket{B}) where {B<:Basis} = - state.data' * (op * state).data -expect( - op::AbstractOperator{B1,B2}, - state::AbstractOperator{B2,B2}, -) where {B1<:Basis,B2<:Basis} = tr(op * state) +expect(op::AbstractOperator{B,B}, state::Ket{B}) where B<:Basis = state.data' * (op * state).data +expect(op::AbstractOperator{B1,B2}, state::AbstractOperator{B2,B2}) where {B1<:Basis,B2<:Basis} = 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::AbstractOperator{B1,B2}, - state::AbstractOperator{B3,B3}, -) where {B1<:Basis,B2<:Basis,B3<:CompositeBasis} +function expect(indices::Vector{Int}, op::AbstractOperator{B1,B2}, state::AbstractOperator{B3,B3}) where {B1<:Basis,B2<:Basis,B3<:CompositeBasis} N = length(state.basis_l.shape) indices_ = complement(N, indices) expect(op, ptrace(state, indices_)) end -function expect( - indices::Vector{Int}, - op::AbstractOperator{B,B}, - state::Ket{B2}, -) where {B<:Basis,B2<:CompositeBasis} +function expect(indices::Vector{Int}, op::AbstractOperator{B,B}, state::Ket{B2}) where {B<:Basis,B2<:CompositeBasis} N = length(state.basis.shape) indices_ = complement(N, indices) expect(op, ptrace(state, indices_)) end expect(index::Int, op::AbstractOperator, state) = expect([index], op, state) -expect(op::AbstractOperator, states::Vector) = [expect(op, state) for state in states] -expect(indices::Vector{Int}, op::AbstractOperator, states::Vector) = - [expect(indices, op, state) for state in states] +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) @@ -340,12 +275,12 @@ Variance of the given operator `op` for the specified `state`. `state` can either be a (density) operator or a ket. """ -function variance(op::AbstractOperator{B,B}, state::Ket{B}) where {B<:Basis} - x = op * state - state.data' * (op * x).data - (state.data' * x.data)^2 +function variance(op::AbstractOperator{B,B}, state::Ket{B}) where B<:Basis + x = op*state + state.data'*(op*x).data - (state.data'*x.data)^2 end -function variance(op::AbstractOperator{B,B}, state::AbstractOperator{B,B}) where {B<:Basis} - expect(op * op, state) - expect(op, state)^2 +function variance(op::AbstractOperator{B,B}, state::AbstractOperator{B,B}) where B<:Basis + expect(op*op, state) - expect(op, state)^2 end """ @@ -353,28 +288,19 @@ 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::AbstractOperator{B,B}, - state::AbstractOperator{BC,BC}, -) where {B<:Basis,BC<:CompositeBasis} +function variance(indices::Vector{Int}, op::AbstractOperator{B,B}, state::AbstractOperator{BC,BC}) where {B<:Basis,BC<:CompositeBasis} N = length(state.basis_l.shape) indices_ = complement(N, indices) variance(op, ptrace(state, indices_)) end -function variance( - indices::Vector{Int}, - op::AbstractOperator{B,B}, - state::Ket{BC}, -) where {B<:Basis,BC<:CompositeBasis} +function variance(indices::Vector{Int}, op::AbstractOperator{B,B}, state::Ket{BC}) where {B<:Basis,BC<:CompositeBasis} N = length(state.basis.shape) indices_ = complement(N, indices) variance(op, ptrace(state, indices_)) end variance(index::Int, op::AbstractOperator, state) = variance([index], op, state) -variance(op::AbstractOperator, states::Vector) = [variance(op, state) for state in states] -variance(indices::Vector{Int}, op::AbstractOperator, states::Vector) = - [variance(indices, op, state) for state in states] +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] """ @@ -382,23 +308,18 @@ variance(indices::Vector{Int}, op::AbstractOperator, states::Vector) = Operator exponential. """ -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().")) +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::AbstractOperator, 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<: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) +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::AbstractOperator) = identityoperator(op) @@ -416,7 +337,7 @@ function check_ptrace_arguments(a::AbstractOperator, indices::Vector{Int}) throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use tr() instead.")) end check_indices(length(a.basis_l.shape), indices) - for i in indices + for i=indices if a.basis_l.shape[i] != a.basis_r.shape[i] throw(ArgumentError("Partial trace can only be applied onto subsystems that have the same left and right dimension.")) end @@ -430,14 +351,12 @@ function check_ptrace_arguments(a::StateVector, indices::Vector{Int}) end 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 +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::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) +multiplicable(a::AbstractOperator, b::AbstractOperator) = multiplicable(a.basis_r, b.basis_l) -Base.size(op::AbstractOperator) = prod(length(op.basis_l), length(op.basis_r)) -Base.size(op::AbstractOperator, i::Int) = (i == 1 ? length(op.basis_l) : length(op.basis_r)) +Base.size(op::AbstractOperator) = prod(length(op.basis_l),length(op.basis_r)) +Base.size(op::AbstractOperator, i::Int) = (i==1 ? length(op.basis_l) : length(op.basis_r)) diff --git a/src/operators_dense.jl b/src/operators_dense.jl index d62ccd1..ba9663d 100644 --- a/src/operators_dense.jl +++ b/src/operators_dense.jl @@ -12,30 +12,22 @@ mutable struct Operator{BL<:Basis,BR<:Basis,T} <: DataOperator{BL,BR} basis_l::BL basis_r::BR data::T - function Operator{BL,BR,T}( - basis_l::BL, - basis_r::BR, - data::T, - ) where {BL<:Basis,BR<:Basis,T} - (length.((basis_l, basis_r)) == size(data)) || - throw(DimensionMismatch("Tried to assign data of size $(size(data)) to bases of length $(length(basis_l)) and $(length(basis_r))!")) - new(basis_l, basis_r, data) + function Operator{BL,BR,T}(basis_l::BL,basis_r::BR,data::T) where {BL<:Basis,BR<:Basis,T} + (length.((basis_l,basis_r))==size(data)) || throw(DimensionMismatch("Tried to assign data of size $(size(data)) to bases of length $(length(basis_l)) and $(length(basis_r))!")) + new(basis_l,basis_r,data) end end -Operator{BL,BR}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T} = - Operator{BL,BR,T}(basis_l, basis_r, data) -Operator(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T} = - Operator{BL,BR,T}(basis_l, basis_r, data) -Operator(b::Basis, data) = Operator(b, b, data) +Operator{BL,BR}(basis_l::BL,basis_r::BR,data::T) where {BL,BR,T} = Operator{BL,BR,T}(basis_l,basis_r,data) +Operator(basis_l::BL,basis_r::BR,data::T) where {BL,BR,T} = Operator{BL,BR,T}(basis_l,basis_r,data) +Operator(b::Basis,data) = Operator(b,b,data) -Base.zero(op::Operator) = Operator(op.basis_l, op.basis_r, zero(op.data)) +Base.zero(op::Operator) = Operator(op.basis_l,op.basis_r,zero(op.data)) Base.eltype(op::Operator) = eltype(op.data) Base.size(op::Operator) = size(op.data) Base.size(op::Operator, d::Int) = size(op.data, d) # Convert data to CuArray with cu(::Operator) -Adapt.adapt_structure(to, x::Operator) = - Operator(x.basis_l, x.basis_r, Adapt.adapt(to, x.data)) +Adapt.adapt_structure(to, x::Operator) = Operator(x.basis_l, x.basis_r, Adapt.adapt(to, x.data)) const DenseOpPureType{BL<:Basis,BR<:Basis} = Operator{BL,BR,<:Matrix} const DenseOpAdjType{BL<:Basis,BR<:Basis} = Operator{BL,BR,<:Adjoint{<:Number,<:Matrix}} @@ -45,18 +37,14 @@ const AdjointOperator{BL<:Basis,BR<:Basis} = Operator{BL,BR,<:Adjoint} """ DenseOperator(b1[, b2, data]) -Dense array implementation of Operator. Converts any given data to `Matrix{ComplexF64}`. +Dense array implementation of Operator. Converts any given data to a dense `Matrix`. """ -DenseOperator(basis_l::Basis, basis_r::Basis, data::T) where {T} = - Operator(basis_l, basis_r, Matrix{ComplexF64}(data)) -DenseOperator(basis_l::Basis, basis_r::Basis, data::Matrix{ComplexF64}) = - Operator(basis_l, basis_r, data) +DenseOperator(basis_l::Basis,basis_r::Basis,data::T) where T = Operator(basis_l,basis_r,Matrix(data)) +DenseOperator(basis_l::Basis,basis_r::Basis,data::Matrix) = Operator(basis_l,basis_r,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::Basis, b2::Basis) = DenseOperator(b1, b2, zeros(ComplexF64, length(b1), length(b2))) DenseOperator(b::Basis) = DenseOperator(b, b) -DenseOperator(op::DataOperator) = - DenseOperator(op.basis_l, op.basis_r, Matrix{ComplexF64}(op.data)) +DenseOperator(op::DataOperator) = DenseOperator(op.basis_l,op.basis_r,Matrix(op.data)) Base.copy(x::Operator) = Operator(x.basis_l, x.basis_r, copy(x.data)) @@ -67,85 +55,57 @@ Convert an arbitrary Operator into a [`DenseOperator`](@ref). """ dense(x::AbstractOperator) = DenseOperator(x) -==(x::DataOperator{BL,BR}, y::DataOperator{BL,BR}) where {BL<:Basis,BR<:Basis} = - (samebases(x, y) && x.data == y.data) +==(x::DataOperator{BL,BR}, y::DataOperator{BL,BR}) where {BL<:Basis,BR<:Basis} = (samebases(x,y) && x.data==y.data) ==(x::DataOperator, y::DataOperator) = false -Base.isapprox( - x::DataOperator{BL,BR}, - y::DataOperator{BL,BR}; - kwargs..., -) where {BL<:Basis,BR<:Basis} = (samebases(x, y) && isapprox(x.data, y.data; kwargs...)) +Base.isapprox(x::DataOperator{BL,BR}, y::DataOperator{BL,BR}; kwargs...) where {BL<:Basis,BR<:Basis} = (samebases(x,y) && isapprox(x.data, y.data; kwargs...)) Base.isapprox(x::DataOperator, y::DataOperator; kwargs...) = false # Arithmetic operations -+(a::Operator{BL,BR}, b::Operator{BL,BR}) where {BL<:Basis,BR<:Basis} = - Operator(a.basis_l, a.basis_r, a.data + b.data) ++(a::Operator{BL,BR}, b::Operator{BL,BR}) where {BL<:Basis,BR<:Basis} = Operator(a.basis_l, a.basis_r, a.data+b.data) +(a::Operator, b::Operator) = throw(IncompatibleBases()) -(a::Operator) = Operator(a.basis_l, a.basis_r, -a.data) --(a::Operator{BL,BR}, b::Operator{BL,BR}) where {BL<:Basis,BR<:Basis} = - Operator(a.basis_l, a.basis_r, a.data - b.data) +-(a::Operator{BL,BR}, b::Operator{BL,BR}) where {BL<:Basis,BR<:Basis} = Operator(a.basis_l, a.basis_r, a.data-b.data) -(a::Operator, b::Operator) = throw(IncompatibleBases()) -*(a::Operator{BL,BR}, b::Ket{BR}) where {BL<:Basis,BR<:Basis} = - Ket{BL}(a.basis_l, a.data * b.data) +*(a::Operator{BL,BR}, b::Ket{BR}) where {BL<:Basis,BR<:Basis} = Ket{BL}(a.basis_l, a.data*b.data) *(a::DataOperator, b::Ket) = throw(IncompatibleBases()) -*(a::Bra{BL}, b::Operator{BL,BR}) where {BL<:Basis,BR<:Basis} = - Bra{BR}(b.basis_r, transpose(b.data) * a.data) +*(a::Bra{BL}, b::Operator{BL,BR}) where {BL<:Basis,BR<:Basis} = Bra{BR}(b.basis_r, transpose(b.data)*a.data) *(a::Bra, b::DataOperator) = throw(IncompatibleBases()) -*(a::Operator{B1,B2}, b::Operator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = - Operator(a.basis_l, b.basis_r, a.data * b.data) +*(a::Operator{B1,B2}, b::Operator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = Operator(a.basis_l, b.basis_r, a.data*b.data) *(a::DataOperator, b::DataOperator) = throw(IncompatibleBases()) -*(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, b * a.data) -*(a::Number, b::Operator) = Operator(b.basis_l, b.basis_r, a * b.data) -function *( - op1::AbstractOperator{B1,B2}, - op2::Operator{B2,B3,T}, -) where {B1<:Basis,B2<:Basis,B3<:Basis,T} - result = Operator{B1,B3,T}( - op1.basis_l, - op2.basis_r, - similar(op2.data, length(op1.basis_l), length(op2.basis_r)), - ) - mul!(result, op1, op2) +*(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, b*a.data) +*(a::Number, b::Operator) = Operator(b.basis_l, b.basis_r, a*b.data) +function *(op1::AbstractOperator{B1,B2}, op2::Operator{B2,B3,T}) where {B1<:Basis,B2<:Basis,B3<:Basis,T} + result = Operator{B1,B3,T}(op1.basis_l, op2.basis_r, similar(op2.data,length(op1.basis_l),length(op2.basis_r))) + mul!(result,op1,op2) return result end -function *( - op1::Operator{B1,B2,T}, - op2::AbstractOperator{B2,B3}, -) where {B1<:Basis,B2<:Basis,B3<:Basis,T} - result = Operator{B1,B3,T}( - op1.basis_l, - op2.basis_r, - similar(op1.data, length(op1.basis_l), length(op2.basis_r)), - ) - mul!(result, op1, op2) +function *(op1::Operator{B1,B2,T}, op2::AbstractOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis,T} + result = Operator{B1,B3,T}(op1.basis_l, op2.basis_r, similar(op1.data,length(op1.basis_l),length(op2.basis_r))) + mul!(result,op1,op2) return result end function *(op::AbstractOperator{BL,BR}, psi::Ket{BR,T}) where {BL<:Basis,BR<:Basis,T} - result = Ket{BL,T}(op.basis_l, similar(psi.data, length(op.basis_l))) - mul!(result, op, psi) + result = Ket{BL,T}(op.basis_l,similar(psi.data,length(op.basis_l))) + mul!(result,op,psi) return result end function *(psi::Bra{BL,T}, op::AbstractOperator{BL,BR}) where {BL<:Basis,BR<:Basis,T} - result = Bra{BR,T}(op.basis_r, similar(psi.data, length(op.basis_r))) - mul!(result, psi, op) + result = Bra{BR,T}(op.basis_r, similar(psi.data,length(op.basis_r))) + mul!(result,psi,op) return result end /(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, a.data ./ b) -dagger(x::Operator) = Operator(x.basis_r, x.basis_l, adjoint(x.data)) -transpose(x::Operator) = Operator(x.basis_r, x.basis_l, transpose(x.data)) +dagger(x::Operator) = Operator(x.basis_r,x.basis_l,adjoint(x.data)) +transpose(x::Operator) = Operator(x.basis_r,x.basis_l,transpose(x.data)) ishermitian(A::DataOperator) = false -ishermitian(A::DataOperator{B,B}) where {B<:Basis} = ishermitian(A.data) +ishermitian(A::DataOperator{B,B}) where B<:Basis = ishermitian(A.data) Base.collect(A::Operator) = Operator(A.basis_l, A.basis_r, collect(A.data)) -tensor(a::Operator, b::Operator) = Operator( - tensor(a.basis_l, b.basis_l), - tensor(a.basis_r, b.basis_r), - kron(b.data, a.data), -) +tensor(a::Operator, b::Operator) = Operator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data)) conj(a::Operator) = Operator(a.basis_l, a.basis_r, conj(a.data)) conj!(a::Operator) = (conj!(a.data); a) @@ -156,13 +116,9 @@ conj!(a::Operator) = (conj!(a.data); a) Outer product ``|x⟩⟨y|`` of the given states. """ -tensor(a::Ket, b::Bra) = Operator( - a.basis, - b.basis, - reshape(kron(b.data, a.data), length(a.basis), length(b.basis)), -) +tensor(a::Ket, b::Bra) = Operator(a.basis, b.basis, reshape(kron(b.data, a.data), length(a.basis), length(b.basis))) -tr(op::Operator{B,B}) where {B<:Basis} = tr(op.data) +tr(op::Operator{B,B}) where B<:Basis = tr(op.data) function ptrace(a::DataOperator, indices::Vector{Int}) check_ptrace_arguments(a, indices) @@ -189,19 +145,17 @@ function ptrace(psi::Bra, indices::Vector{Int}) return Operator(b_, b_, result) end -normalize!(op::Operator) = (rmul!(op.data, 1.0 / tr(op)); op) +normalize!(op::Operator) = (rmul!(op.data, 1.0/tr(op)); op) -function expect(op::DataOperator{B,B}, state::Ket{B}) where {B<:Basis} +function expect(op::DataOperator{B,B}, state::Ket{B}) where B<:Basis state.data' * op.data * state.data end -function expect( - op::DataOperator{B1,B2}, - state::DataOperator{B2,B2}, -) where {B1<:Basis,B2<:Basis} - result = zero(promote_type(eltype(op), eltype(state))) - @inbounds for i = 1:size(op.data, 1), j = 1:size(op.data, 2) - result += op.data[i, j] * state.data[j, i] +function expect(op::DataOperator{B1,B2}, state::DataOperator{B2,B2}) where {B1<:Basis,B2<:Basis} + check_samebases(op, state) + result = zero(promote_type(eltype(op),eltype(state))) + @inbounds for i=1:size(op.data, 1), j=1:size(op.data,2) + result += op.data[i,j]*state.data[j,i] end result end @@ -210,10 +164,7 @@ function exp(op::T) where {B<:Basis,T<:DenseOpType{B,B}} return DenseOperator(op.basis_l, op.basis_r, exp(op.data)) end -function permutesystems( - a::Operator{B1,B2}, - perm::Vector{Int}, -) where {B1<:CompositeBasis,B2<:CompositeBasis} +function permutesystems(a::Operator{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]...) @@ -221,13 +172,9 @@ function permutesystems( data = reshape(data, length(a.basis_l), length(a.basis_r)) return Operator(permutesystems(a.basis_l, perm), permutesystems(a.basis_r, perm), data) end -permutesystems( - a::AdjointOperator{B1,B2}, - perm::Vector{Int}, -) where {B1<:CompositeBasis,B2<:CompositeBasis} = dagger(permutesystems(dagger(a), perm)) +permutesystems(a::AdjointOperator{B1,B2}, perm::Vector{Int}) where {B1<:CompositeBasis,B2<:CompositeBasis} = dagger(permutesystems(dagger(a),perm)) -identityoperator(::Type{T}, b1::Basis, b2::Basis) where {BL,BR,dType,T<:DenseOpType} = - Operator(b1, b2, Matrix{ComplexF64}(I, length(b1), length(b2))) +identityoperator(::Type{T}, b1::Basis, b2::Basis) where {BL,BR,dType,T<:DenseOpType} = Operator(b1, b2, Matrix{ComplexF64}(I, length(b1), length(b2))) """ projector(a::Ket, b::Bra) @@ -240,7 +187,7 @@ projector(a::Ket, b::Bra) = tensor(a, b) Projection operator ``|a⟩⟨a|``. """ -projector(a::Ket) = Operator(a.basis, a.data * a.data') +projector(a::Ket) = Operator(a.basis, a.data*a.data') """ projector(a::Bra) @@ -262,20 +209,16 @@ function _strides(shape::Vector{Int}) N = length(shape) S = zeros(Int, N) S[1] = 1 - for m = 2:N - S[m] = S[m-1] * shape[m-1] + for m=2:N + S[m] = S[m-1]*shape[m-1] end return S end # Dense operator version -@generated function _ptrace( - ::Type{Val{RANK}}, - a, - shape_l::Vector{Int}, - shape_r::Vector{Int}, - indices::Vector{Int}, -) where {RANK} +@generated function _ptrace(::Type{Val{RANK}}, a, + shape_l::Vector{Int}, shape_r::Vector{Int}, + indices::Vector{Int}) where RANK return quote a_strides_l = _strides(shape_l) result_shape_l = copy(shape_l) @@ -288,23 +231,10 @@ end N_result_l = prod(result_shape_l) N_result_r = prod(result_shape_r) result = zeros(eltype(a), N_result_l, N_result_r) - @nexprs 1 (d -> (Jr_{$RANK} = 1; Ir_{$RANK} = 1)) - @nloops $RANK ir (d -> 1:shape_r[d]) (d -> (Ir_{d - 1} = Ir_d; Jr_{d - 1} = Jr_d)) ( - d -> (Ir_d += a_strides_r[d]; if !(d in indices) - Jr_d += result_strides_r[d] - end) - ) begin - @nexprs 1 (d -> (Jl_{$RANK} = 1; Il_{$RANK} = 1)) - @nloops $RANK il (k -> 1:shape_l[k]) ( - k -> ( - Il_{k - 1} = Il_k; Jl_{k - 1} = Jl_k; if (k in indices && il_k != ir_k) - Il_k += a_strides_l[k] - continue - end - ) - ) (k -> (Il_k += a_strides_l[k]; if !(k in indices) - Jl_k += result_strides_l[k] - end)) begin + @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1)) + @nloops $RANK ir (d->1:shape_r[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides_r[d]; if !(d in indices) Jr_d+=result_strides_r[d] end)) begin + @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1)) + @nloops $RANK il (k->1:shape_l[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides_l[k]; continue end)) (k->(Il_k+=a_strides_l[k]; if !(k in indices) Jl_k+=result_strides_l[k] end)) begin result[Jl_0, Jr_0] += a[Il_0, Ir_0] end end @@ -312,12 +242,8 @@ end end end -@generated function _ptrace_ket( - ::Type{Val{RANK}}, - a::Vector, - shape::Vector{Int}, - indices::Vector{Int}, -) where {RANK} +@generated function _ptrace_ket(::Type{Val{RANK}}, a::Vector, + shape::Vector{Int}, indices::Vector{Int}) where RANK return quote a_strides = _strides(shape) result_shape = copy(shape) @@ -325,36 +251,19 @@ end result_strides = _strides(result_shape) N_result = prod(result_shape) result = zeros(eltype(a), N_result, N_result) - @nexprs 1 (d -> (Jr_{$RANK} = 1; Ir_{$RANK} = 1)) - @nloops $RANK ir (d -> 1:shape[d]) (d -> (Ir_{d - 1} = Ir_d; Jr_{d - 1} = Jr_d)) ( - d -> (Ir_d += a_strides[d]; if !(d in indices) - Jr_d += result_strides[d] - end) - ) begin - @nexprs 1 (d -> (Jl_{$RANK} = 1; Il_{$RANK} = 1)) - @nloops $RANK il (k -> 1:shape[k]) ( - k -> ( - Il_{k - 1} = Il_k; Jl_{k - 1} = Jl_k; if (k in indices && il_k != ir_k) - Il_k += a_strides[k] - continue - end - ) - ) (k -> (Il_k += a_strides[k]; if !(k in indices) - Jl_k += result_strides[k] - end)) begin - result[Jl_0, Jr_0] += a[Il_0] * conj(a[Ir_0]) + @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1)) + @nloops $RANK ir (d->1:shape[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides[d]; if !(d in indices) Jr_d+=result_strides[d] end)) begin + @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1)) + @nloops $RANK il (k->1:shape[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides[k]; continue end)) (k->(Il_k+=a_strides[k]; if !(k in indices) Jl_k+=result_strides[k] end)) begin + result[Jl_0, Jr_0] += a[Il_0]*conj(a[Ir_0]) end end return result end end -@generated function _ptrace_bra( - ::Type{Val{RANK}}, - a::Vector, - shape::Vector{Int}, - indices::Vector{Int}, -) where {RANK} +@generated function _ptrace_bra(::Type{Val{RANK}}, a::Vector, + shape::Vector{Int}, indices::Vector{Int}) where RANK return quote a_strides = _strides(shape) result_shape = copy(shape) @@ -362,24 +271,11 @@ end result_strides = _strides(result_shape) N_result = prod(result_shape) result = zeros(eltype(a), N_result, N_result) - @nexprs 1 (d -> (Jr_{$RANK} = 1; Ir_{$RANK} = 1)) - @nloops $RANK ir (d -> 1:shape[d]) (d -> (Ir_{d - 1} = Ir_d; Jr_{d - 1} = Jr_d)) ( - d -> (Ir_d += a_strides[d]; if !(d in indices) - Jr_d += result_strides[d] - end) - ) begin - @nexprs 1 (d -> (Jl_{$RANK} = 1; Il_{$RANK} = 1)) - @nloops $RANK il (k -> 1:shape[k]) ( - k -> ( - Il_{k - 1} = Il_k; Jl_{k - 1} = Jl_k; if (k in indices && il_k != ir_k) - Il_k += a_strides[k] - continue - end - ) - ) (k -> (Il_k += a_strides[k]; if !(k in indices) - Jl_k += result_strides[k] - end)) begin - result[Jl_0, Jr_0] += conj(a[Il_0]) * a[Ir_0] + @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1)) + @nloops $RANK ir (d->1:shape[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides[d]; if !(d in indices) Jr_d+=result_strides[d] end)) begin + @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1)) + @nloops $RANK il (k->1:shape[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides[k]; continue end)) (k->(Il_k+=a_strides[k]; if !(k in indices) Jl_k+=result_strides[k] end)) begin + result[Jl_0, Jr_0] += conj(a[Il_0])*a[Ir_0] end end return result @@ -395,61 +291,28 @@ Fast in-place multiplication of operators/state vectors. Updates `Y` as Julia's 5-arg mul! implementation on the underlying data. See also [`LinearAlgebra.mul!`](@ref). """ -mul!( - result::Operator{B1,B3}, - a::Operator{B1,B2}, - b::Operator{B2,B3}, - alpha, - beta, -) where {B1<:Basis,B2<:Basis,B3<:Basis} = - (LinearAlgebra.mul!(result.data, a.data, b.data, alpha, beta); result) -mul!( - result::Ket{B1}, - a::Operator{B1,B2}, - b::Ket{B2}, - alpha, - beta, -) where {B1<:Basis,B2<:Basis} = - (LinearAlgebra.mul!(result.data, a.data, b.data, alpha, beta); result) -mul!( - result::Bra{B2}, - a::Bra{B1}, - b::Operator{B1,B2}, - alpha, - beta, -) where {B1<:Basis,B2<:Basis} = - (LinearAlgebra.mul!(result.data, transpose(b.data), a.data, alpha, beta); result) +mul!(result::Operator{B1,B3},a::Operator{B1,B2},b::Operator{B2,B3},alpha,beta) where {B1<:Basis,B2<:Basis,B3<:Basis} = (LinearAlgebra.mul!(result.data,a.data,b.data,alpha,beta); result) +mul!(result::Ket{B1},a::Operator{B1,B2},b::Ket{B2},alpha,beta) where {B1<:Basis,B2<:Basis} = (LinearAlgebra.mul!(result.data,a.data,b.data,alpha,beta); result) +mul!(result::Bra{B2},a::Bra{B1},b::Operator{B1,B2},alpha,beta) where {B1<:Basis,B2<:Basis} = (LinearAlgebra.mul!(result.data,transpose(b.data),a.data,alpha,beta); result) rmul!(op::Operator, x) = (rmul!(op.data, x); op) # Multiplication for Operators in terms of their gemv! implementation -function mul!( - result::Operator{B1,B3}, - M::AbstractOperator{B1,B2}, - b::Operator{B2,B3}, - alpha, - beta, -) 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]) - mul!(resultket, M, bket, alpha, beta) - result.data[:, i] = resultket.data +function mul!(result::Operator{B1,B3},M::AbstractOperator{B1,B2},b::Operator{B2,B3},alpha,beta) 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]) + mul!(resultket,M,bket,alpha,beta) + result.data[:,i] = resultket.data end return result end -function mul!( - result::Operator{B1,B3}, - b::Operator{B1,B2}, - M::AbstractOperator{B2,B3}, - alpha, - beta, -) 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, :])) - mul!(resultbra, bbra, M, alpha, beta) - result.data[i, :] = resultbra.data +function mul!(result::Operator{B1,B3},b::Operator{B1,B2},M::AbstractOperator{B2,B3},alpha,beta) 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,:])) + mul!(resultbra,bbra,M,alpha,beta) + result.data[i,:] = resultbra.data end return result end @@ -465,19 +328,13 @@ abstract type DataOperatorStyle{BL<:Basis,BR<:Basis} <: Broadcast.BroadcastStyle struct OperatorStyle{BL<:Basis,BR<:Basis} <: DataOperatorStyle{BL,BR} end # Style precedence rules -Broadcast.BroadcastStyle(::Type{<:Operator{BL,BR}}) where {BL<:Basis,BR<:Basis} = - OperatorStyle{BL,BR}() -Broadcast.BroadcastStyle( - ::OperatorStyle{B1,B2}, - ::OperatorStyle{B3,B4}, -) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(IncompatibleBases()) +Broadcast.BroadcastStyle(::Type{<:Operator{BL,BR}}) where {BL<:Basis,BR<:Basis} = OperatorStyle{BL,BR}() +Broadcast.BroadcastStyle(::OperatorStyle{B1,B2}, ::OperatorStyle{B3,B4}) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(IncompatibleBases()) # Out-of-place broadcasting -@inline function Base.copy( - bc::Broadcast.Broadcasted{Style,Axes,F,Args}, -) where {BL<:Basis,BR<:Basis,Style<:OperatorStyle{BL,BR},Axes,F,Args<:Tuple} +@inline function Base.copy(bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:OperatorStyle{BL,BR},Axes,F,Args<:Tuple} bcf = Broadcast.flatten(bc) - bl, br = find_basis(bcf.args) + bl,br = find_basis(bcf.args) bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf)) return Operator{BL,BR}(bl, br, copy(bc_)) end @@ -485,7 +342,7 @@ find_basis(a::DataOperator, rest) = (a.basis_l, a.basis_r) const BasicMathFunc = Union{typeof(+),typeof(-),typeof(*)} function Broadcasted_restrict_f(f::BasicMathFunc, args::Tuple{Vararg{<:DataOperator}}, axes) - args_ = Tuple(a.data for a in args) + args_ = Tuple(a.data for a=args) return Broadcast.Broadcasted(f, args_, axes) end function Broadcasted_restrict_f(f, args::Tuple{Vararg{<:DataOperator}}, axes) @@ -493,10 +350,7 @@ function Broadcasted_restrict_f(f, args::Tuple{Vararg{<:DataOperator}}, axes) end # In-place broadcasting -@inline function Base.copyto!( - dest::DataOperator{BL,BR}, - bc::Broadcast.Broadcasted{Style,Axes,F,Args}, -) where {BL<:Basis,BR<:Basis,Style<:DataOperatorStyle{BL,BR},Axes,F,Args} +@inline function Base.copyto!(dest::DataOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:DataOperatorStyle{BL,BR},Axes,F,Args} axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc)) # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match if bc.f === identity && isa(bc.args, Tuple{<:DataOperator{BL,BR}}) # only a single input argument to broadcast! @@ -511,12 +365,6 @@ end copyto!(dest.data, bc_) return dest end -@inline Base.copyto!( - A::DataOperator{BL,BR}, - B::DataOperator{BL,BR}, -) where {BL<:Basis,BR<:Basis} = (copyto!(A.data, B.data); A) -@inline Base.copyto!( - dest::DataOperator{BL,BR}, - bc::Broadcast.Broadcasted{Style,Axes,F,Args}, -) where {BL<:Basis,BR<:Basis,Style<:DataOperatorStyle,Axes,F,Args} = +@inline Base.copyto!(A::DataOperator{BL,BR},B::DataOperator{BL,BR}) where {BL<:Basis,BR<:Basis} = (copyto!(A.data,B.data); A) +@inline Base.copyto!(dest::DataOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:DataOperatorStyle,Axes,F,Args} = throw(IncompatibleBases()) diff --git a/src/states.jl b/src/states.jl index 0729a55..4f5ab6d 100644 --- a/src/states.jl +++ b/src/states.jl @@ -120,6 +120,7 @@ Norm of the given bra or ket state. """ norm(x::StateVector) = norm(x.data) + """ normalize(x::StateVector) @@ -127,6 +128,7 @@ Return the normalized state so that `norm(x)` is one. """ normalize(x::StateVector) = x / norm(x) + """ normalize!(x::StateVector) @@ -134,6 +136,7 @@ In-place normalization of the given bra or ket so that `norm(x)` is one. """ normalize!(x::StateVector) = (normalize!(x.data); x) + function permutesystems(state::T, perm::Vector{Int}) where {T<:Ket} @assert length(state.basis.bases) == length(perm) @assert isperm(perm) diff --git a/test/runtests.jl b/test/runtests.jl index 264b9ba..0204ba8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,7 +12,9 @@ names = [ "test_nlevel.jl", "test_state_definitions.jl", "test_metrics.jl", + "test_embed.jl", "test_superoperators.jl", + "test_abstractdata.jl", "test_pauli.jl", ] @@ -29,6 +31,8 @@ if length(unavailable_tests) != 0 error("The following tests could not be found:\n", join(unavailable_tests, "\n")) end +# names = ["test_abstractdata.jl"] + for name in names if startswith(name, "test_") && endswith(name, ".jl") include(name) diff --git a/test/test_abstractdata.jl b/test/test_abstractdata.jl new file mode 100644 index 0000000..af0ceee --- /dev/null +++ b/test/test_abstractdata.jl @@ -0,0 +1,412 @@ +using OpenQuantumSystems +using Test +using LinearAlgebra +import LinearAlgebra: mul! +using Random +import Base: ==, - + +# Implement custom type with AbstractArray interface +mutable struct TestData{T,N,X} <: AbstractArray{T,N} + x::X + function TestData(x::X) where X + x_ = convert(Matrix{ComplexF64}, x) + new{ComplexF64,length(axes(x)),typeof(x_)}(x_) + end +end +Base.size(A::TestData) = size(A.x) +Base.getindex(A::TestData, inds...) = getindex(A.x, inds...) +Base.setindex!(A::TestData, val, inds...) = setindex!(A.x, val, inds...) +Base.isapprox(A::TestData,B::TestData,args...;kwargs...) = isapprox(A.x,B.x,args...;kwargs...) + +# Additional methods +Base.copy(A::TestData) = TestData(copy(A.x)) +Base.kron(A::TestData,B::TestData) = TestData(kron(A.x,B.x)) +# -(A::TestData) = TestData(-A.x) + +# Auxiliary functions +function randtestoperator(args...) + op_ = randoperator(args...) + return Operator(op_.basis_l,op_.basis_r,TestData(copy(op_.data))) +end + +D(a::DataOperator,b::DataOperator,tol=1e-14) = isapprox(a.data,b.data,atol=tol) +D(a::StateVector,b::StateVector,tol=1e-14) = isapprox(a.data,b.data,atol=tol) +D(a::AbstractOperator,b::AbstractOperator,tol=1e-12) = (tol > abs(tracedistance_nh(dense(a), dense(b)))) + +@testset "abstract-data" begin + + Random.seed!(0) + + b1a = GenericBasis(2) + b1b = GenericBasis(3) + b2a = GenericBasis(1) + b2b = GenericBasis(4) + b3a = GenericBasis(1) + b3b = GenericBasis(5) + + b_l = b1a⊗b2a⊗b3a + b_r = b1b⊗b2b⊗b3b + + # Test creation + @test_throws DimensionMismatch Operator(b1a, TestData([1 1 1; 1 1 1])) + @test_throws DimensionMismatch Operator(b1a, b1b, TestData([1 1; 1 1; 1 1])) + op1 = Operator(b1a, b1b, TestData([1 1 1; 1 1 1])) + op2 = Operator(b1b, b1a, TestData([1 1; 1 1; 1 1])) + op3 = op1*op2 + @test op1 == dagger(op2) + + # Test ' shorthand + @test dagger(op2) == op2' + @test transpose(op2) == conj(op2') + + # Test copy + op1 = randtestoperator(b1a) + op2 = copy(op1) + @test op1.data == op2.data + @test !(op1.data === op2.data) + op2.data[1,1] = complex(10.) + @test op1.data[1,1] != op2.data[1,1] + + + # Arithmetic operations + # ===================== + op_zero = DenseOperator(b_l, b_r) + op1 = randtestoperator(b_l, b_r) + op2 = randtestoperator(b_l, b_r) + op3 = randtestoperator(b_l, b_r) + + x1 = Ket(b_r, rand(ComplexF64, length(b_r))) + x2 = Ket(b_r, rand(ComplexF64, length(b_r))) + + xbra1 = Bra(b_l, rand(ComplexF64, length(b_l))) + xbra2 = Bra(b_l, rand(ComplexF64, length(b_l))) + + # Addition + @test_throws DimensionMismatch op1 + dagger(op2) + @test D(op1 + op_zero, op1) + @test D(op1 + op2, op2 + op1) + @test D(op1 + (op2 + op3), (op1 + op2) + op3) + + # Subtraction + @test_throws DimensionMismatch op1 - dagger(op2) + @test D(op1-op_zero, op1) + @test D(op1-op2, op1 + (-op2)) + @test D(op1-op2, op1 + (-1*op2)) + @test D(op1-op2-op3, op1-(op2+op3)) + + # Test multiplication + @test_throws DimensionMismatch op1*op2 + @test D(op1*(x1 + 0.3*x2), op1*x1 + 0.3*op1*x2, 1e-12) + @test D((op1+op2)*(x1+0.3*x2), op1*x1 + 0.3*op1*x2 + op2*x1 + 0.3*op2*x2, 1e-12) + + @test D((xbra1+0.3*xbra2)*op1, xbra1*op1 + 0.3*xbra2*op1) + @test D((xbra1+0.3*xbra2)*(op1+op2), xbra1*op1 + 0.3*xbra2*op1 + xbra1*op2 + 0.3*xbra2*op2) + + @test D(op1*dagger(0.3*op2), 0.3*dagger(op2*dagger(op1))) + @test D((op1 + op2)*dagger(0.3*op3), 0.3*op1*dagger(op3) + 0.3*op2*dagger(op3), 1e-12) + @test D(0.3*(op1*dagger(op2)), op1*(0.3*dagger(op2))) + + tmp = copy(op1) + conj!(tmp) + @test tmp == conj(op1) && conj(tmp.data) == op1.data + + # Internal layout + b1 = GenericBasis(2) + b2 = GenericBasis(3) + b3 = GenericBasis(4) + op1 = randtestoperator(b1, b2) + op2 = randtestoperator(b2, b3) + x1 = randstate(b2) + d1 = op1.data + d2 = op2.data + v = x1.data + @test (op1*x1).data ≈ [d1[1,1]*v[1] + d1[1,2]*v[2] + d1[1,3]*v[3], d1[2,1]*v[1] + d1[2,2]*v[2] + d1[2,3]*v[3]] + @test (op1*op2).data[2,3] ≈ d1[2,1]*d2[1,3] + d1[2,2]*d2[2,3] + d1[2,3]*d2[3,3] + + # Test division + @test D(op1/7, (1/7)*op1) + + # Tensor product + # ============== + op1a = randtestoperator(b1a, b1b) + op1b = randtestoperator(b1a, b1b) + op2a = randtestoperator(b2a, b2b) + op2b = randtestoperator(b2a, b2b) + op3a = randtestoperator(b3a, b3b) + op123 = op1a ⊗ op2a ⊗ op3a + @test isa(op123.data,TestData) + @test op123.basis_l == b_l + @test op123.basis_r == b_r + + # Associativity + @test D((op1a ⊗ op2a) ⊗ op3a, op1a ⊗ (op2a ⊗ op3a)) + + # Linearity + @test D(op1a ⊗ (0.3*op2a), 0.3*(op1a ⊗ op2a)) + @test D((0.3*op1a) ⊗ op2a, 0.3*(op1a ⊗ op2a)) + + # Distributivity + @test D(op1a ⊗ (op2a + op2b), op1a ⊗ op2a + op1a ⊗ op2b) + @test D((op2a + op2b) ⊗ op3a, op2a ⊗ op3a + op2b ⊗ op3a) + + # Mixed-product property + @test D((op1a ⊗ op2a) * dagger(op1b ⊗ op2b), (op1a*dagger(op1b)) ⊗ (op2a*dagger(op2b))) + + # Transpose + @test D(dagger(op1a ⊗ op2a), dagger(op1a) ⊗ dagger(op2a)) + @test D(dagger(op1a ⊗ op2a), dagger(op1a) ⊗ dagger(op2a)) + + # Internal layout + a = Ket(b1a, rand(ComplexF64, length(b1a))) + b = Ket(b2b, rand(ComplexF64, length(b2b))) + ab = a ⊗ dagger(b) + @test ab.data[2,3] == a.data[2]*conj(b.data[3]) + @test ab.data[2,1] == a.data[2]*conj(b.data[1]) + + shape = tuple(op123.basis_l.shape..., op123.basis_r.shape...) + idx = LinearIndices(shape)[2, 1, 1, 3, 4, 5] + @test op123.data[idx] == op1a.data[2,3]*op2a.data[1,4]*op3a.data[1,5] + @test reshape(op123.data, shape...)[2, 1, 1, 3, 4, 5] == op1a.data[2,3]*op2a.data[1,4]*op3a.data[1,5] + + idx = LinearIndices(shape)[2, 1, 1, 1, 3, 4] + @test op123.data[idx] == op1a.data[2,1]*op2a.data[1,3]*op3a.data[1,4] + @test reshape(op123.data, shape...)[2, 1, 1, 1, 3, 4] == op1a.data[2,1]*op2a.data[1,3]*op3a.data[1,4] + + + # Test tr and normalize + op = Operator(GenericBasis(3), TestData([1 3 2;5 2 2;-1 2 5])) + @test 8 == tr(op) + op_normalized = normalize(op) + @test 8 == tr(op) + @test 1 == tr(op_normalized) + op_copy = deepcopy(op) + normalize!(op_copy) + @test tr(op) != tr(op_copy) + @test 1 ≈ tr(op_copy) + @test op === normalize!(op) + + # Test partial tr of operators + b1 = GenericBasis(3) + b2 = GenericBasis(5) + b3 = GenericBasis(7) + op1 = randtestoperator(b1) + op2 = randtestoperator(b2) + op3 = randtestoperator(b3) + op123 = op1 ⊗ op2 ⊗ op3 + + @test D(op1⊗op2*tr(op3), ptrace(op123, 3)) + @test D(op1⊗op3*tr(op2), ptrace(op123, 2)) + @test D(op2⊗op3*tr(op1), ptrace(op123, 1)) + + @test D(op1*tr(op2)*tr(op3), ptrace(op123, [2,3])) + @test D(op2*tr(op1)*tr(op3), ptrace(op123, [1,3])) + @test D(op3*tr(op1)*tr(op2), ptrace(op123, [1,2])) + + @test_throws ArgumentError ptrace(op123, [1,2,3]) + x = randtestoperator(b1, b1⊗b2) + @test_throws ArgumentError ptrace(x, [1]) + x = randtestoperator(b1⊗b1⊗b2, b1⊗b2) + @test_throws ArgumentError ptrace(x, [1, 2]) + x = randtestoperator(b1⊗b2) + @test_throws ArgumentError ptrace(x, [1, 2]) + x = randtestoperator(b1⊗b2, b2⊗b1) + @test_throws ArgumentError ptrace(x, [1]) + + op1 = randtestoperator(b1, b2) + op2 = randtestoperator(b3) + + @test D(op1*tr(op2), ptrace(op1⊗op2, 2)) + + # Test expect + b1 = GenericBasis(3) + b2 = GenericBasis(5) + b3 = GenericBasis(7) + op1 = randtestoperator(b1) + op2 = randtestoperator(b2) + op3 = randtestoperator(b3) + op123 = op1 ⊗ op2 ⊗ op3 + b_l = b1 ⊗ b2 ⊗ b3 + + state = randstate(b_l) + @test expect(op123, state) ≈ dagger(state)*op123*state + @test expect(1, op1, state) ≈ expect(op1, ptrace(state, [2, 3])) + @test expect(2, op2, state) ≈ expect(op2, ptrace(state, [1, 3])) + @test expect(3, op3, state) ≈ expect(op3, ptrace(state, [1, 2])) + + state = randtestoperator(b_l) + @test expect(op123, state) ≈ tr(op123*state) + @test expect(1, op1, state) ≈ expect(op1, ptrace(state, [2, 3])) + @test expect(2, op2, state) ≈ expect(op2, ptrace(state, [1, 3])) + @test expect(3, op3, state) ≈ expect(op3, ptrace(state, [1, 2])) + + # Permute systems + op1 = randtestoperator(b1a, b1b) + op2 = randtestoperator(b2a, b2b) + op3 = randtestoperator(b3a, b3b) + op123 = op1⊗op2⊗op3 + + op132 = op1⊗op3⊗op2 + @test D(permutesystems(op123, [1, 3, 2]), op132) + + op213 = op2⊗op1⊗op3 + @test D(permutesystems(op123, [2, 1, 3]), op213) + + op231 = op2⊗op3⊗op1 + @test D(permutesystems(op123, [2, 3, 1]), op231) + + op312 = op3⊗op1⊗op2 + @test D(permutesystems(op123, [3, 1, 2]), op312) + + op321 = op3⊗op2⊗op1 + @test D(permutesystems(op123, [3, 2, 1]), op321) + + + # Test projector + xket = normalize(Ket(b_l, rand(ComplexF64, length(b_l)))) + yket = normalize(Ket(b_l, rand(ComplexF64, length(b_l)))) + xbra = dagger(xket) + ybra = dagger(yket) + + @test D(projector(xket)*xket, xket) + @test D(xbra*projector(xket), xbra) + @test D(projector(xbra)*xket, xket) + @test D(xbra*projector(xbra), xbra) + @test D(ybra*projector(yket, xbra), xbra) + @test D(projector(yket, xbra)*xket, yket) + + # Test gemv + b1 = GenericBasis(3) + b2 = GenericBasis(5) + op = randtestoperator(b1, b2) + xket = randstate(b2) + xbra = dagger(randstate(b1)) + rket = randstate(b1) + rbra = dagger(randstate(b2)) + alpha = complex(0.7, 1.5) + beta = complex(0.3, 2.1) + @test isa(op.data, TestData) + + rket_ = deepcopy(rket) + OpenQuantumSystems.mul!(rket_,op,xket,complex(1.0),complex(0.)) + @test D(rket_, op*xket) + + @test isa(op.data, TestData) + rket_ = deepcopy(rket) + OpenQuantumSystems.mul!(rket_,op,xket,alpha,beta) + @test D(rket_, alpha*op*xket + beta*rket) + + rbra_ = deepcopy(rbra) + OpenQuantumSystems.mul!(rbra_,xbra,op,complex(1.0),complex(0.)) + @test D(rbra_, xbra*op) + + rbra_ = deepcopy(rbra) + OpenQuantumSystems.mul!(rbra_,xbra,op,alpha,beta) + @test D(rbra_, alpha*xbra*op + beta*rbra) + + # Test gemm + b1 = GenericBasis(37) + b2 = GenericBasis(53) + b3 = GenericBasis(41) + op1 = randtestoperator(b1, b2) + op2 = randtestoperator(b2, b3) + r = randtestoperator(b1, b3) + alpha = complex(0.7, 1.5) + beta = complex(0.3, 2.1) + + r_ = deepcopy(r) + OpenQuantumSystems.mul!(r_,op1,op2,complex(1.0),complex(0.)) + @test D(r_, op1*op2) + + r_ = deepcopy(r) + OpenQuantumSystems.mul!(r_,op1,op2,alpha,beta) + @test D(r_, alpha*op1*op2 + beta*r, 1e-12) + + # Test Hermitian + bspin = SpinBasis(1//2) + bnlevel = NLevelBasis(2) + @test ishermitian(Operator(bspin, bspin, TestData([1.0 im; -im 2.0]))) == true + @test ishermitian(Operator(bspin, bnlevel, TestData([1.0 im; -im 2.0]))) == false + + # Test broadcasting + op1_ = copy(op1) + @test isa(op1_.data, TestData) + op1 .= 2*op1 + @test op1 == op1_ .+ op1_ + op1 .= op1_ + @test op1 == op1_ + op1 .= op1_ .+ 3 * op1_ + @test op1 == 4*op1_ + @test_throws DimensionMismatch op1 .= op2 + bf = FockBasis(3) + op3 = randtestoperator(bf) + @test_throws OpenQuantumSystems.IncompatibleBases op1 .+ op3 + @test_throws ErrorException cos.(op1) + + #################### + + # Test gemv + op1 = randtestoperator(b_l, b_r) + op2 = randtestoperator(b_r, b_l) + op3 = randtestoperator(b_l, b_r) + op = op1 * sparse(op2) * op3 * 0.2 + op_ = 0.2*op1*op2*op3 + tmp = 0.2*prod(dense.([op1*op2*op3])) + + state = Ket(b_r, rand(ComplexF64, length(b_r))) + result_ = Ket(b_l, rand(ComplexF64, length(b_l))) + result = copy(result_) + OpenQuantumSystems.mul!(result,op,state,complex(1.),complex(0.)) + @test D(result, op_*state, 1e-10) + + result = copy(result_) + alpha = complex(1.5) + beta = complex(2.1) + OpenQuantumSystems.mul!(result,op,state,alpha,beta) + @test D(result, alpha*op_*state + beta*result_, 1e-9) + + state = Bra(b_l, rand(ComplexF64, length(b_l))) + result_ = Bra(b_r, rand(ComplexF64, length(b_r))) + result = copy(result_) + OpenQuantumSystems.mul!(result,state,op,complex(1.),complex(0.)) + @test D(result, state*op_, 1e-9) + + result = copy(result_) + alpha = complex(1.5) + beta = complex(2.1) + OpenQuantumSystems.mul!(result,state,op,alpha,beta) + @test D(result, alpha*state*op_ + beta*result_, 1e-9) + + # Test gemm + op1 = randtestoperator(b_l, b_r) + op2 = randtestoperator(b_r, b_l) + op3 = randtestoperator(b_l, b_r) + op = op1 * sparse(op2) * op3 * 0.2 + op_ = 0.2*op1*op2*op3 + + state = randtestoperator(b_r, b_r) + result_ = randtestoperator(b_l, b_r) + result = copy(result_) + result.data[1] = 1 + result.data[1] != result_.data[1] || error("") + OpenQuantumSystems.mul!(result,op,state,complex(1.),complex(0.)) + @test D(result, op_*state, 1e-9) + + result = copy(result_) + alpha = complex(1.5) + beta = complex(2.1) + OpenQuantumSystems.mul!(result,op,state,alpha,beta) + @test D(result, alpha*op_*state + beta*result_, 1e-9) + + state = randtestoperator(b_l, b_l) + result_ = randtestoperator(b_l, b_r) + result = copy(result_) + OpenQuantumSystems.mul!(result,state,op,complex(1.),complex(0.)) + @test D(result, state*op_, 1e-9) + + result = copy(result_) + alpha = complex(1.5) + beta = complex(2.1) + OpenQuantumSystems.mul!(result,state,op,alpha,beta) + @test D(result, alpha*state*op_ + beta*result_, 1e-8) + +end # testset diff --git a/test/test_embed.jl b/test/test_embed.jl new file mode 100644 index 0000000..c8934f9 --- /dev/null +++ b/test/test_embed.jl @@ -0,0 +1,74 @@ +using Test +using OpenQuantumSystems +using Random, SparseArrays, LinearAlgebra + +@testset "embed" begin + +Random.seed!(0) + + # Set up operators + spinbasis = SpinBasis(1//2) + + b1 = NLevelBasis(3) + b2 = SpinBasis(1//2) + b3 = FockBasis(2) + + I1 = dense(identityoperator(b1)) + I2 = dense(identityoperator(b2)) + I3 = dense(identityoperator(b3)) + + b = b1 ⊗ b2 ⊗ b3 + + op1 = DenseOperator(b1, b1, rand(ComplexF64, length(b1), length(b1))) + op2 = DenseOperator(b2, b2, rand(ComplexF64, length(b2), length(b2))) + op3 = DenseOperator(b3, b3, rand(ComplexF64, length(b3), length(b3))) + + + # Test Vector{Int}, Vector{AbstractOperator} + x = embed(b, [1,2], [op1, op2]) + y = op1 ⊗ op2 ⊗ I3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + + x = embed(b, [1,2], [sparse(op1), sparse(op2)]) + y = op1 ⊗ op2 ⊗ I3 + @test 0 ≈ abs(tracedistance_nh(dense(x), y)) + + x = embed(b, 1, op1) + y = op1 ⊗ I2 ⊗ I3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + + x = embed(b, 2, op2) + y = I1 ⊗ op2 ⊗ I3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + + x = embed(b, 3, op3) + y = I1 ⊗ I2 ⊗ op3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + + + # 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)) + + x = embed(b, Dict(1 => op1, 2 => op2)) + y = op1 ⊗ op2 ⊗ I3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + + x = embed(b, Dict([1,3] => sparse(op1⊗op3))) + y = op1 ⊗ I2 ⊗ op3 + @test 0 ≈ abs(tracedistance_nh(dense(x), y)) + + x = embed(b, Dict([1,3] => op1⊗op3)) + y = op1 ⊗ I2 ⊗ op3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + + x = embed(b, Dict([3,1] => sparse(op3⊗op1))) + y = op1 ⊗ I2 ⊗ op3 + @test 0 ≈ abs(tracedistance_nh(dense(x), y)) + + x = embed(b, Dict([3,1] => op3⊗op1)) + y = op1 ⊗ I2 ⊗ op3 + @test 0 ≈ abs(tracedistance_nh(x, y)) + +end # testset diff --git a/test/test_operators.jl b/test/test_operators.jl index ce8745a..8999a8e 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -3,12 +3,10 @@ using OpenQuantumSystems using LinearAlgebra, SparseArrays, Random 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{typeof(b1),typeof(b2)}(b1, b2, data) : throw(DimensionMismatch()) + 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{typeof(b1),typeof(b2)}(b1, b2, data) : throw(DimensionMismatch()) end @testset "operators" begin @@ -30,7 +28,7 @@ end @test basis(op1) == b1 @test length(op1) == length(op1.data) == length(b1)^2 - @test_throws ArgumentError op_test * op_test + @test_throws ArgumentError op_test*op_test @test_throws ArgumentError -op_test @test_throws ArgumentError 1 + op_test @@ -50,27 +48,27 @@ end @test expect(1, op1, ρ) ≈ expect(embed(b, 1, op1), ρ) @test expect(1, op1, ψ) ≈ expect(embed(b, 1, op1), ψ) - @test expect(op, [ρ, ρ]) == [expect(op, ρ) for i = 1:2] + @test expect(op, [ρ, ρ]) == [expect(op, ρ) for i=1:2] @test expect(1, op1, [ρ, ψ]) == [expect(1, op1, ρ), expect(1, op1, ψ)] @test variance(1, op1, ρ) ≈ variance(embed(b, 1, op1), ρ) @test variance(1, op1, ψ) ≈ variance(embed(b, 1, op1), ψ) - @test variance(op, [ρ, ρ]) == [variance(op, ρ) for i = 1:2] + @test variance(op, [ρ, ρ]) == [variance(op, ρ) for i=1:2] @test variance(1, op1, [ρ, ψ]) == [variance(1, op1, ρ), variance(1, op1, ψ)] @test tensor(op_test) === op_test @test_throws ArgumentError tensor(op_test, op_test) @test_throws ArgumentError permutesystems(op_test, [1, 2]) - @test embed(b, b, [1, 2], op) == embed(b, [1, 2], op) - @test embed(b, Dict{Vector{Int},SparseOpType}()) == identityoperator(b) - @test_throws OpenQuantumSystems.IncompatibleBases embed(b1 ⊗ b2, [2], [op1]) + @test embed(b, b, [1,2], op) == embed(b, [1,2], op) + @test embed(b, Dict{Vector{Int}, SparseOpType}()) == identityoperator(b) + @test_throws OpenQuantumSystems.IncompatibleBases embed(b1⊗b2, [2], [op1]) - b_comp = b ⊗ b - @test isapprox(embed(b_comp, [1, [3, 4]], [op1, op]), dense(op1 ⊗ one(b2) ⊗ op)) - @test isapprox(embed(b_comp, [[1, 2], 4], [op, op2]), dense(op ⊗ one(b1) ⊗ op2)) - @test_throws OpenQuantumSystems.IncompatibleBases embed(b_comp, [[1, 2], 3], [op, op2]) - @test_throws OpenQuantumSystems.IncompatibleBases embed(b_comp, [[1, 3], 4], [op, op2]) + b_comp = b⊗b + @test isapprox(embed(b_comp, [1,[3,4]], [op1,op]), dense(op1 ⊗ one(b2) ⊗ op)) + @test isapprox(embed(b_comp, [[1,2],4], [op,op2]), dense(op ⊗ one(b1) ⊗ op2)) + @test_throws OpenQuantumSystems.IncompatibleBases embed(b_comp, [[1,2],3], [op,op2]) + @test_throws OpenQuantumSystems.IncompatibleBases embed(b_comp, [[1,3],4], [op,op2]) function basis_vec(n, N) x = zeros(Complex{Float64}, N) @@ -84,49 +82,34 @@ end end end - embed_op = embed(b_comp, [1, 4], op) - bv = basis_maker(3, 2, 3, 2) + embed_op = embed(b_comp, [1,4], op) + bv = basis_maker(3,2,3,2) all_idxs = [(idx, jdx) for (idx, jdx) in [Iterators.product(0:1, 0:2)...]] - m11 = reshape( - [ - Bra(b_comp, bv(0, idx, jdx, 0)) * embed_op * Ket(b_comp, bv(0, kdx, ldx, 0)) - for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs) - ], - (6, 6), - ) - @test isapprox(m11 / op.data[1, 1], diagm(0 => ones(Complex{Float64}, 6))) - - m21 = reshape( - [ - Bra(b_comp, bv(1, idx, jdx, 0)) * embed_op * Ket(b_comp, bv(0, kdx, ldx, 0)) - for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs) - ], - (6, 6), - ) - @test isapprox(m21 / op.data[2, 1], diagm(0 => ones(Complex{Float64}, 6))) - - m12 = reshape( - [ - Bra(b_comp, bv(0, idx, jdx, 0)) * embed_op * Ket(b_comp, bv(1, kdx, ldx, 0)) - for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs) - ], - (6, 6), - ) - @test isapprox(m12 / op.data[1, 2], diagm(0 => ones(Complex{Float64}, 6))) - - - b_comp = b_comp ⊗ b_comp - OP_test1 = dense(tensor([op1, one(b2), op, one(b1), one(b2), op1, one(b2)]...)) - OP_test2 = embed(b_comp, [1, [3, 4], 7], [op1, op, op1]) - @test isapprox(OP_test1.data, OP_test2.data) - - b8 = b2 ⊗ b2 ⊗ b2 + m11 = reshape([Bra(b_comp, bv(0,idx,jdx,0)) * embed_op * Ket(b_comp, bv(0,kdx,ldx,0)) + for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs)], (6,6)) + @test isapprox(m11 / op.data[1, 1], diagm(0=>ones(Complex{Float64}, 6))) + + m21 = reshape([Bra(b_comp, bv(1,idx,jdx,0)) * embed_op * Ket(b_comp, bv(0,kdx,ldx,0)) + for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs)], (6,6)) + @test isapprox(m21 / op.data[2,1], diagm(0=>ones(Complex{Float64}, 6))) + + m12 = reshape([Bra(b_comp, bv(0,idx,jdx,0)) * embed_op * Ket(b_comp, bv(1,kdx,ldx,0)) + for ((idx, jdx), (kdx, ldx)) in Iterators.product(all_idxs, all_idxs)], (6,6)) + @test isapprox(m12 / op.data[1,2], diagm(0=>ones(Complex{Float64}, 6))) + + + b_comp = b_comp⊗b_comp + OP_test1 = dense(tensor([op1,one(b2),op,one(b1),one(b2),op1,one(b2)]...)) + OP_test2 = embed(b_comp, [1,[3,4],7], [op1,op,op1]) + @test isapprox(OP_test1, OP_test2) + + b8 = b2⊗b2⊗b2 cnot = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0] - op_cnot = DenseOperator(b2 ⊗ b2, cnot) - OP_cnot = embed(b8, [1, 3], op_cnot) - @test ptrace(OP_cnot, [2]) / 2.0 == op_cnot - @test_throws AssertionError embed(b2 ⊗ b2, [1, 1], op_cnot) + op_cnot = DenseOperator(b2⊗b2, cnot) + OP_cnot = embed(b8, [1,3], op_cnot) + @test ptrace(OP_cnot, [2])/2. == op_cnot + @test_throws AssertionError embed(b2⊗b2, [1,1], op_cnot) @test_throws ArgumentError exp(sparse(op1)) diff --git a/test/test_operators_dense.jl b/test/test_operators_dense.jl index 2fc8639..5dbabb2 100644 --- a/test/test_operators_dense.jl +++ b/test/test_operators_dense.jl @@ -6,9 +6,8 @@ using Random, SparseArrays, LinearAlgebra Random.seed!(0) - D(op1::AbstractOperator, op2::AbstractOperator) = - abs(tracedistance_nh(dense(op1), dense(op2))) - D(x1::StateVector, x2::StateVector) = norm(x2 - x1) + D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) + D(x1::StateVector, x2::StateVector) = norm(x2-x1) b1a = GenericBasis(2) b1b = GenericBasis(3) @@ -17,8 +16,8 @@ using Random, SparseArrays, LinearAlgebra b3a = GenericBasis(1) b3b = GenericBasis(5) - b_l = b1a ⊗ b2a ⊗ b3a - b_r = b1b ⊗ b2b ⊗ b3b + b_l = b1a⊗b2a⊗b3a + b_r = b1b⊗b2b⊗b3b # Test creation @test_throws DimensionMismatch DenseOperator(b1a, [1 1 1; 1 1 1]) @@ -36,8 +35,8 @@ using Random, SparseArrays, LinearAlgebra op2 = copy(op1) @test op1.data == op2.data @test !(op1.data === op2.data) - op2.data[1, 1] = complex(10.0) - @test op1.data[1, 1] != op2.data[1, 1] + op2.data[1,1] = complex(10.) + @test op1.data[1,1] != op2.data[1,1] # Arithmetic operations @@ -61,31 +60,22 @@ using Random, SparseArrays, LinearAlgebra # Subtraction @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 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 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, - ) - - @test 1e-11 > D((xbra1 + 0.3 * xbra2) * op1, xbra1 * op1 + 0.3 * xbra2 * op1) - @test 1e-11 > D( - (xbra1 + 0.3 * xbra2) * (op1 + op2), - xbra1 * op1 + 0.3 * xbra2 * op1 + xbra1 * op2 + 0.3 * xbra2 * op2, - ) - - @test 1e-12 > D(op1 * dagger(0.3 * op2), 0.3 * dagger(op2 * dagger(op1))) - @test 1e-12 > D( - (op1 + op2) * dagger(0.3 * op3), - 0.3 * op1 * dagger(op3) + 0.3 * op2 * dagger(op3), - ) - @test 1e-12 > D(0.3 * (op1 * dagger(op2)), op1 * (0.3 * dagger(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) + + @test 1e-11 > D((xbra1+0.3*xbra2)*op1, xbra1*op1 + 0.3*xbra2*op1) + @test 1e-11 > D((xbra1+0.3*xbra2)*(op1+op2), xbra1*op1 + 0.3*xbra2*op1 + xbra1*op2 + 0.3*xbra2*op2) + + @test 1e-12 > D(op1*dagger(0.3*op2), 0.3*dagger(op2*dagger(op1))) + @test 1e-12 > D((op1 + op2)*dagger(0.3*op3), 0.3*op1*dagger(op3) + 0.3*op2*dagger(op3)) + @test 1e-12 > D(0.3*(op1*dagger(op2)), op1*(0.3*dagger(op2))) tmp = copy(op1) conj!(tmp) @@ -101,15 +91,11 @@ using Random, SparseArrays, LinearAlgebra d1 = op1.data d2 = op2.data v = x1.data - @test (op1 * x1).data ≈ [ - d1[1, 1] * v[1] + d1[1, 2] * v[2] + d1[1, 3] * v[3], - d1[2, 1] * v[1] + d1[2, 2] * v[2] + d1[2, 3] * v[3], - ] - @test (op1*op2).data[2, 3] ≈ - d1[2, 1] * d2[1, 3] + d1[2, 2] * d2[2, 3] + d1[2, 3] * d2[3, 3] + @test (op1*x1).data ≈ [d1[1,1]*v[1] + d1[1,2]*v[2] + d1[1,3]*v[3], d1[2,1]*v[1] + d1[2,2]*v[2] + d1[2,3]*v[3]] + @test (op1*op2).data[2,3] ≈ d1[2,1]*d2[1,3] + d1[2,2]*d2[2,3] + d1[2,3]*d2[3,3] # Test division - @test 1e-14 > D(op1 / 7, (1 / 7) * op1) + @test 1e-14 > D(op1/7, (1/7)*op1) # Tensor product # ============== @@ -126,18 +112,15 @@ using Random, SparseArrays, LinearAlgebra @test 1e-13 > D((op1a ⊗ op2a) ⊗ op3a, op1a ⊗ (op2a ⊗ op3a)) # Linearity - @test 1e-13 > D(op1a ⊗ (0.3 * op2a), 0.3 * (op1a ⊗ op2a)) - @test 1e-13 > D((0.3 * op1a) ⊗ op2a, 0.3 * (op1a ⊗ op2a)) + @test 1e-13 > D(op1a ⊗ (0.3*op2a), 0.3*(op1a ⊗ op2a)) + @test 1e-13 > D((0.3*op1a) ⊗ op2a, 0.3*(op1a ⊗ op2a)) # Distributivity @test 1e-13 > D(op1a ⊗ (op2a + op2b), op1a ⊗ op2a + op1a ⊗ op2b) @test 1e-13 > D((op2a + op2b) ⊗ op3a, op2a ⊗ op3a + op2b ⊗ op3a) # Mixed-product property - @test 1e-13 > D( - (op1a ⊗ op2a) * dagger(op1b ⊗ op2b), - (op1a * dagger(op1b)) ⊗ (op2a * dagger(op2b)), - ) + @test 1e-13 > D((op1a ⊗ op2a) * dagger(op1b ⊗ op2b), (op1a*dagger(op1b)) ⊗ (op2a*dagger(op2b))) # Transpose @test 1e-13 > D(dagger(op1a ⊗ op2a), dagger(op1a) ⊗ dagger(op2a)) @@ -147,19 +130,17 @@ using Random, SparseArrays, LinearAlgebra a = Ket(b1a, rand(ComplexF64, length(b1a))) b = Ket(b2b, rand(ComplexF64, length(b2b))) ab = a ⊗ dagger(b) - @test ab.data[2, 3] == a.data[2] * conj(b.data[3]) - @test ab.data[2, 1] == a.data[2] * conj(b.data[1]) + @test ab.data[2,3] == a.data[2]*conj(b.data[3]) + @test ab.data[2,1] == a.data[2]*conj(b.data[1]) shape = tuple(op123.basis_l.shape..., op123.basis_r.shape...) idx = LinearIndices(shape)[2, 1, 1, 3, 4, 5] - @test op123.data[idx] == op1a.data[2, 3] * op2a.data[1, 4] * op3a.data[1, 5] - @test reshape(op123.data, shape...)[2, 1, 1, 3, 4, 5] == - op1a.data[2, 3] * op2a.data[1, 4] * op3a.data[1, 5] + @test op123.data[idx] == op1a.data[2,3]*op2a.data[1,4]*op3a.data[1,5] + @test reshape(op123.data, shape...)[2, 1, 1, 3, 4, 5] == op1a.data[2,3]*op2a.data[1,4]*op3a.data[1,5] idx = LinearIndices(shape)[2, 1, 1, 1, 3, 4] - @test op123.data[idx] == op1a.data[2, 1] * op2a.data[1, 3] * op3a.data[1, 4] - @test reshape(op123.data, shape...)[2, 1, 1, 1, 3, 4] == - op1a.data[2, 1] * op2a.data[1, 3] * op3a.data[1, 4] + @test op123.data[idx] == op1a.data[2,1]*op2a.data[1,3]*op3a.data[1,4] + @test reshape(op123.data, shape...)[2, 1, 1, 1, 3, 4] == op1a.data[2,1]*op2a.data[1,3]*op3a.data[1,4] # Test identityoperator @@ -171,21 +152,17 @@ using Random, SparseArrays, LinearAlgebra I = identityoperator(DenseOpType, b_r) @test isa(I, DenseOpType) @test identityoperator(SparseOpType, b_r) == sparse(I) - @test 1e-11 > D(I * x1, x1) - @test I == - identityoperator(DenseOpType, b1b) ⊗ identityoperator(DenseOpType, b2b) ⊗ - identityoperator(DenseOpType, b3b) + @test 1e-11 > D(I*x1, x1) + @test I == identityoperator(DenseOpType, b1b) ⊗ identityoperator(DenseOpType, b2b) ⊗ identityoperator(DenseOpType, b3b) I = identityoperator(DenseOpType, b_l) @test isa(I, DenseOpType) @test identityoperator(SparseOpType, b_l) == sparse(I) - @test 1e-11 > D(xbra1 * I, xbra1) - @test I == - identityoperator(DenseOpType, b1a) ⊗ identityoperator(DenseOpType, b2a) ⊗ - identityoperator(DenseOpType, b3a) + @test 1e-11 > D(xbra1*I, xbra1) + @test I == identityoperator(DenseOpType, b1a) ⊗ identityoperator(DenseOpType, b2a) ⊗ identityoperator(DenseOpType, b3a) # Test tr and normalize - op = DenseOperator(GenericBasis(3), [1 3 2; 5 2 2; -1 2 5]) + op = DenseOperator(GenericBasis(3), float.([1 3 2;5 2 2;-1 2 5])) @test 8 == tr(op) op_normalized = normalize(op) @test 8 == tr(op) @@ -197,21 +174,21 @@ using Random, SparseArrays, LinearAlgebra @test op === normalize!(op) # Test partial tr of state vectors - psi1 = 0.1 * randstate(b1a) - psi2 = 0.3 * randstate(b2a) - psi3 = 0.7 * randstate(b3a) + psi1 = 0.1*randstate(b1a) + psi2 = 0.3*randstate(b2a) + psi3 = 0.7*randstate(b3a) psi12 = psi1 ⊗ psi2 psi13 = psi1 ⊗ psi3 psi23 = psi2 ⊗ psi3 psi123 = psi1 ⊗ psi2 ⊗ psi3 - @test 1e-13 > D(0.1^2 * 0.3^2 * psi3 ⊗ dagger(psi3), ptrace(psi123, [1, 2])) - @test 1e-13 > D(0.1^2 * 0.7^2 * psi2 ⊗ dagger(psi2), ptrace(psi123, [1, 3])) - @test 1e-13 > D(0.3^2 * 0.7^2 * psi1 ⊗ dagger(psi1), ptrace(psi123, [2, 3])) + @test 1e-13 > D(0.1^2*0.3^2*psi3 ⊗ dagger(psi3), ptrace(psi123, [1, 2])) + @test 1e-13 > D(0.1^2*0.7^2*psi2 ⊗ dagger(psi2), ptrace(psi123, [1, 3])) + @test 1e-13 > D(0.3^2*0.7^2*psi1 ⊗ dagger(psi1), ptrace(psi123, [2, 3])) - @test 1e-13 > D(0.1^2 * psi23 ⊗ dagger(psi23), ptrace(psi123, 1)) - @test 1e-13 > D(0.3^2 * psi13 ⊗ dagger(psi13), ptrace(psi123, 2)) - @test 1e-13 > D(0.7^2 * psi12 ⊗ dagger(psi12), ptrace(psi123, 3)) + @test 1e-13 > D(0.1^2*psi23 ⊗ dagger(psi23), ptrace(psi123, 1)) + @test 1e-13 > D(0.3^2*psi13 ⊗ dagger(psi13), ptrace(psi123, 2)) + @test 1e-13 > D(0.7^2*psi12 ⊗ dagger(psi12), ptrace(psi123, 3)) @test 1e-13 > D(ptrace(psi123, [1, 2]), dagger(ptrace(dagger(psi123), [1, 2]))) @test 1e-13 > D(ptrace(psi123, 3), dagger(ptrace(dagger(psi123), 3))) @@ -227,28 +204,28 @@ using Random, SparseArrays, LinearAlgebra op3 = randoperator(b3) op123 = op1 ⊗ op2 ⊗ op3 - @test 1e-13 > D(op1 ⊗ op2 * tr(op3), ptrace(op123, 3)) - @test 1e-13 > D(op1 ⊗ op3 * tr(op2), ptrace(op123, 2)) - @test 1e-13 > D(op2 ⊗ op3 * tr(op1), ptrace(op123, 1)) + @test 1e-13 > D(op1⊗op2*tr(op3), ptrace(op123, 3)) + @test 1e-13 > D(op1⊗op3*tr(op2), ptrace(op123, 2)) + @test 1e-13 > D(op2⊗op3*tr(op1), ptrace(op123, 1)) - @test 1e-13 > D(op1 * tr(op2) * tr(op3), ptrace(op123, [2, 3])) - @test 1e-13 > D(op2 * tr(op1) * tr(op3), ptrace(op123, [1, 3])) - @test 1e-13 > D(op3 * tr(op1) * tr(op2), ptrace(op123, [1, 2])) + @test 1e-13 > D(op1*tr(op2)*tr(op3), ptrace(op123, [2,3])) + @test 1e-13 > D(op2*tr(op1)*tr(op3), ptrace(op123, [1,3])) + @test 1e-13 > D(op3*tr(op1)*tr(op2), ptrace(op123, [1,2])) - @test_throws ArgumentError ptrace(op123, [1, 2, 3]) - x = randoperator(b1, b1 ⊗ b2) + @test_throws ArgumentError ptrace(op123, [1,2,3]) + x = randoperator(b1, b1⊗b2) @test_throws ArgumentError ptrace(x, [1]) - x = randoperator(b1 ⊗ b1 ⊗ b2, b1 ⊗ b2) + x = randoperator(b1⊗b1⊗b2, b1⊗b2) @test_throws ArgumentError ptrace(x, [1, 2]) - x = randoperator(b1 ⊗ b2) + x = randoperator(b1⊗b2) @test_throws ArgumentError ptrace(x, [1, 2]) - x = randoperator(b1 ⊗ b2, b2 ⊗ b1) + x = randoperator(b1⊗b2, b2⊗b1) @test_throws ArgumentError ptrace(x, [1]) op1 = randoperator(b1, b2) op2 = randoperator(b3) - @test 1e-13 > D(op1 * tr(op2), ptrace(op1 ⊗ op2, 2)) + @test 1e-13 > D(op1*tr(op2), ptrace(op1⊗op2, 2)) # Test expect b1 = GenericBasis(3) @@ -261,36 +238,38 @@ using Random, SparseArrays, LinearAlgebra b_l = b1 ⊗ b2 ⊗ b3 state = randstate(b_l) - @test expect(op123, state) ≈ dagger(state) * op123 * state + @test expect(op123, state) ≈ dagger(state)*op123*state @test expect(1, op1, state) ≈ expect(op1, ptrace(state, [2, 3])) @test expect(2, op2, state) ≈ expect(op2, ptrace(state, [1, 3])) @test expect(3, op3, state) ≈ expect(op3, ptrace(state, [1, 2])) state = randoperator(b_l) - @test expect(op123, state) ≈ tr(op123 * state) + @test expect(op123, state) ≈ tr(op123*state) @test expect(1, op1, state) ≈ expect(op1, ptrace(state, [2, 3])) @test expect(2, op2, state) ≈ expect(op2, ptrace(state, [1, 3])) @test expect(3, op3, state) ≈ expect(op3, ptrace(state, [1, 2])) + @test_throws OpenQuantumSystems.IncompatibleBases expect(op2, op1) + # Permute systems op1 = randoperator(b1a, b1b) op2 = randoperator(b2a, b2b) op3 = randoperator(b3a, b3b) - op123 = op1 ⊗ op2 ⊗ op3 + op123 = op1⊗op2⊗op3 - op132 = op1 ⊗ op3 ⊗ op2 + op132 = op1⊗op3⊗op2 @test 1e-14 > D(permutesystems(op123, [1, 3, 2]), op132) - op213 = op2 ⊗ op1 ⊗ op3 + op213 = op2⊗op1⊗op3 @test 1e-14 > D(permutesystems(op123, [2, 1, 3]), op213) - op231 = op2 ⊗ op3 ⊗ op1 + op231 = op2⊗op3⊗op1 @test 1e-14 > D(permutesystems(op123, [2, 3, 1]), op231) - op312 = op3 ⊗ op1 ⊗ op2 + op312 = op3⊗op1⊗op2 @test 1e-14 > D(permutesystems(op123, [3, 1, 2]), op312) - op321 = op3 ⊗ op2 ⊗ op1 + op321 = op3⊗op2⊗op1 @test 1e-14 > D(permutesystems(op123, [3, 2, 1]), op321) @@ -300,18 +279,18 @@ using Random, SparseArrays, LinearAlgebra xbra = dagger(xket) ybra = dagger(yket) - @test 1e-13 > D(projector(xket) * xket, xket) - @test 1e-13 > D(xbra * projector(xket), xbra) - @test 1e-13 > D(projector(xbra) * xket, xket) - @test 1e-13 > D(xbra * projector(xbra), xbra) - @test 1e-13 > D(ybra * projector(yket, xbra), xbra) - @test 1e-13 > D(projector(yket, xbra) * xket, yket) + @test 1e-13 > D(projector(xket)*xket, xket) + @test 1e-13 > D(xbra*projector(xket), xbra) + @test 1e-13 > D(projector(xbra)*xket, xket) + @test 1e-13 > D(xbra*projector(xbra), xbra) + @test 1e-13 > D(ybra*projector(yket, xbra), xbra) + @test 1e-13 > D(projector(yket, xbra)*xket, yket) # Test operator exponential op = randoperator(b1a) - @test 1e-13 > D(op^2, op * op) - @test 1e-13 > D(op^3, op * op * op) - @test 1e-13 > D(op^4, op * op * op * op) + @test 1e-13 > D(op^2, op*op) + @test 1e-13 > D(op^3, op*op*op) + @test 1e-13 > D(op^4, op*op*op*op) # Test gemv b1 = GenericBasis(3) @@ -325,20 +304,20 @@ using Random, SparseArrays, LinearAlgebra beta = complex(0.3, 2.1) rket_ = deepcopy(rket) - OpenQuantumSystems.mul!(rket_, op, xket, complex(1.0), complex(0.0)) - @test 0 ≈ D(rket_, op * xket) + OpenQuantumSystems.mul!(rket_,op,xket,complex(1.0),complex(0.)) + @test 0 ≈ D(rket_, op*xket) rket_ = deepcopy(rket) - OpenQuantumSystems.mul!(rket_, op, xket, alpha, beta) - @test 1e-13 > D(rket_, alpha * op * xket + beta * rket) + OpenQuantumSystems.mul!(rket_,op,xket,alpha,beta) + @test 1e-13 > D(rket_, alpha*op*xket + beta*rket) rbra_ = deepcopy(rbra) - OpenQuantumSystems.mul!(rbra_, xbra, op, complex(1.0), complex(0.0)) - @test 0 ≈ D(rbra_, xbra * op) + OpenQuantumSystems.mul!(rbra_,xbra,op,complex(1.0),complex(0.)) + @test 0 ≈ D(rbra_, xbra*op) rbra_ = deepcopy(rbra) - OpenQuantumSystems.mul!(rbra_, xbra, op, alpha, beta) - @test 1e-13 > D(rbra_, alpha * xbra * op + beta * rbra) + OpenQuantumSystems.mul!(rbra_,xbra,op,alpha,beta) + @test 1e-13 > D(rbra_, alpha*xbra*op + beta*rbra) # # Test gemm b1 = GenericBasis(37) @@ -351,12 +330,12 @@ using Random, SparseArrays, LinearAlgebra beta = complex(0.3, 2.1) r_ = deepcopy(r) - OpenQuantumSystems.mul!(r_, op1, op2, complex(1.0), complex(0.0)) - @test 1e-13 > D(r_, op1 * op2) + OpenQuantumSystems.mul!(r_,op1,op2,complex(1.0),complex(0.)) + @test 1e-13 > D(r_, op1*op2) r_ = deepcopy(r) - OpenQuantumSystems.mul!(r_, op1, op2, alpha, beta) - @test 1e-10 > D(r_, alpha * op1 * op2 + beta * r) + OpenQuantumSystems.mul!(r_,op1,op2,alpha,beta) + @test 1e-10 > D(r_, alpha*op1*op2 + beta*r) dat = rand(prod(b_r.shape)) x = Ket(b_r, dat) @@ -364,19 +343,19 @@ using Random, SparseArrays, LinearAlgebra @test dm(x) == dm(y) # Test Hermitian - bspin = SpinBasis(1 // 2) + bspin = SpinBasis(1//2) bnlevel = NLevelBasis(2) @test ishermitian(DenseOperator(bspin, bspin, [1.0 im; -im 2.0])) == true @test ishermitian(DenseOperator(bspin, bnlevel, [1.0 im; -im 2.0])) == false # Test broadcasting op1_ = copy(op1) - op1 .= 2 * op1 + op1 .= 2*op1 @test op1 == op1_ .+ op1_ op1 .= op1_ @test op1 == op1_ op1 .= op1_ .+ 3 * op1_ - @test op1 == 4 * op1_ + @test op1 == 4*op1_ @test_throws DimensionMismatch op1 .= op2 bf = FockBasis(3) op3 = randoperator(bf)