Визуализация на Python

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

Библиотек для построения графиков в Python много, однако физику обычно достаточно трёх. Matplotlib даёт статическое изображение, вставляемое затем в статью. Plotly добавляет интерактивность, ценную при разборе сырых данных: наведение курсора с точными числами, масштабирование, много измерений сразу на одном полотне. HoloViews предлагает описывать данные вместо того, чтобы строить график.

Все данные, показанные ниже, вычислены непосредственно в листингах, так что код главы выполняется подряд, от первой строки до последней.

Matplotlib

Matplotlib имеет два интерфейса: первый из них, pyplot, унаследованный от MATLAB, хранит скрытую «текущую фигуру», в которую и добавляют команды наподобие plt.plot. Второй представляет фигуру и оси явными объектами, названными в коде своими именами. Следует сразу использовать второй: в нём видно, к каким осям относится каждый вызов, и код, растущий вместе с задачей, не нарушается, когда фигур становится несколько.

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

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

t = np.linspace(0, 10, 2000)                        # время, мкс
current = 4.2 * np.exp(-((t - 3.5) / 0.45) ** 2)    # импульс тока пучка, мА
current += rng.normal(0, 0.10, t.size)              # шум оцифровки

fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(t, current, lw=0.8)
ax.set_xlabel('Время, мкс')
ax.set_ylabel('Ток пучка, мА')
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig('waveform.png', dpi=300)

subplots возвращает пару, где фигура fig представляет собой весь холст с полями и общим заголовком, а ax — прямоугольник с координатной сеткой, отведённый под сам график. Построение выполняется в осях, сохранение — для фигуры. Вызов tight_layout подбирает поля так, чтобы подписи, поставленные у осей, не срезались краем. savefig записывает файл, а plt.show() открывает окно, блокирующее скрипт до закрытия, и в программе требуется либо одно, либо другое.

Сырой сигнал редко показывают отдельно: рядом помещают сглаженную кривую, а найденный максимум подписывают стрелкой.

from scipy.signal import savgol_filter

smooth = savgol_filter(current, 101, 3)   # окно 101 отсчёт, полином 3-й степени
top = smooth.argmax()

fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(t, current, lw=0.6, alpha=0.4, label='сырой сигнал')
ax.plot(t, smooth, lw=1.6, color='crimson', label='после Савицкого—Голея')
ax.annotate(f'{smooth[top]:.2f} мА при {t[top]:.2f} мкс',
            xy=(t[top], smooth[top]), xytext=(t[top] + 1.0, smooth[top] - 0.5),
            arrowprops={'arrowstyle': '->'})
ax.set_xlabel('Время, мкс')
ax.set_ylabel('Ток пучка, мА')
ax.legend(loc='upper right')
fig.tight_layout()
fig.savefig('waveform-smooth.png', dpi=300)

Две кривые, наложенные в одних осях, различают толщиной, прозрачностью и цветом, а метод legend собирает подписи, переданные параметром label. Метод annotate размещает надпись со стрелкой, протянутой от точки xytext к точке xy на самом графике.

Следующий пример — скан по параметру. Изменяя силу квадрупольной линзы перед экраном, снимают размер пучка, а по параболе, восстановленной подгонкой, вычисляют эмиттанс; сам метод рассмотрен в главе про SciPy. Измерения, снабжённые погрешностями, изображает errorbar.

l = 1.2                                        # расстояние «линза — экран», м
eps, beta0, alpha0 = 0.12e-6, 6.94, -1.69      # эмиттанс, м·рад, и Твисс
gamma0 = (1 + alpha0 ** 2) / beta0

k = np.linspace(-1.5, 1.5, 25)                 # сила линзы, 1/м²
R11, R12 = 1 - k * l, l
sigma = np.sqrt(eps * (R11**2 * beta0 - 2*R11*R12*alpha0 + R12**2*gamma0))
meas = sigma + rng.normal(0, 20e-6, k.size)    # измерение с ошибкой 20 мкм

