Saya memiliki beberapa data eksperimental yang ada seperti:

x = array([1, 1.12, 1.109, 2.1, 3, 4.104, 3.1, ...])
y = array([-9, -0.1, -9.2, -8.7, -5, -4, -8.75, ...])
z = array([10, 4, 1, 4, 5, 0, 1, ...])

Jika nyaman, kita dapat mengasumsikan bahwa data ada sebagai array 3D atau bahkan panda DataFrame:

df = pd.DataFrame({'x': x, 'y': y, 'z': z})

Interpretasinya adalah, untuk setiap posisi x[i], y[i], nilai beberapa variabel adalah z[i]. Ini tidak disampel secara merata, jadi akan ada beberapa bagian yang "diambil sampelnya secara padat" (misalnya antara 1 dan 1,2 dalam x) dan bagian lain yang sangat jarang (misalnya antara 2 dan 3 di x). Karena itu, saya tidak bisa membuangnya begitu saja ke pcolormesh atau contourf.

Yang ingin saya lakukan adalah membuat sampel ulang x dan y secara merata pada interval tertentu dan kemudian menggabungkan nilai z. Untuk kebutuhan saya, z dapat dijumlahkan atau dirata-ratakan untuk mendapatkan nilai yang berarti, jadi ini tidak menjadi masalah. Upaya naif saya adalah seperti ini:

X = np.arange(min(x), max(x), 0.1)  
Y = np.arange(min(y), max(y), 0.1)
x_g, y_g = np.meshgrid(X, Y)
nx, ny = x_g.shape
z_g = np.full(x_g.shape, np.nan)

for ix in range(nx - 1):
    for jx in range(ny - 1):
        x_min = x_g[ix, jx]
        x_max = x_g[ix + 1, jx + 1]
        y_min = y_g[ix, jx]
        y_max = y_g[ix + 1, jx + 1]
        vals = df[(df.x >= x_min) & (df.x < x_max) & 
                  (df.y >= y_min) & (df.y < y_max)].z.values
        if vals.any():
            z_g[ix, jx] = sum(vals)

Ini berfungsi dan saya mendapatkan hasil yang saya inginkan, dengan plt.contourf(x_g, y_g, z_g) tetapi lambat! Saya memiliki ~20k sampel, yang kemudian saya subsampel menjadi ~800 sampel dalam x dan ~500 dalam y, artinya loop for memiliki panjang 400k.

Apakah ada cara untuk membuat vektor/mengoptimalkan ini? Lebih baik lagi jika ada beberapa fungsi yang sudah melakukan ini!

(Juga menandai ini sebagai MATLAB karena sintaks antara numpy/MATLAB sangat mirip dan saya memiliki akses ke kedua perangkat lunak.)

4
Maro K 20 Agustus 2017, 05:31

2 jawaban

Jawaban Terbaik

Berikut adalah solusi Python dalam vektor menggunakan NumPy broadcasting dan matrix multiplication dengan np.dot untuk bagian pengurangan jumlah -

x_mask = ((x >= X[:-1,None]) & (x < X[1:,None]))
y_mask = ((y >= Y[:-1,None]) & (y < Y[1:,None]))

z_g_out = np.dot(y_mask*z[None].astype(np.float32), x_mask.T)

# If needed to fill invalid places with NaNs
z_g_out[y_mask.dot(x_mask.T.astype(np.float32))==0] = np.nan

Perhatikan bahwa kami menghindari penggunaan meshgrid di sana. Dengan demikian, menghemat memori di sana sebagai jerat yang dibuat dengan meshgrid akan sangat besar dan dalam prosesnya diharapkan mendapatkan peningkatan kinerja.

Pembandingan

# Original app
def org_app(x,y,z):    
    X = np.arange(min(x), max(x), 0.1)  
    Y = np.arange(min(y), max(y), 0.1)
    x_g, y_g = np.meshgrid(X, Y)
    nx, ny = x_g.shape
    z_g = np.full(np.asarray(x_g.shape)-1, np.nan)

    for ix in range(nx - 1):
        for jx in range(ny - 1):
            x_min = x_g[ix, jx]
            x_max = x_g[ix + 1, jx + 1]
            y_min = y_g[ix, jx]
            y_max = y_g[ix + 1, jx + 1]
            vals = z[(x >= x_min) & (x < x_max) & 
                      (y >= y_min) & (y < y_max)]
            if vals.any():
                z_g[ix, jx] = sum(vals)
    return z_g

# Proposed app
def app1(x,y,z):
    X = np.arange(min(x), max(x), 0.1)  
    Y = np.arange(min(y), max(y), 0.1)
    x_mask = ((x >= X[:-1,None]) & (x < X[1:,None]))
    y_mask = ((y >= Y[:-1,None]) & (y < Y[1:,None]))

    z_g_out = np.dot(y_mask*z[None].astype(np.float32), x_mask.T)

    # If needed to fill invalid places with NaNs
    z_g_out[y_mask.dot(x_mask.T.astype(np.float32))==0] = np.nan
    return z_g_out

Seperti yang terlihat, untuk pembandingan yang adil, saya menggunakan nilai array dengan pendekatan asli, karena mengambil nilai dari kerangka data dapat memperlambat segalanya.

Waktu dan verifikasi -

In [143]: x = np.array([1, 1.12, 1.109, 2.1, 3, 4.104, 3.1])
     ...: y = np.array([-9, -0.1, -9.2, -8.7, -5, -4, -8.75])
     ...: z = np.array([10, 4, 1, 4, 5, 0, 1])
     ...: 

# Verify outputs
In [150]: np.nansum(np.abs(org_app(x,y,z) - app1(x,y,z)))
Out[150]: 0.0

In [145]: %timeit org_app(x,y,z)
10 loops, best of 3: 19.9 ms per loop

In [146]: %timeit app1(x,y,z)
10000 loops, best of 3: 39.1 µs per loop

In [147]: 19900/39.1  # Speedup figure
Out[147]: 508.95140664961633
2
Divakar 20 Agustus 2017, 06:17

Berikut ini adalah solusi MATLAB:

X = min(x)-1 :.1:max(x)+1; % the grid needs to be expanded slightly beyond the min and max
Y = min(y)-1 :.1:max(y)+1;
x_o = interp1(X, 1:numel(X), x, 'nearest');
y_o = interp1(Y, 1:numel(Y), y, 'nearest');
z_g = accumarray([x_o(:) y_o(:)], z(:),[numel(X) numel(Y)]);
1
rahnema1 20 Agustus 2017, 07:55