Generazione di una matrice numpy con tutte le combinazioni di numeri che sommano a meno di un dato numero

Question

Esistono diversi eleganti esempi di utilizzo di numpy in Python per generare matrici di tutte le combinazioni. Ad esempio la risposta qui: Using numpy to build an array of all combinations of two arrays.

Supponiamo ora che ci sia un vincolo in più, vale a dire che la somma di tutti i numeri non può sommare più di una determinata costante K. Utilizzando un generatore e itertools.product, per un esempio con K=3 dove vogliamo le combinazioni di tre variabili con campi di 0-1,0-3, 0-2 e lo possiamo fare una segue:

from itertools import product 
K = 3 
maxRange = np.array([1,3,2]) 
states = np.array([i for i in product(*(range(i+1) for i in maxRange)) if sum(i)<=K]) 

che restituisce

array([[0, 0, 0], 
     [0, 0, 1], 
     [0, 0, 2], 
     [0, 1, 0], 
     [0, 1, 1], 
     [0, 1, 2], 
     [0, 2, 0], 
     [0, 2, 1], 
     [0, 3, 0], 
     [1, 0, 0], 
     [1, 0, 1], 
     [1, 0, 2], 
     [1, 1, 0], 
     [1, 1, 1], 
     [1, 2, 0]]) 

In linea di principio, l'approccio da https://stackoverflow.com/a/25655090/1479342 può essere utilizzato per generare tutte le possibili combinazioni senza il vincolo e quindi selezionando il sottoinsieme di combinazioni che riassumono a meno di K. Tuttavia, questo approccio genera molte più combinazioni del necessario, specialmente se K è relativamente piccolo rispetto a sum(maxRange).

Ci deve essere un modo per farlo più velocemente e con un minor utilizzo di memoria. Come può essere realizzato utilizzando un approccio vettoriale (ad esempio utilizzando np.indices)?

Answer 1

cura

1. Per completezza, sto aggiungendo qui il codice del PO:

   def partition0(max_range, S): 
       K = len(max_range) 
       return np.array([i for i in itertools.product(*(range(i+1) for i in max_range)) if sum(i)<=S]) 
   

2. Il primo approccio è pura np.indices. È veloce per piccoli input ma consuma molta memoria (OP ha già sottolineato che non è quello che intendeva).

   def partition1(max_range, S): 
       max_range = np.asarray(max_range, dtype = int) 
       a = np.indices(max_range + 1) 
       b = a.sum(axis = 0) <= S 
       return (a[:,b].T) 
   

3. approccio ricorrente sembra essere molto meglio di quelli di cui sopra:

   def partition2(max_range, max_sum): 
       max_range = np.asarray(max_range, dtype = int).ravel()   
       if(max_range.size == 1): 
        return np.arange(min(max_range[0],max_sum) + 1, dtype = int).reshape(-1,1) 
       P = partition2(max_range[1:], max_sum) 
       # S[i] is the largest summand we can place in front of P[i]    
       S = np.minimum(max_sum - P.sum(axis = 1), max_range[0]) 
       offset, sz = 0, S.size 
       out = np.empty(shape = (sz + S.sum(), P.shape[1]+1), dtype = int) 
       out[:sz,0] = 0 
       out[:sz,1:] = P 
       for i in range(1, max_range[0]+1): 
        ind, = np.nonzero(S) 
        offset, sz = offset + sz, ind.size 
        out[offset:offset+sz, 0] = i 
        out[offset:offset+sz, 1:] = P[ind] 
        S[ind] -= 1 
       return out 
   

