Merge branch 'master' into manybody

.github/workflows/breakage.yml

@@ -20,7 +20,7 @@ jobs:
         pkgversion: [latest]
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       # Install Julia
       - uses: julia-actions/setup-julia@v1

.github/workflows/ci-jet.yml

@@ -17,7 +17,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}

.github/workflows/ci.yml

@@ -22,8 +22,12 @@ jobs:
         exclude:
           - os: macOS-latest
             arch: x86
+        include:
+          - os: ubuntu-latest
+            arch: x64
+            version: '1.6'
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}

.github/workflows/docs.yml

@@ -13,7 +13,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'

.github/workflows/invalidations.yml

@@ -19,12 +19,12 @@ jobs:
     - uses: julia-actions/setup-julia@v1
       with:
         version: '1'
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_pr
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         ref: ${{ github.event.repository.default_branch }}
     - uses: julia-actions/julia-buildpkg@v1

Project.toml

@@ -6,11 +6,13 @@ version = "0.4.18"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RandomMatrices = "2576dda1-a324-5b11-aa66-c48ed7e3c618"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
 UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
@@ -19,12 +21,14 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 Adapt = "1, 2, 3.3"
 FFTW = "1.2"
 FastExpm = "1.1.0"
+FastGaussQuadrature = "0.5"
 FillArrays = "0.13, 1"
 LRUCache = "1"
-QuantumInterface = "0.3.0"
+QuantumInterface = "0.3.3"
+RandomMatrices = "0.5"
 Strided = "1, 2"
 UnsafeArrays = "1"
-julia = "1.3"
+julia = "1.6"
 
 [extras]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -37,9 +41,11 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 
 [targets]
-test = ["LinearAlgebra", "SparseArrays", "Random", "Test", "Aqua", "JET", "Adapt", "Dates", "FFTW", "LRUCache", "Strided", "UnsafeArrays", "FillArrays"]
+test = ["LinearAlgebra", "SparseArrays", "Random", "Test", "Aqua", "JET", "Adapt", "Dates", "FFTW", "LRUCache", "Strided", "StridedViews", "UnsafeArrays", "FillArrays", "StatsBase"]

docs/src/api.md

@@ -594,3 +594,21 @@ DenseChiMatrix
 ```@docs
 QuantumOpticsBase.set_printing
 ```
+
+## [LazyTensor functions](@id API: LazyTensor functions)
+
+```@docs
+lazytensor_enable_cache
+```
+
+```@docs
+lazytensor_disable_cache
+```
+
+```@docs
+lazytensor_cachesize
+```
+
+```@docs
+lazytensor_clear_cache
+```

src/QuantumOpticsBase.jl

@@ -3,7 +3,7 @@ module QuantumOpticsBase
 using SparseArrays, LinearAlgebra, LRUCache, Strided, UnsafeArrays, FillArrays
 import LinearAlgebra: mul!, rmul!
 
-import QuantumInterface: dagger, directsum, ⊕, dm, embed, expect, identityoperator, identitysuperoperator,
+import QuantumInterface: dagger, directsum, ⊕, dm, embed, nsubsystems, expect, identityoperator, identitysuperoperator,
         permutesystems, projector, ptrace, reduced, tensor, ⊗, variance, apply!, basis, AbstractSuperOperator
 
 # index helpers
@@ -41,6 +41,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
                 displace, displace_analytical, displace_analytical!,
                 squeeze,
         randstate, randoperator, thermalstate, coherentthermalstate, phase_average, passive_state,
+        randstate_haar, randunitary_haar,
         #spin
                 SpinBasis, sigmax, sigmay, sigmaz, sigmap, sigmam, spinup, spindown,
         #subspace
@@ -65,7 +66,6 @@ export Basis, GenericBasis, CompositeBasis, basis,
         #apply
                 apply!
 
-include("polynomials.jl")
 include("bases.jl")
 include("states.jl")
 include("operators.jl")

src/apply.jl

@@ -1,11 +1,24 @@
+function is_apply_shortcircuit(state, indices, operation)
+    if nsubsystems(state) == 1
+        basis(state)==basis(operation) || throw(ArgumentError("`apply!` failed due to incompatible bases of the state and the operation attempted to be applied on it"))
+    end
+    basis(state)==basis(operation) || return false
+    j = 1
+    for i in indices
+        i == j || return false
+        j+=1
+    end
+    return j-1 == length(indices)
+end
+
 function apply!(state::Ket, indices, operation::Operator)
-    op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation)
+    op = is_apply_shortcircuit(state, indices, operation) ? operation : embed(basis(state), indices, operation)
     state.data = (op*state).data
     state
 end
 
 function apply!(state::Operator, indices, operation::Operator)
