У меня есть программа, которая чередует 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
; возможно, последний кусок будет короче. Поскольку на самом деле это не аспект оптимизации, назовем это бонусом к этой проблеме.