2009-11-08 20 views
11

Ho due array (a e b) con n elementi interi nell'intervallo (0, N).Numpy, problema con i long array

battitura: array con 2^n interi dove il massimo intero assume il valore N = 3^n

voglio calcolare la somma di ogni combinazione di elementi in ae b (sum_ij_ = a_i_ + b_j_ per tutti i, j). Quindi prendi il modulo N (sum_ij_ = sum_ij_% N) e infine calcola la frequenza delle diverse somme.

Per fare ciò velocemente con NumPy, senza alcun loop, ho provato a utilizzare la meshgrid e la funzione bincount.

A,B = numpy.meshgrid(a,b) 
A = A + B 
A = A % N 
A = numpy.reshape(A,A.size) 
result = numpy.bincount(A) 

Ora, il problema è che i miei array di input sono lunghi. E meshgrid mi dà MemoryError quando uso input con elementi 2^13. Vorrei calcolare questo per gli array con 2^15-2^20 elementi.

che è n nel range da 15 a 20

C'è qualche trucchi intelligenti per fare questo con numpy?

Qualsiasi aiuto sarà molto apprezzato.

- jon

+0

E quanto è grande N? – unutbu

+0

Numpy sarà davvero così efficiente? Immagino che staresti meglio in C++, scrivendo le tue funzioni e ottimizzando il più possibile. Da quello che suona come numpy non può gestire array così grandi. Anche se devo dire se hai due array con 2^15 a 2^20 elementi, allora se guardi tutte le loro diverse somme, ti ritroverai con una serie di elementi da 2^30 a 2^40. Il che è molto .. – JSchlather

+0

@unutbu: N ~ 3^n @liberalkid: Penso che tu abbia ragione. Le mie abilità di C++ non sono poi così buone. – jonalm

risposta

7

provare a sistemarlo.il tuo meshgrid è una matrice NxN, blocca fino a 10x10 N/10xN/10 e calcola solo 100 bin, aggiungili alla fine. questo usa solo il ~ 1% di memoria quanto tutto il resto.

+0

Immagino che questa sia la strada da percorrere, ma c'è un modo intelligente per farlo con gli array numpy. Riduzione al minimo dell'uso di cicli for. – jonalm

+0

Ehi, c'è una dimensione ottimale per un blocco? – jonalm

+0

probabilmente il più grande è possibile fare un blocco e mantenerlo comunque al sicuro in ram. – Autoplectic

1

Modifica in risposta al commento di jonalm:

jonalm: N ~ 3^n non n ~ 3^N. N è l'elemento massimo in a e n è il numero di elementi in a.

