Как я с 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 более основательно, с примером решения каких-то типовых задачек. Зачем? Больше языков, хороших и разных, и пусть каждый желающий найдет что-то, что подходит лично ему.

Автор: kirtsar

Источник