Улучшение мин / макс понижающей дискретизации - PullRequest
0 голосов
/ 31 января 2019

У меня есть несколько больших массивов (~ 100 миллионов точек), которые мне нужны для интерактивного построения.Я использую Matplotlib.Построение массивов «как есть» становится очень медленным и бесполезным, так как вы все равно не можете визуализировать столько точек.

Так что я сделал функцию прореживания min / max, которую я привязал к обратному вызову xlim_changed изось.Я выбрал подход мин / макс, потому что данные содержат быстрые всплески, которые я не хочу пропустить, просто просматривая данные.Есть еще обертки, которые обрезают до x-пределов и пропускают обработку при определенных условиях, но соответствующая часть ниже:

def min_max_downsample(x,y,num_bins):
    """ Break the data into num_bins and returns min/max for each bin"""
    pts_per_bin = x.size // num_bins    

    #Create temp to hold the reshaped & slightly cropped y
    y_temp = y[:num_bins*pts_per_bin].reshape((num_bins, pts_per_bin))
    y_out      = np.empty((num_bins,2))
    #Take the min/max by rows.
    y_out[:,0] = y_temp.max(axis=1)
    y_out[:,1] = y_temp.min(axis=1)
    y_out = y_out.ravel()

    #This duplicates the x-value for each min/max y-pair
    x_out = np.empty((num_bins,2))
    x_out[:] = x[:num_bins*pts_per_bin:pts_per_bin,np.newaxis]
    x_out = x_out.ravel()
    return x_out, y_out

Это работает довольно хорошо и достаточно быстро (~ 80 мс на 1e8 баллов и 2kбункера).Отставание очень незначительное, поскольку он периодически пересчитывает и обновляет x & y-данные строки.

Однако моя единственная жалоба - в x-data.Этот код дублирует x-значение левого края каждого бина и не возвращает истинное x-положение пар y min / max.Я обычно устанавливаю количество бинов, чтобы удвоить ширину оси.Таким образом, вы не можете увидеть разницу, потому что мусорные ведра такие маленькие ... но я знаю, что они есть ... и это меня раздражает.

Итак, попытка номер 2, которая возвращает действительные значения x для каждой пары мин / макс.Однако это примерно в 5 раз медленнее.

def min_max_downsample_v2(x,y,num_bins):
    pts_per_bin = x.size // num_bins
    #Create temp to hold the reshaped & slightly cropped y
    y_temp = y[:num_bins*pts_per_bin].reshape((num_bins, pts_per_bin))
    #use argmax/min to get column locations
    cc_max = y_temp.argmax(axis=1)
    cc_min = y_temp.argmin(axis=1)    
    rr = np.arange(0,num_bins)
    #compute the flat index to where these are
    flat_max = cc_max + rr*pts_per_bin
    flat_min = cc_min + rr*pts_per_bin
    #Create a boolean mask of these locations
    mm_mask  = np.full((x.size,), False)
    mm_mask[flat_max] = True
    mm_mask[flat_min] = True  
    x_out = x[mm_mask]    
    y_out = y[mm_mask]  
    return x_out, y_out

На моей машине это занимает примерно 400+ мс, что становится довольно заметно.Итак, мой вопрос в основном, есть ли способ пойти быстрее и обеспечить такие же результаты?Узкое место в основном в функциях numpy.argmin и numpy.argmax, которые немного медленнее, чем numpy.min и numpy.max.

. Возможно, ответ заключается в том, чтобы просто жить с версией №1, поскольку визуальноэто не имеет значения.Или, может быть, попытаться ускорить что-то вроде Cython (который я никогда не использовал).

FYI, использующий Python 3.6.4 в Windows ... пример использования будет выглядеть примерно так:

x_big = np.linspace(0,10,100000000)
y_big = np.cos(x_big )
x_small, y_small = min_max_downsample(x_big ,y_big ,2000) #Fast but not exactly correct.
x_small, y_small = min_max_downsample_v2(x_big ,y_big ,2000) #correct but not exactly fast.

Ответы [ 4 ]

0 голосов
/ 03 февраля 2019

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

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

0 голосов
/ 01 февраля 2019

Так что это не относится к ускорению конкретной функции, о которой идет речь, но показывает несколько способов построения линии с большим количеством точек, достаточно эффективно.Это предполагает, что точки 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)

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

0 голосов
/ 01 февраля 2019

РЕДАКТИРОВАТЬ: добавлено параллельное = True для numba ... еще быстрее

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

Эта версия возвращает правильные значения 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)
0 голосов
/ 31 января 2019

Мне удалось улучшить производительность, используя вывод 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 мс ± 5 мс
  • min_max_downsample_v2: 240 мс ± 8,01 мс
  • min_max_downsample_v3: 164 мс ± 1,23 мс

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

...