-    op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation)
+    op = is_apply_shortcircuit(state, indices, operation) ? operation : embed(basis(state), indices, operation)
     state.data = (op*state*op').data
     state
 end
@@ -15,7 +28,10 @@ function apply!(state::Ket, indices, operation::T) where {T<:AbstractSuperOperat
 end
 
 function apply!(state::Operator, indices, operation::T) where {T<:AbstractSuperOperator}
-    op = basis(state)==basis(operation) ? operation : embed(basis(state), indices, operation)
-    state.data = (op*state).data
-    state
+    if is_apply_shortcircuit(state, indices, operation)
+        state.data = (operation*state).data
+        return state
+    else
+        error("`apply!` does not yet support embedding superoperators acting only on a subsystem of the given state")
+    end
 end

src/operators_lazytensor.jl

@@ -96,7 +96,7 @@ isequal(x::LazyTensor, y::LazyTensor) = samebases(x,y) && isequal(x.indices, y.i
 # Arithmetic operations
 -(a::LazyTensor) = LazyTensor(a, -a.factor)
 
-const single_dataoperator{B1,B2} = LazyTensor{B1,B2,F,I,Tuple{T}} where {B1,B2,F,I,T<:DataOperator} 
+const single_dataoperator{B1,B2} = LazyTensor{B1,B2,F,I,Tuple{T}} where {B1,B2,F,I,T<:DataOperator}
 function +(a::T1,b::T2) where {T1 <: single_dataoperator{B1,B2},T2 <: single_dataoperator{B1,B2}} where {B1,B2}
     if a.indices == b.indices
         op = a.operators[1] * a.factor + b.operators[1] * b.factor
@@ -349,7 +349,7 @@ function _tp_matmul_mid!(result, a::AbstractMatrix, loc::Integer, b, α::Number,
     result_r = Base.ReshapedArray(result, (sz_b_1, size(a, 1), sz_b_3), ())
 
     if a isa FillArrays.Eye
-        # Square Eyes are skipped higher up. This handles the non-square case. 
+        # Square Eyes are skipped higher up. This handles the non-square case.
         size(b, loc) == size(a, 2) && size(result, loc) == size(a, 1) || throw(DimensionMismatch("Dimensions of Eye matrix do not match subspace dimensions."))
         d = min(size(a)...)
 
@@ -505,7 +505,7 @@ function _explicit_isometries(::Type{T}, used_indices, bl::Basis, br::Basis, shi
                 iso_inds = [i + shift]
             else
                 push!(isos, Eye{T}(sl, sr))
-                push!(iso_inds, i + shift)
+                push!(iso_inds::Vector{Int}, i + shift) # Help with inference (detected by JET)
             end
         end
     end
@@ -752,4 +752,3 @@ _mul_puresparse!(result::DenseOpType{B1,B3},h::LazyTensor{B1,B2,F,I,T},op::Dense
 _mul_puresparse!(result::DenseOpType{B1,B3},op::DenseOpType{B1,B2},h::LazyTensor{B2,B3,F,I,T},alpha,beta) where {B1,B2,B3,F,I,T} = (_gemm_puresparse(alpha, op.data, h, beta, result.data); result)
 _mul_puresparse!(result::Ket{B1},a::LazyTensor{B1,B2,F,I,T},b::Ket{B2},alpha,beta) where {B1,B2,F,I,T} = (_gemm_puresparse(alpha, a, b.data, beta, result.data); result)
 _mul_puresparse!(result::Bra{B2},a::Bra{B1},b::LazyTensor{B1,B2,F,I,T},alpha,beta) where {B1,B2,F,I,T} = (_gemm_puresparse(alpha, a.data, b, beta, result.data); result)
-

src/operators_sparse.jl

@@ -6,14 +6,12 @@ const SparseOpPureType{BL,BR} = Operator{BL,BR,<:SparseMatrixCSC}
 const SparseOpAdjType{BL,BR} = Operator{BL,BR,<:Adjoint{<:Number,<:SparseMatrixCSC}}
 const SparseOpType{BL,BR} = Union{SparseOpPureType{BL,BR},SparseOpAdjType{BL,BR}}
 
-
 """
     SparseOperator(b1[, b2, data])
 
 Sparse array implementation of Operator.
 
-The matrix is stored as the julia built-in type `SparseMatrixCSC`
-in the `data` field.
+The matrix is stored as the julia built-in type `SparseMatrixCSC` in the `data` field.
 """
 SparseOperator(b1::Basis, b2::Basis, data) = Operator(b1, b2, SparseMatrixCSC(data))
 SparseOperator(b1::Basis, b2::Basis, data::SparseMatrixCSC) = Operator(b1, b2, data)
@@ -58,7 +56,11 @@ If you only need the result of the exponential acting on a vector,
 consider using much faster implicit methods that do not calculate the entire exponential.
 """
 function exp(op::T; opts...) where {B,T<:SparseOpType{B,B}}
-    return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+    if iszero(op)
+        return identityoperator(op)
+    else
+        return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+    end
 end
 
 function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis}

src/sparsematrix.jl

@@ -246,7 +246,7 @@ function _permutedims(x::AbstractSparseMatrix, shape, perm) # TODO upstream as `
     shape = (shape...,)
     shape_perm = ([shape[i] for i in perm]...,)
     y = spzeros(eltype(x), x.m, x.n)
-    for I in eachindex(x)
+    for I in eachindex(x)::CartesianIndices # Help with inference (detected by JET)
         linear_index = LinearIndices((x.m, x.n))[I.I...]
         cartesian_index = CartesianIndices(shape)[linear_index]
         cartesian_index_perm = [cartesian_index[i] for i=perm]

src/spin.jl

@@ -83,3 +83,20 @@ Spin down state for the given Spin basis.
 """
 spindown(::Type{T}, b::SpinBasis) where T = basisstate(T, b, b.shape[1])
 spindown(b::SpinBasis) = spindown(ComplexF64, b)
+
+
+"""
+    squeeze([T=ComplexF64,] b::SpinBasis, z)
+
+Squeezing operator ``S(z)=\\exp{\\left(\\frac{z^*\\hat{J_-}^2 - z\\hat{J}_+}{2 N}\\right)}`` for the 
+specified Spin-``N/2`` basis with optional data type `T`, computed as the matrix exponential. Too
+large squeezing (``|z| > \\sqrt{N}``) will create an oversqueezed state.
+"""
+function squeeze(::Type{T}, b::SpinBasis, z::Number) where T
+    N = Int(length(b)-1)
+    z = T(z)/2N
+    Jm = sigmam(b)/2
+    Jp = sigmap(b)/2
+    exp(conj(z)*dense(Jm)^2 - z*dense(Jp)^2)
+end
+squeeze(b::SpinBasis, z::T) where {T <: Number} = squeeze(ComplexF64, b, z)