Alternatif cepat untuk numpy.median.reduceat

Berkaitan dengan jawaban ini , apakah ada cara cepat untuk menghitung median atas array yang memiliki grup dengan jumlah elemen yang tidak sama ?

Misalnya:

data =  [1.00, 1.05, 1.30, 1.20, 1.06, 1.54, 1.33, 1.87, 1.67, ... ]
index = [0,    0,    1,    1,    1,    1,    2,    3,    3,    ... ]

Dan kemudian saya ingin menghitung perbedaan antara jumlah dan median per kelompok (misalnya median dari kelompok 0adalah 1.025sehingga hasil pertama adalah 1.00 - 1.025 = -0.025). Jadi untuk array di atas, hasilnya akan muncul sebagai:

result = [-0.025, 0.025, 0.05, -0.05, -0.19, 0.29, 0.00, 0.10, -0.10, ...]

Karena np.median.reduceatbelum ada, apakah ada cara cepat lain untuk mencapai ini? Array saya akan berisi jutaan baris sehingga kecepatan sangat penting!

Indeks dapat dianggap berdampingan dan dipesan (mudah untuk mengubahnya jika tidak).

Contoh data untuk perbandingan kinerja:

import numpy as np

np.random.seed(0)
rows = 10000
cols = 500
ngroup = 100

# Create random data and groups (unique per column)
data = np.random.rand(rows,cols)
groups = np.random.randint(ngroup, size=(rows,cols)) + 10*np.tile(np.arange(cols),(rows,1))

# Flatten
data = data.ravel()
groups = groups.ravel()

# Sort by group
idx_sort = groups.argsort()
data = data[idx_sort]
groups = groups[idx_sort]

Terkadang Anda perlu menulis kode numpy non-idiomatis jika Anda benar - benar ingin mempercepat perhitungan Anda yang tidak dapat Anda lakukan dengan numpy asli.

numbamengkompilasi kode python Anda ke level C. rendah. Karena banyak numpy itu sendiri biasanya secepat C, ini sebagian besar akhirnya berguna jika masalah Anda tidak cocok dengan vektorisasi asli dengan numpy. Ini adalah salah satu contoh (di mana saya berasumsi bahwa indeks berdekatan dan diurutkan, yang juga tercermin dalam contoh data):

import numpy as np
import numba

# use the inflated example of roganjosh https://stackoverflow.com/a/58788534
data =  [1.00, 1.05, 1.30, 1.20, 1.06, 1.54, 1.33, 1.87, 1.67]
index = [0,    0,    1,    1,    1,    1,    2,    3,    3] 

data = np.array(data * 500) # using arrays is important for numba!
index = np.sort(np.random.randint(0, 30, 4500))               

# jit-decorate; original is available as .py_func attribute
@numba.njit('f8[:](f8[:], i8[:])') # explicit signature implies ahead-of-time compile
def diffmedian_jit(data, index): 
    res = np.empty_like(data) 
    i_start = 0 
    for i in range(1, index.size): 
        if index[i] == index[i_start]: 
            continue 

        # here: i is the first _next_ index 
        inds = slice(i_start, i)  # i_start:i slice 
        res[inds] = data[inds] - np.median(data[inds]) 

        i_start = i 

    # also fix last label 
    res[i_start:] = data[i_start:] - np.median(data[i_start:])

    return res

Dan berikut adalah beberapa timing menggunakan %timeitkeajaiban IPython :

>>> %timeit diffmedian_jit.py_func(data, index)  # non-jitted function
... %timeit diffmedian_jit(data, index)  # jitted function
...
4.27 ms ± 109 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
65.2 µs ± 1.01 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Menggunakan contoh data yang diperbarui dalam pertanyaan angka-angka ini (yaitu runtime dari fungsi python vs runtime dari JT-accelerated functio) adalah

>>> %timeit diffmedian_jit.py_func(data, groups) 
... %timeit diffmedian_jit(data, groups)
2.45 s ± 34.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
93.6 ms ± 518 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Ini berarti 65x speedup dalam case yang lebih kecil dan 26x speedup dalam case yang lebih besar (dibandingkan dengan kode gila, tentu saja) menggunakan kode akselerasi. Kelebihan lainnya adalah bahwa (tidak seperti vektorisasi khas dengan numpy asli) kami tidak memerlukan memori tambahan untuk mencapai kecepatan ini, ini semua tentang kode tingkat rendah yang dioptimalkan dan disusun yang akhirnya dijalankan.

Fungsi di atas mengasumsikan bahwa array int numpy secara int64default, yang sebenarnya tidak terjadi pada Windows. Jadi alternatifnya adalah menghapus tanda tangan dari panggilan ke numba.njit, memicu kompilasi tepat waktu yang tepat. Tetapi ini berarti bahwa fungsi tersebut akan dikompilasi selama eksekusi pertama, yang dapat mencampuri hasil pengaturan waktu (kita dapat mengeksekusi fungsi tersebut secara manual, menggunakan tipe data yang representatif, atau hanya menerima bahwa eksekusi pengaturan waktu pertama akan jauh lebih lambat, yang seharusnya diabaikan). Ini persis apa yang saya coba cegah dengan menentukan tanda tangan, yang memicu kompilasi sebelumnya.

Bagaimanapun, dalam kasus JIT yang tepat dekorator yang kita butuhkan adalah adil

