Быстрый тест производительности Python для вычислительных задач

• отказаться от итераций и работать только с матричными операциями (NumPy)
• свести задачу к набору стандартных базовых операций (SciPy)
• спуститься к уровню C/C++ (и только здесь возникает С++ код)
• добавить информацию о типах для JIT компилятора

Думаю, с транспонированием 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)

    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(*,*)