fig, ax = plt.subplots(figsize=(6, 3.4))
ax.errorbar(k, meas * 1e3, yerr=0.02, fmt='o', ms=4, capsize=3, label='экран')
ax.plot(k, sigma * 1e3, color='crimson', label='модель')
ax.set_xlabel(r'Сила линзы $k$, м$^{-2}$')
ax.set_ylabel(r'Размер пучка $\sigma$, мм')
ax.legend()
fig.tight_layout()
fig.savefig('quadscan.png', dpi=300)

Подписи осей поддерживают подмножество TeX, заключённое в знаки доллара, поэтому индексы и греческие буквы записываются привычным образом. Строку, передаваемую в подпись, необходимо помечать префиксом r, иначе обратный слеш будет обработан Python раньше matplotlib.

Двумерный массив отображает imshow, окрашивающий каждую ячейку по её значению; рассмотрим для него матрицу отклика замкнутой орбиты, вычисленную в главе про коррекцию.

n_bpm, n_cor, nu = 64, 48, 8.42
phi_b = np.sort(rng.uniform(0, 2 * np.pi * nu, n_bpm))   # фазы мониторов
phi_c = np.sort(rng.uniform(0, 2 * np.pi * nu, n_cor))   # фазы корректоров
beta_b = rng.uniform(4.0, 24.0, n_bpm)                   # бета-функция, м
beta_c = rng.uniform(4.0, 24.0, n_cor)

dphi = np.abs(phi_b[:, None] - phi_c[None, :])
R = (np.sqrt(beta_b[:, None] * beta_c[None, :]) / (2 * np.sin(np.pi * nu))
     * np.cos(dphi - np.pi * nu))
print(R.shape)

lim = np.abs(R).max()
fig, ax = plt.subplots(figsize=(5.5, 4.5))
im = ax.imshow(R, cmap='RdBu_r', vmin=-lim, vmax=lim, aspect='auto')
fig.colorbar(im, ax=ax, label=r'$R_{ij}$, мм/мрад')
ax.set_xlabel('Номер корректора $j$')
ax.set_ylabel('Номер монитора $i$')
fig.tight_layout()
fig.savefig('response.png', dpi=300)

(64, 48)

Палитру выбирают по смыслу величины: матрице отклика, где знак меняется от элемента к элементу, соответствует расходящаяся RdBu_r с белым нулём, и границы vmin и vmax, заданные вручную, обязаны быть симметричными. Иначе белый цвет сместится с нуля и изображение будет вводить в заблуждение. Для знакопостоянной величины лучше брать viridis, не меняющий светлоту скачками и потому различимый даже в чёрно-белой печати. Радужный jet, унаследованный от старых пакетов, создаёт ложные границы там, где данные меняются плавно.

Ячейки у imshow равны между собой, а физическая сетка часто неравномерна, и для неё применяют pcolormesh, принимающий сами координаты узлов. Физические границы изображения задаются параметром extent, и номера пикселей на осях сменяются метрами.

Настройки, повторяющиеся от графика к графику, выносятся в rcParams один раз на весь скрипт.

plt.rcParams.update({
    'figure.figsize': (6.5, 3.6),   # ширина колонки статьи
    'font.size': 11,
    'axes.grid': True,
    'grid.alpha': 0.3,
    'savefig.dpi': 300,
    'savefig.bbox': 'tight',
})

fig, (up, down) = plt.subplots(2, 1, sharex=True, figsize=(6.5, 5))
up.plot(t, current, lw=0.6)
up.set_ylabel('Ток, мА')
down.plot(t, smooth, color='crimson')
down.set_ylabel('Сглажено, мА')
down.set_xlabel('Время, мкс')
fig.align_ylabels()
fig.savefig('report.pdf')          # вектор, не рассыпающийся при печати

