4371106380

Build Type

github

Committed by GitHub

Commit Message

Merge cdfafba80 into 881920951

Pull Request Pull Request #196: Apple silicon workaround

Run Details

15 of 15 new or added lines in 4 files covered. (100.0%)

4459 of 4841 relevant lines covered (92.11%)

9410575.87 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

87.16

/src/BitStringAddresses/sortedparticlelist.jl

"""
    select_int_type(M)

Select unsigned integer type that can hold values up to `M`.
"""
function select_int_type(M)
    if M ≤ 0
        throw(ArgumentError("`M` must be positive!"))
    elseif M ≤ typemax(UInt8)
        return UInt8
    elseif M ≤ typemax(UInt16)
        return UInt16
    elseif M ≤ typemax(UInt32)
        return UInt32
    else
        return UInt64
    end
end

"""
    SortedParticleList{N,M,T<:Unsigned}

Type for storing sparse fock states. Stores the mode number of each particle as an entry
with only its mode stored. The entries are always kept sorted.

Iterating over `SortedParticleList`s yields occupied modes as a tuple of occupation number,
mode number, and position in list.

# Constructors

* `SortedParticleList{N,M,T}(::SVector{N,T})`: unsafe constructor. Does not sort input.

* `SortedParticleList(arr::AbstractVector)`: convert ONR to `SortedParticleList`

"""
struct SortedParticleList{N,M,T<:Unsigned}
    storage::SVector{N,T}
end
function Base.show(io::IO, ss::SortedParticleList{N,M,T}) where {N,M,T}
    print(io, "SortedParticleList{$N,$M,$T}(", Int.(ss.storage), ")")
end

function SortedParticleList{N,M}() where {N,M}
    T = select_int_type(M)
    return SortedParticleList{N,M,T}(ones(SVector{N,T}))
end
function from_onr(::Type{S}, onr) where {N,M,T,S<:SortedParticleList{N,M,T}}
    spl = zeros(MVector{N,T})
    curr = 1
    for (n, v) in enumerate(onr)
        for _ in 1:v
            spl[curr] = n
            curr += 1
        end
    end
    return SortedParticleList{N,M,T}(SVector(spl))
end
function from_onr(::Type{S}, onr::SparseVector) where {N,M,T,S<:SortedParticleList{N,M,T}}
    spl = zeros(MVector{N,T})
    curr = 1
    for (n, v) in zip(rowvals(onr), nonzeros(onr))
        for _ in 1:v
            spl[curr] = n
            curr += 1
        end
    end
    return SortedParticleList{N,M,T}(SVector(spl))
end

function Base.isless(ss1::SortedParticleList, ss2::SortedParticleList)
    return isless(ss1.storage, ss2.storage)
end

num_particles(::Type{<:SortedParticleList{N}}) where {N} = N
num_modes(::Type{<:SortedParticleList{<:Any,M}}) where {M} = M

###
### General functions
###
Base.eltype(::SortedParticleList) = Tuple{Int,Int,Int}

function Base.length(ss::SortedParticleList{<:Any,<:Any,T}) where {T}
    curr = zero(T)
    res = 0
    for i in ss.storage
        res += (i != curr)
        curr = i
    end
    return res
end

function Base.iterate(ss::SortedParticleList{N}, i=1) where {N}
    @inbounds if i > N
        return nothing
    else
        occnum = 1
        mode = ss.storage[i]
        offset = i
        i += 1
        while i ≤ N && ss.storage[i] == mode
            occnum += 1
            i += 1
        end
        return (occnum, Int(mode), i - occnum - 1), i
    end
end

function Base.reverse(ss::SortedParticleList{N,M,T}) where {N,M,T}
    new_storage = map(reverse(ss.storage)) do i
        T(M) - i + one(T)
    end
    return SortedParticleList{N,M,T}(new_storage)
end

# Somehow this is faster than the default method.
Base.hash(ss::SortedParticleList, u::UInt) = hash(ss.storage, u)

# In this case, getting the ONR is the same for bosons and fermions, and assumes the address
# is not malformed.
function onr(ss::SortedParticleList{<:Any,M}) where {M}
    mvec = zeros(MVector{M,Int})
    @inbounds for (occnum, mode, _) in ss
        mvec[mode] = occnum
    end
    return SVector(mvec)
end

# Same as above.
function find_mode(ss::SortedParticleList, n)
    offset = 0
    for (occnum, mode, _) in ss
        if mode == n
            return (occnum, mode, offset)
        elseif mode > n
            return (0, n, offset)
        end
        offset += occnum
    end
    return (0, n, offset)
end

"""
    move_particles(ss::SortedParticleList, dsts, srcs)

Move several particles at once. Moves `srcs[i]` to `dsts[i]`.  `dsts` and `srcs` should be
tuples of [`BoseFSIndex`](@ref) or [`FermiFSIndex`](@ref). The legality of the moves is not
checked - the result of an illegal move is undefined!
"""
function move_particles(ss::SortedParticleList{N,M,T}, dsts, srcs) where {N,M,T}
    new_storage = ss.storage
    for (dst, src) in zip(dsts, srcs)
        src_pos = 1
        @inbounds for i in 1:N
            src_pos = max(src_pos, i * (new_storage[i] == src.mode % T))
        end
        @boundscheck 0 < dst.mode ≤ M || throw(BoundsError(ss, dst.mode))
        new_storage = setindex(new_storage, dst.mode % T, src_pos)
    end
    return SortedParticleList{N,M,T}(sort(new_storage))
end

###
### Bose interface
###
function from_bose_onr(::Type{S}, onr) where{S<:SortedParticleList}
    from_onr(S, onr)
end
to_bose_onr(ss::SortedParticleList, _) = onr(ss)

bose_num_occupied_modes(ss::SortedParticleList) = length(ss)

bose_find_mode(ss::SortedParticleList, n) = find_mode(ss, n)

@inline function bose_excitation(
    ss::SortedParticleList{N,M,T}, creations, destructions
) where {N,M,T}
    creations_rev = reverse(creations)
    value = bose_excitation_value(creations_rev, destructions)
    if iszero(value)
        return ss, 0.0
    else
        return move_particles(ss, creations_rev, destructions), √value
    end
end

Base.length(bom::BoseOccupiedModes{<:Any,<:Any,<:SortedParticleList}) = length(bom.storage)
function Base.iterate(bom::BoseOccupiedModes{<:Any,<:Any,<:SortedParticleList}, i=1)
    it = iterate(bom.storage, i)
    if isnothing(it)
        return nothing
    else
        res, i = it
        return BoseFSIndex(res...), i
    end
end

###
### Fermi interface
###
function from_fermi_onr(::Type{S}, onr) where {S<:SortedParticleList}
    from_onr(S, onr)
end

# Fix offsets and occupation numbers after creation/destruction operator is applied.
# The idea behind these is to allow computing the value from the indices alone.
function _fix_pos_create(c, index)
    index = @set index.offset += (c.mode < index.mode)
    index = @set index.occnum += (c.mode == index.mode)
    return index
end
_fix_pos_create(c) = Base.Fix1(_fix_pos_create, c)
function _fix_pos_destroy(d, index)
    index = @set index.offset -= (d.mode < index.mode)
    index = @set index.occnum -= (d.mode == index.mode)
    return index
end
_fix_pos_destroy(d) = Base.Fix1(_fix_pos_destroy, d)

"""
    fermi_excitation_value_spl(
        creations::NTuple{_,FermiFSIndex}, destructions::NTuple{_,::FermiFSIndex}
    ) -> {-1,0,1}

Compute the value of an excitation from indices. Starts by applying all destruction
operators, and then applying all creation operators. The operators must be given in reverse
order. Will return 0 if move is illegal.

Note that this function only works on indices obtained from a [`SortedParticleList`](@ref).
"""
@inline fermi_excitation_value_spl(::Tuple{}, ::Tuple{}) = 1.0
@inline function fermi_excitation_value_spl((c, cs...), ::Tuple{})
    cs = map(_fix_pos_create(c), cs)
    return fermi_excitation_value_spl(cs, ()) * ifelse(isodd(c.offset), -1, 1) * (c.occnum == 0)
end
@inline function fermi_excitation_value_spl(cs, (d, ds...))
    cs = map(_fix_pos_destroy(d), cs)
    ds = map(_fix_pos_destroy(d), ds)
    return fermi_excitation_value_spl(cs, ds) * ifelse(isodd(d.offset), -1, 1) * (d.occnum == 1)
end

fermi_find_mode(ss::SortedParticleList, n) = FermiFSIndex(find_mode(ss, n))
function fermi_find_mode(ss::SortedParticleList, ns::Tuple)
    # It's OK to do that instead of the fancy method used with bosons, because the
    # assumption is that `N` is small.
    return map(n -> FermiFSIndex(find_mode(ss, n)), ns)
end

@inline function fermi_excitation(
    ss::SortedParticleList{N,M,T}, creations::NTuple{K}, destructions::NTuple{K}
) where {N,M,T,K}
    creations_rev = reverse(creations)
    destructions_rev = reverse(destructions)
    value = fermi_excitation_value_spl(creations_rev, destructions_rev)
    if iszero(value)
        return ss, 0.0
    else
        return move_particles(ss, creations_rev, destructions), float(value)
    end
end

Base.length(fom::FermiOccupiedModes{N,<:SortedParticleList}) where {N} = length(fom.storage)
function Base.iterate(fom::FermiOccupiedModes{<:Any,<:SortedParticleList}, i=1)
    itr = iterate(fom.storage, i)
    if isnothing(itr)
        return nothing
    else
        res, i = itr
        return FermiFSIndex(res...), i
    end
end

1	"""
2	select_int_type(M)
3
4	Select unsigned integer type that can hold values up to `M`.
5	"""
6	function select_int_type(M)	324✔
7	if M ≤ 0	324✔
8	throw(ArgumentError("`M` must be positive!"))	×
9	elseif M ≤ typemax(UInt8)	324✔
10	return UInt8	324✔
11	elseif M ≤ typemax(UInt16)	×
12	return UInt16	×
13	elseif M ≤ typemax(UInt32)	×
14	return UInt32	×
15	else
16	return UInt64	×
17	end
18	end
19
20	"""
21	SortedParticleList{N,M,T<:Unsigned}
22
23	Type for storing sparse fock states. Stores the mode number of each particle as an entry
24	with only its mode stored. The entries are always kept sorted.
25
26	Iterating over `SortedParticleList`s yields occupied modes as a tuple of occupation number,
27	mode number, and position in list.
28
29	# Constructors
30
31	* `SortedParticleList{N,M,T}(::SVector{N,T})`: unsafe constructor. Does not sort input.
32
33	* `SortedParticleList(arr::AbstractVector)`: convert ONR to `SortedParticleList`
34
35	"""
36	struct SortedParticleList{N,M,T<:Unsigned}
37	storage::SVector{N,T}	13,997✔
38	end
39	function Base.show(io::IO, ss::SortedParticleList{N,M,T}) where {N,M,T}	×
40	print(io, "SortedParticleList{$N,$M,$T}(", Int.(ss.storage), ")")	×
41	end
42
43	function SortedParticleList{N,M}() where {N,M}	×
44	T = select_int_type(M)	×
45	return SortedParticleList{N,M,T}(ones(SVector{N,T}))	×
46	end
47	function from_onr(::Type{S}, onr) where {N,M,T,S<:SortedParticleList{N,M,T}}	6,875✔
48	spl = zeros(MVector{N,T})	6,875✔
49	curr = 1	6,875✔
50	for (n, v) in enumerate(onr)	6,877✔
51	for _ in 1:v	1,487,559✔
52	spl[curr] = n	116,890✔
53	curr += 1	116,890✔
54	end	120,830✔
55	end	2,742,343✔
56	return SortedParticleList{N,M,T}(SVector(spl))	6,875✔
57	end
58	function from_onr(::Type{S}, onr::SparseVector) where {N,M,T,S<:SortedParticleList{N,M,T}}	61✔
59	spl = zeros(MVector{N,T})	61✔
60	curr = 1	61✔
61	for (n, v) in zip(rowvals(onr), nonzeros(onr))	122✔
62	for _ in 1:v	1,570✔
63	spl[curr] = n	804✔
64	curr += 1	804✔
65	end	823✔
66	end	1,509✔
67	return SortedParticleList{N,M,T}(SVector(spl))	61✔
68	end
69
70	function Base.isless(ss1::SortedParticleList, ss2::SortedParticleList)	×
71	return isless(ss1.storage, ss2.storage)	×
72	end
73
74	num_particles(::Type{<:SortedParticleList{N}}) where {N} = N	1✔
75	num_modes(::Type{<:SortedParticleList{<:Any,M}}) where {M} = M	1✔
76
77	###
78	### General functions
79	###
80	Base.eltype(::SortedParticleList) = Tuple{Int,Int,Int}	×
81
82	function Base.length(ss::SortedParticleList{<:Any,<:Any,T}) where {T}	131✔
83	curr = zero(T)	131✔
84	res = 0	131✔
85	for i in ss.storage	131✔
86	res += (i != curr)	701✔
87	curr = i	701✔
88	end	1,271✔
89	return res	131✔
90	end
91
92	function Base.iterate(ss::SortedParticleList{N}, i=1) where {N}	12,691,123✔
93	@inbounds if i > N	12,691,123✔
94	return nothing	314,985✔
95	else
96	occnum = 1	11,761,283✔
97	mode = ss.storage[i]	11,761,283✔
98	offset = i	11,761,283✔
99	i += 1	11,761,283✔
100	while i ≤ N && ss.storage[i] == mode	12,081,089✔
101	occnum += 1	319,806✔
102	i += 1	319,806✔
103	end	319,806✔
104	return (occnum, Int(mode), i - occnum - 1), i	11,761,283✔
105	end
106	end
107
108	function Base.reverse(ss::SortedParticleList{N,M,T}) where {N,M,T}	20✔
109	new_storage = map(reverse(ss.storage)) do i	20✔
110	T(M) - i + one(T)	300✔
111	end
112	return SortedParticleList{N,M,T}(new_storage)	20✔
113	end
114
115	# Somehow this is faster than the default method.
116	Base.hash(ss::SortedParticleList, u::UInt) = hash(ss.storage, u)	354✔
117
118	# In this case, getting the ONR is the same for bosons and fermions, and assumes the address
119	# is not malformed.
120	function onr(ss::SortedParticleList{<:Any,M}) where {M}	123,335✔
121	mvec = zeros(MVector{M,Int})	123,335✔
122	@inbounds for (occnum, mode, _) in ss	123,335✔
123	mvec[mode] = occnum	2,368,029✔
124	end	4,612,723✔
125	return SVector(mvec)	123,335✔
126	end
127
128	# Same as above.
129	function find_mode(ss::SortedParticleList, n)	491,520✔
130	offset = 0	491,520✔
131	for (occnum, mode, _) in ss	491,520✔
132	if mode == n	2,968,766✔
133	return (occnum, mode, offset)	23,443✔
134	elseif mode > n	2,945,323✔
135	return (0, n, offset)	424,214✔
136	end
137	offset += occnum	2,521,109✔
138	end	4,998,355✔
139	return (0, n, offset)	43,863✔
140	end
141
142	"""
143	move_particles(ss::SortedParticleList, dsts, srcs)
144
145	Move several particles at once. Moves `srcs[i]` to `dsts[i]`. `dsts` and `srcs` should be
146	tuples of [`BoseFSIndex`](@ref) or [`FermiFSIndex`](@ref). The legality of the moves is not
147	checked - the result of an illegal move is undefined!
148	"""
149	function move_particles(ss::SortedParticleList{N,M,T}, dsts, srcs) where {N,M,T}	7,041✔
150	new_storage = ss.storage	7,041✔
151	for (dst, src) in zip(dsts, srcs)	7,041✔
152	src_pos = 1	7,561✔
153	@inbounds for i in 1:N	7,561✔
154	src_pos = max(src_pos, i * (new_storage[i] == src.mode % T))	126,200✔
155	end	244,839✔
156	@boundscheck 0 < dst.mode ≤ M \|\| throw(BoundsError(ss, dst.mode))	7,561✔
157	new_storage = setindex(new_storage, dst.mode % T, src_pos)	7,561✔
158	end	8,081✔
159	return SortedParticleList{N,M,T}(sort(new_storage))	7,041✔
160	end
161
162	###
163	### Bose interface
164	###
165	function from_bose_onr(::Type{S}, onr) where{S<:SortedParticleList}	4,839✔
166	from_onr(S, onr)	4,839✔
167	end
168	to_bose_onr(ss::SortedParticleList, _) = onr(ss)	123,335✔
169
170	bose_num_occupied_modes(ss::SortedParticleList) = length(ss)	701✔
171
172	bose_find_mode(ss::SortedParticleList, n) = find_mode(ss, n)	×
173
174	@inline function bose_excitation(	123,464✔
175	ss::SortedParticleList{N,M,T}, creations, destructions
176	) where {N,M,T}
177	creations_rev = reverse(creations)	123,464✔
178	value = bose_excitation_value(creations_rev, destructions)	123,464✔
179	if iszero(value)	123,464✔
180	return ss, 0.0	118,457✔
181	else
182	return move_particles(ss, creations_rev, destructions), √value	5,007✔
183	end
184	end
185
186	Base.length(bom::BoseOccupiedModes{<:Any,<:Any,<:SortedParticleList}) = length(bom.storage)	×
187	function Base.iterate(bom::BoseOccupiedModes{<:Any,<:Any,<:SortedParticleList}, i=1)	5,712,598✔
188	it = iterate(bom.storage, i)	10,414,379✔
189	if isnothing(it)	9,921,378✔
190	return nothing	24,815✔
191	else
192	res, i = it	5,194,782✔
193	return BoseFSIndex(res...), i	5,194,782✔
194	end
195	end
196
197	###
198	### Fermi interface
199	###
200	function from_fermi_onr(::Type{S}, onr) where {S<:SortedParticleList}	2,097✔
201	from_onr(S, onr)	2,097✔
202	end
203
204	# Fix offsets and occupation numbers after creation/destruction operator is applied.
205	# The idea behind these is to allow computing the value from the indices alone.
206	function _fix_pos_create(c, index)	163,840✔
207	index = @set index.offset += (c.mode < index.mode)	163,840✔
208	index = @set index.occnum += (c.mode == index.mode)	163,840✔
209	return index	163,840✔
210	end
211	_fix_pos_create(c) = Base.Fix1(_fix_pos_create, c)	245,760✔
212	function _fix_pos_destroy(d, index)	737,280✔
213	index = @set index.offset -= (d.mode < index.mode)	737,280✔
214	index = @set index.occnum -= (d.mode == index.mode)	737,280✔
215	return index	737,280✔
216	end
217	_fix_pos_destroy(d) = Base.Fix1(_fix_pos_destroy, d)	491,520✔
218
219	"""
220	fermi_excitation_value_spl(
221	creations::NTuple{_,FermiFSIndex}, destructions::NTuple{_,::FermiFSIndex}
222	) -> {-1,0,1}
223
224	Compute the value of an excitation from indices. Starts by applying all destruction
225	operators, and then applying all creation operators. The operators must be given in reverse
226	order. Will return 0 if move is illegal.
227
228	Note that this function only works on indices obtained from a [`SortedParticleList`](@ref).
229	"""
230	@inline fermi_excitation_value_spl(::Tuple{}, ::Tuple{}) = 1.0	1✔
231	@inline function fermi_excitation_value_spl((c, cs...), ::Tuple{})	245,760✔
232	cs = map(_fix_pos_create(c), cs)	245,760✔
233	return fermi_excitation_value_spl(cs, ()) * ifelse(isodd(c.offset), -1, 1) * (c.occnum == 0)	245,760✔
234	end
235	@inline function fermi_excitation_value_spl(cs, (d, ds...))	245,760✔
236	cs = map(_fix_pos_destroy(d), cs)	245,760✔
237	ds = map(_fix_pos_destroy(d), ds)	245,760✔
238	return fermi_excitation_value_spl(cs, ds) * ifelse(isodd(d.offset), -1, 1) * (d.occnum == 1)	245,760✔
239	end
240
241	fermi_find_mode(ss::SortedParticleList, n) = FermiFSIndex(find_mode(ss, n))	491,520✔
242	function fermi_find_mode(ss::SortedParticleList, ns::Tuple)	×
243	# It's OK to do that instead of the fancy method used with bosons, because the
244	# assumption is that `N` is small.
245	return map(n -> FermiFSIndex(find_mode(ss, n)), ns)	×
246	end
247
248	@inline function fermi_excitation(	122,880✔
249	ss::SortedParticleList{N,M,T}, creations::NTuple{K}, destructions::NTuple{K}
250	) where {N,M,T,K}
251	creations_rev = reverse(creations)	122,880✔
252	destructions_rev = reverse(destructions)	122,880✔
253	value = fermi_excitation_value_spl(creations_rev, destructions_rev)	122,880✔
254	if iszero(value)	122,880✔
255	return ss, 0.0	120,846✔
256	else
257	return move_particles(ss, creations_rev, destructions), float(value)	2,034✔
258	end
259	end
260
261	Base.length(fom::FermiOccupiedModes{N,<:SortedParticleList}) where {N} = length(fom.storage)	×
262	function Base.iterate(fom::FermiOccupiedModes{<:Any,<:SortedParticleList}, i=1)	1,475,650✔
263	itr = iterate(fom.storage, i)	2,582,384✔
264	if isnothing(itr)	2,459,412✔
265	return nothing	122,972✔
266	else
267	res, i = itr	1,229,706✔
268	return FermiFSIndex(res...), i	1,229,706✔
269	end
270	end

joachimbrand / Rimu.jl / 4371106380

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous