Хабрахабр

[Из песочницы] Как я с Python на Julia переходил (и зачем)

Немного предыстории о Python

Python — замечательный язык. Несколько языков я и до него пробовал: Pascal в школе; Си, Си с классами, Си++ — в университете. Последние два (три) привили стойкое отвращение к программированию: вместо решения задачи возишься с аллокациями и деструкторами (страшные слова из прошлого), мыслишь в терминах низкоуровневых примитивов. Мое мнение — Си не подходит для решения учебных и научных задач (во всяком случае, в области математики). Уверен, что мне возразят, но я никому не пытаюсь ничего навязать, просто высказываю своё мнение.

Впервые в жизни я писал на несколько уровней абстракции выше, чем это принято в Си. Python стал в своё время откровением. Расстояние между задачей и кодом сократилось как никогда ранее.

Казалось бы, задача очень простая: есть массив длины несколько (>= 10) мегабайт, есть набор тестов, которые надо применить к данному массиву.
Я бы так наверное всю жизнь и писал бы на Python, если бы не пришлось внезапно реализовывать статистические тесты NIST.

Чем хорош numpy?

В python для работы с массивами есть хороший пакет numpy, который по сути является высокоуровневым интерфейсом над быстрыми библиотеками наподобие LAPACK-а. Казалось бы — весь джентельменский набор для научных вычислений имеется (Sympy, Numpy, Scipy, Matplotlib), чего ещё желать?

Нужно ещё постараться, чтобы операции были векторизованные (как в R), т.е. Каждый, кто имел дело с numpy, знает, что он хорош, но не во всём. Если вдруг ваша задачка по каким-то причинам не вписывается в эту парадигму, то у вас возникают проблемы. в некотором смысле единообразны для всего массива.

Да хотя бы взять тот же NIST: вычислить длину линейного регистра сдвига по алгоритму Берлекэмпа-Месси; подсчитать длину самой длинной подпоследовательности из подряд идущих единиц и так далее. О каких таких не-векторизуемых задачах я толкую? Я не исключаю того, что возможно есть какое-то хитроумное решение, которое сведет задачу к векторизованной.

Хитроумно?

Как пример из того же NIST: необходимо было подсчитать число «переключений» последовательности, где под «переключением» я имею в виду смену подряд идущих единиц (...1111...) на подряд идущие нули (...000...), и наоборот. Для этого можно взять исходную последовательность без последнего элемента (x[: -1]) и вычесть из неё последовательность, сдвинутую на 1 (x[2: ]), а затем рассчитать число ненулевых элементов в полученной новой последовательности. Всё вместе будет выглядеть как:

np.count_nonzero(x[:-1] - x[1:])

Это может выглядеть как занимательная разминка для ума, но по сути здесь происходит нечто неестественное, некий трюк, который будет неочевиден читателю через некоторое непродолжительное время. Не говоря уже о том, что это всё равно медленно — аллокацию памяти никто не отменял.

Не-векторизованные операции в Python — это долго. Как с ними бороться, если numpy становится недостаточно? Обычно предлагают несколько решений:

  1. Numba JIT. Если бы она работала так, как это описано на заглавной страничке Numba, то я думаю, что здесь стоило бы закончить рассказ. За давностью я уже немного подзабыл, что с ней у меня пошло не так; суть в том, что ускорение было не такое впечатляющее, на какое я рассчитывал, увы.
  2. Cython. ОК, поднимите руки те, кто считает, что cython — это действительно красивое, элегантное решение, которое не разрушает сам смысл и дух Python? Я так не считаю; если вы дошли до Cython, то можете уже перестать себя обманывать и перейти на нечто менее изощрённое вроде C++ и Си.
  3. Перепиши узкие места на Си и дергай их из твоего любимого Питона. ОК, а что, если у меня вся программа — это сплошные вычисления и узкие места? Си не предлагать! Я уже не говорю о том, что в этом решении нужно хорошо знать уже не один, но два языка — и Python, и Си.

Here comes the JULIA!

Намаявшись с существующими решениями и не найдя (не сумев запрограммировать) ничего хорошего, решил переписать на чем-то «более быстром». По сути, если вы пишете «молотилку чисел» в 21 веке с нормальной поддержкой массивов, векторизированными операциями «из коробки» и т.д. и т.п., то выбор у вас не особо густой:

  1. Fortran. И не надо смеяться, кто из нас не дёргал BLAS/LAPACK из своих любимых языков? Fortran — это действительно хороший (лучше СИ!) язык для НАУЧНЫХ вычислений. Не взял я его по той причине, что с времен Фортрана уже много чего придумали и добавили в языки программирования; я надеялся на нечто более современное.
  2. R страдает от тех же проблем, что и Python (векторизация).
  3. Matlab — наверное да, у меня к сожалению денег нет, чтобы проверить.
  4. Julia — лошадка тёмная, взлетит-не взлетит (а было это естественно до stable-версии 1.0)