Параметр sharex=True связывает оси подграфиков, и подписи по времени печатаются только под нижним. Для печати следует сохранять векторный формат, .pdf или .svg, не размываемый при увеличении; растровый .png с плотностью 300 точек на дюйм предназначен для слайдов и веба. Все изображения, собранные в одну статью, полезно строить одним скриптом, поскольку правка стиля тогда обновляет сразу весь набор.

Plotly

При самостоятельном разборе данных требуется подвести курсор и прочитать точное число, выделить рамкой фрагмент и увеличить его, скрыть лишнюю кривую щелчком по легенде. Эти возможности предоставляет Plotly, строящий график в браузере поверх библиотеки plotly.js.

pip install plotly

Рассмотрим огибающую пучка, вычисленную вдоль тракта. Данные даёт уравнение Капчинского—Владимирского, выведенное в отдельной главе, а тракт состоит из пяти соленоидов, расставленных через полметра.

import numpy as np
from scipy.integrate import solve_ivp

gamma = 1 + 2.0 / 0.511                            # пучок 2 МэВ
beta = np.sqrt(1 - gamma ** -2)
perv = 2 * 100.0 / (17045 * (beta * gamma) ** 3)   # первеанс при токе 100 А
eps = 1e-6                                         # эмиттанс, м·рад

z_lens = np.arange(0.5, 3.0, 0.5)                  # пять соленоидов через 0.5 м
z = np.linspace(0, 3, 400)                         # сетка вдоль тракта, м

def envelope(k):
    """Радиус пучка a(z), м, при жёсткостях соленоидов k."""
    def rhs(s, y):
        ks = (k * np.exp(-((s - z_lens) / 0.06) ** 2)).sum()
        a, ap = y
        return [ap, -ks * a + perv / a + eps ** 2 / a ** 3]
    return solve_ivp(rhs, (0, 3), [5e-3, 0.0], t_eval=z, max_step=0.005).y[0]

a = envelope(np.full(5, 25.0)) * 1e3               # огибающая, мм
print(f'{a.min():.2f} {a.max():.2f}')

2.36 7.13

Первое изображение соберём вручную из объектов graph_objects, дающих полный контроль над каждой линией.

import plotly.graph_objects as go

fig = go.Figure()
fig.add_scatter(x=z, y=a, name='верх', line={'color': 'crimson'},
                hovertemplate='z = %{x:.2f} м<br>a = %{y:.2f} мм<extra></extra>')
fig.add_scatter(x=z, y=-a, name='низ', line={'color': 'crimson'},
                fill='tonexty', fillcolor='rgba(220,20,60,0.15)',
                hoverinfo='skip')
for zi in z_lens:
    fig.add_vline(x=zi, line_dash='dot', line_color='gray')
fig.update_layout(xaxis_title='z, м', yaxis_title='Радиус пучка, мм',
                  showlegend=False, height=350)
fig.write_html('envelope.html', include_plotlyjs='cdn')

Параметр hovertemplate задаёт содержимое всплывающей подсказки: %{x} и %{y} подставляют координаты точки, а <extra></extra> убирает служебную рамку с именем кривой. Метод write_html создаёт рядом самодостаточный файл, открываемый в любом браузере без Python. С ключом include_plotlyjs='cdn' он загружает plotly.js из сети, а со значением True встраивает библиотеку внутрь, и файл разрастается до нескольких мегабайт.

Скан по параметру удобно показывать картой, поэтому выполним расчёт для сорока одной жёсткости соленоидов и объединим полученные огибающие в один двумерный массив.

k_grid = np.linspace(15.0, 35.0, 41)
scan = np.array([envelope(np.full(5, k)) for k in k_grid]) * 1e3

heat = go.Figure(go.Heatmap(
    x=z, y=k_grid, z=scan, colorscale='Viridis',
    colorbar={'title': 'a, мм'},
    hovertemplate='z = %{x:.2f} м<br>k = %{y:.1f} 1/м²'
                  '<br>a = %{z:.2f} мм<extra></extra>'))
heat.update_layout(xaxis_title='z, м',
                   yaxis_title='Жёсткость соленоидов k, 1/м²', height=380)
