Вычислить интеграл сплайна по кускам более эффективно

У меня есть программа, которая чередует 2 фазы: REST и REG. На каждом этапе по 2 положительный количества рассчитываются a и d. К концу REG фазе, я вычисляю оценку на основе предыдущего REG и REST фазы. Это вычисление очень неуклюже, и я хотел бы более быструю реализацию (особенно, я бы хотел избавиться от цикла for).

Сгенерируйте образцы данных:

#%% Imports
from scipy.interpolate import InterpolatedUnivariateSpline
import numpy as np
from matplotlib import pyplot as plt

#%% Generate random data
timestamps_REST = np.linspace(0, 10, 30, endpoint=True)
a_REST = np.random.rand(30)
d_REST = np.random.rand(30)

timestamps_REG = np.linspace(11, 45, 80, endpoint=True)
a_REG = np.random.rand(80)
d_REG = np.random.rand(80)

#%% Spline interpolation
a_REST_spline = InterpolatedUnivariateSpline(timestamps_REST, a_REST, k=1)
a_REG_spline = InterpolatedUnivariateSpline(timestamps_REG, a_REG, k=1)
d_REST_spline = InterpolatedUnivariateSpline(timestamps_REST, d_REST, k=1)
d_REG_spline = InterpolatedUnivariateSpline(timestamps_REG, d_REG, k=1)

В образцах данных, созданных ниже, количества a и d оцениваются в:

  • 30 точек с равным интервалом для REST фаза
  • 80 точек с равным интервалом для REG фаза

В REST и REG продолжительность фазы / продолжительность разные, 0 to 10 для REST фаза и 11 to 45 для фазы REG.

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

Визуализировать / Графики

Для визуализации данных вы можете использовать приведенный ниже фрагмент:

#%% Plot
f, ax = plt.subplots(2, 2, figsize=(10, 5))
ax[0, 0].set_title('REST', fontsize=20)
ax[0, 1].set_title('REG', fontsize=20)
for a in ax.flatten():
    a.set_yticks([])
ax[0, 0].set_ylabel('a', fontsize=20)
ax[1, 0].set_ylabel('d', fontsize=20)

# Spline x-axis
plot_spline_REST_timestamps = np.linspace(0, 10, 2000, endpoint=True)
plot_spline_REG_timestamps = np.linspace(11, 45, 6800, endpoint=True)

# a REST
ax[0, 0].plot(timestamps_REST, a_REST, 'ro', ms=5)
ax[0, 0].plot(plot_spline_REST_timestamps, a_REST_spline(plot_spline_REST_timestamps), 'b', lw=2)
ax[0, 0].fill_between(plot_spline_REST_timestamps, a_REST_spline(plot_spline_REST_timestamps), color="b", alpha=0.3)

# d Rest
ax[1, 0].plot(timestamps_REST, d_REST, 'ro', ms=5)
ax[1, 0].plot(plot_spline_REST_timestamps, d_REST_spline(plot_spline_REST_timestamps), 'teal', lw=2)
ax[1, 0].fill_between(plot_spline_REST_timestamps, d_REST_spline(plot_spline_REST_timestamps), color="teal", alpha=0.3)

# a REG
ax[0, 1].plot(timestamps_REG, a_REG, 'ro', ms=5)
ax[0, 1].plot(plot_spline_REG_timestamps, a_REG_spline(plot_spline_REG_timestamps), 'b', lw=2)
ax[0, 1].fill_between(plot_spline_REG_timestamps, a_REG_spline(plot_spline_REG_timestamps), color="b", alpha=0.3)

# d REG
ax[1, 1].plot(timestamps_REG, d_REG, 'ro', ms=5)
ax[1, 1].plot(plot_spline_REG_timestamps, d_REG_spline(plot_spline_REG_timestamps), 'teal', lw=2)
ax[1, 1].fill_between(plot_spline_REG_timestamps, d_REG_spline(plot_spline_REG_timestamps), color="teal", alpha=0.3)

Пример данных

Вычисление баллов: код для оптимизации

Как уже было сказано, я вычисляю оценку на основе интегралов. Идея состоит в том, чтобы вычислить оценку как:

max(0, Area(a_REG) - Area(a_REST)) + max(0, Area(d_REST) - Area(d_REG))

где Площадь представляет собой окрашенную область на графике выше, то есть площадь между кривой и 0, то есть интегралом. Проще говоря, оценка увеличивается, если:

  • область в синий в REG фаза больше, чем площадь в синий в REST фаза
  • область в бирюзовый в REG фаза меньше площади в бирюзовый в REG фаза

NB: при необходимости интегральное вычисление можно изменить на другой метод.

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

Версия фрагмента ниже нуждается в некоторой оптимизации.

#%% Score
score = 0

# compute baseline for REST phase
rest_chunks = np.linspace(timestamps_REST[0], timestamps_REST[-1], num=11, endpoint=True)
rest_baselines = [[], []]
for k, _ in enumerate(rest_chunks):
    if k == 0:
        continue
    
    rest_baselines[0].append(a_REST_spline.integral(rest_chunks[k-1], rest_chunks[k]) / (rest_chunks[k] - rest_chunks[k-1]))
    rest_baselines[1].append(d_REST_spline.integral(rest_chunks[k-1], rest_chunks[k]) / (rest_chunks[k] - rest_chunks[k-1]))

rest_baseline = (np.median(rest_baselines[0]), np.median(rest_baselines[1]))

# Score per chunks
reg_chunks = np.linspace(timestamps_REG[0], timestamps_REG[-1], num=21, endpoint=True)
for k, _ in enumerate(reg_chunks):
    if k == 0:
        continue
    
    # alpha
    score += np.max((a_REG_spline.integral(reg_chunks[k-1], reg_chunks[k]) / (reg_chunks[k] - reg_chunks[k-1])) - rest_baseline[0], 0) * 50 * 10**13
    
    # delta
    score += np.max(rest_baseline[0] - (d_REG_spline.integral(reg_chunks[k-1], reg_chunks[k]) / (reg_chunks[k] - reg_chunks[k-1])), 0) * 50 * 10**13

Наконец, вы можете заметить выше, что размер блока зависит от продолжительности фазы и от num параметр в np.linspace() функция. Вместо этого я бы предпочел определить настройку chunk_size, а затем, если длительность больше, чем chunk_size алгоритмы разрезают сигнал на куски размера chunk_size; возможно, последний кусок будет короче. Поскольку на самом деле это не аспект оптимизации, назовем это бонусом к этой проблеме.

0

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *