Numpy, problema con arrays largos

Tengo dos matrices (ayb) con n elementos enteros en el rango (0, N).

typo: arrays con 2 ^ n enteros donde el entero más grande toma el valor N = 3 ^ n

Quiero calcular la sum de cada combinación de elementos en a y b (sum_ij_ = a_i_ + b_j_ para todos i, j ). Luego tome el módulo N (sum_ij_ = sum_ij_% N), y finalmente calcule la frecuencia de las diferentes sums.

Para hacer esto rápido con numpy, sin ningún tipo de bucles, traté de usar el meshgrid y la función bincount.

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

Ahora, el problema es que mis matrices de entrada son largas. Y meshgrid me da MemoryError cuando uso entradas con 2 ^ 13 elementos. Me gustaría calcular esto para matrices con 2 ^ 15-2 ^ 20 elementos.

    eso es n en el rango de 15 a 20

    ¿Hay algún truco inteligente para hacer esto con numpy?

    Cualquier ayuda será muy apreciada.

    – jon

    Intenta triturarlo. su meshgrid es una matriz NxN, bloque que hasta 10×10 N / 10xN / 10 y solo calcula 100 contenedores, súmelos al final. Esto solo usa ~ 1% de la memoria que todo el proceso.

    Revisa tus matemáticas, eso es mucho espacio que estás pidiendo:

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

    Si cada uno de sus elementos era solo un byte, eso ya es un terabyte de memoria.

    Agrega un bucle o dos. Este problema no es adecuado para maximizar su memoria y minimizar su cálculo.

    Editar en respuesta al comentario de jonalm:

    jonalm: N ~ 3 ^ n no n ~ 3 ^ N. N es el elemento máximo en a y n es el número de elementos en a.

    n es ~ 2 ^ 20. Si N es ~ 3 ^ n, entonces N es ~ 3 ^ (2 ^ 20)> 10 ^ (500207). Los científicos estiman ( http://www.stormloader.com/ajy/reallife.html ) que solo hay alrededor de 10 ^ 87 partículas en el universo. Por lo tanto, no hay una manera (ingenua) de que una computadora pueda manejar un int de tamaño 10 ^ (500207).

    jonalm: Sin embargo, soy un poco curioso acerca de la función pv () que define. (No logro ejecutarlo ya que text.find () no está definido (supongo que está en otro módulo)). ¿Cómo funciona esta función y cuál es su ventaja?

    pv es una pequeña función auxiliar que escribí para depurar el valor de las variables. Funciona como print (), excepto cuando dices pv (x) imprime el nombre de la variable literal (o cadena de expresión), dos puntos y luego el valor de la variable.

    Si pones

     #!/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) 

    en un guión debes obtener

     x: 1 

    La modesta ventaja de utilizar pv over print es que le ahorra escribir. En lugar de tener que escribir

     print('x: %s'%x) 

    usted puede simplemente abofetear

     pv(x) 

    Cuando hay varias variables para rastrear, es útil etiquetar las variables. Simplemente me cansé de escribirlo todo.

    La función pv funciona utilizando el módulo de rastreo para ver la línea de código utilizada para llamar a la función pv. (Consulte http://docs.python.org/library/traceback.html#module-traceback ) Esa línea de código se almacena como una cadena en el texto variable. text.find () es una llamada al método de cadena habitual find (). Por ejemplo, si

     text='pv(x)' 

    entonces

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

    Estoy asumiendo que n ~ 3 ^ N, y n ~ 2 ** 20

    La idea es trabajar en el módulo N. Esto reduce el tamaño de los arreglos. La segunda idea (importante cuando n es enorme) es usar una cantidad de ndarrays de tipo “objeto” porque si usa un tipo de entero, corre el riesgo de desbordar el tamaño del máximo entero permitido.

     #!/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)) 

    Puede cambiar n para que sea 2 ** 20, pero a continuación muestro lo que sucede con small n para que la salida sea más fácil de leer.

     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 el número de 0s, 1s, 2s, 3s en wb tiene el número de 0s, 1s, 2s, 3s en 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') 

    Piense en un 0 como un token o chip. Del mismo modo para 1,2,3.

    Piense que wa = [24 28 28 20] significa que hay una bolsa con 24 fichas, 28 fichas 1, 28 fichas 2, 20 fichas.

    Tienes un wa-bag y un wb-bag. Cuando robas un chip de cada bolsa, las “sums” y creas un nuevo chip. Usted “mod” la respuesta (módulo N).

    Imagine tomar un chip de la bolsa de wb y agregarlo con cada chip en la bolsa de wa.

     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) 

    Ya que hay 34 fichas 1 en la bolsa wb, cuando las agregas contra todas las fichas en la bolsa wa = [24 28 28 20], obtienes

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

    Este es solo el conteo parcial debido a los 34 1-chips. También tienes que manejar los otros tipos de chips en el wb-bag, pero esto te muestra el método utilizado a continuación:

     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] 

    Este es el resultado final: 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]