Подтвердить что ты не робот

Numpy, проблема с длинными массивами

У меня есть два массива (a и b) с n целыми элементами в диапазоне (0, N).

typo: массивы с 2 ^ n целыми числами, где наибольшее целое число принимает значение N = 3 ^ n

Я хочу рассчитать сумму каждой комбинации элементов в и b (sum_ij_ = a_i_ + b_j_ для всех i, j). Затем возьмем модуль N (sum_ij_ = sum_ij_% N) и, наконец, вычислим частоту различных сумм.

Чтобы сделать это быстро с numpy, без каких-либо циклов, я попытался использовать функцию meshgrid и bincount.

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

Теперь проблема в том, что мои входные массивы длинны. И meshgrid дает мне MemoryError, когда я использую входы с 2 ^ 13 элементами. Я хотел бы рассчитать это для массивов с элементами 2 ^ 15-2 ^ 20.

то есть n в диапазоне от 15 до 20

Есть ли какие-нибудь умные трюки, чтобы сделать это с помощью numpy?

Любая помощь будет высоко оценена.

- джон

4b9b3361

Ответ 1

попробуйте сделать это. ваш meshgrid является матрицей NxN, блокируйте до 10x10 N/10xN/10 и просто вычисляйте 100 бит, добавьте их в конец. это использует только 1% столько памяти, сколько делает все это.

Ответ 2

Изменить в ответ на комментарий jonalm:

jonalm: N ~ 3 ^ n не n ~ 3 ^ N. N - максимальный элемент в a, а n - число элементов в a.

n равно ~ 2 ^ 20. Если N равно ~ 3 ^ n, то N равно ~ 3 ^ (2 ^ 20) > 10 ^ (500207). Ученые оценивают (http://www.stormloader.com/ajy/reallife.html), что во Вселенной насчитывается всего около 10 ^ 87 частиц. Таким образом, нет (наивного) способа, которым компьютер может обрабатывать int размером 10 ^ (500207).

jonalm: Я, тем не менее, немного интересуюсь функцией pv(), которую вы определяете. (Я не удается запустить его, поскольку text.find() не определен (угадать его в другом модуль)). Как работает эта функция и каково ее преимущество?

pv - небольшая вспомогательная функция, которую я написал для отладки значения переменных. Он работает как print(), за исключением случаев, когда вы произносите pv (x), он печатает как литеральное имя переменной (или строку выражения), так и двоеточие, а затем значение переменной.

Если вы положили

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

в script вы должны получить

x: 1

Скромное преимущество использования pv over print заключается в том, что он экономит ваше печатание. Вместо того, чтобы написать

print('x: %s'%x)

вы можете просто пощекотать

pv(x)

Если отслеживать несколько переменных, полезно отметить эти переменные. Я просто устал писать все это.

Функция pv работает, используя модуль трассировки, чтобы заглянуть в строку кода используется для вызова самой функции pv. (См. http://docs.python.org/library/traceback.html#module-traceback). Эта строка кода хранится как строка в переменном тексте. text.find() - вызов обычного метода string find(). Например, если

text='pv(x)'

затем

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

Я предполагаю, что n ~ 3 ^ N и n ~ 2 ** 20

Идея заключается в работе модуля N. Это сокращает размер массивов. Вторая идея (важная, когда n огромна) заключается в использовании numpy ndarrays типа 'object', потому что если вы используете целочисленный dtype, вы рискуете переполнять размер максимально допустимого целого числа.

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

Вы можете изменить n на 2 ** 20, но ниже я покажу, что происходит с небольшим n поэтому выход легче читать.

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 содержит число 0s, 1s, 2s, 3s в wb содержит число 0s, 1s, 2s, 3s в 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')

Подумайте о 0 как токен или чип. Аналогично для 1,2,3.

Подумайте о wa = [24 28 28 20], так как есть сумка с 24 0-чипами, 28 1-стружкой, 28 2-стружками, 20 3-струбцинами.

У тебя есть сумка и сумка. Когда вы рисуете чип из каждого мешка, вы "добавляете" их вместе и формируете новый чип. Вы "mod" ответ (по модулю N).

Представьте, что вы взяли 1-чип из wb-мешка и добавили его с каждым чипом в сумке.

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)

Поскольку в сумке wb есть 34 1-чипа, когда вы добавляете их против всех фишек в сумке wa = [24 28 28 20], вы получаете

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

Это всего лишь частичный счет из-за 34 1-фишек. Вы также должны обращаться с другим типы чипов в wb-bag, но это показывает вам метод, используемый ниже:

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]

Это конечный результат: 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]

Ответ 3

Проверьте свою математику, что много места вы просите:

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

Если каждый из ваших элементов был всего одним байтом, это уже один терабайт памяти.

Добавьте цикл или два. Эта проблема не подходит для максимизации вашей памяти и минимизации ваших вычислений.