heat.write_html('scan.html', include_plotlyjs='cdn')

При наведении курсора на точку карты читаются три числа: положение вдоль тракта, жёсткость линз и радиус пучка. В matplotlib для того же пришлось бы обращаться к массиву по индексам. Изображение, сохранённое в HTML, остаётся интерактивным и после закрытия Python: масштаб меняется колесом, а прямоугольник, выделенный курсором, растягивается на всё поле, и такими картами удобно искать рабочую точку установки.

Миллион точек plotly отрисовывает уже медленно, поскольку каждая линия, переданная в браузер, существует отдельным объектом JavaScript. Для больших массивов предусмотрен Scattergl, у которого отрисовка переложена на видеокарту.

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

import pandas as pd
import plotly.express as px

rng = np.random.default_rng(1)
trials = rng.uniform(18.0, 32.0, size=(300, 5))     # 300 случайных настроек
loss = np.array([np.sqrt(np.mean((envelope(k) - 4e-3) ** 2)) * 1e3
                 for k in trials])                  # отклонение от 4 мм, мм

runs = pd.DataFrame(trials, columns=[f'k{i}' for i in range(1, 6)])
runs['loss'] = loss

par = px.parallel_coordinates(runs, color='loss',
                              color_continuous_scale='Blues_r')
par.write_html('parallel.html', include_plotlyjs='cdn')
print(f'лучшая настройка: {loss.min():.2f} мм, худшая: {loss.max():.2f} мм')

лучшая настройка: 1.18 мм, худшая: 3.18 мм

Тысяча вариантов настройки оптики линака: каждая линия — один расчёт

Изображение выше построено тем же вызовом, только на настоящей задаче: тысяча прогонов оптимизации оптики линейного ускорителя инжектора ЦКП «СКИФ». Осей восемь, и крайняя левая несёт значение целевой функции, заключённое между 0.461 и 17.752, а остальные отведены под добавки к токам квадруполей, то есть каналы MG-LA1:QLD1-Iadd, MG-LA1:QLF1-Iadd, MG-LA2:QLD2-Iadd и так далее. Цветом закодирована та же целевая функция, и тёмные линии отмечают удачные прогоны.

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

HoloViews

Matplotlib строит то, что ему указано, и потому требует расписывать каждый шаг, от подписи оси до заданного вручную цвета.

HoloViews построен на принципе «аннотировать данные, а не строить график». Пользователь описывает семантику данных, то есть что является аргументом, что значением и в каких они единицах, а библиотека сама выбирает представление. Один и тот же объект, собранный однажды, отрисовывается разными бэкендами: интерактивным bokeh, matplotlib или plotly.

pip install holoviews bokeh

Основой является элемент, контейнер с данными, знающий свой тип визуализации. hv.Curve задаёт непрерывную зависимость, hv.Scatter показывает дискретные измерения, hv.Image отображает величину на регулярной сетке. Измерениям задают удобочитаемые подписи и единицы, попадающие затем на оси и во всплывающие подсказки.

import numpy as np
import holoviews as hv

hv.extension('bokeh')     # можно 'matplotlib' или 'plotly'

rng = np.random.default_rng(42)
t = np.linspace(0, 10, 2000)
current = 4.2 * np.exp(-((t - 3.5) / 0.45) ** 2) + rng.normal(0, 0.10, t.size)

dim_t = hv.Dimension('t', label='Время', unit='мкс')
dim_i = hv.Dimension('I', label='Ток пучка', unit='мА')

wave = hv.Curve((t, current), dim_t, dim_i)
wave

В Jupyter объект, возвращённый последней строкой ячейки, отображается автоматически, так что аналог %matplotlib inline не требуется. С бэкендом bokeh панорамирование, масштабирование и подсказки доступны сразу, и подозрительный выброс, замеченный на графике, рассматривается вблизи без перестроения изображения.

