Improve `manybodyoperator` performance once more

aryavorskiy commented 1 year ago

state_transition! instead of coefficient allows to loop over occupations only once. The state index can be found by searching in the occupations list - but it is still way faster.
The occupation vectors are usually sorted lexicographically (and we can sort them if not). We can make use of it to make searching in the occupations list more performant. These two tricks yield yet another ~300x performance boost :)

The OccupationsIterator interface can make it even faster, because for fermion states one can avoid storing all occupation vectors by calculating the state index using a formula - but I'd better leave it for another pull request.

Before (36 modes, 3 fermions):

BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max):  3.925 s …   3.952 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.939 s              ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.939 s ± 19.174 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                                                       █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  3.93 s         Histogram: frequency by time        3.95 s <

 Memory estimate: 2.63 MiB, allocs estimate: 69.

After:

BenchmarkTools.Trial: 341 samples with 1 evaluation.
 Range (min … max):  12.618 ms … 48.976 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     14.262 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   14.672 ms ±  2.320 ms  ┊ GC (mean ± σ):  1.43% ± 4.52%

       ▁    █ ▄▁▃ ▂▁
  ▅▂▄▅▇█▇█▇█████████▇▆▅▄▄▆▃▃▂▃▃▂▄▃▁▃▃▄▂▃▃▃▃▂▂▂▁▁▂▂▂▂▁▁▁▁▁▂▁▃▂ ▃
  12.6 ms         Histogram: frequency by time        19.7 ms <

 Memory estimate: 2.98 MiB, allocs estimate: 14377.

UPD: Important detail is that the default state ordering in manybody bases is now slightly changed if Nparticles is a vector. Hope that does not break anything.

codecov[bot] commented 1 year ago

Codecov Report

Merging #128 (ac60054) into master (06809f7) will decrease coverage by 0.18%. The diff coverage is 97.51%.

@@            Coverage Diff             @@
##           master     #128      +/-   ##
==========================================
- Coverage   93.55%   93.37%   -0.18%     
==========================================
  Files          24       24              
  Lines        3055     3018      -37     
==========================================
- Hits         2858     2818      -40     
- Misses        197      200       +3

Files	Coverage Δ
src/manybody.jl	`98.40% <97.51%> (-1.25%)`	:arrow_down:

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

Krastanov commented 1 year ago

This looks great, thank you as always!

I am not too worried about the change of ordering in the returned vector. However, I think this has at least one change that probably should be considered breaking. Previously the return type for many of these functions was Vector, but now it is Occupations, and Occupations does not permit operations that users might already be relying on.

master:

julia> fermionstates(2,2)
1-element Vector{Vector{Int64}}:
 [1, 1]

julia> fermionstates(2,2)[1]
2-element Vector{Int64}:
 1
 1

this pr:

julia> fermionstates(2,2)
QuantumOpticsBase.Occupations{Vector{Int64}}([[1, 1]])

julia> fermionstates(2,2)[1]
ERROR: MethodError: no method matching getindex(::QuantumOpticsBase.Occupations{Vector{Int64}}, ::Int64)
Stacktrace:
 [1] top-level scope
   @ REPL[2]:1

I think we need to support a bit more of the Vector API before being comfortable merging this.

I also pushed an empty commit to this branch to rerun the CI. Feel free to overwrite it.

aryavorskiy commented 1 year ago

Done. The best way to support the vector API is implementing the AbstractVector interface :) I removed the OccupationsIterator interface completely - one can extend this functionality by subtyping AbstractVector directly. The ManyBodyBasis also accepts any type of containers now.

Krastanov commented 1 year ago

I do not think that Occupations{T} <: AbstractVector{Vector{T}} actually helps here. Subtyping in julia does not guarantee that any particular interface is being implemented or followed. I think this is an object-oriented idiom that does not really apply to Julia. Unless you are using these for something, may we just do struct Occupations{T} without any supertypes? Actually, if ManyBodyBasis can take any container, what is the reason to create Occupations?

Naturally, I might be wrong about these suggestions of mine, thus please run me through your reasoning to make sure I am not missing something.

Implementing formal interfaces in Julia is itself a very interesting topic with many threads on Discourse and experiments in Github repos, but no proper solution yet.

aryavorskiy commented 1 year ago

I'm sorry I disappeared for a while :) Implementing the AbstractVector interface may be useful for non-trivial use cases like applying map or find to an Occupations object. As for the Occupations structure itself, its main and only feature is that it wraps a sorted array which makes searching much faster. Is's probably better to call it in some other way (maybe SortedContainer?), but I think the structure itself should be kept as it is now. I can attach some benchmarks to prove there is a significant difference in performance

Krastanov commented 1 year ago

No worries!

I agree that implementing the AbstractVector interface would be useful, but subclassing AbstractVector does not contribute to that goal, unless there are some specific methods that are already defined for AbstractVector and that rely on methods defined for Occupations.

Are you saying that your definitions of size and getindex are enough for "deriving" map, find, length, collect, etc, thanks to subtyping AbstractVector? If that is the case, I would agree with the subtyping (but I would avoid calling it an interface as I understand "interface" to have a different more formal meaning).

If these indeed work now, could you add a few "interface" unit tests? Checking that getting an index, getting a length, collecting, etc all work? collect is probably the main one as it would be the way to fix broken user code (if it turns out that we actually broke something)

aryavorskiy commented 1 year ago

Sure, I'll rename Occupations to SortedVector and add some unit tests. The old name feels kinda bad because there is no functionality designed specially for occupation vectors.

Speaking of interfaces, I used this manual page which says that redefining size and getindex along with setting IndexStyle to Linear() implements the default read-only vector semantics.

Krastanov commented 1 year ago

Speaking of interfaces, I used [https://docs.julialang.org/en/v1/manual/interfaces/](this manual page) which says that redefining size and getindex along with setting IndexStyle to Linear() implements the default read-only vector semantics.

Oh, I understand now! I am all on board with this. I do not envision other points to discuss and would be eager to merge this after the addition of a couple of interface tests. Thank you so much for this contribution!!!

There are a couple of "correctness of interface implementation testing" packages in the ecosystem. I will add an issue to track the possible use of such a package in the future (not specifically for this contribution, rather in general). It can be a good google summer of code onboarding tasks too.

aryavorskiy commented 1 year ago

Done. By the way, this PR also seems to cover the issue #147 - I added a test to ensure.

qojulia / QuantumOpticsBase.jl

Improve `manybodyoperator` performance once more #128

Codecov Report