Add lazy superoperator types and measurement utility functions

src/QuantumOpticsBase.jl

@@ -21,6 +21,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
         #operators
                 AbstractOperator, DataOperator, expect, variance,
                 identityoperator, ptrace, reduced, embed, dense, tr, sparse,
+                _drop_singular_bases, _branch_prob, _overlap, _project_and_drop,
         #operators_dense
                 Operator, DenseOperator, DenseOpType, projector, dm,
         #operators_sparse
@@ -41,6 +42,8 @@ export Basis, GenericBasis, CompositeBasis, basis,
                 SuperOperator, DenseSuperOperator, DenseSuperOpType,
                 SparseSuperOperator, SparseSuperOpType, spre, spost, sprepost, liouvillian,
                 identitysuperoperator,
+        #operators_lazysuperoperator
+                AbstractLazySuperOperator, LazyPrePost, LazySuperSum, LazySuperTensor,
         #fock
                 FockBasis, number, destroy, create,
                 fockstate, coherentstate, coherentstate!,
@@ -78,6 +81,7 @@ include("bases.jl")
 include("states.jl")
 include("operators.jl")
 include("operators_dense.jl")
+include("operators_project_drop.jl")
 include("sparsematrix.jl")
 include("operators_sparse.jl")
 include("operators_lazysum.jl")
@@ -86,6 +90,7 @@ include("operators_lazytensor.jl")
 include("time_dependent_operator.jl")
 include("states_lazyket.jl")
 include("superoperators.jl")
+include("operators_lazysuperoperator.jl")
 include("spin.jl")
 include("fock.jl")
 include("charge.jl")

src/operators.jl

