Implement new basis interface

Author: akirakyle

CHANGELOG.md

@@ -2,5 +2,5 @@
 
 ## v1.0.0
 
-- First release, implementing necessary changes from reworking abstract types in `QuantumInterface`.
+- First release, implementing necessary changes to be compatible with the breaking release of `QuantumInterface` 0.4.0.
 

Project.toml

@@ -1,6 +1,6 @@
 name = "QuantumOpticsBase"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "0.6.0"
+version = "1.0.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/src/api.md

@@ -141,7 +141,7 @@ QuantumOpticsBase.check_samebases
 ```
 
 ```@docs
-@samebases
+@compatiblebases
 ```
 
 ```@docs
@@ -640,4 +640,4 @@ lazytensor_cachesize
 
 ```@docs
 lazytensor_clear_cache
-```
+```

src/QuantumOpticsBase.jl

@@ -10,8 +10,8 @@ import QuantumInterface: dagger, directsum, ⊕, dm, embed, nsubsystems, expect,
 # index helpers
 import QuantumInterface: complement, remove, shiftremove, reducedindices!, check_indices, check_sortedindices, check_embed_indices
 
-export Basis, GenericBasis, CompositeBasis, basis,
-        tensor, ⊗, permutesystems, @samebases,
+export Basis, GenericBasis, CompositeBasis, basis, basis_l, basis_r,
+        tensor, ⊗, permutesystems, @compatiblebases,
         #states
                 StateVector, Bra, Ket, basisstate, sparsebasisstate, norm,
                 dagger, normalize, normalize!,

src/bases.jl

@@ -1,3 +1,3 @@
 import QuantumInterface: Basis, basis, GenericBasis, CompositeBasis,
-    equal_shape, IncompatibleBases, @samebases, samebases, check_samebases,
-    multiplicable, check_multiplicable, reduced, ptrace, permutesystems
+    equal_shape, IncompatibleBases, @compatiblebases, samebases, check_samebases,
+    addible, check_addible, multiplicable, check_multiplicable, reduced, ptrace, permutesystems

src/charge.jl

@@ -28,6 +28,7 @@ struct ChargeBasis{T} <: Basis
 end
 
 Base.:(==)(b1::ChargeBasis, b2::ChargeBasis) = (b1.ncut == b2.ncut)
+Base.length(b::ChargeBasis) = b.dim
 
 """
     ShiftedChargeBasis(nmin, nmax) <: Basis
@@ -50,6 +51,7 @@ end
 
 Base.:(==)(b1::ShiftedChargeBasis, b2::ShiftedChargeBasis) =
     (b1.nmin == b2.nmin && b1.nmax == b2.nmax)
+Base.length(b::ShiftedChargeBasis) = b.dim
 
 """
     chargestate([T=ComplexF64,] b::ChargeBasis, n)

src/manybody.jl

@@ -50,6 +50,10 @@ end
 ManyBodyBasis(onebodybasis::B, occupations::O) where {B,O} = ManyBodyBasis{B,O}(onebodybasis, occupations)
 ManyBodyBasis(onebodybasis::B, occupations::Vector{T}) where {B,T} = ManyBodyBasis(onebodybasis, SortedVector(occupations))
 
+==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash == b2.occupations_hash && b1.onebodybasis == b2.onebodybasis
+Base.length(b::ManyBodyBasis) = length(b.occupations)
+
+
 allocate_buffer(occ) = similar(occ)
 allocate_buffer(mb::ManyBodyBasis) = allocate_buffer(first(mb.occupations))
 
@@ -89,8 +93,6 @@ bosonstates(T::Type, Nmodes::Int, Nparticles::Vector{Int}) = union((bosonstates(
 bosonstates(T::Type, onebodybasis::Basis, Nparticles) = bosonstates(T, length(onebodybasis), Nparticles)
 bosonstates(arg1, arg2) = bosonstates(OccupationNumbers{BosonStatistics,Int}, arg1, arg2)
 
-==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash == b2.occupations_hash && b1.onebodybasis == b2.onebodybasis
-
 """
     basisstate([T=ComplexF64,] mb::ManyBodyBasis, occupation::Vector)
 

src/metrics.jl

@@ -199,15 +199,15 @@ The `indices` argument can be a single integer or a collection of integers.
 function ptranspose(rho::DenseOpType{B,B}, indices=1) where B<:CompositeBasis
     # adapted from qutip.partial_transpose (https://qutip.org/docs/4.0.2/modules/qutip/partial_transpose.html)
     # works as long as QuantumOptics.jl doesn't change the implementation of `tensor`, i.e. tensor(a,b).data = kron(b.data,a.data)
-    nsys = length(rho.basis_l.shape)
+    nsys = nsubsystems(basis_l(rho))
     mask = ones(Int, nsys)
     mask[collect(indices)] .+= 1
     pt_dims = reshape(1:2*nsys, (nsys,2)) # indices of the operator viewed as a tensor with 2nsys legs
     pt_idx = [[pt_dims[i,mask[i]] for i = 1 : nsys]; [pt_dims[i,3-mask[i]] for i = 1 : nsys] ] # permute the legs on the subsystem of `indices`
     # reshape the operator data into a 2nsys-legged tensor and shape it back with the legs permuted
-    data = reshape(permutedims(reshape(rho.data, Tuple([rho.basis_l.shape; rho.basis_r.shape])), pt_idx), size(rho.data))
+    data = reshape(permutedims(reshape(rho.data, Tuple([size(basis_l(rho)); size(basis_r(rho))])), pt_idx), size(rho.data))
 
-    return DenseOperator(rho.basis_l,data)
+    return DenseOperator(basis_l(rho),data)
 
 end
 
@@ -276,7 +276,7 @@ function can be provided (for example to compute the entanglement-renyi entropy)
 """
 function entanglement_entropy(psi::Ket{B}, partition, entropy_fun=entropy_vn) where B<:CompositeBasis
     # check that sites are within the range
-    @assert all(partition .<= length(psi.basis.bases))
+    @assert all(partition .<= nsubsystems(psi.basis))
 
     rho = ptrace(psi, partition)
     return entropy_fun(rho)
@@ -285,17 +285,18 @@ end
 function entanglement_entropy(rho::DenseOpType{B,B}, partition, args...) where {B<:CompositeBasis}
     # check that sites is within the range
     hilb = rho.basis_l
-    all(partition .<= length(hilb.bases)) || throw(ArgumentError("Indices in partition must be within the bounds of the composite basis."))
-    length(partition) <= length(hilb.bases) || throw(ArgumentError("Partition cannot include the whole system."))
+    N = nsubsystems(hilb)
+    all(partition .<= N) || throw(ArgumentError("Indices in partition must be within the bounds of the composite basis."))
+    length(partition) <= N || throw(ArgumentError("Partition cannot include the whole system."))
 
     # build the doubled hilbert space for the vectorised dm, normalized like a Ket.
     b_doubled = hilb^2
     rho_vec = normalize!(Ket(b_doubled, vec(rho.data)))
 
     if partition isa Tuple
-        partition_ = tuple(partition..., (partition.+length(hilb.bases))...)
+        partition_ = tuple(partition..., (partition.+N)...)
     else
-        partition_ = vcat(partition, partition.+length(hilb.bases))
+        partition_ = vcat(partition, partition.+N)
     end
 
     return entanglement_entropy(rho_vec,partition_,args...)

src/operators.jl

@@ -9,8 +9,7 @@ Abstract type for operators with a data field.
 This is an abstract type for operators that have a direct matrix representation
 stored in their `.data` field.
 """
-abstract type BLROperator{BL,BR} <: AbstractOperator end
-abstract type DataOperator{BL,BR} <: BLROperator{BL,BR} end
+abstract type DataOperator{BL,BR} <: AbstractOperator end
 
 
 # Common error messages
@@ -24,14 +23,14 @@ Embed operator acting on a joint Hilbert space where missing indices are filled
 """
 function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
                indices, op::T) where T<:DataOperator
-    N = length(basis_l.bases)
-    @assert length(basis_r.bases) == N
+    N = nsubsystems(basis_l)
+    @assert nsubsystems(basis_r) == N
 
-    reduce(tensor, basis_l.bases[indices]) == op.basis_l || throw(IncompatibleBases())
-    reduce(tensor, basis_r.bases[indices]) == op.basis_r || throw(IncompatibleBases())
+    reduce(tensor, basis_l[indices]) == op.basis_l || throw(IncompatibleBases())
+    reduce(tensor, basis_r[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, eltype(op), basis_l.bases[i], basis_r.bases[i]) for i in index_order]
+    index_order = [idx for idx in 1:N if idx ∉ indices]
+    all_operators = AbstractOperator[identityoperator(T, eltype(op), basis_l[i], basis_r[i]) for i in index_order]
 
     for idx in indices
         pushfirst!(index_order, idx)
@@ -45,8 +44,8 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
 
     # Reorient the matrix to act in the correctly ordered basis.
     # Get the dimensions necessary for index permuting.
-    dims_l = [b.shape[1] for b in basis_l.bases]
-    dims_r = [b.shape[1] for b in basis_r.bases]
+    dims_l = size(basis_l)
+    dims_r = size(basis_r)
 
     # Get the order of indices to use in the first reshape. Julia indices go in
     # reverse order.
@@ -74,12 +73,12 @@ end
 function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
                 index::Integer, op::T) where T<:DataOperator
 
-    N = length(basis_l.bases)
+    N = nsubsystems(basis_l)
 
     # Check stuff
-    @assert N==length(basis_r.bases)
-    basis_l.bases[index] == op.basis_l || throw(IncompatibleBases())
-    basis_r.bases[index] == op.basis_r || throw(IncompatibleBases())
+    @assert nsubsystems(basis_r) == N
+    basis_l[index] == op.basis_l || throw(IncompatibleBases())
+    basis_r[index] == op.basis_r || throw(IncompatibleBases())
     check_indices(N, index)
 
     # Build data
@@ -92,8 +91,8 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
             data = kron(data, op.data)
             i -= length(index)
         else
-            bl = basis_l.bases[i]
-            br = basis_r.bases[i]
+            bl = basis_l[i]
+            br = basis_r[i]
             id = SparseMatrixCSC{Tnum}(I, length(bl), length(br))
             data = kron(data, id)
             i -= 1
@@ -110,18 +109,21 @@ Expectation value of the given operator `op` for the specified `state`.
 
 `state` can either be a (density) operator or a ket.
 """
-expect(op::BLROperator{B,B}, state::Ket{B}) where B = dot(state.data, (op * state).data)
+function expect(op::AbstractOperator, state::Ket)
+    check_multiplicable(op,op); check_multiplicable(op,state)
+    dot(state.data, (op * state).data)
+end
 
 # TODO upstream this one
 # expect(op::AbstractOperator{B,B}, state::AbstractKet{B}) where B = norm(op * state) ^ 2
 
-function expect(indices, op::BLROperator{B,B}, state::Ket{B2}) where {B,B2<:CompositeBasis}
-    N = length(state.basis.shape)
+function expect(indices, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis}
+    N = nsubsystems(basis(state))
     indices_ = complement(N, indices)
     expect(op, ptrace(state, indices_))
 end
 
-expect(index::Integer, op::BLROperator{B,B}, state::Ket{B2}) where {B,B2<:CompositeBasis} = expect([index], op, state)
+expect(index::Integer, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis} = expect([index], op, state)
 
 """
     variance(op, state)
@@ -130,44 +132,42 @@ Variance of the given operator `op` for the specified `state`.
 
 `state` can either be a (density) operator or a ket.
 """
-function variance(op::BLROperator{B,B}, state::Ket{B}) where B
+function variance(op::AbstractOperator, state::Ket)
+    check_multiplicable(op,op); check_multiplicable(op,state)
     x = op*state
     state.data'*(op*x).data - (state.data'*x.data)^2
 end
 
-function variance(indices, op::BLROperator{B,B}, state::Ket{BC}) where {B,BC<:CompositeBasis}
-    N = length(state.basis.shape)
+function variance(indices, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis}
+    N = nsubsystems(basis(state))
     indices_ = complement(N, indices)
     variance(op, ptrace(state, indices_))
 end
 
-variance(index::Integer, op::BLROperator{B,B}, state::Ket{BC}) where {B,BC<:CompositeBasis} = variance([index], op, state)
+variance(index::Integer, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis} = variance([index], op, state)
 
 # Helper functions to check validity of arguments
 function check_ptrace_arguments(a::AbstractOperator, indices)
     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
-    rank = length(a.basis_l.shape)
-    if rank != length(a.basis_r.shape)
+    rank = nsubsystems(basis_l(a))
+    if rank != nsubsystems(basis_r(a))
         throw(ArgumentError("Partial trace can only be applied onto operators wich have the same number of subsystems in the left basis and right basis."))
     end
     if rank == length(indices)
         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)
+    check_indices(nsubsystems(basis_l(a)), indices)
     for i=indices
-        if a.basis_l.shape[i] != a.basis_r.shape[i]
+        if size(basis_l(a))[i] != size(basis_r(a))[i]
             throw(ArgumentError("Partial trace can only be applied onto subsystems that have the same left and right dimension."))
         end
     end
 end
 function check_ptrace_arguments(a::StateVector, indices)
-    if length(basis(a).shape) == length(indices)
+    if nsubsystems(basis(a)) == length(indices)
         throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use tr() instead."))
     end
-    check_indices(length(basis(a).shape), indices)
+    check_indices(nsubsystems(basis(a)), indices)
 end
-
-multiplicable(a::AbstractOperator, b::Ket) = multiplicable(a.basis_r, b.basis)
-multiplicable(a::Bra, b::AbstractOperator) = multiplicable(a.basis, b.basis_l)