hv.Points внешне похож на Scatter, однако семантика у него иная: обе координаты выступают независимыми измерениями, а точки просто располагаются на плоскости, как частицы на фазовой плоскости. Готовые элементы комбинируются операторами, из которых + размещает графики рядом, а * накладывает их в одних осях.

x = rng.normal(0, 1.0, 2000)        # координата частиц, мм
xp = rng.normal(0, 0.3, 2000)       # угол, мрад
phase = hv.Points((x, xp), ['x', "x'"])

grid = np.linspace(-3, 3, 200)
xx, yy = np.meshgrid(grid, grid)
spot = hv.Image(np.exp(-2 * (xx ** 2 + yy ** 2) / 1.2 ** 2),   # пятно на экране
                bounds=(-3, -3, 3, 3), kdims=['x', 'y'], vdims='I')

sparse = hv.Scatter((t[::200], current[::200]), dim_t, dim_i)
overlay = wave * sparse                       # в одних осях
layout = (wave + phase + spot).cols(3)        # рядом, в три колонки

Совмещение работает благодаря тому, что каждый элемент самодостаточен: HoloViews определяет, что измерения у wave и sparse совпали, и сам выравнивает оси, подписи и легенду. Число колонок в layout задаёт метод .cols(), вызываемый на готовой компоновке.

В matplotlib скан по параметру превращается в цикл, объединяющий десяток подграфиков в одну фигуру, а HoloMap превращает тот же параметр в слайдер.

from scipy.integrate import solve_ivp

# модель тракта (z_lens, perv, eps, z) взята из раздела про Plotly
dim_z = hv.Dimension('z', label='Положение', unit='м')
dim_a = hv.Dimension('a', label='Радиус пучка', unit='мм')

def curve(k):
    """Огибающая пучка при жёсткости всех соленоидов k."""
    def rhs(s, y):
        ks = (k * np.exp(-((s - z_lens) / 0.06) ** 2)).sum()
        a, ap = y
        return [ap, -ks * a + perv / a + eps ** 2 / a ** 3]
    a = solve_ivp(rhs, (0, 3), [5e-3, 0.0], t_eval=z, max_step=0.005).y[0]
    return hv.Curve((z, a * 1e3), dim_z, dim_a)

hmap = hv.HoloMap({k: curve(k) for k in [20, 25, 30, 35]}, kdims='k')
hv.save(hmap, 'envelope-scan.html')     # слайдер останется живым и без Python

dmap = hv.DynamicMap(curve, kdims='k').redim.range(k=(15.0, 35.0))

HoloMap вычисляет все кадры заранее и хранит их в памяти, зато сохраняется в самодостаточный HTML, который можно передать коллеге. DynamicMap вычисляет кадр по движению слайдера, что необходимо, когда параметр непрерывен или каждый расчёт дорог, однако без запущенного ядра Python объект не отображается, и в статический файл он не сохраняется.

Оформление отделено от семантики и задаётся методом .opts() поверх готового объекта.

overlay.opts(
    hv.opts.Curve(width=650, height=280, color='black'),
    hv.opts.Scatter(size=7, color='crimson'),
)

# logz — логарифмическая шкала цвета, нужная пучкам и спектрам
# с большим динамическим диапазоном
spot.opts(width=380, height=340, cmap='inferno', colorbar=True, logz=True)

Часть опций зависит от бэкенда: ширина width и высота height заданы в пикселях у bokeh, а у matplotlib вместо них используются fig_inches и aspect. Полный список опций элемента выводит hv.help(hv.Curve).

Вокруг HoloViews сформировалась экосистема HoloViz. Пакет hvPlot после одного импорта добавляет к DataFrame аксессор .hvplot, заменяющий встроенный df.plot() из pandas; возвращает он обычные объекты HoloViews со всеми слайдерами и оверлеями. Библиотека Panel связывает графики, виджеты и текст в единый интерфейс, превращаемый командой panel serve app.py в веб-приложение. Мониторинг установки становится страницей в браузере без единой строки JavaScript.

Полезные ссылки