@@ -111,8 +111,7 @@ Expectation value of the given operator `op` for the specified `state`.
 """
 expect(op::AbstractOperator{B,B}, state::Ket{B}) where B = dot(state.data, (op * state).data)
 
-# TODO upstream this one
-# expect(op::AbstractOperator{B,B}, state::AbstractKet{B}) where B = norm(op * state) ^ 2
+expect(op::AbstractOperator{B,B}, state::AbstractKet{B}) where B = norm(op * state) ^ 2
 
 function expect(indices, op::AbstractOperator{B,B}, state::Ket{B2}) where {B,B2<:CompositeBasis}
     N = length(state.basis.shape)
@@ -170,3 +169,4 @@ end
 
 multiplicable(a::AbstractOperator, b::Ket) = multiplicable(a.basis_r, b.basis)
 multiplicable(a::Bra, b::AbstractOperator) = multiplicable(a.basis, b.basis_l)
+

src/operators_lazysuperoperator.jl

@@ -0,0 +1,120 @@
+using QuantumInterface: AbstractSuperOperator
+
+"""
+    AbstractLazySuperOperator{B1,B2} <: AbstractSuperOperator{B1,B2}
+
+Base type for lazy superoperator implementations that compute their action
+on operators without explicitly storing the full superoperator matrix.
+"""
+abstract type AbstractLazySuperOperator{B1,B2} <: AbstractSuperOperator{B1,B2} end
+
+"""
+    LazyPrePost{B,DT} <: AbstractLazySuperOperator{Tuple{B,B},Tuple{B,B}}
+
+Lazy superoperator that applies a pre-operator and post-operator to a quantum state/operator.
+For an operator `ρ`, this computes `preop * ρ * postop'`.
+
+# Arguments
+- `preop::Operator{B,B,DT}`: Operator applied from the left
+- `postop::Operator{B,B,DT}`: Operator applied from the right (conjugate transposed)
+
+# Example
+```julia
+# Create Lindblad-type superoperator: L ρ L† - (1/2){L†L, ρ}
+jump_op = sigmax()
+dissipator = LazyPrePost(jump_op, jump_op)
+evolved_state = dissipator * initial_state
+```
+"""
+struct LazyPrePost{B,DT} <: AbstractLazySuperOperator{Tuple{B,B},Tuple{B,B}}
+    preop::Operator{B,B,DT}
+    postop::Operator{B,B,DT}
+end
+
+function LazyPrePost(preop::T,postop::T) where {B,DT,T<:Operator{B,B,DT}}
+    LazyPrePost{B,DT}(preop,postop)
+end
+
+"""
+    LazySuperSum{B,F,T} <: AbstractLazySuperOperator{Tuple{B,B},Tuple{B,B}}
+
+Lazy sum of superoperators with corresponding factors.
+Computes `sum(factor[i] * sop[i] * ρ for i in 1:length(sops))`.
+
+# Fields
+- `basis::B`: Quantum basis for the superoperator
+- `factors::F`: Vector of scaling factors
+- `sops::T`: Vector of superoperators to sum
+"""
+struct LazySuperSum{B,F,T} <: AbstractLazySuperOperator{Tuple{B,B},Tuple{B,B}}
+    basis::B
+    factors::F
+    sops::T
+end
+
+"""
+    LazySuperTensor{B,T} <: AbstractLazySuperOperator{Tuple{B,B},Tuple{B,B}}
+
+Lazy tensor product of superoperators.
+For composite quantum systems, applies each superoperator to its respective subsystem.
+
+# Fields  
+- `basis::B`: Composite quantum basis
+- `sops::T`: Vector of superoperators for each subsystem
+"""
+struct LazySuperTensor{B,T} <: AbstractLazySuperOperator{Tuple{B,B},Tuple{B,B}}
+    basis::B
+    sops::T
+end
+
+# Basis functions
+QuantumOpticsBase.basis(sop::LazyPrePost) = basis(sop.preop)
+QuantumOpticsBase.basis(sop::LazySuperSum) = sop.basis
+
+# Embed functions  
+QuantumOpticsBase.embed(bl,br,index,op::LazyPrePost) = 
+    LazyPrePost(embed(bl,br,index,op.preop), embed(bl,br,index,op.postop))
+
+QuantumOpticsBase.embed(bl,br,index,op::LazySuperSum) = 
+    LazySuperSum(bl, op.factors, [embed(bl,br,index,o) for o in op.sops])
+
+# Application to operators
+function Base.:(*)(sop::LazyPrePost, op::Operator)
+    # Apply preop from left and postop from right (conjugated)
+    # TODO: Optimize to avoid creating intermediate spre/spost objects
+    # TODO: Implement in-place version with pre-allocated buffers
+    result = op
+    result = spre(sop.preop) * result  # Left multiplication
+    result = spost(sop.postop) * result  # Right multiplication
+    result
+end
+
+function Base.:(*)(ssop::LazySuperSum, op::Operator)
+    result = zero(op)
+    for (factor, sop) in zip(ssop.factors, ssop.sops)
+        result += factor * (sop * op)
+    end
+    result
+end
+
+function Base.:(*)(ssop::LazySuperTensor, op::Operator)
+    result = op
+    for sop in ssop.sops
+        result = sop * result
+    end
+    result
+end
+
+# Composition of lazy superoperators
+Base.:(*)(l::LazyPrePost, r::LazyPrePost) = 
+    LazyPrePost(l.preop * r.preop, r.postop * l.postop)
+
+Base.:(+)(ops::LazyPrePost...) = 
+    LazySuperSum(basis(first(ops)), fill(1, length(ops)), ops)
+
+# Tensor product of superoperators
+function QuantumInterface.tensor(sops::AbstractSuperOperator...)
+    b = QuantumInterface.tensor(basis.(sops)...)
+    @assert length(sops) == length(b.bases) "tensor products of superoperators over composite bases are not implemented yet"
+    LazySuperTensor(b, [embed(b, b, i, s) for (i,s) in enumerate(sops)])
+end

src/operators_project_drop.jl

@@ -0,0 +1,83 @@
+# Project and drop functions for quantum measurements and post-selection
+
+"""
+    _drop_singular_bases(state)
+
+Remove singular (dimension-1) bases from a quantum state or operator.
+This is useful after measurement or post-selection operations that reduce
+certain subsystems to classical states.
+
+# Examples
+```julia
+# After measuring a subsystem in a definite state
+reduced_ket = _drop_singular_bases(post_measurement_ket)
+reduced_op = _drop_singular_bases(post_measurement_operator)
+```
+"""
+function _drop_singular_bases(ket::Ket)
+    b = tensor([b for b in basis(ket).bases if length(b)>1]...)
+    return Ket(b, ket.data)
+end
+
+function _drop_singular_bases(op::Operator)
+    b = tensor([b for b in basis(op).bases if length(b)>1]...)
+    return Operator(b, op.data)
+end
+
+"""
+    _branch_prob(state)
+
+Calculate the probability of a quantum branch/measurement outcome.
+Returns the square of the norm for kets, or the trace for density operators.
+"""
+_branch_prob(psi::Ket) = norm(psi)^2
+_branch_prob(op::Operator) = real(sum((op.data[i, i] for i in 1:size(op.data,1))))
+
+"""
+    _overlap(left, right)
+
+Calculate the overlap between two quantum states/operators.
+For kets: |⟨left|right⟩|²
+For ket-operator pairs: ⟨left|right|left⟩
+"""
+_overlap(l::Ket, r::Ket) = abs2(l'*r)
+_overlap(l::Ket, op::Operator) = real(l'*op*l)
+
+"""
+    _project_and_drop(state, project_on, basis_index)
+
+Project a quantum state onto a specific state in one subsystem and remove
+that subsystem from the composite state. This is useful for conditional
+evolution after measurement.
+
+# Arguments
+- `state`: Quantum state (Ket or Operator) to project
+- `project_on`: State to project onto
+- `basis_index`: Index of the subsystem to project and remove
+
+# Returns
+Reduced state with the projected subsystem removed.
+"""
+function _project_and_drop(state::Ket, project_on, basis_index)
+    singularbasis = GenericBasis(1)
+    singularket = basisstate(singularbasis,1)
+    proj = projector(singularket, project_on')
+    basis_r = collect(Any,basis(state).bases)
+    basis_l = copy(basis_r)
+    basis_l[basis_index] = singularbasis
+    emproj = embed(tensor(basis_l...),tensor(basis_r...),basis_index,proj)
+    result = emproj*state
+    return _drop_singular_bases(result)
+end
+
+function _project_and_drop(state::Operator, project_on, basis_index)
+    singularbasis = GenericBasis(1)
+    singularket = basisstate(singularbasis,1)
+    proj = projector(singularket, project_on')
+    basis_r = collect(Any,basis(state).bases)
+    basis_l = copy(basis_r)
+    basis_l[basis_index] = singularbasis
+    emproj = embed(tensor(basis_l...),tensor(basis_r...),basis_index,proj)
+    result = emproj*state*emproj'
+    return _drop_singular_bases(result)
+end