Skip to content

Commit

Permalink
[WIP] Improve generation performance in manybody.jl (#121)
Browse files Browse the repository at this point in the history 
* Improve generation performance in manybody.jl

---------

Co-authored-by: Stefan Krastanov <[email protected]>
  • Loading branch information
@aryavorskiy @Krastanov
aryavorskiy and Krastanov authored Jul 16, 2023
1 parent 1d4ad57 commit 8f16def
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 56 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "QuantumOpticsBase"
uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
version = "0.4.9"
version = "0.4.10"

[deps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
Expand Down
110 changes: 55 additions & 55 deletions src/manybody.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,9 @@ end
Creation operator for the i-th mode of the many-body basis `b`.
"""
function create(::Type{T}, b::ManyBodyBasis, index) where T
result = SparseOperator(T, b)
Is = Int[]
Js = Int[]
Vs = T[]
# <{m}_i| at |{m}_j>
for i=1:length(b)
occ_i = b.occupations[i]
Expand All @@ -90,11 +92,13 @@ function create(::Type{T}, b::ManyBodyBasis, index) where T
end
for j=1:length(b)
if isnonzero(occ_i, b.occupations[j], index)
result.data[i, j] = sqrt(occ_i[index])
push!(Is, i)
push!(Js, j)
push!(Vs, sqrt(occ_i[index]))
end
end
end
result
return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
end
create(b::ManyBodyBasis, index) = create(ComplexF64, b, index)

Expand All @@ -104,7 +108,9 @@ create(b::ManyBodyBasis, index) = create(ComplexF64, b, index)
Annihilation operator for the i-th mode of the many-body basis `b`.
"""
function destroy(::Type{T}, b::ManyBodyBasis, index) where T
result = SparseOperator(T, b)
Is = Int[]
Js = Int[]
Vs = T[]
# <{m}_j| a |{m}_i>
for i=1:length(b)
occ_i = b.occupations[i]
Expand All @@ -113,11 +119,13 @@ function destroy(::Type{T}, b::ManyBodyBasis, index) where T
end
for j=1:length(b)
if isnonzero(occ_i, b.occupations[j], index)
result.data[j, i] = sqrt(occ_i[index])
push!(Is, j)
push!(Js, i)
push!(Vs, sqrt(occ_i[index]))
end
end
end
result
return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
end
destroy(b::ManyBodyBasis, index) = destroy(ComplexF64, b, index)

Expand All @@ -127,11 +135,7 @@ destroy(b::ManyBodyBasis, index) = destroy(ComplexF64, b, index)
Particle number operator for the i-th mode of the many-body basis `b`.
"""
function number(::Type{T}, b::ManyBodyBasis, index) where T
result = SparseOperator(T, b)
for i=1:length(b)
result.data[i, i] = b.occupations[i][index]
end
result
diagonaloperator(b, T[b.occupations[i][index] for i in 1:length(b)])
end
number(b::ManyBodyBasis, index) = number(ComplexF64, b, index)

Expand All @@ -141,11 +145,7 @@ number(b::ManyBodyBasis, index) = number(ComplexF64, b, index)
Total particle number operator.
"""
function number(::Type{T}, b::ManyBodyBasis) where T
result = SparseOperator(T, b)
for i=1:length(b)
result.data[i, i] = sum(b.occupations[i])
end
result
diagonaloperator(b, T[sum(b.occupations[i]) for i in 1:length(b)])
end
number(b::ManyBodyBasis) = number(ComplexF64, b)

Expand Down Expand Up @@ -178,7 +178,9 @@ end
Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes.
"""
function transition(::Type{T}, b::ManyBodyBasis, to, from) where T
result = SparseOperator(T, b)
Is = Int[]
Js = Int[]
Vs = T[]
# <{m}_j| at_to a_from |{m}_i>
for i=1:length(b)
occ_i = b.occupations[i]
Expand All @@ -188,11 +190,13 @@ function transition(::Type{T}, b::ManyBodyBasis, to, from) where T
for j=1:length(b)
occ_j = b.occupations[j]
if isnonzero(occ_j, occ_i, to, from)
result.data[j, i] = sqrt(occ_i[from])*sqrt(occ_j[to])
push!(Is, j)
push!(Js, i)
push!(Vs, sqrt(occ_i[from]) * sqrt(occ_j[to]))
end
end
end
result
return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
end
transition(b::ManyBodyBasis, to, from) = transition(ComplexF64, b, to, from)

Expand Down Expand Up @@ -240,7 +244,7 @@ function manybodyoperator_1(basis::ManyBodyBasis, op::Operator)
result = DenseOperator(basis)
@inbounds for n=1:N, m=1:N
for j=1:S, i=1:S
C = coefficient(basis.occupations[m], basis.occupations[n], [i], [j])
C = coefficient(basis.occupations[m], basis.occupations[n], i, j)
if C != 0.
result.data[m,n] += C*op.data[i,j]
end
Expand All @@ -252,22 +256,25 @@ manybodyoperator_1(basis::ManyBodyBasis, op::AdjointOperator) = dagger(manybodyo

function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOpPureType)
N = length(basis)
S = length(basis.onebodybasis)
result = SparseOperator(basis)
M = op.data
Is = Int[]
Js = Int[]
Vs = ComplexF64[]
@inbounds for colindex = 1:M.n
for i=M.colptr[colindex]:M.colptr[colindex+1]-1
row = M.rowval[i]
value = M.nzval[i]
for m=1:N, n=1:N
C = coefficient(basis.occupations[m], basis.occupations[n], [row], [colindex])
C = coefficient(basis.occupations[m], basis.occupations[n], row, colindex)
if C != 0.
result.data[m, n] += C*value
push!(Is, m)
push!(Js, n)
push!(Vs, C * value)
end
end
end
end
return result
return SparseOperator(basis, sparse(Is, Js, Vs, N, N))
end

function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
Expand All @@ -280,7 +287,7 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::Operator)
occupations = basis.occupations
@inbounds for m=1:N, n=1:N
for l=1:S, k=1:S, j=1:S, i=1:S
C = coefficient(occupations[m], occupations[n], [i, j], [k, l])
C = coefficient(occupations[m], occupations[n], (i, j), (k, l))
result.data[m,n] += C*op_data[i, j, k, l]
end
end
Expand All @@ -290,7 +297,9 @@ end
function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType)
N = length(basis)
S = length(basis.onebodybasis)
result = SparseOperator(basis)
Is = Int[]
Js = Int[]
Vs = ComplexF64[]
occupations = basis.occupations
rows = rowvals(op.data)
values = nonzeros(op.data)
Expand All @@ -302,11 +311,13 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType)
index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2 + row])
C = coefficient(occupations[m], occupations[n], index[1:2], index[3:4])
if C!=0.
result.data[m,n] += C*value
push!(Is, m)
push!(Js, n)
push!(Vs, C * value)
end
end
end
return result
return SparseOperator(basis, sparse(Is, Js, Vs, N, N))
end


Expand Down Expand Up @@ -354,7 +365,7 @@ function onebodyexpect_1(op::Operator, state::Ket)
for m=1:N, n=1:N
value = conj(state.data[m])*state.data[n]
for i=1:S, j=1:S
C = coefficient(occupations[m], occupations[n], [i], [j])
C = coefficient(occupations[m], occupations[n], i, j)
if C != 0.
result += C*op.data[i,j]*value
end
Expand All @@ -371,7 +382,7 @@ function onebodyexpect_1(op::Operator, state::Operator)
@inbounds for s=1:N, t=1:N
value = state.data[t,s]
for i=1:S, j=1:S
C = coefficient(occupations[s], occupations[t], [i], [j])
C = coefficient(occupations[s], occupations[t], i, j)
if !iszero(C)
result += C*op.data[i,j]*value
end
Expand All @@ -391,7 +402,7 @@ function onebodyexpect_1(op::SparseOpPureType, state::Ket)
row = M.rowval[i]
value = M.nzval[i]
for m=1:N, n=1:N
C = coefficient(occupations[m], occupations[n], [row], [colindex])
C = coefficient(occupations[m], occupations[n], row, colindex)
if C != 0.
result += C*value*conj(state.data[m])*state.data[n]
end
Expand All @@ -412,7 +423,7 @@ function onebodyexpect_1(op::SparseOpPureType, state::Operator)
row = M.rowval[i]
value = M.nzval[i]
for s=1:N, t=1:N
C = coefficient(occupations[s], occupations[t], [row], [colindex])
C = coefficient(occupations[s], occupations[t], row, colindex)
if C != 0.
result += C*value*state.data[t,s]
end
Expand All @@ -426,29 +437,18 @@ end
"""
Calculate the matrix element <{m}|at_1...at_n a_1...a_n|{n}>.
"""
function coefficient(occ_m, occ_n, at_indices, a_indices)
occ_m = copy(occ_m)
occ_n = copy(occ_n)
C = 1.
for i=at_indices
if occ_m[i] == 0
return 0.
end
C *= sqrt(occ_m[i])
occ_m[i] -= 1
end
for i=a_indices
if occ_n[i] == 0
return 0.
end
C *= sqrt(occ_n[i])
occ_n[i] -= 1
end
if occ_m == occ_n
return C
else
return 0.
Base.@propagate_inbounds function coefficient(occ_m, occ_n, at_indices, a_indices)
any(==(0), (occ_m[m] for m in at_indices)) && return 0.
any(==(0), (occ_n[n] for n in a_indices)) && return 0.
C = prod(√, (occ_m[m] for m in at_indices)) * prod(√, (occ_n[n] for n in a_indices))
for i in 1:length(occ_m)
vm = occ_m[i]
vn = occ_n[i]
i in at_indices && (vm -= 1)
i in a_indices && (vn -= 1)
vm != vn && return zero(C)
end
return C
end

function _distribute_bosons(Nparticles, Nmodes, index, occupations, results)
Expand Down

2 comments on commit 8f16def

@Krastanov
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/87598

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.4.10 -m "<description of version>" 8f16def2df329249ffe08958b23fff9211e6ff4d
git push origin v0.4.10

Please sign in to comment.