Некоторые очевидные плюсы Julia

  1. Похож на Python, во всяком случае такой же «высокоуровневый», с возможностью, если надо, спуститься до битов в числах.
  2. Нет возни с аллокациями памяти и тому подобными вещами.
  3. Мощная система типов. Типы прописываются опционально и очень дозированно. Программу можно написать, не указав ни единого типа — и если делать это УМЕЮЧИ, то программа будет быстрая. Но есть нюансы.
  4. Легко писать пользовательские типы, которые будут такими же быстрыми, как и встроенные.
  5. Multiple dispatch. Об этом можно говорить часами, но на мой взгляд — это самое лучшее, что есть в Julia, она есть основа дизайна всех программ и вообще основа философии языка.
    Благодаря multiple dispatch многие вещи пишутся очень единообразно.

    Примеры с multiple dispatch

    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.

  6. 6.Массивы родные, встроенные. Векторизация из коробки.

Минусы, умолчать о которых было бы преступлением!

  1. Новый язык. С одной стороны, конечно, минус. В чем-то, возможно, незрелый. С другой стороны — учёл грабли многих прошлых языков, стоит «на плечах гигантов», вобрал в себя много интересного и нового.
  2. Сразу же быстрые программы вряд ли получится писать. Есть некоторые неочевидные вещи, с которыми бороться очень легко: ты наступаешь на грабли, спрашиваешь помощи у сообщества, опять наступаешь… и т.д. В основном это type instability, нестабильность типов и global variables. В целом, насколько я могу судить по себе, программист, который хочет быстро писать на Julia, проходит несколько этапов: а) пишет как на Python. Это здорово, и так можно, но иногда будут проблемы с производительностью. б) пишет как на Си: маниакально прописывая типы везде, где только можно. в) постепенно понимает, где нужно (очень дозированно) прописывать типы, а где они мешают.
  3. Экосистема. Некоторые пакеты — сырые, в том смысле, что где-то постоянно что-то отваливается; некоторые зрелые, но несостыкуются друг с другом (необходима переконвертация в другие типы, к примеру). С одной стороны это очевидно плохо; с другой, в Julia есть много интересных и смелых идей как раз таки потому, что «стоим на плечах гигантов» — накоплен колоссальный опыт наступания на грабли и «как делать не надо», и это учитывается (частично) разработчиками пакетов.
  4. Нумерация массивов с 1. Ха, шучу, это конечно же плюс! Нет, ну серьезно, какое дело, с какого индекса начинается нумерация? К этому привыкаешь за 5 минут. Никто же не жалуется, что целые числа называются integer, а не whole. Это всё вопросы вкуса. И да, возьмите хотя бы того же Кормена по алгоритмам — там всё нумеруют с единицы, и переделывать на Python иногда наоборот бывает неудобно.

Ну а чего по скорости-то?

Пора бы и вспомнить, ради чего всё затевалось.

Примечание: Python я благополучно забыл, поэтому пишите в комментарии свои улучшения, я постараюсь их запустить на своём ноутбуке и посмотреть время исполнения.

Итак, два маленьких примера с микробенчмарками.

Нечто векторизованное

Задача: на вход функции подаётся вектор X из 0 и 1. Необходимо преобразовать его в вектор из 1 и (-1) (1 ->1, 0 -> -1) и подсчитать, сколько коэф-тов из преобразования Фурье данного вектора по абсолютному значению превосходят границу boundary.


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

Ничего удивительного мы здесь и не увидим — все всё равно считают не сами, а передают в хорошо оптимизированные фортран-программы.

Нечто не-векторизуемое

На вход подается массив из 0 и 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 секунд. Вполне может быть, что это я дурак, и было пространство для оптимизаций и т.д. и т.п. Но это я, каков я есть, со всеми своими достоинствами и недостатками, и на мой взгляд надо считать не сферическую скорость программ в вакууме на языках программирования, а как я лично могу запрограммировать это (без совсем уж грубых ошибок).

Зачем я всё это написал?

Людям здесь стало интересно; я решил написать, когда время появится. В дальнейшем думаю пройтись по Julia более основательно, с примером решения каких-то типовых задачек. Зачем? Больше языков, хороших и разных, и пусть каждый желающий найдет что-то, что подходит лично ему.

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть