Оптимизация предложений для части кода Джулии и Python

Question

У меня есть программа для симуляции, и внутри программы у меня есть функция.Я понял, что функция потребляет большую часть времени моделирования.Итак, я пытаюсь оптимизировать функцию в первую очередь.Функция выглядит следующим образом:

Julia version 1.1:

function fun_jul(M,ksi,xi,x)

  F(n,x) = sin(n*pi*(x+1)/2)*cos(n*pi*(x+1)/2);
  K = length(ksi);
  Z = zeros(length(x),K);
  for n in 1:M
    for k in 1:K
        for l in 1:length(x)
          Z[l,k] +=  (1-(n/(M+1))^2)^xi*F(n,ksi[k])*F(n,x[l]);
        end
    end
  end

return Z

end

Я также переписываю вышеуказанную функцию в python + numba для сравнения следующим образом

Python + numba

import numpy as np
from numba import prange, jit    
@jit(nopython=True, parallel=True)
def fun_py(M,ksi,xi,x):

    K = len(ksi);
    F = lambda nn,xx: np.sin(nn*np.pi*(xx+1)/2)*np.cos(nn*np.pi*(xx+1)/2);

    Z = np.zeros((len(x),K));
    for n in range(1,M+1):
        for k in prange(0,K):
            Z[:,k] += (1-(n/(M+1))**2)**xi*F(n,ksi[k])*F(n,x);


    return Z

Но коды Джулии очень медленные, вот мои результаты:

Результаты Джулии:

using BenchmarkTools
N=400; a=-0.5; b=0.5; x=range(a,b,length=N); cc=x; M = 2*N+100; xi = M/40;
@benchmark fun_jul(M,cc,xi,x)

BenchmarkTools.Trial: 
  memory estimate:  1.22 MiB
  allocs estimate:  2
  --------------
  minimum time:     25.039 s (0.00% GC)
  median time:      25.039 s (0.00% GC)
  mean time:        25.039 s (0.00% GC)
  maximum time:     25.039 s (0.00% GC)
  --------------
  samples:          1
  evals/sample:     1

Результаты Python:

N=400;a = -0.5;b = 0.5;x = np.linspace(a,b,N);cc = x;M = 2*N + 100;xi = M/40;

%timeit fun_py(M,cc,xi,x);
1.2 s ± 10.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Любая помощь поприветствуется улучшение кодов как для julia, так и для python + numba.

Обновлено

На основе ответа @Przemyslaw Szufel и других постов я улучшил коды numba и julia,Теперь оба распараллелены.Вот времена

Время Python + Numba:

@jit(nopython=True, parallel=True)
def fun_py(M,ksi,xi,x):

    K = len(ksi);
    F = lambda nn,xx: np.sin(nn*np.pi*(xx+1)/2)*np.cos(nn*np.pi*(xx+1)/2);

Z = np.zeros((K,len(x)));
for n in range(1,M+1):
    pw = (1-(n/(M+1))**2)**xi; f=F(n,x)
    for k in prange(0,K):
        Z[k,:] = Z[k,:] + pw*F(n,ksi[k])*f;


    return Z


N=1000; a=-0.5; b=0.5; x=np.linspace(a,b,N); cc=x; M = 2*N+100; xi = M/40;

%timeit fun_py(M,cc,xi,x);
733 ms ± 13.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Время Джулии

N=1000; a=-0.5; b=0.5; x=range(a,b,length=N); cc=x; M = 2*N+100; xi = M/40;
@benchmark fun_jul2(M,cc,xi,x)

BenchmarkTools.Trial: 
  memory estimate:  40.31 MiB
  allocs estimate:  6302
  --------------
  minimum time:     705.470 ms (0.17% GC)
  median time:      726.403 ms (0.17% GC)
  mean time:        729.032 ms (1.68% GC)
  maximum time:     765.426 ms (5.27% GC)
  --------------
  samples:          7
  evals/sample:     1

Answer 1

Я получил 300 мс в одном потоке (вместо 28 с на моей машине) со следующим кодом.

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

Вот код:

function fun_jul2(M,ksi,xi,x)
  F(n,x) = sin(n*pi*(x+1))/2;
  K = length(ksi);
  L = length(x);
  Z = zeros(length(x),K);
  for n in 1:M
    F_im1= [F(n,ksi[k]) for k in 1:K]
    F_im2 = [F(n,x[l]) for l in 1:L]
    pow = (1-(n/(M+1))^2)^xi    
    for k in 1:K
       for l in 1:L    
          Z[l,k] += pow*F_im1[k]*F_im2[l];
      end
    end
  end
  Z
end

julia> fun_jul2(M,cc,xi,x) ≈ fun_jul(M,cc,xi,x)
true

julia> @time fun_jul2(M,cc,xi,x);
  0.305269 seconds (1.81 k allocations: 6.934 MiB, 1.60% gc time)

** EDIT: с многопоточностью и границами, предложенными DNF:

function fun_jul3(M,ksi,xi,x)
  F(n,x) = sin(n*pi*(x+1))/2;
  K = length(ksi);
  L = length(x);
  Z = zeros(length(x),K);
  for n in 1:M
    F_im1= [F(n,ksi[k]) for k in 1:K]
    F_im2 = [F(n,x[l]) for l in 1:L]
    pow = (1-(n/(M+1))^2)^xi    
    Threads.@threads for k in 1:K
       for l in 1:L    
          @inbounds Z[l,k] += pow*F_im1[k]*F_im2[l];
      end
    end
  end
  Z
end

А теперь время выполнения (не забудьте запустить set JULIA_NUM_THREADS=4 или эквивалент Linux перед запуском Julia):

julia>  fun_jul2(M,cc,xi,x) ≈ fun_jul3(M,cc,xi,x)
true

julia> @time fun_jul3(M,cc,xi,x);
  0.051470 seconds (2.71 k allocations: 6.989 MiB)

Вы также можете попробовать поэкспериментировать с распараллеливанием вычислений F_im1 и F_im2.

Answer 2

Для производительности просто предварительно вычислите тригонометрическую часть. Действительно, sin является дорогостоящей операцией:

%timeit np.sin(1.)
712 ns ± 2.22 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%timeit 1.2*3.4
5.88 ns ± 0.016 ns per loop (mean ± std. dev. of 7 runs, 100000000 loops each)

В питоне:

@jit
def fun_py2(M,ksi,xi,x):

    NN = np.arange(1,M+1)
    Fksi = np.sin(np.pi*np.outer(NN,ksi+1))/2  # sin(a)cos(a) is sin(2a)/2
    Fx = np.sin(np.pi*np.outer(NN,x+1))/2
    U = (1-(NN/(M+1))**2)**xi
    Z = np.zeros((len(x),len(ksi)))

    for n in range(len(NN)):
        for k in range(len(ksi)):
            for l in range(len(x)):
                Z[k,l] += U[n] * Fksi[n,k] * Fx[n,l];
    return Z

Для улучшения в 30 раз:

np.allclose(fun_py(M,cc,xi,x),fun_py2(M,cc,xi,x))
True

%timeit fun_py(M,cc,xi,x)
1.14 s ± 4.47 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit fun_py2(M,cc,xi,x)
29.5 ms ± 375 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Это не вызывает параллелизма. Полагаю, то же самое произойдет и с Юлией.

Answer 3

Вы можете выполнять или не выполнять оптимизацию циклов на любом языке, в котором есть циклы.Основное отличие состоит в том, что код numba векторизован для внутреннего цикла, а код Julia - нет.Чтобы векторизовать версию Julia, иногда необходимо изменить операторы на их векторизованные версии с помощью., Так что + становится, например,. +.

Так как я не могу заставить Numba правильно установить на мой старый компьютер с Windows 10Я запускал приведенные ниже версии кода для бесплатных версий Linux в Интернете.Это означает, что мне пришлось использовать интерфейс Python для timeit (), а не для командной строки.

Запускать в Jupyter на mybinder, возможно, с 1 потоком, поскольку он не указан.:

import timeit
timeit.timeit("""
@jit(nopython=True, parallel=True)
def fun_py(M,ksi,xi,x):

    K = len(ksi);
    F = lambda nn,xx: np.sin(nn*np.pi*(xx+1)/2)*np.cos(nn*np.pi*(xx+1)/2);

    Z = np.zeros((len(x),K));
    for n in range(1,M+1):
        for k in prange(0,K):
            Z[:,k] += (1-(n/(M+1))**2)**xi*F(n,ksi[k])*F(n,x);


    return Z


N=400; a = -0.5; b = 0.5; x = np.linspace(a,b,N); cc = x;M = 2*N + 100; xi = M/40;
fun_py(M,cc,xi,x)
""", setup ="import numpy as np; from numba import prange, jit", number=5)

Out [1]: 61.07768889795989

Ваша машина должна быть намного быстрее, чем Jupyter, ForBonder.

Я запустил эту оптимизированную версию кода Джулии ниже, вJupyter на JuliaBox, указано одноядерное ядро:

using BenchmarkTools

F(n, x) = sinpi.(n * (x .+ 1) / 2) .* cospi.(n * (x .+ 1) / 2)

function fun_jul2(M, ksi, xi, x)
    K = length(ksi)
    Z = zeros(length(x), K)
    for n in 1:M, k in 1:K
        Z[:, k] .+= (1 - (n / (M + 1))^2)^xi * F(n, ksi[k]) * F(n, x)
    end
    return Z
end


const N=400; const a=-0.5; const b=0.5; const x=range(a,b,length=N); 
const cc=x; const M = 2*N+100; const xi = M/40;

@btime fun_jul2(M, cc, xi, x)

8,076 с (1080002 выделения: 3,35 ГиБ)