4. Dopo un breve pensiero, sono stato in grado di prendere un po 'più lontano. Se conosciamo in anticipo il numero di possibili partizioni, possiamo allocare memoria sufficiente in una volta. (È in qualche modo simile a cartesian in an already linked thread.)

   Innanzitutto, abbiamo bisogno di una funzione che contenga le partizioni.

   def number_of_partitions(max_range, max_sum): 
       ''' 
       Returns an array arr of the same shape as max_range, where 
       arr[j] = number of admissible partitions for 
         j summands bounded by max_range[j:] and with sum <= max_sum 
       ''' 
       M = max_sum + 1 
       N = len(max_range) 
       arr = np.zeros(shape=(M,N), dtype = int)  
       arr[:,-1] = np.where(np.arange(M) <= min(max_range[-1], max_sum), 1, 0) 
       for i in range(N-2,-1,-1): 
        for j in range(max_range[i]+1): 
         arr[j:,i] += arr[:M-j,i+1] 
       return arr.sum(axis = 0) 
   

   La funzione principale:

   def partition3(max_range, max_sum, out = None, n_part = None): 
       if out is None: 
        max_range = np.asarray(max_range, dtype = int).ravel() 
        n_part = number_of_partitions(max_range, max_sum) 
        out = np.zeros(shape = (n_part[0], max_range.size), dtype = int) 
   
       if(max_range.size == 1): 
        out[:] = np.arange(min(max_range[0],max_sum) + 1, dtype = int).reshape(-1,1) 
        return out 
   
       P = partition3(max_range[1:], max_sum, out=out[:n_part[1],1:], n_part = n_part[1:])   
       # P is now a useful reference 
   
       S = np.minimum(max_sum - P.sum(axis = 1), max_range[0]) 
       offset, sz = 0, S.size 
       out[:sz,0] = 0 
       for i in range(1, max_range[0]+1): 
        ind, = np.nonzero(S) 
        offset, sz = offset + sz, ind.size 
        out[offset:offset+sz, 0] = i 
        out[offset:offset+sz, 1:] = P[ind] 
        S[ind] -= 1 
       return out 
   

5. Alcuni test:

   max_range = [3, 4, 6, 3, 4, 6, 3, 4, 6] 
   for f in [partition0, partition1, partition2, partition3]: 
       print(f.__name__ + ':') 
       for max_sum in [5, 15, 25]: 
        print('Sum %2d: ' % max_sum, end = '') 
        %timeit f(max_range, max_sum) 
       print() 
   
   partition0: 
   Sum 5: 1 loops, best of 3: 859 ms per loop 
   Sum 15: 1 loops, best of 3: 1.39 s per loop 
   Sum 25: 1 loops, best of 3: 3.18 s per loop 
   
   partition1: 
   Sum 5: 10 loops, best of 3: 176 ms per loop 
   Sum 15: 1 loops, best of 3: 224 ms per loop 
   Sum 25: 1 loops, best of 3: 403 ms per loop 
   
   partition2: 
   Sum 5: 1000 loops, best of 3: 809 µs per loop 
   Sum 15: 10 loops, best of 3: 62.5 ms per loop 
   Sum 25: 1 loops, best of 3: 262 ms per loop 
   
   partition3: 
   Sum 5: 1000 loops, best of 3: 853 µs per loop 
   Sum 15: 10 loops, best of 3: 59.1 ms per loop 
   Sum 25: 1 loops, best of 3: 249 ms per loop 
   

   E qualcosa di più grande:

   %timeit partition0([3,6] * 5, 20) 
   1 loops, best of 3: 11.9 s per loop 
   
   %timeit partition1([3,6] * 5, 20) 
   The slowest run took 12.68 times longer than the fastest. This could mean that an intermediate result is being cached 
   1 loops, best of 3: 2.33 s per loop 
   # MemoryError in another test 
   
   %timeit partition2([3,6] * 5, 20) 
   1 loops, best of 3: 877 ms per loop 
   
   %timeit partition3([3,6] * 5, 20) 
   1 loops, best of 3: 739 ms per loop 
   

Answer 2

Non so che cos'è un approccio numpy, ma ecco una soluzione abbastanza pulita. Sia A un array di numeri interi e sia k un numero che viene fornito come input.

Iniziare con un array vuoto B; mantenere la somma dell'array B in una variabile s (inizialmente impostata su zero). Applicare la seguente procedura:

• se la somma s dell'array B è inferiore k, quindi (i) aggiungere alla raccolta, (ii) e per ogni elemento dall'array originale A , aggiungere tale elemento a B e aggiornare s, (iii) cancellarlo da A e (iv) applicare la procedura in modo ricorsivo; (iv) quando la chiamata ritorna, aggiungere l'elemento a A e aggiornare s; altro non fare nulla.

Questo approccio bottom-up prugne rami non validi nella fase iniziale e visite solo i sottoinsiemi necessari (vale a dire quasi solo i sottoinsiemi che somma a meno di k).