Threading over permutations

General Usage
1

I am trying to parallelize a loop over the Combinatorics permutations(1:n) without using collect() since the permutations I am working with are too big and they cause memory issues. I think there is a problem with just using Threads.@threads for perm in permutations(1:n) since I can’t get the index of a permutation without using collect. Is there any way to use threads but over the permutations(1:n) without using collect? For example, my function using collect does the following:

function find_lamda_coef(r, indep_M, groundset_size, lambda, t)
    """
    Computes c_λ where [M] = Σ_λ^* [Ω_λ^*] using the characters t

    Parameters:
    - r: rank of matroid
    - idep_M: independent subsets of matroid
    -groundset_size: size of groundset of the matroid
    -lambda: partition
    -t: list of integers 

    Returns:
    - c_λ
    """
    n = groundset_size - 1 
    perms = collect(permutations(1:n+1))
    summands_list = zeros(length(perms))

    Threads.@threads for i in 1:length(perms)
        perm = perms[i]
        B_perm = find_lex_first_basis(r, indep_M, perm)
        x_sym = [t[i] for i in B_perm]
        y_sym = [-t[i] for i in 1:length(t)]
        s_lam = s_lambda_x_y(lambda, x_sym, y_sym)
        denom = one(parent(t[1]))
            
        for j in 1:n
            denom *= (t[perm[j]] - t[perm[j+1]])
        end
        summands_list[i] = s_lam/denom
    end
return sum(summands_list)
end
2

This doesn’t exactly use Combinatorics.permutations, but if you generate the permutations on 1:n via Heap’s algorithm, then you can easily divide them into n parts according to where n goes. Each of these n parts can be assigned to a separate thread.

EDIT: Well, the problem is that it’s not clear which permutation you have to start with if you want n to be send to some other element k. The example on the Wikipedia page suggests that the first permutation in Heap’s algorithm that sends n to k maps the other elements in increasing or decreasing order, depending on the parity of n-k.

3

Hi, and welcome to the Julia community!

The issue is that permutations returns a generator, i.e. a way to (sequentially) produce permutations without actually computing them yet. Using collect does perform this computation, for all entries at once.

The easiest fix would be to use perm = nthperm(1:n+1, i) inside of the @threads for loop.

I expect this to work mostly fine in the case of permutations (though I’m not necessarily knowledgeable on how to efficiently list them). But in certain cases it might be more efficient to iterate, rather than get each element individually. For example, obtaining the first 10 Fibonacci numbers jointly is more efficient than finding them individually (i.e. without caching/memoisation, which is what the iteration state does). When direct indexing is not possible, you would need to manually perform the iteration in the main thread, creating new tasks when previous ones have finished.

Manual task creation example

In the code below we manually iterate over the permutations. When we have resources available, we create a Task to handle the current permutation. To know when a task is done (without constantly looping over them), we use a Channel which will contain the indices of Tasks which have finished.

using Combinatorics

function handle_permutation(perm, channel, task_idx)
    sleep_time = rand() * 5  # Sleep for a random time up to 5 s, to simulate heterogenous computation times
    sleep(sleep_time)
    println("$task_idx: $perm ($sleep_time)") 
    put!(channel, task_idx)  # Signal our task is done
end

function threaded_handle_permutations(perms, nb_tasks)
    channel = Channel{Int}(nb_tasks)
    tasks = Array{Task}(undef, nb_tasks)
    it = iterate(perms)  # Either a permutation and state, or nothing when we are finished
    for task_idx = 1:nb_tasks  # Setup, where we don't have to wait until Tasks are done
        isnothing(it) && break
        perm, state = it
        let perm = perm
            tasks[task_idx] = Task(() -> handle_permutation(perm, channel, task_idx))
        end
        schedule(tasks[task_idx])  # (or use @spawn)
        it = iterate(perms, state)
    end

    while !isnothing(it)
        perm, state = it
        task_idx = take!(channel)  # tasks[task_idx] is ready
        let perm = perm
            tasks[task_idx] = Task(() -> handle_permutation(perm, channel, task_idx))
        end
        schedule(tasks[task_idx])
        it = iterate(perms, state)
    end
    wait.(tasks)  # Wait until all Tasks have finished
    close(channel)
end

threaded_handle_permutations(permutations(1:4), Threads.nthreads())

with example output

4: [1, 3, 4, 2] (1.2117935089662901)
2: [1, 2, 4, 3] (1.5704595207841299)
5: [1, 4, 2, 3] (2.7562927397700028)
7: [2, 1, 3, 4] (3.427912051023216)
1: [1, 2, 3, 4] (3.9286201986106724)
3: [1, 3, 2, 4] (4.150059037617181)
8: [2, 1, 4, 3] (4.41249434461974)
1: [3, 1, 2, 4] (0.8946810612178713)
6: [1, 4, 3, 2] (4.943705002447017)
8: [3, 2, 1, 4] (0.5524897265787826)
4: [2, 3, 1, 4] (4.140999678913283)
2: [2, 3, 4, 1] (4.170628020371655)
1: [3, 2, 4, 1] (0.9551984725843093)
7: [2, 4, 3, 1] (2.4064765159373414)
6: [3, 4, 1, 2] (1.0672378759488854)
5: [2, 4, 1, 3] (3.5553860548631944)
5: [4, 3, 2, 1] (0.5065772298310989)
3: [3, 1, 4, 2] (3.9094692535870177)
6: [4, 3, 1, 2] (2.789216187966278)
2: [4, 1, 3, 2] (3.504362858170555)
8: [3, 4, 2, 1] (4.593180548236061)
7: [4, 2, 3, 1] (3.935449181248672)
4: [4, 1, 2, 3] (4.61345171679344)
1: [4, 2, 1, 3] (4.4173984501241454)