Создание тепловой карты путем выборки и группировки из трехмерного массива

У меня есть некоторые экспериментальные данные, которые существуют так:

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, ...])

Если это удобно, мы можем предположить, что данные существуют в виде 3D-массива или даже панды DataFrame:

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

Интерпретация заключается в том, что для каждой позиции x[i], y[i] значение некоторой переменной равно z[i]. Они отбираются неравномерно, поэтому некоторые части будут иметь "плотную выборку" (например, от 1 до 1,2 в x), а другие - очень разреженные (например, от 2 до 3 в x). Из-за этого я не могу просто бросить их в pcolormesh или contourf.

Вместо этого я хотел бы сделать повторную выборку x и y равномерно с некоторым фиксированным интервалом, а затем агрегировать значения z. Для моих нужд 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(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)

Это работает, и я получаю желаемый результат с plt.contourf(x_g, y_g, z_g), но это МЕДЛЕННО! У меня есть ~20 тыс. отсчетов, которые я затем разделяю на ~800 отсчетов по x и ~500 отсчетов по y, что означает, что цикл for имеет длину 400 тыс.

Есть ли способ векторизовать/оптимизировать это? Еще лучше, если есть какая-то функция, которая уже это делает!

(Также помечаю это как MATLAB, потому что синтаксис между numpy/MATLAB очень похож, и у меня есть доступ к обоим программам.)


python numpy matplotlib matlab contourf
person Maro K    schedule 20.08.2017    source источник
comment
Возможное решение в pandas: stackoverflow.com/questions/42689070/ (хотя, вероятно, не так эффективно, как решение numpy ниже).   -  person ImportanceOfBeingErnest    schedule 20.08.2017

Ответы (2)


arrow_upward
2
arrow_downward

Вот векторизованное решение Python, использующее NumPy broadcasting и matrix multiplication с np.dot для части уменьшения суммы -

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

Обратите внимание, что здесь мы избегаем использования meshgrid. Таким образом, экономия памяти там, поскольку меши, созданные с помощью meshgrid, будет огромной, и в процессе мы надеемся получить улучшение производительности.

Бенчмаркинг

# 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

Как видно, для справедливого бенчмаркинга я использую значения массива с оригинальным подходом, поскольку выборка значений из фрейма данных может замедлить работу.

Сроки и проверка -

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
person Divakar    schedule 20.08.2017

arrow_upward
1
arrow_downward

Вот решение 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)]);
person rahnema1    schedule 20.08.2017