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

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

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

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

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

Чем хорош numpy?

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

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

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

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 я благополучно забыл, поэтому пишите в комментарии свои улучшения, я постараюсь их запустить на своём ноутбуке и посмотреть время исполнения.

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

 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 

Хочу заметить также, что представленные версии не вошли в итоговую программу, а были далее оптимизированы; здесь же они приведены для примера и без ненужных усложнений (пробрасывание дополнительной памяти в Julia, чтобы делать rfft in-place и т.д.).

Что вышло в итоге?

Итоговый код для Python и для Julia я не покажу: это тайна (по крайней мере пока что). На момент, когда я бросил допиливать python-версию, она отрабатывала все NIST-тесты на массиве длины 16 мегабайт двоичных знаков за ~ 50 секунд. На Julia на данный момент все тесты прогоняются на том же объеме ~ 20 секунд. Вполне может быть, что это я дурак, и было пространство для оптимизаций и т.д. и т.п. Но это я, каков я есть, со всеми своими достоинствами и недостатками, и на мой взгляд надо считать не сферическую скорость программ в вакууме на языках программирования, а как я лично могу запрограммировать это (без совсем уж грубых ошибок).

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

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