NetOnZeroDXC：多变量时间序列网络识别软件

软件X 10（2019）100316原始软件出版物NetOnZeroDXC：通过零延迟互相关从多变量时间序列中识别网络的软件包Alessio Perinellia，Leonardo Riccia，b，a意大利特伦托38123特伦托大学物理系b意大利罗韦雷托大学心智/脑科学中心CIMeC，邮编：38068ar t i cl e i nf o文章历史记录：收到2019年收到修订版，2019年7月8日接受，2019年保留字：互相关网络多元时间序列分析a b st ra ct网络的识别是一个与许多研究领域相关的问题。通过分析与网络节点相关的时间序列来评估这些节点之间的链接是一个至关重要的课题在神经科学和气候研究等不同领域。在这项工作中，我们提出NetOnZeroDXC，一个软件包，实现了最近发布的方法，用于评估多变量时间序列的网络。该方法依赖于时间序列对之间的零延迟互相关的计算该软件以通用的方式实现了算法的各个阶段，从而构成了一种可用于不同领域的网络分析的新工具©2019作者由爱思唯尔公司出版这是CC BY-NC-ND下的开放获取文章许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。代码元数据C1当前代码版本v1.0C2指向此代码版本所用 代码/存储库的永久链接https：//github。Com/ElsevierSoftwareX/SOFTX_2019_131C3 Code Ocean compute capsuleC4法律代码许可证GNU GPL v3使用的C5代码版本控制系统C6使用C++、GNU Scientific Libraries、wxWidgets的软件代码语言、工具和服务C7编译要求，操作环境依赖性Windows 64位。Linux（dependencies：GNU Scientific Libraries，wxWidgets library）C8 If available开发人员文档/手册C9支持电子邮件的问题leonardo。ricci@unitn. It1. 介绍从网络角度分析复杂系统是许多研究领域的共同问题，从神经科学[1]到气候研究[2，3]和经济学[4]。可以通过首先评估可能的网络元件（此后称为节点）之间的链路来识别网络结构。通常通过分析为两个节点候选者中的每一个记录的可观测量的时间序列来确定链接的存在。在最近发表的工作[5]，提出了一种新的方法来识别一组时间序列中的网络。该方法依赖于时间序列对之间的零延迟互相关的计算：相关性的显著性通过以下方式评估：通讯作者：Department of Physics，University of Trento，38123 Trento，Italy。电子邮件地址：leonardo. unitn.it（L. Ricci）。https://doi.org/10.1016/j.softx.2019.100316替代数据生成方法。该算法的主要特点是它不仅提供了网络结构，在节点之间的链接，但它也估计的时间尺度，在链接出现，从而允许研究网络的动态该方法最初是在神经科学的背景下开发的，适用于任何实验情况下，节点的特点是标量时间序列。本软件包实现了上述方法的所有不同阶段，最终从一组输入时间序列中产生连接矩阵。该软件包提供了一个图形用户界面，允许预处理数据集，并以简单的方式调整分析参数。此外，该软件包提供了一个图形界面，允许在颜色尺度连接矩阵的结果的演示文稿。与现有的计算互相关的工具相比，我们的软件包实现了一个完整的分析管道，其中包括生成用于估计2352-7110/©2019作者。由爱思唯尔公司出版。这是一篇开放获取的文章，使用CC BY-NC-ND许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx2A. Perinelli和L.Ricci/SoftwareX 10（2019）100316−−−日本最大+−;; −;; −Fig. 1. 分析方法的流水线在阶段1（输入时间序列的集合）到阶段5（最终连接矩阵）之间。红色图标对应于必须由用户提供的输入数据或参数。算法步骤以蓝色突出显示。所有步骤都在正文中详细描述。(For解释在这个图例中，颜色的引用，读者可以参考本文的网络版本意义这种方法允许处理存在非高斯动态2. 网络的零时延互相关分析在本包中实现的方法在图中以图形方式总结。1.一、算法的输入是对应于N个候选节点的N个时间序列（阶段1）。该方法自动执行平均减法和normalization上分析的时间序列，没有额外的preprocessing是必需的。然而，为了适当地选择算法所覆盖的窗口宽度范围，应该通过常规的时间序列分析技术[6]和非线性技术[7，8]进行频谱分析。在下面解释的主处理流水线之前，并且为了可能地突出潜在的平稳过程，可以计算和可视化整个时间序列的互相关，或者可替代地相关的p值。第一步包括计算N（N1）/2对节点中的每一对的相关图。相关性计算为分析中的两个时间序列的两个窗口之间的零延迟皮尔逊相关系数r（k，w）窗口的宽度w（以样本数表示）从基本宽度w0到Nww0，其中w0和Nw都是用户设置的整数参数;此外，w0必须要公平。最大窗口宽度必须小于时间序列Λ的长度。具有最大长度Nww0的第一窗口跨越具有0和Nww01之间的索引的时间序列的点。两个中间点中最左边的点（以下简称为窗口的中间点）的索引对于较短的窗口长度，该起始值也保持相同。因此，参数k对应于任何运行的中点的位置，窗口，然后以w0的步长逐渐递增，直到仍然容纳时间序列内的最大窗口的值：k=kmin+m·w0，其中0mm<$.（一）NWW0因此，相关图，即特定节点对的所有相关系数r（k，w）的集合，具有Nw行，并且 列数mmax1，取决于时间序列的长度。相关图的示例可见于图1的右侧部分，阶段2。用于计算相关系数的窗口方案的图如图2所示。在某些研究领域中，零延迟互相关的估计会受到噪声的严重影响。在这些情况下，更方便的是将r（k，w）（在两个时间序列之间具有零延迟）评估为r（k，w τ）和r（k，w τ）的平均值，其中r（k，w τ）（r（k，w τ））是延迟τ的第一（第二）时间序列与第二（第一）时间序列之间的互相关。延迟τ必须指定为正整数，因此以采样时间单位表示，必须远小于要研究的时间尺度。MEG分析提供了一个例子，其中源泄漏现象破坏了r（k，w）的直接评估[9]。必须注意的是，延迟的使用将在上面的m_max定义中使用的时间序列长度的值从Λ改变为Λ τ。第二步包括评估每个相关系数的显著性并产生p值图（阶段3），即一组p值p（k，w），每个p值与相应的相关系数r（k，w）相关联.给定k，w，为了评估p值，生成一组M对代理时间序列通过迭代幅度调整傅立叶变换（IAAFT）算法[10，11]从两个原始时间序列中提取。A. Perinelli和L.Ricci/SoftwareX 10（2019）1003163ˆ=.ˆˆˆˆˆˆˆˆ图二. 在节点对i，j的情况下评估相关系数r（k，w）。该算法既保留了值的分布，又（近似地）保留了原始时间序列的自相关性。 有关IAAFT算法的限制和警告的详细讨论可以在参考文献中找到。[11、12]。对于每对替代时间序列，生成相关图。然后，通过对M个替代值的集合内的原始r（k，w）相关系数进行排名，并通过M对排名进行归一化，来计算p值p（k，w）。要生成的替代对的数量M根据p评估的期望分辨率来设置。还可以根据以下表达式，使用基于两个时间序列p=1−F（f，ν1=1，ν2=w− 2），其中F（x，ν1，ν2）是具有ν1，ν2的累积F分布自由度和W1r2（k，w）然而，必须强调的是，Fisher F检验的应用设定了正态性要求，但不一定符合待分析的时间序列通过计算效率η（w）进一步处理p值图中包含的信息。在给定的窗口宽度w处，效率η（w）被评估为宽度w的窗口的分数，其中，显著性阈值< α是由用户设置的参数。该步骤为每对节点i， j产生效率函数η（w）（阶段4）。效率也可以通过使用非重叠窗口来计算，然而，这种操作增加了参数η上的统计噪声。最后，在最后一个流水线步骤中，评估链路的存在性和相应的可观测性时间尺度。如果效率超过给定的效率阈值η0，则认为节点i和节点j之间的链路存在。最小窗口宽度，如果有的话，其中η（W）<$η0确定节点i和节点j之间的链路的可观测性的时间尺度。在对称矩阵W（阶段5）中收集关于链路的存在和相应的可观测性的时间尺度的信息： 是时间尺度，如果有的话，在节点i和j之间的链路是可观察的。对角线W的元素等于零，因为每个节点在任何时间尺度上都与自身相关。时标矩阵W是算法第一部分的最终结果。但是，同一系统可能有多组记录，对应于每个节点的多个时间序列。例如，如果W1，A是建立在通过多通道设备对受试者A进行的第一次脑电图记录上的矩阵，则可以在以后的时间复制该记录，从而产生W2，A，W3，A，. . . . 此外，多个系统具有相同的预期网络结构可能正在调查中。作为示例，同样的EEG分析也可以在对象B，C，. . . .然后，有必要首先合并包含在同一系统内的时间尺度矩阵中的信息，然后可能合并不同系统之间的结果。这两个后处理步骤总结在图3中。在同一系统内，包含在时间尺度的矩阵中的信息被合并以产生矩阵WA，mire，其第i， j个元素通过对单个矩阵的所有第i， j个元素进行排序并挑选第i，j个元素内的R来评估，其中R是要由用户设置的参数。不存在时间尺度的矩阵元素被认为比任何有限元都大如果给出不同的系统，则通过应用相同的排序过程（尽管借助于不同的排序参数R）将具有其WS，merge，information的每个系统合并在单个矩阵中。通过排序合并信息是合理的，因为时间尺度的分布是未知的，因此样本均值不一定是可靠的统计量。3. 封装结构该软件包由两个图形用户界面（GUI）应用程序和两个命令行可执行文件组成其中一个应用程序netOnZe- roDXC_analysis执行从输入时间序列到时间尺度矩阵的整个分析管道，如图1所示。前两个分析步骤，即计算相 关 图 和 p 值 图 ， 也 可 以 通 过 命 令 行 程 序netOnZeroDXC_diagram执行。第三步，即计算效率，由命令行程序netOnZeroDXC_efficiency执行。这两个命令行程序与GUI分开提供，以便于用户在远程计算集群上执行部分分析最后，第二个GUI应用程序netOnZeroDXC_merge，-1F4A. Perinelli和L.Ricci/SoftwareX 10（2019）100316图三. 合并来自不同记录和系统的结果。两个参数Rwithin和Rmerge必须由用户提供实现图3中所示的流水线的步骤，即来自同一系统内和不同系统之间的不同记录的结果的合并。代码是用C++编写的。 除了标准的C++库，GNU科学库（GSL）[13，14]被用来提供随机数生成和快速傅立叶变换例程。图形界面构建在wxWidgets库[15]上。源文件netOnZeroDXC_algorithm.cpp包含算法的核心实现，由上述四个程序共享。代码的其余部分输入/输出功能（*_io.cpp）和图形交互面临 （*_gui_*.cpp）。 的 ''主'' 源 文件 是 标记的*_main.cpp。图1显示了软件包组件及其在分析管道中的作用。四、提供了详细介绍软件功能的用户手册，以及Windows和Linux平台的设置说明。包中包含一组示例3.1. 时间尺度所有程序都以普通ASCII文件的形式提供输入数据，其中列由有效的分隔符（制表符，空格或逗号）分隔，并在应用程序中分配。同样，所有结果都保存为具有相同列分隔符的ASCII文件对输入文件的扩展名不设置任何约束。生成的输出文件扩展名为通过专用输入面板在netOnZeroDXC_analysis应用程序中加载 多 变 量 时 间 序 列 （ 或 等 效 地 将 其 输 入 命 令 行 程 序netOnZeroDXC_diagram也可以加载中间分析结果，而不是原始的多变量时间序列。通过这种方式，可以重复一些分析步骤，例如调整相关参数，同时避免从一开始就重新启动分析过程的必要性出于同样的原因，可以选择保存所有中间步骤，并且可以从这些保存的步骤重新启动当一组文件（每个文件包含一个时间序列）被馈送到应用程序时，每个文件名必须包含一个共享部分和一个单一的运行标签，它是由应用程序通过查看分隔符（下划线或破折号）并丢弃共享部分（例如文件扩展名）来排序的。 有效文件名的示例是“sequence_a. csv”、“sequence_b. csv”等。另一方面，在包含所有时间序列的单个文件的情况下，节点被自动地标记为渐进整数。同样地，如果一组文件（每个文件包含单个p值图或单个效率函数）被馈送到应用程序，则相应的文件名需要包含共享部分和两个运行标签。同样，这些标签是由应用程序以类似的方式整理出来的。有效文件名的示例有“pdag_12_15.dat”、“pdag_12_16.dat”等。运行标签按字典顺序排序。加载数据后，配置面板允许用户设置分析的目标阶段以及所有必需的参数（图中以红色标记）。①的人。“运行"按钮触发计算。有关netOnZeroDXC_analysis应用程序功能的更多详细信息，请参见软件包附带的用户手册。一个主要的选择涉及利用并行计算的可能性，如果可用的话，以加速p值图的评估，这是管道中计算要求最高的步骤。为此，命令行程序netOnZeroDXC_diagram提供，它执行一个单一的p值图的一对时间序列的计算。必须通过命令行标志向程序最后，命令行程序netOnZeroDXC_efficiency根据单个p值图计算效率函数.对于这两个命令行程序，用户手册以及可以通过-h，-help调用的软件版本中描述了必要的标志netOnZeroDXC_analysis应用程序可以通过使用颜色刻度以图形方式显示时间刻度矩阵在预览窗口内，两个控件允许更改显著性阈值α和效率阈值η0。预览矩阵相应地实时更新。A. Perinelli和L.Ricci/SoftwareX 10（2019）1003165图四、 包组件及其在分析管道中的作用。图五. netOnZeroDXC_analysis应用程序的布局。步骤参考示例1，并用数字突出显示：（1）输入面板;（2）选择分析的目标输出;（3）分析参数;（4）有关输出文件的设置;（5）通过并行计算进行分析（启用）;（6）“预 览 "按 钮， 用于 显 示 时标 矩 阵 的预 览。3.2. 矩阵合并合并来自同一系统内不同记录的矩阵的步骤，可能在不同系统之间（见图1）。3）通过使用netOnZeroDXC_merge应用程序进行。在应用程序中，可以适当地设置记录/系统处理好了该应用程序可以加载时间尺度矩阵或效率函数。后一种可能性是在用户通过以下方式执行先前的分析步骤的情况下提供的：命令行程序，这确实产生效率的功能。在大型数据集的情况下，为了避免手动输入数据，可以使用6A. Perinelli和L.Ricci/SoftwareX 10（2019）100316见图6。netOnZeroDXC_merge应用程序的布局必须向应用程序提供适当的配置文件，其中列出了记录和系统以及要加载的文件名。关于配置的负载和示例配置文件的详细信息可以分别在用户手册和examples文件夹中找到。4. 说明性实例在本节中，提供了四个示例来说明该软件包的功能。这些示例以及相应的数据可以在随软件包提供的“ex-amples”文件夹的子文件夹中找到4.1. 一组序列在文件夹“example_1”中包含的示例中，执行针对单个记录的完整分析流水线。分析的数据集由15个时间序列组成，每个时间序列500个点，存储为名为“se-quences.dat”的单个文件中的制表符分隔列。该文件加载在应用程序netOnZeroDXC_通过勾选序列：单个 参数根据图1设置。 五、在标准笔记本电脑上计算需要几分钟最后，显示了时间尺度计算矩阵的图形预览：由于分析是从原始多变量时间序列开始进行的，因此应用程序计算p因此，可以修整显著性阈值α和效率阈值η0两者。图5示出了示例的概要以及两个不同阈值配置的时间尺度的矩阵。请注意，p值图也可以通过使用命令行程序netOnZeroDXC_diagram来计算。对于每一对i，j，调用如下所示netOnZeroDXC_diagram-isequences.dat-osample_pdiag_i_j.dat-n i> j>-W50 -L 4-M500其中从1运行到14，从+1运行到15。4.2. 一组p值图在文件夹“example_2”中包含的示例中这种情况发生在p值图已经可用时.通过勾选“起始数据”控制框至“p值图表：N（N-1）/2个文件”，可在应用程序netOnZeroDXC_analysis参数根据中A. Perinelli和L.Ricci/SoftwareX 10（2019）1003167−==示例文件夹。再次，示出了作为可修整显著性阈值α和效率阈值η0的函数的时间尺度的计算矩阵的图形预览。4.3. 合并多个记录在文件夹“example_3”中包含的示例中，执行多个记录和系统的合并。数据集由包含效率函数的文件组成。因为N10节点被调查，N（N1）/2有45个文件，每个记录。数据集的结构在配置文件“EQUIP_efficiency.cfg”中描述。在这个文件中，对应于一个记录的 每 一 行 都 包 含 系 统 标 签 、 记 录 标 签 、 设 置 为 “e” （ 代 表“efficiencies”）的模式标签和文件名。后一个字符串标识“data”文件夹中的一个文件，该文件夹包含要为给定记录加载的效率文件的列表。 配置文件通过‘‘Configured 输入控件根据以下内容设置：见图6。一旦数据被加载，数据集必须被验证和预处理。合并矩阵的图形预览最终可用，作为表征合并操作的可修整排名参数的函数。由于效率函数是可用的，因此效率阈值η0也可以被修整。 图6显示了示例和合并矩阵的摘要的时间尺度。示例3的文件夹还包含一个配置文件4.4. 从多个记录中合并矩阵在文件夹“example_4”中包含的示例中，执行了来自同一系统的六个该数据集由六个时间尺度矩阵组成，每个矩阵包含24个节点。同样，数据集经过验证和预处理，合并矩阵的图形预览显示。在这种情况下，因为存在单个系统，所以可以修剪的唯一参数是与系统内的记录相关联的排名参数。5. 结论在这项工作中，NetOnZeroDXC包。该软件包实现了一种方法，用于识别网络的多变量时间序列，并被设计为用户和相关的多阶段算法之间的一个易于使用的接口。该软件包提供了两个图形界面应用程序，分别用于分析来自单一来源的多变量时间序列数据和合并来自不同来源的分析。软件包中还包括两个命令行程序，每个程序都执行第一个分析应用程序中实现的关键作为分析应用程序的替代方案，这些命令行程序适用于在计算机集群上运行的分析。该软件包可以应用于所有研究领域，其中复杂系统的每个节点由记录的时间序列描述，并且研究底层网络结构本文提出的方法的一个发展方向是用互信息代替互相关，以研究非线性相关。或者，可以研究部分相关性，以便可能区分直接链接和中介链接。该方法的一个主要改进将涉及非零延迟互相关的评估，以评估不仅网络结构，而且网络内的层次关系竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作引用[1]Buckner RL，Andrews-Hanna JR，Schacter DL.大脑的默认网络：解剖学，功能和与疾病的相关性。Ann New York Acad Sci 2008;1124（ 1）：1-38.http：//dx. 多岛或g/10。1196/annals. 1440. 011。[2]Tsonis AA，Swanson KL，Roebber PJ.网络与气候有什么关系？Bull AmMeteorol Soc 2006;87 （ 5 ） ： 585-95. http://dx.doi.org/10 的 网 站 。1175/BAMS-87-5-585。[3]Gelbrecht M ， Boers N ， Kurths J.风 流 动 的 复 杂 网 络 表 示 。 混 沌 2017;27（3）：035808。 http：//dx. 多岛或g/10。1063/1。四九七七六九九。[4]Mantegna RN.金融市场的等级结构。欧洲物理学杂志1999;11（1）：193-7.http://dx.doi.org/10.1007/s100510050929。[5]放大图片作者：A.不同时间尺度分辨率下大脑网络的相关性。混沌2018;28（6）：063127. http：//dx. 多岛或g/10。1063/1。5025242.[6]按WH，Teukolsky SA，Vetterling WT，Flannery BP。Numerical RecipesinC，3rd ed.剑桥大学出版社; ISBN：978-0-521-88068-8，2007.[7]Franchi M，Ricci L.通过发散率方法计算的最大李雅普诺夫指数的统计特性。Phys Rev E 2014;90（6）：062920。 http：//dx. 多岛或g/10。1103/PhysRevE.九点零。062920.[8] 作者：A.混沌有限维系统产生的时间序列的合适嵌入维数和滞后的识别Phys Rev E 2018;98 （ 5 ） ： 052226 。http ： //dx. 多 岛 或 g/10 。1103/PhysRevE. 九点八。052226.[9]Brookes MJ，Woolrich MW，Price D. In：Supek S，Aine C，editors.磁-脑电图。Berlin：Springer-Verlag; 2014 ，p. 321 -58，Ch. An Introduction toMEG Connectivity Measurements，h t t p：/ / d x. 多岛或g/10。1007/978-3-642-33045-216.[10] Schreiber T，Schmitz A.改进了非线性测试的替代数据物理学评论快报1996;77（4）：635-8. http：//dx. 多岛或g/10。1103/PhysRevLett. 七点七。635.[11]Schreiber T ， Schmitz A. 替 代 时 间 序 列 。 Physica D 2000;142 （ 3 http ：//dx.doi.org/10.1016/S0167-2789（00）00043-9.[12][10] Räth C，Gliozzi M，Papadakis IE，Brinkmann W.重访生成代理时间序列的算法。Phys Rev Lett 2012;109（14）：144101。http：//dx. 多岛或g/10。1103/PhysRevLett. 109. 144101.[13]Galassi M，et al. GNU Scientific Library Reference Manual，3rd ed.. ISBN：0954612078，2009年。[14]GSL网站：https：//www. gnu. 或g/software/gsl/.[15]wxWidgets库网站：https：//www. wxwidgets. 或g/。