Cython¶

In [2]:
%load_ext Cython
In [16]:
%%cython

from cython.parallel import prange # параллельный range

from libc.math cimport sin # С-функция


ctypedef double (*func)(double x) nogil # объявление типа функции, передеваемой в integrate
# func - любая функция, принимающая аргумент x (double) и возвращающая double (число двойной точности с плавающей запятой)
# GIL выключен, т.к. из nogil функции (integrate) нельзя вызвать функцию с GIL


# Cython-функция, GIL выключен, т.е. вычисление осуществляется на самом деле параллельно
# все переменные имеют тип double (с float получается неточный ответ)
# f не может быть Python-функцией (т.к. её нельзя вызвать из Cython-функции без GIL и при этом она возвращает значение типа Python, а не Cython)
cdef double integrate(func f, float a, float b, int n_iter = 1000) nogil:
    if b == a:
        return 0

    cdef double h = (b-a)/n_iter
    cdef double z = 0
    cdef double x = a + h

    # исходный вариант (однопоточный)
    # while x <= b - h:
    #     z = z + f(x)
    #     x = x + h
    
    # вариант с параллельным range
    cdef int i
    for i in prange(n_iter):
        z += + f(x + i*h) # к начальному значению x прибавляется смещение по h

    cdef double y = (f(a)+ f(b)) /2
    z = h*(z+y)

    return z


# для timeit
# из другой ячейки (и вообще из Python-функции - def) нельзя вызвать Cython-функцию (cdef)
# cpdef можно вызвать как из Python, так и из Cython-функций
# аргумента f нет, т.к. его нельзя передать из Python (нужна Cython-функция)
# функция объявлена в этой же ячейке из-за области видимости (Cython-функцию можно вызывать только в той же ячейке)
cpdef integrate_f(a, b, n_iter):
    return integrate(sin, a, b, n_iter)
In [9]:
integrate_f(1, 2.5, n_iter=10**6)
Out[9]:
1.341446850875759

Joblib¶

In [19]:
def integrate(f, a: float, b: float, *, n_iter: int = 1000):
    if b == a:
        return 0

    h = (b-a)/float(n_iter)
    z = 0
    x = a + h

    while x <= b - h:
        z = z + f(x)
        x = x + h

    y = (f(a)+ f(b)) /2
    z = h*(z+y)

    return round(z, 8)
In [20]:
from joblib import Parallel, delayed

def integrate_async(f, a: float, b: float, *, n_jobs: int = 2, n_iter: int = 1000, backend='threading'):
    step = (b - a) / n_jobs # интервал делится на части, кол-во частей равно числу потоков
    # каждый поток вычисляет значение интеграла на части интервала от a до b

    with Parallel(n_jobs=n_jobs, backend=backend) as p:
        # каждый поток вычисляет значение интеграла на части интервала от a до b
        fs = (delayed(integrate)(f, a + i * step, a + (i + 1) * step, n_iter=n_iter // n_jobs)
              for i in range(n_jobs))
        return sum(p(fs)) # значение интеграла равно сумме значений, вычисленных для каждой части в потоках
In [21]:
import math
integrate_async(math.sin, 1, 2.5, n_iter=10**6)
Out[21]:
1.34144502