Comments 38
А если вместо векторов использовать массивы, то думаю, что C++ будет работать ещё быстрее.
И тут как раз возникает ряд возможных решений с той иной степенью красоты:
- отказаться от итераций и работать только с матричными операциями (NumPy)
- свести задачу к набору стандартных базовых операций (SciPy)
- спуститься к уровню C/C++ (и только здесь возникает С++ код)
- добавить информацию о типах для JIT компилятора
И я как раз пробую насколько хорошо работает каждое из этих решений.
Python интерпретируемый язык, Python использует внутри себя GIL, что делает невозможным параллелизацию на уровне самого кода
И по всей статье ни одного (даже плохого примера) с параллельностью. Есть и multiprocessing, и concurrent futures…
Да и тема «чистого» питона не раскрыта, можно генератор написать, можно map функцию использовать.
В общем посыл не понятен, что Вы хотели этой статьей показать, что питон медленный и не «параллельный»? Так мы это и так знаем)
"… заменить 'range' на 'prange' для параллельного выполнения внешнего цикла"?
А можно ли на чистом Python писать на самом деле быстрый и параллельный код?
Вот с этого места.
не увидел «чистого» питона в статье, уж простите.
Возможно у нас просто разные взгляды на «чистый питон», тогда я уточню. Для меня это пайтон, который ставится из коробки и используется только функционал его встроенных модулей.
Что для вас наиболее чистое решение этой задачи?
Если же требуется решение в реальном времени, то задачи с недетерминированным временем в рантайм выносить нельзя. Сегодня у Вас 5000 точек, завтра будет 10к, после -1млн. Очевидно что человек, выполнивший такую операцию повесит сервер для всех.
Что касается реализации, о которой вы спросили, то лично я бы таки выносил это в очередь, и делал на concurrent futures, с очередью, что касается чистого пайтона. Если бы скорость меня не устроила, я бы написал расширение DLL/SO на С++/С и просто бы его использовал и пайтона. Возможно это не самый простой способ, ноя не люблю полагаться на сторонний код, если реализация подобного занимает немного времени.
Второй момент, я не придирался к вашей реализации, возможно в жизни она действительно хороша и эффективна, я даже не сомневаюсь в этом. Единственное что меня смутило, что Вы обозвали питон не параллельным из коробки, указали примеры на сторонних библиотеках, при этом не привели решения, который таки дает «коробочный» питон.
Думаю, с транспонированием q что-то не так. Зачем там транспонирование — не понятно. cdist принимает на вход данные в виде такой матрицы:
[(x1, y1, z1, ..., n1),
(x2, y2, z2, ..., n2),
...
(xK, yK, zK, ..., nK)
]
Поэтому, если у автора данные хранятся именно в такой форме и если транспонировать один из наборов данных, то будет ошибка:
ValueError: XA and XB must have the same number of columns (i.e. feature dimension.)
import numpy as np
N = 5000
p = np.random.rand(N, 1)
q = np.random.rand(N, 1)
# no broadcasting
print((p - q).shape)
>>> (5000, 1)
# with broadcasting
print((p - q.T).shape)
>>> (5000, 5000)
Сейчас я попробую сделать тест без транспонирования и погляжу на результаты.
А не пробовали запустить тест на PyPy? Интересно посмотреть, что он может продемонстрировать в сравнении с другими участниками тестирования.
И если есть возможность, поделитесь кодом, который использовался для тестирования, я бы с удовольствием поигрался у себя на компьютере.
Спасибо за тест!
Код я выложил:
github.com/alkalinin/open-nuance/tree/master/src-cpp/tests
github.com/alkalinin/open-nuance/tree/master/src-py/open_nuance/tests
Как насчёт того, чтобы попробовать сравнить еще и Intel Distribution for Python?
USE IFPORT
REAL(4) p(5000,3),q(5000,3),R(5000,5000)
integer i,o
CALL SRAND(100)
do j=1,1000
N = 5000
do i=1,N
do o=1,3
p(i,o)=RAND()
q(i,o)=RAND()
end do
end do
CALL CPU_time(t1)
DO CONCURRENT (i=1:size(p,1))
DO CONCURRENT (o=1:size(q,1))
rx=p(i,1)-q(o,1)
ry=p(i,2)-q(o,2)
rz=p(i,3)-q(o,3)
R(i,o)=1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
end do
end do
CALL CPU_time(t2)
print *, t1, t2, t2-t1
dt=dt+(t2-t1)
print *, R(1,1) !чтобы компилятор не халтурил и выполнял всю работу
end do
print *, 'Avg time:', dt/1000.
read(*,*)
1) real(8) будет много шустрее работать
2) print обязательно выносить из-под цикла
Можно еще поиграть с опциями, я как-то «оптимизировал» код методом подбора параметров компиляции с 80+ секунд до 29…
- real(8) — плавающая точка с двойной точностью, такие операции выполняются медленнее. У меня вышло 96 мс.
- принт в данном случае выносить из цикла не надо, так как затраты на вывод информации не учитываются (считается время исполнения кода между двумя вызовами cpu_time)
А оптимизацией я заниматься не буду, так как оптимизируя опции компилятора вкупе с использованием старого кода, досконально не разбираясь в каждом привнесенном или унесенном эффекте, легко что-нибудь сломать (например, изменив опции qsave и qzero можно получить совсем неадекватные результаты, если в коде все переменные «инициализируются» implicit'ом), а в моем случае скорость работы кода не критична, но все должно считаться правильно)
А вот про замедление на real(8) это неожиданно, у меня он самый быстрый…
Надо смотреть условия, код, настройки компилятора, при которых real(8) самый быстрый, там точно что-то не так и это значит, что не реал8 быстрый, а реал4 слишком медленный
Upd. все расчёты велись для double, float выдаёт 041мс
Я конечно тоже считаю numpy неотъемлемой частью питона и не буду городить списки списков для двумерных массивов чисел, но всё-таки код на чистом питоне будет выглядеть как-то так:
def get_R_pure_python(p, q):
R = [[0 for x in range(len(q))] for y in range(len(p))]
for i in range(len(p)):
for j in range(len(q)):
rx = p[i][0] - q[j][0]
ry = p[i][1] - q[j][1]
rz = p[i][2] - q[j][2]
R[i][j] = 1 / (1 + math.sqrt(rx * rx + ry * ry + rz * rz))
return R
Кроме того, этот код быстрее чем numpy без broadcasting и его уже намного проще подсунуть pypy. В итоге на i5-2500 я получил следующие результаты:
| numpy 1.13.1 (без broadcasting) | 42532 мс | 477.89x |
| python 3.5.2 | 22247 мс | 249.97x |
| python 2.7.12 | 17773 мс | 199.70x |
| pypy 5.1.2 | 621 мс | 6.98x |
| numpy 1.13.1 (с broadcasting) | 511 мс | 5.74x |
| scipy 0.19.1 | 280 мс | 3.15x |
| numba 0.35.0rc1 | 89 мс | 1x |
Для python 2.7.12 функцию range я заменил на xrange там где это надо.
P.S. У кого-нибудь есть возможность поправить комментарий, чтобы заработал Markdown, потому что у меня, к сожалению, не работает.
Начнем тест скорости с кода на чистом Python:
R = np.empty((p.shape[0], q.shape[1]))
это разве не Numpy методы?
т.е. это уже не совсем чистый python?
/занудамодеоф
Быстрый тест производительности Python для вычислительных задач