A differentiable simulation package for performing inference of synchrotron-radiation-based diagnostics

Watt, R.; O'Shea, B.

doi:10.1107/S1600577524000663

计算机程序

的日志
同步加速器
辐射

国际标准编号：1600-5775

第31卷| 第2部分| 2024年3月| 第409-419页

https://doi.org/10.1107/S16057752400663

打开

访问

用于执行同步辐射诊断推理的可微分模拟包

罗比·瓦特 ^一 ^*和布伦丹·奥谢 ^一

^一SLAC国家加速器实验室，2575 Sand Hill Road，Menlo Park，CA，USA
^*通信电子邮件：rwatt@slac.stanford.edu

编辑：A.Stevenson，澳大利亚同步加速器，澳大利亚(收到日期：2023年11月3日； 2024年1月18日接受；在线2024年2月16日)

粒子加速器的发展方向是不断提高束流质量、电流和重复率。由于光束和诊断信号的相互破坏，这对直接拦截光束的传统诊断提出了挑战。另一种方法是从弯曲磁铁中发射的同步辐射中非侵入性地推断光束参数。然而，从测量的辐射模式推断光束分布是一项复杂且计算昂贵的任务。为了应对我们提出的挑战SYRIPY公司(PYthon中的同步加速器辐射)，一个软件包，用作执行基于同步电子辐射的诊断推断的工具。SYRIPY公司已使用开发PyTorch公司这使得它既可区分又能够利用GPU的高性能，这是执行统计推断的两个重要特征。该软件包由三个模块组成：粒子跟踪器、Lienard–Wiechart解算器和傅里叶光学传播器，可以对同步辐射检测进行从头到脚的模拟。SYRIPY公司已对其进行基准测试SRW公司，该领域中流行的数值包，显示出良好的一致性，速度提高了50倍。最后，我们演示了如何SYRIPY公司可用于使用随机变分推断执行波束参数的贝叶斯推断。

关键词：加速器诊断;可微分模拟;机器学习.

类似文章

1.简介

加速器开发领域不断推动边界，以产生越来越明亮的光束，提高峰值电流并提高重复率。这些高强度光束是揭示不同研究领域新见解的独特工具，从超高梯度等离子体尾流场加速器到探测强场量子电动力学（Yakimenko等。, 2016 ). 例如，先进加速器实验测试设施II（FACET-II）将很快拥有以10GeV的能量输送电荷为2nC的光束的能力，其归一化横向发射度小于10µm，最高可达200kA峰值电流（亚基曼科等。, 2019 ). 然而，这些强束对传统加速器诊断提出了相当大的挑战，尤其是那些需要在束流路径中放置材料的诊断。例如，光束尺寸测量依赖于捕获光束穿过薄箔时发射的光学跃迁辐射。由于表面加热，这些箔片会迅速退化（Stupakov，2013 )导致运营困难，包括成本和时间支出不断增加。此外，这些破坏性诊断对下游光束质量产生了负面影响，从而限制了它们与实验同时进行。因此，人们越来越喜欢利用机器学习来克服这些障碍的单次、无创诊断（Emma等。, 2018 , 2021 ).

当光束沿着加速器传播时，它会穿过弯曲的磁铁并发射同步辐射。这种辐射对于直线加速器和环形加速器来说都是常见的，因此它是一种很有希望的单次激发、非侵入性诊断方法。当进入或离开弯曲磁铁时，光束受到快速波动的边缘场的影响。在这种情况下，产生的辐射强度可以超过标准同步辐射的强度，这种现象称为边缘辐射（Titov&Yarov，1991 ; 丘巴，1995年一 ; 博世，1999 ; 杰洛尼等。, 2009 ). 当光束穿过两个连续弯曲的磁铁时，发射的辐射会产生干扰，从而产生振铃强度分布。这些光环对光束的大小和发散具有很高的灵敏度，这使得边缘辐射成为单发、非侵入性诊断的主要候选者。这一潜力在西伯利亚1号电子储存环和FERMI自由电子激光器的应用中已经得到了探索（Chubar，1995b条 ; 菲奥里托等。, 2014 ).

为了通过统计推断从测量的光强分布中提取光束信息，需要一个系统模型。这无法通过分析实现，因此需要进行数值模拟。为此，存在许多公开可用的软件包，包括光谱（田中和北村，2001年 )和同步辐射车间(SRW公司)（Chubar和Elleaume，1998年 ). 这些软件包被广泛使用，并根据实验结果进行了广泛的基准测试。然而，在应用统计推断时，必须进行大量模拟。能够对边缘辐射建模的现有软件包仅限于基于CPU的操作，这对于该目标来说是次优的。因此，我们开发了一个新的包SYRIPY公司(PYthon中的同步加速器辐射)建立在PyTorch公司库（Paszke等。, 2019 ).SYRIPY公司是专门设计用于执行基于同步电子辐射的诊断推断的工具。通过PyTorch公司，代码在图形处理单元（GPU）上以本机方式运行，使我们能够利用大规模并行体系结构实现高数值效率。此外，SYRIPY公司可以利用PyTorch公司的自动微分包，用于计算输出强度剖面相对于模拟输入的梯度。这种高效率和梯度信息是在高维空间中应用推理方案的宝贵工具。

在本文中，我们将回顾和演示我们的新同步辐射工具包叙利亚人我们将从第2节开始通过讨论该软件包求解的方程组，详细介绍了具体的数值实现。在第3节中我们将提供与分析计算和数值计算相比较的基准结果。最后，我们将演示该软件包的一个应用程序，用它对模拟实验数据进行贝叶斯推理。

2.理论和数值实现

