All combinations of a numpy 2d array filled with 0s and 1s-CodePudding

Given K, I need to have all the possibile combinations of K x 2 numpy matrices so that in each matrix there are all 0s except for two 1s in different rows and columns. Something like this for K = 5:

[[1,0],[0,1],[0,0],[0,0][0,0]]
[[1,0],[0,0],[0,1],[0,0][0,0]]
[[1,0],[0,0],[0,0],[0,1][0,0]]
[[1,0],[0,0],[0,0],[0,0][0,1]]
[[0,0],[1,0],[0,1],[0,0][0,0]]
[[0,0],[1,0],[0,0],[0,1][0,0]]
... and so on

So the resulting array should be a K x 2 x (K*(K-1)/2). I want to avoid loops since it's not an efficient way when K is big enough (in my specific case K = 300)

CodePudding user response：

I can't think of an elegant solution but here's a not-so-elegant pure numpy one:

import numpy as np

def combination_matrices(K):
    # get combination indices
    i, j = np.indices((K, K))
    comb_indices = np.transpose((i < j).nonzero())  # (num_combs, 2) array where ones are
    num_combs = comb_indices.shape[0]  # K*(K-1)/2

    # create a matrix of the desired shape, first axis enumerates combinations
    matrices = np.zeros((num_combs, K, 2), dtype=int)
    # broadcasting assignment of ones
    comb_range, col_index = np.ogrid[:num_combs, :2]
    matrices[comb_range, comb_indices, col_index] = 1
    return matrices

This first uses the indices of a (K, K)-shaped array to find the index pairs for every combination (these are indices that encode the upper triangle of the array, excluding the diagonal). Then we use a bit tricky broadcasting assignment (heavy , python's itertools doesn't currently support this. So simplest solution is to use the multiset tools of the sympy library.

The following code took about ~2.5 minutes to run on my machine, so is fairly fast for a single thread. You're looking at 179700 unique permutations for K=300.

(I took inspiration from https://stackoverflow.com/a/40289807/10739860)

from collections import Counter
from math import factorial, prod

import numpy as np
from sympy.utilities.iterables import multiset_permutations
from tqdm import tqdm


def No_multiset_permutations(multiset: list) -> int:
    """Calculates the No. possible permutations given a multiset.
    See: https://en.wikipedia.org/wiki/Permutation#Permutations_of_multisets

    :param multiset: List representing a multiset.
    """
    value_counts = Counter(multiset).values()
    denominator = prod([factorial(val) for val in value_counts])
    return int(factorial(len(multiset)) / denominator)


def multiset_Kx2_permutations(K: int) -> np.ndarray:
    """This will generate all possible unique Kx2 permutations of an array
    withsize K where two values are 1 and the rest are 0.

    :param K: The size of the array.
    """
    # Construct number multiset, e.g. K=5 gives [1, 1, 0, 0, 0, 0, 0, 0, 0, 0]
    numbers = [1, 1]   [0] * (K - 1) * 2

    # Use sympy's multiset_permutations to get a multiset permutation generator
    generator = multiset_permutations(numbers)

    # Calculate the No. possible permutations
    number_of_perms = No_multiset_permutations(numbers)

    # Get all permutations, bonus progress bar is included :)
    unique_perms = [next(generator) for _ in tqdm(range(number_of_perms))]

    # Reshape each permutation to Kx2
    unique_perms = np.array(unique_perms, dtype=np.int8)
    return unique_perms.reshape(-1, K, 2)


if __name__ == "__main__":
    solution = multiset_Kx2_permutations(300)