- PVSM.RU - https://www.pvsm.ru -
Python — замечательный язык. Несколько языков я и до него пробовал: Pascal в школе; Си, Си с классами, Си++ — в университете. Последние два (три) привили стойкое отвращение к программированию: вместо решения задачи возишься с аллокациями и деструкторами (страшные слова из прошлого), мыслишь в терминах низкоуровневых примитивов. Мое мнение — Си не подходит для решения учебных и научных задач (во всяком случае, в области математики). Уверен, что мне возразят, но я никому не пытаюсь ничего навязать, просто высказываю своё мнение.
Python стал в своё время откровением. Впервые в жизни я писал на несколько уровней абстракции выше, чем это принято в Си. Расстояние между задачей и кодом сократилось как никогда ранее.
Я бы так наверное всю жизнь и писал бы на Python, если бы не пришлось внезапно реализовывать статистические тесты NIST. Казалось бы, задача очень простая: есть массив длины несколько (>= 10) мегабайт, есть набор тестов, которые надо применить к данному массиву.
В python для работы с массивами есть хороший пакет numpy, который по сути является высокоуровневым интерфейсом над быстрыми библиотеками наподобие LAPACK-а. Казалось бы — весь джентельменский набор для научных вычислений имеется (Sympy, Numpy, Scipy, Matplotlib), чего ещё желать?
Каждый, кто имел дело с numpy, знает, что он хорош, но не во всём. Нужно ещё постараться, чтобы операции были векторизованные (как в R), т.е. в некотором смысле единообразны для всего массива. Если вдруг ваша задачка по каким-то причинам не вписывается в эту парадигму, то у вас возникают проблемы.
О каких таких не-векторизуемых задачах я толкую? Да хотя бы взять тот же NIST: вычислить длину линейного регистра сдвига по алгоритму Берлекэмпа-Месси; подсчитать длину самой длинной подпоследовательности из подряд идущих единиц и так далее. Я не исключаю того, что возможно есть какое-то хитроумное решение, которое сведет задачу к векторизованной.
np.count_nonzero(x[:-1] - x[1:])
Это может выглядеть как занимательная разминка для ума, но по сути здесь происходит нечто неестественное, некий трюк, который будет неочевиден читателю через некоторое непродолжительное время. Не говоря уже о том, что это всё равно медленно — аллокацию памяти никто не отменял.
Не-векторизованные операции в Python — это долго. Как с ними бороться, если numpy становится недостаточно? Обычно предлагают несколько решений:
Намаявшись с существующими решениями и не найдя (не сумев запрограммировать) ничего хорошего, решил переписать на чем-то «более быстром». По сути, если вы пишете «молотилку чисел» в 21 веке с нормальной поддержкой массивов, векторизированными операциями «из коробки» и т.д. и т.п., то выбор у вас не особо густой:
rand() # выдать случайное число в диапазоне от 0 до 1
rand(10) # массив длины 10 из случайных чисел от 0 до 1
rand(3,5) # матрица размера 3 х 5 из случайных ....
using Distributions
d = Normal() # Нормальное распределение с параметрами 0, 1
rand(d) # случайное нормально распределенное число
rand(d, 10) # массив длины 10 ... и так далее
То есть, первым аргументом может быть любое (одномерное) распределение из Distributions, но вызов функции останется буквально тем же. И не только распределение (можно передавать собственно сам ГСЧ в виде объекта MersenneTwister и многое другое). Ещё один (на мой взгляд показательный) пример — навигация в DataFrames без loc/iloc.
Пора бы и вспомнить, ради чего всё затевалось.
Примечание: Python я благополучно забыл, поэтому пишите в комментарии свои улучшения, я постараюсь их запустить на своём ноутбуке и посмотреть время исполнения.
Итак, два маленьких примера с микробенчмарками.
def process_fft(x, boundary):
return np.count_nonzero(np.abs(np.fft.fft(2*x-1)) > boundary)
%timeit process_fft(x, 2000)
84.1 ms ± 335 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
function process_fft(x, boundary)
count(t -> t > boundary, abs.(fft(2*x.-1)))
end
@benchmark process_fft(x, 2000)
BenchmarkTools.Trial:
memory estimate: 106.81 MiB
allocs estimate: 61
--------------
minimum time: 80.233 ms (4.75% GC)
median time: 80.746 ms (4.70% GC)
mean time: 85.000 ms (8.67% GC)
maximum time: 205.831 ms (52.27% GC)
--------------
samples: 59
evals/sample: 1
Ничего удивительного мы здесь и не увидим — все всё равно считают не сами, а передают в хорошо оптимизированные фортран-программы.
def longest(x):
maxLen = 0
currLen = 0
# This will count the number of ones in the block
for bit in x:
if bit == 1:
currLen += 1
maxLen = max(maxLen, currLen)
else:
maxLen = max(maxLen, currLen)
currLen = 0
return max(maxLen, currLen)
%timeit longest(x)
599 ms ± 639 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
function longest(x)
maxLen = 0
currLen = 0
# This will count the number of ones in the block
for bit in x
if bit == 1
currLen += 1
maxLen = max(maxLen, currLen)
else
maxLen = max(maxLen, currLen)
currLen = 0
end
end
return max(maxLen, currLen)
end
@benchmark longest(x)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 9.094 ms (0.00% GC)
median time: 9.248 ms (0.00% GC)
mean time: 9.319 ms (0.00% GC)
maximum time: 12.177 ms (0.00% GC)
--------------
samples: 536
evals/sample: 1
Разница очевидна «невооруженным» взглядом. Советы по «допиливанию» и/или векторизации numpy-версии приветствуются. Хочу также обратить внимание, что программы практически идентичны. Для примера я не прописал ни одного типа в Julia (сравни с предыдущим) — всё равно всё работает быстро.
Хочу заметить также, что представленные версии не вошли в итоговую программу, а были далее оптимизированы; здесь же они приведены для примера и без ненужных усложнений (пробрасывание дополнительной памяти в Julia, чтобы делать rfft in-place и т.д.).
Итоговый код для Python и для Julia я не покажу: это тайна (по крайней мере пока что). На момент, когда я бросил допиливать python-версию, она отрабатывала все NIST-тесты на массиве длины 16 мегабайт двоичных знаков за ~ 50 секунд. На Julia на данный момент все тесты прогоняются на том же объеме ~ 20 секунд. Вполне может быть, что это я дурак, и было пространство для оптимизаций и т.д. и т.п. Но это я, каков я есть, со всеми своими достоинствами и недостатками, и на мой взгляд надо считать не сферическую скорость программ в вакууме на языках программирования, а как я лично могу запрограммировать это (без совсем уж грубых ошибок).
Людям здесь [1] стало интересно; я решил написать, когда время появится. В дальнейшем думаю пройтись по Julia более основательно, с примером решения каких-то типовых задачек. Зачем? Больше языков, хороших и разных, и пусть каждый желающий найдет что-то, что подходит лично ему.
Автор: kirtsar
Источник [2]
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/python/297126
Ссылки в тексте:
[1] здесь: https://habr.com/post/423811/
[2] Источник: https://habr.com/post/427879/?utm_campaign=427879
Нажмите здесь для печати.