同步辐射产生和检测的从头到尾模拟，从初始电子束参数到探测器的预期强度分布，可以分为三个阶段。首先，通过跟踪穿过感兴趣区域的电子束来获得电子轨迹。其次，利用这些轨迹，计算了初始下游波前的电磁场。最后，场通过光学元件传播到探测器平面。这些阶段如图1所示图中显示了束团压缩器中心边缘辐射的产生和检测示意图。在本节中，我们将详细介绍用于求解此系统的理论和数值实现。该实现在许多模拟参数上高度并行，包括电子、观测点和时间样本。因此，并行化硬件(即GPU）是进行这些计算的理想选择。

图1
图中显示了束流压缩机中央部分边缘辐射的产生和检测。这表明一束300MeV的电子束通过两个0.5T偶极子并发射同步辐射。然后将场传播到探测器的傅里叶平面上（f）=10 cm透镜。

2.1. 粒子跟踪

计算的第一部分包括从束流中采样电子并生成其穿过感兴趣区域的轨迹。如果我们假设束内电子之间的相互作用可以忽略不计，那么它们的运动由洛伦兹运动方程控制，

$[{{\rm{d}}{\bf{p}}\在{{\rm{d}t}}=-e\，c\，{\boldbeta}\次{\bf B}，\eqno（1）]$

哪里第页是电子动量，β是相对论速度，c（c）是光速和e（电子）这个电子电荷。 B类≡B类(第页)表示使用的定义磁场，可以由漂移空间、偶极子和四极子组成。方程式（1）使用四阶龙格-库塔格式求解。由于轨迹的独立性，对束内多个电子进行并行计算是微不足道的。

2.2. 辐射求解器

生成电子轨迹样本后，下一步是计算在下游平面产生的同步辐射(即波前）。对每个采样轨迹单独执行此计算。由单个电子任意运动产生的电磁场由Liénard–Wiechart标量势和矢量势给出（Jackson，1999）; 朗道，2013年 ),

$[\eqalign{\phi（{\bf{r}}，t）&={e}\over{4\pi\epsilon_{0}}\left[{1}\over{（1-{\bv{n}}\cdot{\boldbeta}）r}}\right]_{\rm{ret}}_｛0｝c}}\左[{{{\boldbeta}}\over{（1-{\bf{n}}\cdot{\boldeta}）R}}\right]_{\rm{ret}}，}\eqno（2）]$

其中使用国际单位制，第页= (x个, 年, z（z）)是观察点，t吨是观测时间，R（右）= |第页−第页_{e（电子）}|是电子和观察点之间的距离，第页_{e（电子）}是电子位置，n个是从电子指向观测点的单位矢量(即 n个=R（右）/|R（右）|),ε₀是电常数和[…]_{ret（雷特）}表示括号内的项是在延迟时间计算的

$[t^{\，\prime}=t-{R}\在{c}}.\eqno（3）]$

更常见的电场可以通过以下定义用标量势和矢量势表示，

$[{\bf{E}}=-\nabla\phi-{\partial{\bf{A}}}\over{\paratil{t}}.\eqno（4）]$

大多数诊断只对有限光谱范围内的辐射敏感。因此，与时域相比，在频域中计算电场更有效。将傅里叶变换应用于方程（4）产量

$[{\bf{E}}{\omega}=-\nabla\phi{\omega}-i\omega{\bf{A}}{{\omega}，\eqno（5）]$

哪里

$[\eqalign{\phi_{\omega}&={e}\over{4\pi\epsilon_{0}}\int_{-\infty}^{，\infty}{1}\over{R（t^{\，\prime}）}}\exp\left\{i\omega\left[t^{{4\pi\epsilon上的\cr{\bf{A}}_{\omega}&={{e}\_{0}c}}\int_｛-\infty｝^｛\，\infty｝｛｛｛\boldbeta｝（t^｛\$

和ω是场的角频率。使用方程式（5）和（6），傅里叶域电场由（Chubar，1995）给出一)

${4\pi\epsilon上的[{\bf{E}}_{\omega}={{ie\omeganeneneep_｛0｝c}}\int\limits_{-\infty}^{\infty}{{1}\over{R}}\left[{\boldbeta}-{\bf{n}}\leaft（1+{{ic}\over{\omega{R}{}}\right）\right]\exp\left[i\omega\left（t^{\，\prime}+R/c\right。\等式（7）]$

通过计算方程式（7）在观测点网格上，得到了初始波前。这些观测点中的每一个都是独立的，这再一次使计算易于并行化。

方程式（7）形式为

$[I=\int\limits_{-\infty}^{\infty}{\bf{f}}（t）\exp\big[I\omega\，g（t）\ big]\，{\rm{d}}t，\eqno（8）]$

哪里（f）(t吨)是一个缓慢变化的函数和exp[我ω克(t吨)]快速振荡。这使得求解方程（7）使用标准求积方法在数值上不可行。为了理解原因，我们可以研究图1中的设置，取电子能量和偶极场强分别为～100 MeV和～0.1 T。电子正在靠近c（c），所以需要 $[{\cal O}（10^{-8}\，{\rm{s}}）]$ 遍历设置，设置积分的极限。如果我们在同步加速器的峰值进行计算发射光谱[ω≃ 10¹⁵，使用方程式（24）]积分的振荡部分的周期为～10⁻¹⁵ s.为了避免较大的数值误差，被积函数必须进行密集采样，以解决这些振荡。对于我们的示例，这至少需要10个⁷不实用的样本。

求解方程（7）使用实际计算资源需要专门为高振荡函数设计的求积方法。为了应用这些，我们首先将积分分为三部分，

$[I\，=\int_{-\infty}^{，t_{L}}+\int__{L}{^{$

哪里我_C类在电子通过感兴趣区域的轨迹上积分。感兴趣区域之外(我_{L（左）/R（右）})，我们假设电子以恒定速度向无穷远处移动β_{L（左）/R（右）}.忽略我_{L（左）/R（右）}将导致在t吨_L（左）和t吨_R（右）电子的位置作为时间的函数简单地由下式给出R（右）=c（c）β_{L（左）/R（右）} (t吨−t吨_{L（左）/R（右）}) +R（右）_{L（左）/R（右）}，其中R（右）_{L（左）/R（右）}是积分边界上电子的位置t吨_{L（左）/R（右）}。使用此（f）(t吨)和克(t吨)可以计算为任何顺序，从而允许我们求解我_{L（左）/R（右）}使用渐近展开。这涉及到对等式（8）连续应用分部积分，生成项按顺序递增的序列ω⁻¹第一次订购时（Stein&Murphy，1993 )

$[\eqaligno{I&=\int_{a}^{，b}{\bf{f}}（t）\exp\left[I\omega\，g（t exp\left[I\omega\，g（t）\right]\，{\rm{d}t&（10）\cr&={{1}\over{I\omega}\left\{{{{bf{f}（t）}\over{g^{prime}\在{i\omega}}\int_{a}^{，b}{{{rm{d}}\上在{{rm}d}}\left[{{bf{f}}（t）}\上$

其中，第一项是积分的近似值，第二项是误差。要继续扩展，将对错误应用相同的过程。对于这项工作感兴趣的系统，ω≫1，使扩张迅速收敛。因此，发现一阶展开就足够了。

之间的轨迹t吨_L（左）和t吨_R（右）是用户定义磁场的复杂参数函数。有关的信息（f）(t吨),克(t吨)它们在边界处的导数不足以求解我_C类因此，不可能应用渐近展开式，而是采用菲龙方法（菲龙，1930 ). 该方法与常用的Simpson法则有相似之处，因为二次近似适用于积分的非振荡部分。在应用Filon方法之前，我们首先通过改变变量来去除不规则的非平稳相位x个=克(t吨),

$[I=\int\limits_{g（t_{R}）}^{g（t_{L}）}{｛\bf｝｝\left[g^｛-1}（x）\right]｝\over{g^｛\prime｝\left[g^｛-1}（x）\right]｝｝\exp（I\omega\，x）\，｛\rm｛d｝｝x=\int_｛a｝^｛\，b｝｛\bf｛h｝｝（x）\exp（I\omega\，x）\，｛\rm｛d｝｝x \eqno（11）]$

然后将积分离散为n个间隔和小时在端点和中心点用二次曲线插值(x个₁,x个₂和x个_三)在每个间隔中，即小时(x个) ≃v（v）(x个) =c（c）₁+c（c）₂x个+c（c）_三x个²二次拟合的参数c（c）_我通过求解线性系统得到每个区间内

$[{\bfv}（x{1}$

将此二次近似应用于方程（11）产量

$[\int\limits_{a}^{b}{\bf{h}}（x）\exp（i\omega\，x）\，{\rm{d}}x\simeq\sum_{j\，=\，0}^n-1}\int\limits_{x{2j}}^{x{c}}\左[{\bf{c}{（j）}{1}+{bf{c}^{_{2} x个+{\bf{c}}^{（j）}_{3} x个^{2} \right]\exp（i\omega\，x）\，{\rm{d}}x.\eqno（13）]$

欧拉公式用于用正弦和余弦函数表示复指数，允许使用矩的分析公式来求解每个区间内的积分 $[\textstyle\int\！{x^{\，m}\sin（\omega\，x）}\，{\rm{d}x]$ 和 $[\textstyle\int\！{x^{\，m}\cos（\omega\，x）}\，{\rm{d}}x]$ .对于固定的区间大小，方程（13）的近似误差衰变为 $[{\cal O}（\omega^{-2}）]$ （Stein&Murphy，1993年). 这与用于计算的一阶渐近展开式相同我_{L（左）/R（右）}.

2.2.1. 提高数值效率

方程式（10）和（13）可以使用双精度浮点格式（FP64）轻松求解以获得初始波前场。然而，利用单精度（FP32）是有利的，因为GPU通常针对这种格式进行了优化。基于安培架构的Nvidia GPU（例如，用于进行本工作模拟的Nvidia RTX A6000）的FP32与FP64理论峰值性能比为32:1。我们无法直接使用FP32进行计算，因为灾难性抵消可能会产生数字错误。为了理解原因，我们可以研究方程（11）的主要发射区域当分母克′(t吨) = 1 −n个 · β很小。考虑到电子是高度相对论的，n个 · β≃1当电子向观察者移动时。获取克′(t吨)需要取两个相似数字之间的差值，如果使用FP32，则会导致较大的相对误差。为了避免这种抵消误差，我们可以应用一个小的观测角近似值(即 Z= |z（z）−z（z）_{e（电子）}| ≫X（X）= |x个−x个_{e（电子）}|,Y（Y）= |年−年_{e（电子）}|和β_z（z）≫β_x个, β_年)它允许我们重写相位梯度克′(t吨)和相位克(t吨)作为

$[\eqalign{&g^{\prime}（t）={{\gamma^{-2}+|{\boldbeta}{\perp}|^{2}}\over{2}{-{2Z\，{\bf{R}}{\perp}\cdot{\boldebeta}_{\perp}-\beta{z}|{\bf{R}{{{{\perp2}|^2}}\ over{4Z^{2{+|{{R}}_{\perp}|^{2}}，\cr&g（t）=\int_{-\infty}^{，t}g^{\prime}$

哪里γ是洛伦兹因子R（右）_⊥和β_⊥分别是横向相对位置和速度。用方程（14）定义相位梯度避免使用相似数字之间的差异。

意识到辐射发射在以下区域占主导地位克′(t吨)它很小，可以让我们在数值上做进一步的改进。为了从菲龙方法中获得较小的数值误差，小时(x个)必须在每个间隔内用二次方精确近似。图2(一)显示了小时_x个(t吨)[x个-的组件小时(x个)]对于图1所示的系统当电子通过压缩器时，被积函数经历尖峰和平坦区域。为了用二次曲线近似平坦区域，需要低密度的样品，而尖峰则需要高密度的样品。因此，求解方程（13）使用均匀间隔的时间样本效率很低。

图2
的实部和虚部小时_x个(x个)带有(一)显示针对均匀间隔时间样本绘制的函数，以及(b条)显示绘制的函数u个.

一种更有效的方法将重新分配样本，以便在以下情况下密度更高小时(x个)变化迅速，即何时|[克′(t吨)⁻¹]'|很大。这可以通过逆变换采样实现。这里，目标分布的累积分布函数（CDF）第页(t吨) ∝ |[克′(t吨)⁻¹]‘|已计算，

$[u=C（t）=\textstyle\int\limits_{-\infty}^{\，t}p（t）\，{\rm{d}t，\eqno（15）]$

和它的逆t吨=C类⁻¹(u个)获得。如果我们在u个并根据C类⁻¹(u个)，我们将在t吨具有所需的密度。如图2所示(b条)显示与相同的功能(一)现在密谋反对u个。山峰现在变宽了，这意味着小时(x个)可以使用较少的样本通过二次曲线精确近似。

到目前为止，我们只根据在一个点上观察到的发射重新分配了样品。然而，对于沿x个-尺寸（光束弯曲的横向尺寸）目标分布将发生变化。作为n个和β现在平行于轨迹的不同部分，图2中的峰值(一)将发生变化。更新的时间样本对于这个新的观测点来说不是最佳的。为了解决这个问题|[克_我′(t吨′)⁻¹]'|是在一组M（M）沿线位置x个，新的目标分布由该集合上任意给定值的最大值给出t吨,

$[p（t）\，\propto\，\max_{i\，\in\，M}\left|\left[g_{i}^{prime}（t）^{\，-1}\right]^{prime}\right |。\等式（16）]$

对于单个观测点，重新分配样本需要额外的开销，从而增加了执行计算所需的时间。然而，这一过程仅在有限数量的观测点上进行，远远少于二维波前中的总数(即 $[M\ll N^{\，2}]$ 哪里N个是两者中波前观测点的数量x个和年). 因此，与总仿真时间相比，开销时间很小。最重要的是，正如我们将在下一节中看到的那样，计算具有有限发射度的光束的强度需要对多个电子轨迹求和。鉴于光束较小，发射峰值将出现在所有轨迹的类似位置。因此，这一过程只需要对中心轨迹和用于束流中所有电子的更新时间样本进行。我们发现图3所示的示例计算该过程可以将样本数量减少五倍，而不会显著增加误差。

图3
波前强度分布(一)初始波前和(b条)图1所示设置的探测器平面

.辐射波长为λ=5µm和总样本N个_x个N个_年= 500².

2.2.2. 电子束的发射

当用相机测量辐射剖面时，直接获得的数量是总的光子通量电子束的密度（每单位表面积、每单位相对光谱间隔的光子数）。这可以通过将束内单个电子的电场贡献相加并平方得到，

$[{{\rm{d}}^{2} N个_{\rm{ph}}}\在{\rm}}\Sigma\上，{\rm{d}}\omega/\omega}}={{\epsilon_{0}c}\在{\hbar\pi}}\left|\sum\limits_{i\，=\，1}^{N_{\rm{e}}e_{\omega}（{\bf{r}}_{i}，{\bf{p}}_}i}）\right|^{2}，\eqno（17）]上$

哪里N个_酸碱度是光子数，Σ是表面积，N个_{e（电子）}是光束中的电子数，以及第页_我和第页_我分别是电子的初始位置和动量。这个总和可以分解为时间连贯和非连贯的部分（Hirschmugl等。, 1991 ),

$[\eqaligno{\left|\sum\limits_{i\，=\，1}^{N_{\rm{e}}e_{\omega}}，{\bf{p}}）\，f（{\bf{r}}，}\bf}f（{\bf{r}}，{\bf{p}}）\，{\rm{d}}{\bf{r}}\cr&=I{\rm{CSR}}+I{\rm{ISR}}，&（18）}]$

哪里（f）(第页, 第页)是电子束分布函数。如果光束比辐射波长长，相干项可以忽略，然后通过将单电子强度积分到（f）(第页, 第页). 方程式（18）包含六维积分，这些积分在数值上过于昂贵，无法使用标准求积方法求解。因此，使用蒙特卡罗方法。这里，积分近似为在采样点计算的单电子电场/强度之和{第页_我, 第页_我} ≃（f）(第页, 第页). 图3(一)显示了图1所示系统的强度模式示例，使用此处描述的方法计算。

2.3. 波前传播

计算的最后一部分涉及通过光学系统将波前传播到探测器位置。这是使用标量衍射理论进行的（Goodman，2005 ). 这里进行了近轴近似（小观测角），允许我们将场分量视为独立的，而忽略纵向场。在这种假设下，场在自由空间中的传播由瑞利-索末菲衍射积分给出，

$[E_{\perp}（x_{2}，y_2}）={{z}\over{i\lambda}}\，{\int\！\int}_{\！\Sigma}E_{\perp}}}y{1}，\eqno（19）]$

哪里x个₁, 年₁和x个₂, 年₂分别为传播前后的横坐标，第页₁₂= [(x个₂−x个₁)²+ (年₂−年₁)²+z（z）²]^1/2是两个平面上位置之间的距离，Σ是源平面的面积z（z）是传播距离。方程式（19）是一个卷积积分，可以用卷积定理写成

$[E_{\perp}（x_{2}，y_{2{）={{z}\over{i\lambda}}\，{\cal{F}}^{\，-1}\Big\{{\cal}F}}\left[E_{（x_1}，y_{1}）\right]\，{\ cal{F{}\left[h_{RS}（x_{1{，y_1}）\ right]\Big\}，\eqno（20）]$

哪里小时_RS系列(x个,年)=经验(知识库₁₂)/第页₁₂²和 $[{\cal{F}}]$ 表示傅里叶变换。使用啁啾对该方程进行数值求解z（z）使用Bluestein算法实现的变换（CZT）). CZT类似于更常见的快速傅里叶变换（FFT），但计算效率更高。它允许独立于输入样本设置变换后样本，并内部化波前的零填充（Leutenger等。, 2006 ).

在小观察角度下，计算经验(知识库₁₂)容易出现数值错误。如果是这种情况，我们可以进行菲涅耳近似，从而使用一阶泰勒展开给出第页₁₂≃z（z）+ 0.5(x个₂−x个₁)²+ 0.5(年₂−年₁)²然后将其插入方程式（19）中给予

$[\eqaligno{E_{perp}（x_{2}，y_{2{）={}&{{exp（ikz）}\over{i\lambda{z}}\，int\_{2} -x个_{1} ）^{2}+（y_{2} -年_{1} ）^{2}\右]\右}\，{\rm{d}}x{1}\，}$

这就是菲涅耳衍射解。这也是一个卷积型积分，用与方程（20）相同的方法求解.

为了模拟通过透镜或光圈的传播，只需将场乘以透射函数即可。这些透射函数由下式给出

$[\eqaligno{t_{rm{L}}（x，y）&=\exp\left[-i\，{{k}\over{2f}}}}}\右]&（22）}]$

分别用于透镜和光圈，其中（f）是透镜焦距，Π(x个)是矩形函数w个_一是孔径半径。

为了测量某个入射场的发散度，可以在透镜的傅里叶平面上成像。我们可以通过将场乘以方程（22）中给出的透镜透过率函数来对该系统进行建模和传播z（z）=（f）使用方程式（21）然而，对于低（f）/#镜头这种计算容易产生混叠错误。与透射函数相乘引入了一个二次相位项，它具有增加源带宽的效果。然后需要初始波前中的高密度样本，以避免傅里叶域中的混叠。这增加了传播字段的内存需求。为了避免这种情况，我们可以直接替换E类(x个₁, 年₁) =t吨_L（左）(x个₁, 年₁)E类(x个₁, 年₁)到方程（21）中并重新安排给予

$[\eqaligno{E（x_{2}，y_{2{）={}&{{\exp（ikz）}\over{i\lambda{f}}\exp\left[i\，{{k}\over{2f}}（x_}2}^{2}+y_{2]^{2{}）\right]&（23）\cr&\times\int\！\int_｛\ Sigma｝E（x_｛1｝，y_｛1｝）\exp\left[-i｛2 \pi｝\ over｛\lambda｛f｝｝｝（x_{2} x个_{1} +年_{2} 年_{1} ）\右]\，｛\rm｛d｝｝x_｛1｝\，｛\rm｛d｝｝$

二次项相消。这被称为夫琅和费衍射积分，与方程式（19）不同和（21），它不是卷积型积分，而是简单的缩放傅里叶变换。

图3(b条)图3显示了应用于波前的传播计算示例(一). 场已传播到（f）=10 cm透镜，孔径为2 cm。使用公式（23）中给出的夫琅和费衍射积分进行计算.

3.基准模拟

对第2节中讨论的实现进行基准测试，我们可以将模拟与发射辐射的分析计算进行比较。很少有系统能做到这一点，其中之一是电子在恒定磁场中作圆周运动时的发射。Jackson（1999）给出了完整的计算结果)使用光子通量密度（单位立体角dΣ=R（右）²d日Ω)给定为

$[\eqaligno{{{\rm{d}}^{2} N个_｛\rm｛ph｝｝｝\over｛\rm｛d｝｝\Omega\，｛\rm｛d｝｝\Omega/\Omega｝｝＝｛｝&｛\e ^｛2｝｝\over｛12\pi^{3} c（c）\epsilon_{0}\hbar}}\左（{{omega\rho}\上{c}}\右）^{2}\左^{2}_{2/3}（\psi）+{θ$

哪里ρ是曲率半径，θ是极角，K（K）_2/3和K（K）_1/3是第二类修正贝塞尔函数，其参数为

$[\psi={{omega\rho}\ over{3c}}\ left（{{1}\ over{gamma^{2}}}+theta^{2{right）^{\！3/2}.\eqno（25）]$

方程式（24）的绘图对于固定频率和观察角度，可以在图4中找到(一)和4(b条)分别是。这些显示了光子通量在1T磁场中进行圆周运动的100MeV电子束发出的密度。图中还显示了光子通量使用第2节中讨论的方法计算的密度两种计算结果非常一致。

图4
使用解析（实线红线）和数值（虚线黑线）方法计算的圆周运动光子通量密度分布。面板(一)显示了固定频率下观测角度的分布ω= 10¹⁶ 弧度s⁻¹.面板(b条)用同轴观察点显示频率分布(θ= 0).

对于本工作感兴趣的系统(例如图1)发射的辐射分布不能进行解析计算。因此，为了对这种模拟进行基准测试，我们比较了SYRIPY公司广泛使用的SRW公司包裹。这种比较可以在图5中找到，同时显示通量密度对于单电子，使用与图3相同的模拟参数。再一次SYRIPY公司模拟结果与基准计算结果吻合良好。

图5
使用以下公式计算单个电子的光子通量密度SRW公司（实线红线）和SYRIPY公司（黑色虚线）。面板(一)显示了穿过波前中心的水平线。面板(b条)显示了穿过波前中心的垂直线。

在执行高度并行化的任务时，GPU与CPU相比表现出优越的性能。因此，优化SYRIPY公司对于GPU执行，与SRW公司.性能比较SRW公司和SYRIPY公司如表1所示在这里，我们可以看到执行基准计算所花费的时间由宏电子数归一化。这个SRW公司模拟是使用AMD Ryzen Threadripper Pro 3955WX，16核，32线程CPU执行的，该CPU使用Python的多处理包在多个内核上进行并行化。这个SYRIPY公司在单精度和双精度模式下使用NVIDIA A6000 RTX GPU进行模拟。这使我们能够比较这两种模式的性能；然而，对于这些模拟参数，每像素～0.6%的单精度误差是可以接受的。在双精度模式下，我们发现SRW公司，计算速度增加了～2倍。然而，在单精度模式下，发现了较大的～15×改进。

表1
比较SRW公司和SYRIPY公司模拟运行时间由模拟的宏电子数归一化

波前形状	SRW单处理器	SRW多处理器	FP64（秒）	FP32（秒）
500×500（二维）	1.86	0.086	0.039	0.0057
500×1（1D）	0.014	0.00064	7.0 × 10⁻⁵	1.2 × 10⁻⁵

在上面的SYRIPY公司模拟中，逐个计算了宏观电子。对于较大的2D波前，由于内存限制，无法并行模拟多个宏电子。然而，对于样本数较少的波前，情况并非如此(例如一维模拟N个_年= 1). 通过批量计算单个宏电子波前，减少了开销时间，从而更有效地利用GPU资源。表1还将模拟性能与一维波前（在水平面上，相当于图5）进行了比较(一)]. 这些模拟包括10个⁶宏电子，批量大小为10⁴用于SYRIPY公司在这种情况下，SYRIPY公司与相比，计算速度显著提高SRW公司双精度模式的改进幅度为～9×，单精度模式的提高幅度为～53×。

4.基于梯度的贝叶斯推理

从实验测量推断潜在变量的过程，例如从观察到的强度分布推断光束参数，是一个可以使用统计推断技术解决的反问题的例子。该解决方案需要实现系统的正向模型，通常是模拟。通过调整模拟的输入参数，可以对输出进行操作，直到其与实验观察结果一致。SYRIPY公司主要设计为执行基于同步辐射诊断的统计推断的正向模型。

如果潜在变量的数量很大，求解逆问题可能会非常昂贵。这是由于必须搜索的高维空间，也是维数诅咒的直接后果。直接求解第2节中的方程使用C++/CUDA实现可以消除Python开销，从而提高性能。然而，我们选择使用PyTorch公司而不是CUDA，因为它允许我们使用其自动区分软件包。梯度信息允许我们以智能的方式搜索输入空间，大大减少了所需的资源。这为应用利用高维输入空间的推理方案提供了可能性。

当正向模型不是内射函数时，反问题不成立，即多个不同的输入具有相同的输出（Tarantola，2005 ). 例如，电子束的发散和光斑大小都有助于辐射的展宽。在这种情况下，对潜在参数进行点估计可能会导致错误的结果。贝叶斯推理提供了一种解决不适定逆问题的方法，因为与点估计相反，它推导了输入空间上的概率分布。除此之外，贝叶斯方法为不确定性量化提供了一种稳健的方法。在本节中，我们将演示如何SYRIPY公司可用于从测量的强度分布图中提取光束大小和发散的分布。我们将使用图1所示的相同设置，但我们省略了透镜，使强度分布对光束大小敏感。对于这个简化的例子，我们将通过在波前中心画一条水平线，将测量的强度分布限制为1D。这减少了要推断的参数数量，因为强度分布仅取决于x个-梁尺寸的构件，σ_x个和发散， $[\sigma_{x^{\prime}}]$ 本例中使用的强度分布如图6所示(一)，我们假设σ_x个=300µm， $[\sigma_{x^{\prime}}]$ =150µrad和高斯像素噪声σ_N个=像素计数的2%。

图6
(一)使用合成实验数据，通过贝叶斯推断推断波束参数。(b条)使用最可能的输入光束参数的强度剖面。

我们希望获得所有未知参数的联合概率分布x个= $[（\sigma_{x}，\sigma _{x^{prime}}，\sigma_{N}）]$ 给定测量的强度分布年这是通过应用贝叶斯规则来实现的

$[p（{\bf{x}}|{\bf{y}}）={p$

哪里第页(年|x个)是可能性，第页(x个)是优先分配和第页(年) = $[\textstyle\intp（{\bf{y}}|{\bf{x}}）p（{\ff{x}）\，{\rm{d}}{\bf-x}}]$ 是边际分布。可能性是在给定输入光束参数的情况下获得测量强度分布的概率。像素噪声是高斯的，因此可能性的适当形式是

$[p（{\bf{y}}|{\bf{x}}）=\prod\limits^{无}_{i} ｛｛1｝\over｛\sigma_｛N｝\sqrt｛2\pi｝｝｝\，\exp\left\｛｛｛｛\big[\，f_{i｝（｛\bf｛x｝｝）-y_｛i｝\big]^｛2｝｝\over｛2\sigma_｛N｝^｛2｝｝｝｝\right\｝，\eqno（27）]$

哪里（f）_我(x个)是从正向模型在像素处获得的强度我,年_我是像素处的测量强度我和N个是像素的总数。先验分布编码了我们在观测前对光束参数的先验知识年。我们将仅假设参数的数量级知识，并使用以下无信息的统一先验：第页(σ_x个) = $[｛\cal U｝]$ （100µm、1000µm）， $[p（\sigma_x^{，\prime}）]$ = $[{\cal U}]$ （100µrad，1000µrad）和第页(σ_N个) = $[{\cal U}（0,0.1）]$ .

使用第页(年|x个)和第页(x个)定义，方程式（26）可以解决。然而，直接执行此操作需要计算第页(年)，是所有梁参数的积分。通过我们简化的1D示例，这需要大量计算。然而，对于更现实的高维问题(例如包括年-维数和有限能量扩散），积分很难处理。因此，需要一个近似推理方案。这里，我们将应用随机变分推理（SVI）（Hoffman&Blei，2015）)，使用随机优化器求解方程（26）使用变分推理近似。这包括假设后验值可以通过已知参数形式的分布来近似，即第页(x个|年) ≃q个_θ(x个). 这称为带参数的变分分布θ变分推理的目的是找到θ这样的话q个_θ(x个)最接近真实后方。这是通过最小化两个分布之间的距离来实现的。为此需要一个指标，如Kullback-Leibler（KL）散度，

$[D_{{\rm{KL}}（q||p）=\int\limits_{-\infty}^{，\infty}q_{\boldtheta}}\eqno（28）]$

然而，无法直接计算D类_吉隆坡(q个||第页)因为它需要事先了解第页(x个|年)（我们试图获得的确切分布）。因此，为了继续，我们可以将KL散度改写为

$[\eqalinno｛D_｛｛\rm｛KL｝｝｝｝｝（q||p）&=｛\bf｛E｝｝_｛q_｛\boldtheta｝｝（｛\bf｛x｝）｝\left[\log p（｛\bf｛x｝，｛\bf｛y｝）-\log q_｛\boldtheta｝（｛\bf｛x｝）\right]+\log p（｛\bf｛y｝）｝）\cr&=｛\rm｛ELBO｝｝+\log p（｛\bf｛y｝）｝），&（29）｝]$

哪里 $[{\rm{ELBO}}]$ = $[{\bf{E}}_{q{{\boldtheta}}（{\bf{x}}）}[\log p$ 代表证据下限（霍夫曼和布莱，2015). 方程（29）中的第二项独立于θ; 因此，最大限度地减小相对于θ相当于最小化ELBO，可以直接计算。

通过这个过程，推理方案已从积分问题简化为优化问题，使其在数值上易于处理。在每个更新步骤中_θ需要ELBO，可以通过以下方式高效获得SYRIPY公司是可微的。推理方案的流程图如图7所示为了简化此方案的实施，概率编程库烟火（宾厄姆等。, 2018 )已使用。选择多元高斯分布作为变分分布，以考虑参数之间的相关性。这导致了一个九维优化问题，因为我们正在推断协方差的三个平均值和六个唯一元素。应用SVI的结果如图8所示，显示了所有的2D和1D边缘。正如预期的那样，光束大小和发散度之间存在很强的相关性。从边际分布中，我们可以提取潜在参数的估计值和不确定性σ_x个= 323.6 ± 43.4, $[\sigma_｛x^｛\prime｝｝]$ =145.5±8.1和σ_N个= 0.0213 ± 0.0017. 每个参数的基本真实值在1范围内σ每个估计值的界限，如图8所示带有垂直虚线。最后，使用后向模式进行正向模拟的结果（相当于最大值后部估计）如图6所示(b条). 合成实验数据与最可能的预测结果吻合得很好。

图7
流程图显示了使用可微分模拟的变分贝叶斯推理方案。后验分布由多元高斯近似q个_θ(x个)带参数θ。从该分布中提取样本，并用作模拟的输入。模拟输出和测量的强度分布用于计算ELBO，其中第页(x个, 年) =第页(年|x个)第页(x个). 使用自动微分，ELBO相对于θ已计算。通过最小化ELBOq个_θ和第页(x个|年)也是最小化的。

图8
显示所有2D和1D边缘的后验分布配对图。虚线垂直黑线对应地面真实值。

5.总结

在这份手稿中，我们介绍了一个新颖的Python包SYRIPY公司，专门设计用于促进基于同步电子辐射的诊断的统计推断。SYRIPY公司由三个核心模块组成：粒子跟踪器、Liénard–Wiechart解算器和基于傅里叶光学的传播模块。这使得能够对同步辐射的产生和检测进行从头到脚的模拟。该软件包是使用库开发的PyTorch公司，允许SYRIPY公司在CPU和GPU上本机运行。特别是，Liénard–Wiechart解算器和Fourier光学模块具有高度并行性，使得代码在GPU上运行时非常高效。使用开发包PyTorch公司因为底层库支持自动计算渐变。

我们给出了基准计算结果，表明该软件包与分析和数值结果都吻合良好。对于只需要单个浮点精度的模拟，SYRIPY公司与相比，速度显著提高（～50倍）SRW公司这是GPU与CPU相比更高的指令吞吐量的直接结果。

叙利亚人它既快速又可微分，是执行统计推断的理想工具。我们已经通过使用该软件包对模拟实验数据进行贝叶斯推断来证明了这种能力。随着SVI的应用，求解贝叶斯方程的复杂任务简化为更易于管理的优化问题。即使如此，通过我们简化的1D示例，必须优化九个参数。由于维度诅咒，如果没有梯度信息，这将是很难做到的。此外SYRIPY公司不限于贝叶斯推理。该软件包的其他应用可以包括通过最大似然如Roussel所述，使用神经网络参数化估计或预测光束的完整横向相空间等。(2023 ).

叙利亚人在GitHub上公开提供，网址为https://github.com/robbiewatt1/SYRIPY网站.

资金筹措信息

根据合同DE-AC02-76SF00515，这项工作得到了SLAC国家加速器实验室能源部的支持。

工具书类

Bingham，E.，Chen，J.P.，Jankowiak，M.，Obermeyer，F.，Pradhan，N.，Karaletsos，T.，Singh，R.，Szerlip，P.，Horsfall，P.&Goodman，N.D.（2018）。J.马赫。学习。物件。 20, 1–6. 谷歌学者
 Bluestein，L.（1970）。IEEE传输。音频电声。 18, 451–455. 交叉参考谷歌学者
 Bosch，R.（1999）。编号。仪器。方法物理学。决议A,431, 320–333. 科学网交叉参考中国科学院谷歌学者
 丘巴尔，O.V.（1995年一).科学评论。仪器。 66, 1872–1874. 交叉参考中国科学院谷歌学者
 Chubar，O.（1995年b条).1995年粒子加速器会议记录（PAC1995）1995年5月1日至5日，第4卷，美国德克萨斯州达拉斯，第2402–2404页。谷歌学者
 Chubar，O.和Elleaume，P.（1998年）。第六届欧洲粒子加速器会议记录（EPAC98）1998年6月22日至26日，瑞典斯德哥尔摩，第1177至1179页谷歌学者
 Emma，C.、Edelen，A.、Hanuka，A.、O'Shea，B.和Scheinker，A.（2021年）。问询处,12, 61. 交叉参考谷歌学者
 Emma，C.、Edelen，A.、Hogan，M.J.、O'Shea，B.、White，G.和Yakimenko，V.（2018年）。物理学。修订加速度。横梁,21, 112802. 交叉参考谷歌学者
 费隆，L.N.G.（1930）。程序。R.Soc.爱丁堡。 49, 38–47. 交叉参考谷歌学者
 Fiorito，R.、Shkvarunets，A.、Castronovo，D.、Cornacchia，M.、Di Mitri，S.、Kishek，R.，Tschalaer，C.和Veronese，M.（2014）。物理学。修订版ST加速度。横梁,17, 122803. 交叉参考谷歌学者
 Geloni，G.、Kocharyan，V.、Saldin，E.、Schneidmiller，E.和Yurkov，M.（2009年）。编号。仪器。方法物理学。决议A,605, 409–429. 科学网交叉参考中国科学院谷歌学者
 Goodman，J.W.（2005）。傅立叶光学导论。罗伯茨公司。谷歌学者
 Hirschmugl，C.J.，Sagurton，M.&Williams，G.P.（1991）。物理学。版次A,44, 1316–1320. 交叉参考中国科学院公共医学科学网谷歌学者
 Hoffman，M.D.和Blei，D.M.（2015）。2015年第18届国际人工智能与统计会议（AISTATS）会议记录2015年5月9日至12日，美国加利福尼亚州圣地亚哥，第361-369页谷歌学者
 Jackson，J.D.（1999）。经典电动力学第三版，纽约：Wiley。谷歌学者
 Landau，L.D.（2013）。经典场论，第2卷。爱思唯尔。谷歌学者
 Leutenger，M.、Rao，R.、Leitgeb，R.A.和Lasser，T.（2006年）。选择。快递,14, 11277–11291. 交叉参考公共医学谷歌学者
 Paszke，A.、Gross，S.、Massa，F.、Lerer，A.、Bradbury，J.、Chanan，G.、Killeen，T.、Lin，Z.、Gimelshein，N.、Antiga，L.、Desmaison，A.、Köpf，A.，Yang，E.、DeVito，Z.，Raison，M.、Tejani，A.、Chilamkurthy，S.，Steiner，B.、Fang，L.，Bai，J.和Chintala，S.（2019年）。arXiv公司:1912.01703. 谷歌学者
 Roussel，R.，Edelen，A.，Mayes，C.，Ratner，D.，Gonzalez-Aguilera，J.P.，Kim，S.，Wisniewski，E.&Power，J.（2023年）。物理学。修订稿。 130, 145001. 交叉参考公共医学谷歌学者
 Stein，E.M.和Murphy，T.S.（1993年）。谐波分析：实变量方法、正交性和振荡积分，第3卷。普林斯顿大学出版社。谷歌学者
 Stupakov，G.（2013）。通过入射相对论电子束熔化薄箔。技术报告SLAC-PUB-15729。SLAC国家加速器实验室，美国加利福尼亚州门洛帕克谷歌学者
 Tanaka，T.和Kitamura，H.（2001）。J.同步辐射。 8, 1221–1228. 科学网交叉参考中国科学院 IUCr日志谷歌学者
 Tarantola，A.（2005）。模型参数估计的反问题理论和方法。暹罗。谷歌学者
 Titov，A.和Yarov，A.（1991）。编号。仪器。方法物理学。决议A,308, 117–119. 交叉参考谷歌学者
 Yakimenko，V.，Alsberg，L.，Bong，E.，Bouchard，G.，Clarke，C.，Emma，C.，Green，S.，Hast，C.，Hogan，M.J.，Seabury，J.，Lipkowitz，N.，O'Shea，B.，Storey，D.，White，G.&Yocky，G.（2019年）。物理学。修订加速度。横梁,22, 101301. 交叉参考谷歌学者
 Yakimenko，V.、Cai，Y.、Clarke，C.、Green，S.、Hast，C.、Hogan，M.、Lipkowitz，N.、Phinney，N.，White，G.和Yocky，G.（2016）。第七届国际粒子加速器会议记录（IPAC’16），2016年5月8日至13日，韩国釜山，第1067至1070页。图OBB02谷歌学者