Подход № 1
corr2_coeff_rowwise показывает, как сделать поэлементную корреляцию между строками. Мы могли бы сократить его до случая поэлементной корреляции между двумя столбцами. Таким образом, мы получим цикл, который использует corr2_coeff_rowwise. Затем мы попытаемся векторизовать его и увидеть, что в нем есть фрагменты, которые можно векторизовать:
- Получение средних значений с помощью
mean. Это может быть векторизовано с использованием единого фильтра.
- Затем мы получили различия между этими средними значениями относительно скользящих элементов из входных массивов. Чтобы портировать на векторизованный, мы бы использовали
broadcasting.
Остальное остается прежним, чтобы получить первый из двух выходов от pearsonr.
Чтобы получить второй вывод, мы возвращаемся к source code. Это должно быть просто, учитывая вывод первого коэффициента.
Итак, имея в виду, мы бы получили что-то вроде этого -
import scipy.special as special
from scipy.ndimage import uniform_filter
def sliding_corr1(a,b,W):
# a,b are input arrays; W is window length
am = uniform_filter(a.astype(float),W)
bm = uniform_filter(b.astype(float),W)
amc = am[W//2:-W//2+1]
bmc = bm[W//2:-W//2+1]
da = a[:,None]-amc
db = b[:,None]-bmc
# Get sliding mask of valid windows
m,n = da.shape
mask1 = np.arange(m)[:,None] >= np.arange(n)
mask2 = np.arange(m)[:,None] < np.arange(n)+W
mask = mask1 & mask2
dam = (da*mask)
dbm = (db*mask)
ssAs = np.einsum('ij,ij->j',dam,dam)
ssBs = np.einsum('ij,ij->j',dbm,dbm)
D = np.einsum('ij,ij->j',dam,dbm)
coeff = D/np.sqrt(ssAs*ssBs)
n = W
ab = n/2 - 1
pval = 2*special.btdtr(ab, ab, 0.5*(1 - abs(np.float64(coeff))))
return coeff,pval
Таким образом, чтобы получить окончательный вывод от входов из серии панд -
out = sliding_corr1(df['x'].to_numpy(copy=False),df['y'].to_numpy(copy=False),50)
Подход № 2
Очень похоже на Approach #1, но мы будем использовать numba для повышения эффективности памяти вместо шага № 2 от предыдущего подхода.
from numba import njit
import math
@njit(parallel=True)
def sliding_corr2_coeff(a,b,amc,bmc):
L = len(a)-W+1
out00 = np.empty(L)
for i in range(L):
out_a = 0
out_b = 0
out_D = 0
for j in range(W):
d_a = a[i+j]-amc[i]
d_b = b[i+j]-bmc[i]
out_D += d_a*d_b
out_a += d_a**2
out_b += d_b**2
out00[i] = out_D/math.sqrt(out_a*out_b)
return out00
def sliding_corr2(a,b,W):
am = uniform_filter(a.astype(float),W)
bm = uniform_filter(b.astype(float),W)
amc = am[W//2:-W//2+1]
bmc = bm[W//2:-W//2+1]
coeff = sliding_corr2_coeff(a,b,amc,bmc)
ab = W/2 - 1
pval = 2*special.btdtr(ab, ab, 0.5*(1 - abs(np.float64(coeff))))
return coeff,pval
Подход № 3
Очень похоже на предыдущее, за исключением того, что мы увеличиваем все коэффициенты до numba -
@njit(parallel=True)
def sliding_corr3_coeff(a,b,W):
L = len(a)-W+1
out00 = np.empty(L)
for i in range(L):
a_mean = 0.0
b_mean = 0.0
for j in range(W):
a_mean += a[i+j]
b_mean += b[i+j]
a_mean /= W
b_mean /= W
out_a = 0
out_b = 0
out_D = 0
for j in range(W):
d_a = a[i+j]-a_mean
d_b = b[i+j]-b_mean
out_D += d_a*d_b
out_a += d_a*d_a
out_b += d_b*d_b
out00[i] = out_D/math.sqrt(out_a*out_b)
return out00
def sliding_corr3(a,b,W):
coeff = sliding_corr3_coeff(a,b,W)
ab = W/2 - 1
pval = 2*special.btdtr(ab, ab, 0.5*(1 - np.abs(coeff)))
return coeff,pval
Сроки -
In [181]: df = pd.DataFrame({'x': np.random.rand(10000), 'y': np.random.rand(10000)})
In [182]: %timeit sliding_corr2(df['x'].to_numpy(copy=False),df['y'].to_numpy(copy=False),50)
100 loops, best of 3: 5.05 ms per loop
In [183]: %timeit sliding_corr3(df['x'].to_numpy(copy=False),df['y'].to_numpy(copy=False),50)
100 loops, best of 3: 5.51 ms per loop
Примечание:
sliding_corr1, кажется, занимает много времени на этом наборе данных и, скорее всего, из-за требования к памяти из его шага # 2.
Узкое место после использования функций numba, а затем переходит к вычислению p-val с помощью special.btdtr.