@numba.njit
def diffmedian_jit(...):

Perhatikan bahwa timing di atas yang saya perlihatkan untuk fungsi yang dikompilasi jit hanya berlaku setelah fungsi dikompilasi. Ini bisa terjadi pada definisi (dengan kompilasi bersemangat, ketika tanda tangan eksplisit diteruskan ke numba.njit), atau selama panggilan fungsi pertama (dengan kompilasi malas, ketika tidak ada tanda tangan dilewatkan ke numba.njit). Jika fungsi ini hanya akan dieksekusi sekali maka waktu kompilasi juga harus dipertimbangkan untuk kecepatan metode ini. Biasanya hanya layak mengkompilasi fungsi jika total waktu kompilasi + eksekusi kurang dari runtime yang tidak dikompilasi (yang sebenarnya benar dalam kasus di atas, di mana fungsi python asli sangat lambat). Ini sebagian besar terjadi ketika Anda sering memanggil fungsi yang dikompilasi.

Seperti yang dicatat oleh max9111 dalam komentar, salah satu fitur penting numbaadalah cachekata kunci untuk jit. Passing cache=Trueto numba.jitakan menyimpan fungsi yang dikompilasi ke disk, sehingga selama eksekusi berikutnya dari modul python yang diberikan fungsi akan dimuat dari sana daripada dikompilasi ulang, yang lagi-lagi dapat membantu Anda menjalankan runtime dalam jangka panjang.

Salah satu pendekatan akan digunakan di Pandassini murni untuk memanfaatkan groupby. Saya telah menggembungkan ukuran input sedikit untuk memberikan pemahaman yang lebih baik tentang timing (karena ada overhead dalam membuat DF).

import numpy as np
import pandas as pd

data =  [1.00, 1.05, 1.30, 1.20, 1.06, 1.54, 1.33, 1.87, 1.67]
index = [0,    0,    1,    1,    1,    1,    2,    3,    3]

data = data * 500
index = np.sort(np.random.randint(0, 30, 4500))

def df_approach(data, index):
    df = pd.DataFrame({'data': data, 'label': index})
    df['median'] = df.groupby('label')['data'].transform('median')
    df['result'] = df['data'] - df['median']

Memberi yang berikut timeit:

%timeit df_approach(data, index)
5.38 ms ± 50.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Untuk ukuran sampel yang sama, saya mendapatkan pendekatan dict dari Aryerez menjadi:

%timeit dict_approach(data, index)
8.12 ms ± 3.47 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Namun, jika kita menambah input dengan faktor lain 10, waktunya menjadi:

%timeit df_approach(data, index)
7.72 ms ± 85 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit dict_approach(data, index)
30.2 ms ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Namun, dengan mengorbankan beberapa kemudahan, jawaban oleh Divakar menggunakan numpy murni muncul di:

%timeit bin_median_subtract(data, index)
573 µs ± 7.48 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Mengingat dataset baru (yang seharusnya sudah ditetapkan di awal):

%timeit df_approach(data, groups)
472 ms ± 2.52 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit bin_median_subtract(data, groups) #https://stackoverflow.com/a/58788623/4799172
3.02 s ± 31.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit dict_approach(data, groups) #https://stackoverflow.com/a/58788199/4799172
<I gave up after 1 minute>

# jitted (using @numba.njit('f8[:](f8[:], i4[:]') on Windows) from  https://stackoverflow.com/a/58788635/4799172
%timeit diffmedian_jit(data, groups)
132 ms ± 3.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Mungkin Anda sudah melakukan ini, tetapi jika tidak, lihat apakah itu cukup cepat:

median_dict = {i: np.median(data[index == i]) for i in np.unique(index)}
def myFunc(my_dict, a): 
    return my_dict[a]
vect_func = np.vectorize(myFunc)
median_diff = data - vect_func(median_dict, index)
median_diff

Keluaran:

array([-0.025,  0.025,  0.05 , -0.05 , -0.19 ,  0.29 ,  0.   ,  0.1  ,
   -0.1  ])

Berikut adalah pendekatan berbasis NumPy untuk mendapatkan binned-median untuk nilai nampan / indeks positif -

def bin_median(a, i):
    sidx = np.lexsort((a,i))

    a = a[sidx]
    i = i[sidx]

    c = np.bincount(i)
    c = c[c!=0]

    s1 = c//2

    e = c.cumsum()
    s1[1:] += e[:-1]

    firstval = a[s1-1]
    secondval = a[s1]
    out = np.where(c%2,secondval,(firstval+secondval)/2.0)
    return out

Untuk mengatasi kasus spesifik kami yang dikurangkan -

def bin_median_subtract(a, i):
    sidx = np.lexsort((a,i))

    c = np.bincount(i)

    valid_mask = c!=0
    c = c[valid_mask]    

    e = c.cumsum()
    s1 = c//2
    s1[1:] += e[:-1]
    ssidx = sidx.argsort()
    starts = c%2+s1-1
    ends = s1

    starts_orgindx = sidx[np.searchsorted(sidx,starts,sorter=ssidx)]
    ends_orgindx  = sidx[np.searchsorted(sidx,ends,sorter=ssidx)]
    val = (a[starts_orgindx] + a[ends_orgindx])/2.
    out = a-np.repeat(val,c)
    return out