Улучшить минимальную/максимальную понижающую дискретизацию

Мне удалось повысить производительность, используя вывод arg(min|max) непосредственно для индексации массивов данных. Это происходит за счет дополнительного вызова np.sort, но ось для сортировки имеет только два элемента (минимальные/максимальные индексы), а общий массив довольно мал (количество ячеек):

def min_max_downsample_v3(x, y, num_bins):
    pts_per_bin = x.size // num_bins

    x_view = x[:pts_per_bin*num_bins].reshape(num_bins, pts_per_bin)
    y_view = y[:pts_per_bin*num_bins].reshape(num_bins, pts_per_bin)
    i_min = np.argmin(y_view, axis=1)
    i_max = np.argmax(y_view, axis=1)

    r_index = np.repeat(np.arange(num_bins), 2)
    c_index = np.sort(np.stack((i_min, i_max), axis=1)).ravel()

    return x_view[r_index, c_index], y_view[r_index, c_index]

Я проверил тайминги для вашего примера и получил:

• min_max_downsample_v1: 110 ms ± 5 ms
• min_max_downsample_v2: 240 ms ± 8.01 ms
• min_max_downsample_v3: 164 ms ± 1.23 ms

Я также проверил возврат сразу после вызовов arg(min|max), и результат был равным 164 мс, т.е. после этого больше не было реальных накладных расходов.

Таким образом, это не касается ускорения конкретной рассматриваемой функции, но показывает несколько эффективных способов построения линии с большим количеством точек. Это предполагает, что точки x упорядочены и выбраны равномерно (или близко к равномерному).

Настраивать

from pylab import *

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

def calc_rand(y, factor):
    split = y[:len(y)//factor*factor].reshape(-1, factor)
    idx = randint(0, split.shape[-1], split.shape[0])
    return split[arange(split.shape[0]), idx]

А вот мин и макс чтобы увидеть огибающую сигнала

def calc_env(y, factor):
    """
    y : 1D signal
    factor : amount to reduce y by (actually returns twice this for min and max)
    Calculate envelope (interleaved min and max points) for y
    """
    split = y[:len(y)//factor*factor].reshape(-1, factor)
    upper = split.max(axis=-1)
    lower = split.min(axis=-1)
    return c_[upper, lower].flatten()

Следующая функция может принимать любой из них и использовать их для уменьшения отрисовываемых данных. По умолчанию фактически снято 5000 точек, что должно намного превышать разрешение монитора. Данные кэшируются после их уменьшения. Память может быть проблемой, особенно при больших объемах данных, но она не должна превышать объем, требуемый исходным сигналом.

def plot_bigly(x, y, *, ax=None, M=5000, red=calc_env, **kwargs):
    """
    x : the x data
    y : the y data
    ax : axis to plot on
    M : The maximum number of line points to display at any given time
    kwargs : passed to line
    """
    assert x.shape == y.shape, "x and y data must have same shape!"
    if ax is None:
        ax = gca()

    cached = {}

    # Setup line to be drawn beforehand, note this doesn't increment line properties so
    #  style needs to be passed in explicitly
    line = plt.Line2D([],[], **kwargs)
    def update(xmin, xmax):
        """
        Update line data

        precomputes and caches entire line at each level, so initial
        display may be slow but panning and zooming should speed up after that
        """
        # Find nearest power of two as a factor to downsample by
        imin = max(np.searchsorted(x, xmin)-1, 0)
        imax = min(np.searchsorted(x, xmax) + 1, y.shape[0])
        L = imax - imin + 1
        factor = max(2**int(round(np.log(L/M) / np.log(2))), 1)

        # only calculate reduction if it hasn't been cached, do reduction using nearest cached version if possible
        if factor not in cached:
            cached[factor] = red(y, factor=factor)

        ## Make sure lengths match correctly here, by ensuring at least
        #   "factor" points for each x point, then matching y length
        #  this assumes x has uniform sample spacing - but could be modified
        newx = x[imin:imin + ((imax-imin)//factor)* factor:factor]
        start = imin//factor
        newy = cached[factor][start:start + newx.shape[-1]]
        assert newx.shape == newy.shape, "decimation error {}/{}!".format(newx.shape, newy.shape)

        ## Update line data
        line.set_xdata(newx)
        line.set_ydata(newy)

    update(x[0], x[-1])
    ax.add_line(line)
    ## Manually update limits of axis, as adding line doesn't do this
    #   if drawing multiple lines this can quickly slow things down, and some
    #   sort of check should be included to prevent unnecessary changes in limits
    #   when a line is first drawn.
    ax.set_xlim(min(ax.get_xlim()[0], x[0]), max(ax.get_xlim()[1], x[1]))
    ax.set_ylim(min(ax.get_ylim()[0], np.min(y)), max(ax.get_ylim()[1], np.max(y)))

    def callback(*ignore):
        lims = ax.get_xlim()
        update(*lims)

    ax.callbacks.connect('xlim_changed', callback)

    return [line]

Вот тестовый код

L=int(100e6)
x=linspace(0,1,L)
y=0.1*randn(L)+sin(2*pi*18*x)
plot_bigly(x,y, red=calc_env)

На моей машине это отображается очень быстро. Масштабирование имеет небольшую задержку, особенно когда оно на большую величину. С панорамированием проблем нет. Использование случайного выбора вместо минимума и максимума немного быстрее и имеет проблемы только при очень высоких уровнях масштабирования.

EDIT: добавлено parallel=True для numba... еще быстрее

В итоге я создал гибрид однопроходной процедуры argmin+max и улучшенной индексации из ответа @a_guest и ссылки на это связанный одновременный вопрос min max.

Эта версия возвращает правильные значения x для каждой пары min/max y и благодаря numba на самом деле немного быстрее, чем "быстрая, но не совсем правильная" версия.

from numba import jit, prange
@jit(parallel=True)
def min_max_downsample_v4(x, y, num_bins):
    pts_per_bin = x.size // num_bins
    x_view = x[:pts_per_bin*num_bins].reshape(num_bins, pts_per_bin)
    y_view = y[:pts_per_bin*num_bins].reshape(num_bins, pts_per_bin)    
    i_min = np.zeros(num_bins,dtype='int64')
    i_max = np.zeros(num_bins,dtype='int64')

    for r in prange(num_bins):
        min_val = y_view[r,0]
        max_val = y_view[r,0]
        for c in range(pts_per_bin):
            if y_view[r,c] < min_val:
                min_val = y_view[r,c]
                i_min[r] = c
            elif y_view[r,c] > max_val:
                max_val = y_view[r,c]
                i_max[r] = c                
    r_index = np.repeat(np.arange(num_bins), 2)
    c_index = np.sort(np.stack((i_min, i_max), axis=1)).ravel()        
    return x_view[r_index, c_index], y_view[r_index, c_index]

Сравнение скоростей с использованием timeit показывает, что код numba примерно в 2,6 раза быстрее и обеспечивает лучшие результаты, чем v1. Это чуть более чем в 10 раз быстрее, чем последовательное выполнение argmin и argmax numpy.

%timeit min_max_downsample_v1(x_big ,y_big ,2000)
96 ms ± 2.46 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit min_max_downsample_v2(x_big ,y_big ,2000)
507 ms ± 4.75 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit min_max_downsample_v3(x_big ,y_big ,2000)
365 ms ± 1.27 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit min_max_downsample_v4(x_big ,y_big ,2000)
36.2 ms ± 487 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Вы пробовали pyqtgraph для интерактивного построения графиков? Он более отзывчив, чем matplotlib.

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