Generare tutti i possibili risultati di k sfere in n bin (somma di risultati multinomiali/categoriali)

Question

Supponiamo di disporre di contenitori n in cui stiamo lanciando palle . Che cosa è un veloce (vale a dire utilizzando numpy/scipy anziché codice python) per generare tutti i possibili risultati come una matrice?

Per esempio, se n = 4 e k = 3, vorremmo il seguente numpy.array:

3 0 0 0 
2 1 0 0 
2 0 1 0 
2 0 0 1 
1 2 0 0 
1 1 1 0 
1 1 0 1 
1 0 2 0 
1 0 1 1 
1 0 0 2 
0 3 0 0 
0 2 1 0 
0 2 0 1 
0 1 2 0 
0 1 1 1 
0 1 0 2 
0 0 3 0 
0 0 2 1 
0 0 1 2 
0 0 0 3 

Mi scuso se ogni permutazione è stato mancato, ma questo è l'idea generale. Le permutazioni generate non devono essere in alcun ordine particolare, ma l'elenco precedente era conveniente per iterare categoricamente attraverso di esse mentalmente.

Ancora meglio, c'è un modo per mappare ogni numero intero da 1 a multiset number (la cardinalità di questo elenco) direttamente a una determinata permutazione?

Questa domanda è legata a quelle che seguono, che sono attuate in R con molto diversi servizi:

• Generating all permutations of N balls in M bins
• Generate a matrix of all possible outcomes for throwing n dice (ignoring order)

riferimenti anche correlati:

• https://en.wikipedia.org/wiki/Stars_and_bars_(combinatorics)
• https://en.wikipedia.org/wiki/Multiset#Counting_multisets
• https://en.wikipedia.org/wiki/Combinatorial_number_system

Answer 1

Ecco una soluzione generatore utilizzando itertools.combinations_with_replacement, non so se sarà adatto alle vostre esigenze.

def partitions(n, b): 
    masks = numpy.identity(b, dtype=int) 
    for c in itertools.combinations_with_replacement(masks, n): 
     yield sum(c) 

output = numpy.array(list(partitions(3, 4))) 
# [[3 0 0 0] 
# [2 1 0 0] 
# ... 
# [0 0 1 2] 
# [0 0 0 3]] 

La complessità di questa funzione cresce esponenzialmente, per cui v'è un confine discreto tra ciò che è fattibile e cosa no.

Si noti che mentre gli array di numpy devono conoscere le loro dimensioni durante la costruzione, ciò è facilmente possibile poiché il numero del multiset è facilmente reperibile. Al di sotto di potrebbe essere essere un metodo migliore, non ho avuto tempi.

from math import factorial as fact 
from itertools import combinations_with_replacement as cwr 

nCr = lambda n, r: fact(n)/fact(n-r)/fact(r) 

def partitions(n, b): 
    partition_array = numpy.empty((nCr(n+b-1, b-1), b), dtype=int) 
    masks = numpy.identity(b, dtype=int) 
    for i, c in enumerate(cwr(masks, n)): 
     partition_array[i,:] = sum(c) 
    return partition_array 

Answer 2

ecco un'implementazione naif con comprehensions elenco, non è sicuro su prestazioni rispetto a NumPy

def gen(n,k): 
    if(k==1): 
     return [[n]] 
    if(n==0): 
     return [[0]*k] 
    return [ g2 for x in range(n+1) for g2 in [ u+[n-x] for u in gen(x,k-1) ] ] 

> gen(3,4) 
[[0, 0, 0, 3], 
[0, 0, 1, 2], 
[0, 1, 0, 2], 
[1, 0, 0, 2], 
[0, 0, 2, 1], 
[0, 1, 1, 1], 
[1, 0, 1, 1], 
[0, 2, 0, 1], 
[1, 1, 0, 1], 
[2, 0, 0, 1], 
[0, 0, 3, 0], 
[0, 1, 2, 0], 
[1, 0, 2, 0], 
[0, 2, 1, 0], 
[1, 1, 1, 0], 
[2, 0, 1, 0], 
[0, 3, 0, 0], 
[1, 2, 0, 0], 
[2, 1, 0, 0], 
[3, 0, 0, 0]] 

Answer 3

A scopo di riferimento, il codice seguente utilizza Ehrlich's algorithm per scorrere tutte le possibili combinazioni di un multinsieme in C++ , Javascript, e Python:

https://github.com/ekg/multichoose

Questo può essere convertito nel formato sopra usando this method.In particolare,

for s in multichoose(k, set): 
    row = np.bincount(s, minlength=len(set) + 1) 

Questo non è ancora pura NumPy, ma può essere utilizzato per riempire un preassegnati numpy.array abbastanza rapidamente.

Answer 4

Ecco la soluzione che ho trovato per questo.

import numpy, itertools 
def multinomial_combinations(n, k, max_power=None): 
    """returns a list (2d numpy array) of all length k sequences of 
    non-negative integers n1, ..., nk such that n1 + ... + nk = n.""" 
    bar_placements = itertools.combinations(range(1, n+k), k-1) 
    tmp = [(0,) + x + (n+k,) for x in bar_placements] 
    sequences = numpy.diff(tmp) - 1 
    if max_power: 
     return sequences[numpy.where((sequences<=max_power).all(axis=1))][::-1] 
    else: 
     return sequences[::-1] 

Nota 1: Il [:: - 1] alla fine inverte semplicemente l'ordine in modo che corrisponda all'output di esempio.

Nota 2: Trovare queste sequenze equivale a trovare tutti i modi per disporre n stelle e barre k-1 in (per riempire n + punti k-1) (vedere stars and bars thm 2).

Nota 3: L'argomento max_power consente di restituire solo sequenze in cui tutti gli interi sono al di sotto di un numero massimo.