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


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

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

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

Чем хорош numpy?


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

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

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

Комментарии (21)


  1. denaspireone
    26.10.2018 17:37

    мощно


  1. robert_ayrapetyan
    26.10.2018 18:26
    +1

    То вы сетуете на производительность, то не хотите возиться с аллокациями… Вы уж определитесь. Увеличили с 50 сек до 20, проверьте что ли упрощенный вариант на С, вдруг там в 100 раз быстрее, а вы его ругаете?


    1. kirtsar Автор
      26.10.2018 22:00

      Да, в этом и был смысл: можно ли писать «высокоуровнево», не возясь с низкоуровневыми вещами, при этом не сильно теряя в скорости?
      Это интересный вопрос: а что, если писать на Си? К сожалению, я сам плохо пишу на Си, чтобы сравнивать. Могу привести несколько интересных ссылок со сравнениями (в том числе с Си):
      1. www.ibm.com/developerworks/community/blogs/jfp/entry/A_Comparison_Of_C_Julia_Python_Numba_Cython_Scipy_and_BLAS_on_LU_Factorization?lang=en
      2. modelingguru.nasa.gov/docs/DOC-2625
      3. julialang.org/benchmarks (бенчмарки с официальной страницы Julia).
      Проблемы всех бенчмарков в том, что их пишут либо предвзятые люди, либо специалисты в одном языке. По бенчмаркам видно, что 100х очень вряд ли…


  1. bfDeveloper
    26.10.2018 19:13

    Простите, что встреваю со своими тараканами, Julia действительно хорошо подходит для задач обработки данных и статистики, но не могу не порекомендовать язык D.
    Вы получите: производительность C, можно не возиться с аллокациями (а можно и возиться, если это оправдано), хорошая библиотека алгоритмов (std.algorithm), mir — одна из самых быстрых либ для матриц, алгоритмов и линейной алгебры, во многом вдохновлена numpy, и много-много всего.


    1. gecube
      26.10.2018 19:17
      +3

      А в rust ничего подобного нет? Не совсем понимаю целевую аудиторию D


    1. evocatus
      26.10.2018 19:42

      А почему бы сразу не CUDA? (естественно, через библиотеку типа cuBlas)


  1. DoctorMoriarty
    26.10.2018 19:29
    -1

    вместо решения задачи возишься с аллокациями и деструкторами

    … что, на самом деле, означает вот что:

    "я ниасилил аллокации и деструкторы"

    Никто же не жалуется, что целые числа называются integer, а не whole. Это всё вопросы вкуса

    Нет, это не вопросы вкуса. Это разница между традициями выбора научной лексики в английском и русском языком. Название для целых чисел в английском языке взято еще ~ в 17 веке из латинского языка, от in+tangere, «ненарушенный, нетронутый, чистый», а не образовано от собственно английского слова со значением «целый».


    1. kirtsar Автор
      26.10.2018 20:40

      1. Думаю, вы верно всё прочитали между строк (это я о «не осилил»). Скажем так: я не думаю, что это очень сложная задача, тем более что есть много качественной литературы по языку Си и С++. Вопрос в том, действительно ли необходимо это знать тем, кто не занимается профессионально программированием, а использует его для каких-то своих целей?
      Раньше вы могли сказать: я хочу, чтоб писать было легко, как на python, а летало бы не сильно медленнее Си, и вам бы справедливо ответили: «ты либо крест сними,… ». Я утверждаю, что теперь вы можете достичь скорости, сравнимой со скоростью Си, при этом не опускаясь на несколько уровней абстракции ниже, чем это требуется для решения задачи.

      2. Спасибо за ликбез. Кстати говоря, «whole number» в английском языке есть, и так говорят. В учебниках правда не встречал.


      1. DoctorMoriarty
        26.10.2018 21:26

        между строк

        Скажем так, если бы вы так и написали, что это вызвало затруднения лично у вас и в связи со спецификой вашей деятельности — саркастического комментария бы не было :-) В той же форме, как это написано в тексте — это воспринимается (несмотря на дисклаймер в следующем предложении) как «объективная» претензия в духе «зачем это нужно вообще?». Я бы на вашем месте написал что-то вроде «но вместо того, чтобы просто заняться решением моей задачи, я, увы, должен был бы разбираться со специфической для этого языка необходимостью корректного управления памятью».

        и так говорят. (...) В учебниках

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


  1. Nikles
    26.10.2018 19:38

    Небольшое замечание: я бы убрал лишнюю операцию их исходного алгоритма.

    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
            else:
                maxLen = max(maxLen, currLen)
                currLen = 0
        return max(maxLen, currLen)


  1. iroln
    26.10.2018 23:58
    +1

    50/20 секунд — это не особо впечатляющий прирост производительности (2.5 раза), не находите?


    Вот тут есть некоторые занимательные бенчи
    https://hackernoon.com/performance-analysis-julia-python-c-dd09f03282a3


    И ещё есть такой комментарий на HN


    This has been my experience as well: Python-style Julia code often does not run faster than in Python. Only of you rewrite it in C-style do you get C-like performance.
    This is not necessarily a bad thing. Just an observation from me translating my Python code to Julia.

    https://news.ycombinator.com/item?id=17204750


    1. slovak
      27.10.2018 00:17

      А не имеет ли С style тех же нeдостатков, что и cython?


      1. iroln
        27.10.2018 00:33

        Cython заметно отличается от Python, это отдельные файлы со своим расширением, надо возиться для настройки окружения для его компиляции, интеграции скомпилированных модулей с программой и т. п., а указание типов в Julia не выходит за рамки программирования именно на Julia.


    1. kirtsar Автор
      27.10.2018 00:23

      Ну смотрите, по бенчмаркам из первого примера:
      1. Версия 0.4.5, а сейчас уже 1.0.
      2. (и это гораздо важнее) — замер времени с помощью tic-toc. Это неправильно из-за JIT-компиляции. Надо замерять время после «прогрева» и усреднять по большому количеству запусков; для этого в питоне есть %timeit (точнее в ipython); в Julia для этого есть макрос @benchmark.
      3. В С-программе массив заранее аллокируется; в Julia используется динамический, без sizehint даже.
      В целом, как мне кажется, я честно предупредил в статье, что нужен некоторый навык для написания программ на Julia (как и на любом другом языке).

      2. По второму комментарию. Я привел пример в самой статье когда я буквально написал один и тот же код на Python и на Julia, второй работал гораздо быстрее. Это был python-style или c-style? Я на это могу возразить только своими личными ощущениями от написания кода на julia: на ней можно писать быстро, и код по моим ощущениям ближе (гораздо) к python.

      Согласен, возможно в этой задаче не такой впечатляющий. По мне это довольно ощутимо. Я почти уверен, что на Си эта граница не сильно упадёт, просто из соображений сложности некоторых алгоритмов (дискретное преобразование Фурье, алгоритм Берлекэмпа-Месси и другие затратные тесты).


      1. iroln
        27.10.2018 00:27

        В плане работы JIT было бы интересно сравнить Julia и PyPy. :)


  1. eyeofhell
    27.10.2018 10:19

    ИМХО, ускорение в 2.5 раза не стоит перехода на новый язык. Учитвая то, насколько Python силен экосистемой и батарейками.


    1. kirtsar Автор
      27.10.2018 16:13
      +2

      Возможно и так, про экосистему согласен, но на Julia тоже появились довольно интересные проекты (на мой взгляд), и я не знаю, есть ли аналоги в Питоне, и если есть, то насколько они продвинутые?
      К примеру:
      1. flux (https://github.com/FluxML/Flux.jl) — машинное и глубокое обучение с возможностью писать очень разноуровневый код на pure Julia (от сложений на GPU до описания сетей в несколько строк).
      2. DifferentialEquations (http://docs.juliadiffeq.org/latest/) — решатель всевозможных дифф уравнений, в т.ч. стохастических и прочей экзотики.
      3. JuliaOpt (https://www.juliaopt.org/) — всяческие задачи оптимизации, в том числе DSL для постановки и решения оптимизационных задач — JuMP.
      4. Грамматика графики Gadfly (http://gadflyjl.org/stable/).
      5. Как любитель математики не могу не отметить связку Nemo-Singular-Hecke (http://nemocas.org/) — это конечно не SAGE, но тоже много интересного, от конечных полей до довольно нетривиальных вещей типа алгебраической теории чисел.
      6. Некоторые «маленькие» вещи мне больше нравится в Julia, а именно — задание всяческих вероятностных распределений (пакет Distributions), аналог pandas (DataFrames).

      Тем не менее, есть некоторые вещи, которые есть в Python и нет в Julia, и по которым я скучаю. К примеру, всякие простые модели машинного обучения мне больше нравятся в sklearn. Тут стоит отметить, что питоний код легко запускать из Julia (см., к примеру, на такое: github.com/cstjean/ScikitLearn.jl). Но по мне это как-то немного глупо — писать на одном языке, всё время используя штуки из другого )


  1. alec_kalinin
    27.10.2018 14:22
    +7

    Введение
    Я был одним из тех, кто просил написать статью про переход с Python на Julia. Большое спасибо!


    Прежде всего я хочу согласиться с автором статьи. Да, впихнуть невекторизуемые операции в NumPy порой очень непросто. Cython это страшно, а Numba это непредсказуемо. Именно поэтому я пару лет назад чуть не ушел с Python на что-то другое. Но в то время чего-то другого я не нашел, и остался на Python.


    Сейчас я выступлю в роли защитника Python и все-таки продолжу утверждать, что будущее численных расчетов все-таки за Python


    Улучшаем код на Python
    Я взял код из статьи за основу и провел свои собственные тесты. Длина выборки N=10^6 элементов типа np.int32.


    Тест 1. Чистый Python.


    def get_longest_0(x):
        maxLen = 0
        currLen = 0
    
        for bit in x:
            if bit == 1:
                currLen += 1
            else:
                maxLen = max(maxLen, currLen)
                currLen = 0
        return max(maxLen, currLen)

    Результаты: 3.167 s


    ОК, это исходная точка.


    Тест 2. Все-таки запихиваем все в NumPy
    Автор статьи правильно начал рассуждать о том, что код можно векторизовать путем сдвига и вычитания. Но решил, что это не очень хорошее решение. Но я все-таки довел это способ до конца.


    def get_longest_1(x):
        x_ext = np.hstack(([0], x, [0]))
        x_difs = np.diff(x_ext)
        lengths = np.where(x_difs < 0)[0] - np.where(x_difs > 0)[0]
        max_len = np.max(lengths)
    
        return max_len

    Результаты: 0.189 s


    Быстрее где-то в 16 раз.


    Кроме скорости выполнения по мне так такой способ гораздо нагляднее, чем первый итерационный способ в цикле. Но для этого нужно сменить свой способ мышления. Нужно думать в векторизованном стиле. Даже скорее в стиле потокового вычислительного графа, где мы отделяем операции о данных. А для операций мы строим полусимвольный граф вычислений. Но об этом позднее.


    Тест 3. И все-таки Numba
    Код, как в Tecт 1, но всего одна аннотация.


    @numba.jit
    def get_longest_2(x):
        maxLen = 0
        currLen = 0
    
        for bit in x:    ...

    и...


    Результаты: 0.009 s


    Еще раз в 20 быстрее, чем NumPy.


    Да, Numba порой непредсказуема, но она очень активно развивается. И сейчас все больше случаев, где она просто хорошо работает.


    Рассуждение о будущем научных вычислений
    Язык Julia можно рассматривать как продолжение линии Fortran, где нам нужно написать числодробилку для каких-нибудь сеточных методов в таком императивном стиле, плюс распараллелить все по типу OpenMP. В этом смысле язык Julia идеален.


    Но проблема то в том, что для современного железа это слишком низкий уровень, не позволяющий достичь максимальной производительности. Современное железно это многоядерные (реально многоядерные) процессоры со своими уровня кэша + GPU. Потоки данных между таким железом нужно гонять достаточно нетривиальным образом для максимальной производительности.


    В стиле Fortran эффективный код для этого написать крайне нетривиально. И старый добрый LAPACK уже не справляется.


    Современный ответ на эту проблему это концепция потоковых полусимвольный графов вычислений. Что-то типа TensorFlow. Мы отделяем граф вычислений от потоков данных. Граф вычислений собираем из готовых блоков, которые автоматически распределяются под разнообразное железо.


    На мой взгляд Julia слишком низкоуровневая для этого. А вот Python подходит под эту концепцию на все 100%.


    А теперь поглядим на Cython и Numba с этой точки зрения. В Numba один и тот же код можно автоматически распараллелить для мультиядерных CPU аннотацией '@numba.jit(nopython=True, parallel=True)', а при желании автоматически все отправить на GPU аннотацией '@cuda.jit'.


    Проект, где Cython зашел идеально, это cupy. Это интерфейс NumPy, но реализованный полностью на GPU. Более того, они предлагают вообще писать гетерогенный код, который с точки исходного кода полностью одинаков, а проблемы как эффективно выполнить его на GPU полностью берет на себя CuPy.


    CuPy один из лучших примеров, как нужно писать на Cython. Cython тут выступает таким клеем между Python и низкоуровневым GPU кодом. Например, вот тут
    https://github.com/cupy/cupy/blob/master/cupy/cuda/cublas.pyx один из немногих случаев, когда я читал Cython код без головной боли.


    Вот в рамках такой концепции мне Python кажется более прогрессивным языком, чем Julia. Хотя может быть я не очень знаю куда дальше собирается идти Julia.


    1. kirtsar Автор
      27.10.2018 23:36

      Спасибо, очень дельный комментарий!
      Смотрите, я для полноты картины у себя тоже запустил.

      # without numba
      %timeit get_longest_1(x)
      14.2 ms ± 14.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
      
      @numba.jit
      def get_longest_2(x):
          x_ext = np.hstack(([0], x, [0]))
          x_difs = np.diff(x_ext)
          lengths = np.where(x_difs < 0)[0] - np.where(x_difs > 0)[0]
          max_len = np.max(lengths)
          return max_len
      
      
      In [14]: %timeit get_longest_2(x)
      12.5 ms ± 247 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
      


      По каким-то причинам моя numba старается не так сильно, как ваша! Примерно та же история была, когда я как раз в первый раз на python писал.
      Замечу, что векторизация совершенно уравняла на данной конкретной задаче шансы.

      По поводу GPU: советую всё же обратить немного внимания, например, на это: github.com/JuliaGPU
      В частности: «Because CuArray is an AbstractArray, it doesn't have much of a learning curve; just use your favourite array ops as usual». AbstractArray — это общий тип всех массивов, что-то наподобие классов типов в Haskell.
      По параллельным вычислениям — тоже кое-что реализовано (см. docs.julialang.org/en/v1/manual/parallel-computing). Т. е. это буквально тоже на уровне языка, делается в один макрос.
      Также обратите внимание на проект Celeste (какие-то астрономические расчеты, я деталей не знаю к сожалению):
      At the 2017 JuliaCon conference, Jeffrey Regier, Keno Fischer and others announced that the Celeste project used Julia to achieve «peak performance of 1.54 petaFLOPS using 1.3 million threads» on 9300 Knights Landing (KNL) nodes of the Cori (Cray XC40) supercomputer (the 5th fastest in the world at the time; 8th fastest as of November 2017). Julia thus joins C, C++, and Fortran as high-level languages in which petaFLOPS computations have been written.


    1. sairus777
      28.10.2018 01:46

      Современный ответ на эту проблему это концепция потоковых полусимвольный графов вычислений.

      На мой взгляд Julia слишком низкоуровневая для этого. А вот Python подходит под эту концепцию на все 100%.

      Из чего вообще следует, что Julia низкоуровневая?
      Честно говоря, нигде не встречал подобных аргументов как недостатков Джулии. Насколько я знаю, у неё нет проблем ни с символьными вычислениями, ни с вычислительным графом — см. IR representation. Да и разработана она была с заделом на параллельность и многопоточность, как напрямую следует из ее описания.

      P.S. Для чистоты эксперимента надо бы добавить результаты на Julia на той же архитектуре…


  1. sairus777
    28.10.2018 02:29

    Как оптимизировали, профилировщиком @profile пользовались?
    Добавляли макросы типа, @inbounds, @simd, @fastmath?
    Был какой-нибудь прирост за их счет?

    Мы с коллегой как-то погоняли С++ и Julia на тупом алгоритме простых чисел — Julia оказалась быстрее на 15-20% (правда, для С++ был gcc, когда перешли на llvm, С++ стал быстрее на ~5%).