n è ~ 2^20. Se N è ~ 3^n allora N è ~ 3^(2^20)> 10^(500207). Gli scienziati stimano (http://www.stormloader.com/ajy/reallife.html) che nell'universo ci sono solo circa 10^87 particelle. Quindi non esiste (in modo ingenuo) un computer in grado di gestire un int di dimensione 10^(500207).

jonalm: Sono comunque un po 'curioso sulla funzione pv() che definite. (I non riesco a eseguirlo in quanto text.find() non è definito (indovinalo in un altro modulo )). Come funziona questa funzione e qual è il suo vantaggio?

pv è una piccola funzione di supporto che ho scritto per eseguire il debug del valore delle variabili. Funziona come print() tranne quando si dice pv (x) che stampa sia il nome della variabile letterale (o stringa di espressione), i due punti, e quindi il valore della variabile.

Se mettete

#!/usr/bin/env python 
import traceback 
def pv(var): 
    (filename,line_number,function_name,text)=traceback.extract_stack()[-2] 
    print('%s: %s'%(text[text.find('(')+1:-1],var)) 
x=1 
pv(x) 

in uno script si dovrebbe ottenere

x: 1 

Il modesto vantaggio di utilizzare pv nel corso di stampa è che consente di risparmiare la digitazione.Invece di dover scrivere

print('x: %s'%x) 

si può solo schiaffo giù

pv(x) 

Quando ci sono più variabili per tenere traccia, è utile per etichettare le variabili. Mi sono solo stancato di scrivere tutto.

La funzione pv funziona utilizzando il modulo traceback per visualizzare la riga del codice utilizzata per chiamare la funzione pv stessa. (Vedere http://docs.python.org/library/traceback.html#module-traceback) Tale riga di codice viene memorizzata come stringa nel testo variabile. text.find() è una chiamata al solito metodo di stringa find(). Ad esempio, se

text='pv(x)' 

poi

text.find('(') == 2    # The index of the '(' in string text 
text[text.find('(')+1:-1] == 'x' # Everything in between the parentheses 

sto supponendo n ~ 3^N, e n ~ 2 ** 20

L'idea è di lavorare modulo N. This tagli in basso sulla dimensione degli array. La seconda idea (importante quando n è enorme) consiste nell'utilizzare numpy ndarrays di tipo "oggetto" perché se si utilizza un numero intero di dtype si corre il rischio di traboccare la dimensione del numero intero massimo consentito.

#!/usr/bin/env python 
import traceback 
import numpy as np 

def pv(var): 
    (filename,line_number,function_name,text)=traceback.extract_stack()[-2] 
    print('%s: %s'%(text[text.find('(')+1:-1],var)) 

È possibile modificare n per essere 2 ** 20, ma al di sotto vi mostro quello che succede con piccola n quindi l'uscita è più facile da leggere.

n=100 
N=int(np.exp(1./3*np.log(n))) 
pv(N) 
# N: 4 

a=np.random.randint(N,size=n) 
b=np.random.randint(N,size=n) 
pv(a) 
pv(b) 
# a: [1 0 3 0 1 0 1 2 0 2 1 3 1 0 1 2 2 0 2 3 3 3 1 0 1 1 2 0 1 2 3 1 2 1 0 0 3 
# 1 3 2 3 2 1 1 2 2 0 3 0 2 0 0 2 2 1 3 0 2 1 0 2 3 1 0 1 1 0 1 3 0 2 2 0 2 
# 0 2 3 0 2 0 1 1 3 2 2 3 2 0 3 1 1 1 1 2 3 3 2 2 3 1] 
# b: [1 3 2 1 1 2 1 1 1 3 0 3 0 2 2 3 2 0 1 3 1 0 0 3 3 2 1 1 2 0 1 2 0 3 3 1 0 
# 3 3 3 1 1 3 3 3 1 1 0 2 1 0 0 3 0 2 1 0 2 2 0 0 0 1 1 3 1 1 1 2 1 1 3 2 3 
# 3 1 2 1 0 0 2 3 1 0 2 1 1 1 1 3 3 0 2 2 3 2 0 1 3 1] 

wa contiene il numero di 0, 1s, 2s, 3s in un wb contiene il numero di 0, 1s, 2s, 3s in b

wa=np.bincount(a) 
wb=np.bincount(b) 
pv(wa) 
pv(wb) 
# wa: [24 28 28 20] 
# wb: [21 34 20 25] 
result=np.zeros(N,dtype='object') 

pensare a un 0 come segno o chip. Allo stesso modo per 1,2,3.

Pensate a wa = [24 28 28 20] nel senso che esiste una borsa con 24 chip 0, 28 1-chip, 28 2-chip, 20 3-chip.

Hai un wa-bag e un wb-bag. Quando peschi un gettone da ciascun sacchetto, li "aggiungi" e forma un nuovo chip. "Mod" la risposta (modulo N).

Immagina di prendere un 1-chip dal wb-bag e aggiungerlo con ogni chip nel wa-bag.

1-chip + 0-chip = 1-chip 
1-chip + 1-chip = 2-chip 
1-chip + 2-chip = 3-chip 
1-chip + 3-chip = 4-chip = 0-chip (we are mod'ing by N=4) 

Poiché ci sono 34 1-chips nel sacchetto wb, quando si aggiunge contro tutte le chips nel wa = [24 28 28 20] Borsa, si ottiene

34*24 1-chips 
34*28 2-chips 
34*28 3-chips 
34*20 0-chips 

Questo è solo il conteggio parziale dovuto ai 34 1-chips. Hai anche a gestire le altre tipi di chip nel wb-bag, ma questo vi mostra il metodo utilizzato di seguito:

for i,count in enumerate(wb): 
    partial_count=count*wa 
    pv(partial_count) 
    shifted_partial_count=np.roll(partial_count,i) 
    pv(shifted_partial_count) 
    result+=shifted_partial_count 
# partial_count: [504 588 588 420] 
# shifted_partial_count: [504 588 588 420] 
# partial_count: [816 952 952 680] 
# shifted_partial_count: [680 816 952 952] 
# partial_count: [480 560 560 400] 
# shifted_partial_count: [560 400 480 560] 
# partial_count: [600 700 700 500] 
# shifted_partial_count: [700 700 500 600] 

pv(result)  
# result: [2444 2504 2520 2532] 

Questo è il risultato finale: 2444 0s, 2504 1s, 2520 2s, 2532 3s .

# This is a test to make sure the result is correct. 
# This uses a very memory intensive method. 
# c is too huge when n is large. 
if n>1000: 
    print('n is too large to run the check') 
else: 
    c=(a[:]+b[:,np.newaxis]) 
    c=c.ravel() 
    c=c%N 
    result2=np.bincount(c) 
    pv(result2) 
    assert(all(r1==r2 for r1,r2 in zip(result,result2))) 
# result2: [2444 2504 2520 2532] 
+0

Si noti che 'c% = N' funziona (e può utilizzare il doppio della memoria). – EOL

+0

@EOL, sì, c% = N è migliore. Tuttavia, la definizione di 'c = (a [:] + b [:, np.newaxis])' significa che hai già perso la battaglia, poiché si tratta di un enorme array di forme 2-d (n, n) mentre il precedente la soluzione non utilizza nient'altro che una coppia di matrici 1-d di forma (N). – unutbu

+0

Grazie mille per la risposta, mi piace questo metodo. Ma non penso che questo mi sarà d'aiuto dato che tutti i numeri nella matrice a (e in b) sono diversi (non ho detto questo, il mio male). bincount (a) sarà composto solo da 1 e 0. N ~ 3^n not n ~ 3^N. N è l'elemento massimo in aen è il numero di elementi in a. Sono comunque un po 'curioso riguardo alla funzione pv() che definisci. (Non riesco a eseguirlo in quanto text.find() non è definito (indovinate il suo in un altro modulo)). Come funziona questa funzione e qual è il suo vantaggio? – jonalm

1

Controllare la matematica, che è un sacco di spazio che stai chiedendo:

2^20 * 2^20 = 2^40 = 1 099 511 627 776

Se ognuno di vostra gli elementi erano solo un byte, è già un terabyte di memoria.

Aggiungi un ciclo o due. Questo problema non è adatto per massimizzare la memoria e ridurre al minimo il tuo calcolo.