Performance update for permutations.jl

[depial]

This update contains improvements to the performance of permutations.jl, keeping the underlying algorithm (nextpermutation()) in place (see Issue #204 for various benchmarks, with the benchmarking after the first post directly relevant to this implementation). The main strategy was to move potential overhead away from the performance critical methods of nextpermutation and iterate to the constructors by standardizing input. One step involves separating Permutations from MultiSetPermutations by reverting to a previous version of the iterate method for Permutations (resulting in a large performance boost).

Special attention was paid to reducing the number of allocations made by both permutations() and multiset_permutations(). To this end, one of the larger changes involves now modifying the state in place during iteration (while the data in the structs remains unchanged). I can't currently see this as an issue since the algorithm is serial and can't be parallelized.

In total, these modifications see a cut in allocations to 1/3 and 1/2 their current numbers for permutations() and multiset_permutations(), respectively, and bring their performance into line with Heap's Permutation Algorithm.

Other notes:

• mutlitset_permutations() and permutations() now have the same performance on collections with unique elements (where v.1.1.0 has multiset_permutations() outperforming permutations()).
• The constructors now homogenize the input to the structs, with data and m always being a Vector{T} where T is the element type of the input.
• Attention has been paid to accept any input which is indexable (i.e. no need to be iterable).
• Tested and working with Vectors, Multidimensional Arrays, Sparse Arrays and Offset Arrays (i.e. covering LinearIndex, CartesianIndex and offset indices).
• permutations() is now type safe (in line with multiset_permutations()), always returning a Permutations type.
• multiset_permutations(m, t) is now linear (vs the current $O(n^2)$ version), however, this is likely not terribly important since input size is highly limited.

Heap's Algorithm used in performance comparisons

Note: Below is an implementation of Heap's algorithm which I used to compare performance. It was written to be comparable in structure to nextpermutation() used in this update, but it's not actually correct for permlen < length(data), since it always produces all permutations (i.e. there are potential duplicates).

function heapermutations(data, permlen=length(data), perm=collect(eachindex(data)), datalen=length(data))
    state = ones(Int, datalen)
    output = [data[view(perm, 1:permlen)]]
    i = 1
    while i ≤ datalen
        @inbounds(if state[i] < i
            if isodd(i)
                perm[1], perm[i] = perm[i], perm[1]
            else
                perm[state[i]], perm[i] = perm[i], perm[state[i]]
            end
            push!(output, data[view(perm, 1:permlen)])
            state[i] += 1
            i = 1
        else
            state[i] = 1
            i += 1
        end)
    end
    output
end

[codecov]

Codecov Report

❌ Patch coverage is 98.88889% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 97.19%. Comparing base (b808ce2) to head (153f827).

Files with missing lines	Patch %	Lines
src/permutations.jl	98.88%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #205      +/-   ##
==========================================
+ Coverage   97.17%   97.19%   +0.02%     
==========================================
  Files           8        8              
  Lines         813      857      +44     
==========================================
+ Hits          790      833      +43     
- Misses         23       24       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[depial]

I noticed that derangements() was a filter of multiset_permutations(), so I've coded up a derangement-specific implementation which mirrors that of permutations and multiset permutations. It now has it's own type Derangements and nextderangement() method, which is an iterative version of Rohl's algorithm described here in recursive form.

Some benefits of the new implementation:

• More performant in both time and space.
• Increasingly more performant as multiplicity increases in multisets.
• Includes support for derangements of size t.
• Parallel type implementation to Permutations and MultiSetPermutations

Benchmarking

Unique set elements

Current implementation (v1.1.0) run on collect(derangements(1:10))

BenchmarkTools.Trial: 6 samples with 1 evaluation per sample.
 Range (min … max):  475.889 ms …    1.359 s  ┊ GC (min … max):  0.00% … 58.50%
 Time  (median):        1.009 s               ┊ GC (median):    35.18%
 Time  (mean ± σ):   976.497 ms ± 295.074 ms  ┊ GC (mean ± σ):  41.24% ± 21.38%

  █                            █ █          ██                █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁█▁▁▁▁▁▁▁▁▁▁██▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  476 ms           Histogram: frequency by time          1.36 s <

 Memory estimate: 1.16 GiB, allocs estimate: 19479014.

New Derangements implementation run on collect(derangements(1:10))

BenchmarkTools.Trial: 21 samples with 1 evaluation per sample.
 Range (min … max):  131.592 ms … 551.025 ms  ┊ GC (min … max):  0.00% … 62.50%
 Time  (median):     223.156 ms               ┊ GC (median):    30.68%
 Time  (mean ± σ):   239.835 ms ± 100.070 ms  ┊ GC (mean ± σ):  36.27% ± 21.40%

  █▁█    ▁▁▁ █ █▁▁  ▁▁▁      ▁     ▁   ▁                      ▁  
  ███▁▁▁▁███▁█▁███▁▁███▁▁▁▁▁▁█▁▁▁▁▁█▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  132 ms           Histogram: frequency by time          551 ms <

 Memory estimate: 218.85 MiB, allocs estimate: 2669983.

---

Multiset elements

Current implementation (v1.1.0) run on collect(derangements([1, 1, 2, 4, 5, 5, 6, 7, 7, 9])

BenchmarkTools.Trial: 63 samples with 1 evaluation per sample.
 Range (min … max):  40.422 ms … 154.164 ms  ┊ GC (min … max):  0.00% … 43.45%
 Time  (median):     75.977 ms               ┊ GC (median):    25.91%
 Time  (mean ± σ):   80.700 ms ±  26.078 ms  ┊ GC (mean ± σ):  24.51% ± 18.12%

  █   ▁   ▁▁ ▁  ▄▁     ▁▁█   ▁   ▁     ▁  ▁▁▁            ▁      
  █▁▁▆█▆▁▆██▁█▁▆██▆▆▆▁▆███▆▁▆█▆▁▆█▆▁▆▆▆█▁▆███▁▁▆▁▆▁▆▁▆▆▁▁█▁▁▁▆ ▁
  40.4 ms         Histogram: frequency by time          133 ms <

 Memory estimate: 142.82 MiB, allocs estimate: 2352076.

New Derangements implementation run on collect(derangements([1, 1, 2, 4, 5, 5, 6, 7, 7, 9]))

BenchmarkTools.Trial: 581 samples with 1 evaluation per sample.
 Range (min … max):  5.055 ms … 65.012 ms  ┊ GC (min … max):  0.00% … 81.06%
 Time  (median):     7.474 ms              ┊ GC (median):     0.00%
 Time  (mean ± σ):   8.576 ms ±  5.002 ms  ┊ GC (mean ± σ):  14.50% ± 17.52%

  █▇▃▂▁ ▄▂▂▂▃▅▃▁                                              
  ███████████████▇█▇▇█▆▇▇▇▆▇▇▆▇▅▄▆▁▄▁▅▁▆▁▁▁▁▁▄▄▅▄▆▄▁▁▄▄▁▁▁▄▄ ▇
  5.06 ms      Histogram: log(frequency) by time     27.2 ms <

 Memory estimate: 13.37 MiB, allocs estimate: 168113.

[depial]

Note: I've changed the three argument multiset_permutations(m::Vector, f::Vector{<:Integer}, t::Integer) to be more clearly an outer constructor MultiSetPermutations(m::Vector, f::Vector{<:Integer}, t::Integer) since it appears this method is not actually meant to be exported.

If it is meant to be exported, I believe it would need a docstring which explains how to construct m and f in the way it is done in the two argument method multiset_permutations(a, t::Integer).