Floyd算法:求多源解负权图最短路径问题

思想:动态规划

Floyd-Warshall算法(Floyd-Warshall algorithm)是解决任意两点间的最短路径的一种算法,可以正确处理有向图或负权的最短路径问题,同时也被用于计算有向图的传递闭包。

Floyd-Warshall算法的时间复杂度为O(N^3),空间复杂度为O(N^2)。

算法描述

1) 算法思想原理:

Floyd算法是一个经典的动态规划算法。用通俗的语言来描述的话,首先我们的目标是寻找从点i到点j的最短路径。

从动态规划的角度看问题,我们需要为这个目标重新做一个诠释(这个诠释正是动态规划最富创造力的精华所在)

从任意节点i到任意节点j的最短路径不外乎2种可能,1是直接从i到j,2是从i经过若干个节点k到j。所以,我们假设Dis(i,j)为节点u到节点v的最短路径的距离,对于每一个节点k,我们检查Dis(i,k) + Dis(k,j) < Dis(i,j)是否成立,如果成立,证明从i到k再到j的路径比i直接到j的路径短,我们便设置Dis(i,j) = Dis(i,k) + Dis(k,j),这样一来,当我们遍历完所有节点k,Dis(i,j)中记录的便是i到j的最短路径的距离。

2) 算法描述:

a. 从任意一条单边路径开始。所有两点之间的距离是边的权,如果两点之间没有边相连,则权为无穷大。   

b. 对于每一对顶点 u 和 v,看看是否存在一个顶点 w 使得从 u 到 w 再到 v 比己知的路径更短。如果是更新它。

Floyd算法过程矩阵的计算—-十字交叉法

方法:两条线,从左上角开始计算一直到右下角 如下所示

先来简单分析下,由于矩阵中对角线上的元素始终为0,因此以k为中间点时,从上一个矩阵到下一个矩阵变化时,矩阵的第k行,第k列和对角线上的元素是不发生改变的(对角线上都是0,因为一个顶点到自己的距离就是0,一直不变;而当k为中间点时,k到其他顶点(第k行)和其他顶点到k(第k列)的距离是不变的

因此每一步中我们只需要判断4*4-3*4+2=6个元素是否发生改变即可,也就是要判断既不在第k行第k列又不在对角线上的元素。具体计算步骤如下:以k为中间点(1)“三条线”:划去第k行,第k列,对角线(2)“十字交叉法”:对于任一个不在三条线上的元素x,均可与另外在k行k列上的3个元素构成一个2阶矩阵x是否发生改变与2阶矩阵中不包含x的那条对角线上2个元素的和有关,若二者之和小于x,则用它们的和替换x,对应的Path矩阵中的与x相对应的位置用k来替代。

给出矩阵,其中矩阵A是邻接矩阵,而矩阵Path记录u,v两点之间最短路径所必须经过的点(指的是u->….->该点->v)

分别以每一个点作为K点转接v

image
image
image
image

代码实现

java 版本

核心代码只有 4 行,实在令人赞叹。  [java]

private void floyd() {
    for (int k = 0; k < n; k++) {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                a[i][j] = Math.min(a[i][j], a[i][k] + a[k][j]);
            }
        }
    }

    // 打印
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            System.out.println(i + " " + j + ":" + a[i][j]);
        }
    }
}

正确性分析

核心代码只有四行,三行循环,一行更新,的确十分简洁优雅,可是这四行代码为什么就能求出任意两点的最短路径呢?

看代码的特点,很显然特别像是一种动态规划,确实,之前也说过该算法是基于动态规划的。

我们从动态规划的角度考虑,解动态规划题目的重点就是合理的定义状态,划分阶段,我们定义 f[k][i][j] 为经过前k的节点,从i到j所能得到的最短路径,f[k][i][j]可以从f[k-1][i][j]转移过来,即不经过第k个节点,也可以从f[k-1][i][k]+f[k-1][k][j]转移过来,即经过第k个节点。

观察一下这个状态的定义,满足不满足最优子结构和无后效性原则。

  • 最优子结构

图结构中一个显而易见的定理:

最短路径的子路径仍然是最短路径, 这个定理显而易见,比如一条从a到e的最短路a->b->c->d->e 那么 a->b->c 一定是a到c的最短路c->d->e一定是c到e的最短路,反过来,如果说一条最短路必须要经过点k,那么i->k的最短路加上k->j的最短路一定是i->j 经过k的最短路,因此,最优子结构可以保证。

  • 无后效性

状态f[k][i][j]由f[k-1][i][j],f[k-1][i,k] 以及f[k-1][k][j]转移过来,很容易可以看出,f[k] 的状态完全由f[k-1]转移过来,只要我们把k放倒最外层循环中,就能保证无后效性。

满足以上两个要求,我们即证明了Floyd算法是正确的。

我们最后得出一个状态转移方程:f[k][i][j] = min(f[k-1][i][k],f[k-1][i][k],f[k-1][k][j]) ,可以观察到,这个三维的状态转移方程可以使用滚动数组进行优化。

K 为什么要放在最外层

采用动态规划思想,f[k][i][j] 表示 i 和 j 之间可以通过编号为 1…k 的节点的最短路径。

f[0][i][j] 初值为原图的邻接矩阵。

f[k][i][j]则可以从f[k-1][i][j]转移来,表示 i 到 j 不经过 k 这个节点。

也可以 f[k-1][i][k]+f[k-1][k][j] 从转移过来,表示经过 k 这个点。

意思即 f[k][i][j]=min(f[k-1][i][j], f[k-1][i][k]+f[k-1][k][j])

然后你就会发现 f 最外层一维空间可以省略,因为 f[k] 只与 f[k-1] 有关。

虽然这个算法非常简单,但也需要找点时间理解这个算法,就不会再有这种问题啦。

Floyd算法的本质是DP,而k是DP的阶段,因此要写最外面。

想象一个图,讨论的是要从1点到达3点,是直接走还是经过中间点2,从而确定两点之间的最短路径

滚动数组优化

滚动数组是一种动态规划中常见的降维优化的方式,应用广泛(背包dp等),可以极大的减少空间复杂度。

在我们求出的状态转移方程中,我们在更新f[k]层状态的时候,用到f[k-1]层的值,f[k-2] f[k-3]层的值就直接废弃了。

所以我们大可让第一层的大小从n变成2,再进一步,我们在f[k]层更新f[k][i][j]的时候,f[k-1][i][k] 和 f[k-1][k][j] 我们如果能保证,在更新k层另外一组路径m->n的时候,不受前面更新过的f[k][i][j]的影响,就可以把第一维度去掉了。

我们现在去掉第一个维度,写成我们在代码中的那样,就是f[i][j] 依赖 f[i][k] + f[k][j] 我们在更新f[m][n]的时候,用到了f[m][k] + f[k][n] 假设f[i][j]的更新影响到了f[m][k] 或者 f[k][m] 即 {m=i,k=j} 或者 {k=i,n=j} 这时候有两种情况,j和k是同一个点,或者i和k是同一个点,那么 f[i][j] = f[i][j] + f[j][j],或者f[i][j] = f[i][i]+f[i][j] 这时候,我们所谓的“前面更新的点对”还是这两个点本来的路径,也就是说,唯一两种在某一层先更新的点,影响到后更新的点的情况,是完全合理的,所以说,我们即时把第一维去掉,也满足无后效性原则。

因此可以用滚动数组优化。

优化之后的状态转移方程即为:f[i][j] = min(f[i][j],f[i][k]+f[k][j])

求具体路径:

我们上面仅仅是知道了最短路径的长度,实际应用中我们常常是需要知道具体的路径,在Floyd算法中怎么求具体路径呢,很简单,我们只需要记录下来在更新f[i][j]的时候,用到的中间节点是哪个就可以了。

假设我们用path[i][j]记录从i到j松弛的节点k,那么从i到j,肯定是先从i到k,然后再从k到j, 那么我们在找出path[i][k] , path[k][j]即可。

即 i到k的最短路是 i -> path[i][k] -> k -> path[k][j] -> k q` 然后求path[i][k]和path[k][j] ,一直到某两个节点没有中间节点为止,代码如下:

在更新路径的时候:  [java]

if(a[i][k]>temp){
    a[i][j] = temp;
    path[i][j] = k;
}

求路径的时候:  [java]

public String getPath(int[][] path, int i, int j) {
    if (path[i][j] == -1) {
        return " "+i+" "+j;
    } else {
        int k = path[i][j];
        return getPath(path, i, k)+" "+getPath(path, k, j)+" ";
    }
}

总结

Floyd算法只能在不存在负权环的情况下使用,因为其并不能判断负权环,上面也说过,如果有负权环,那么最短路将无意义,因为我们可以不断走负权环,这样最短路径值便成为了负无穷。

但是其可以处理带负权边但是无负权环的情况。

以上就是求多源最短路的Floyd算法,基于动态规划,十分优雅。

但是其复杂度确实是不敢恭维,因为要维护一个路径矩阵,因此空间复杂度达到了O(n^2),时间复杂度达到了O(n^3),只有在数据规模很小的时候,才适合使用这个算法,但是在实际的应用中,求单源最短路是最常见的,比如在路由算法中,某个节点仅需要知道从这个节点出发到达其他节点的最短路径,而不需要知道其他点的情况,因此在路由算法中最适合的应该是单源最短路算法,Dijkstra算法。

Hierholzer 算法 —有向图求解欧拉路径(环路)

  • 欧拉路径
    如果在一张图中,可以从一点出发遍历所有的边,那么遍历过程中的这条路径就叫做欧拉路径。如果这条路径是闭合的,那就称为欧拉回路
    简单地说,如果一张图可以被“一笔画”,那么“一笔画”的那个轨迹就叫做欧拉路径。
  • 欧拉图判定定理
    含有欧拉回路的图称为欧拉图,含有欧拉路径的图称为半欧拉图。
    无向图中,如果所有顶点的度数都为偶数,则为欧拉图;如果有两个顶点的度数为奇数,其他的为偶数,则为半欧拉图。
    有向图中,如果所有顶点的入度等于出度,那么就是欧拉图。判定半欧拉图的方法也很简单,大家可以自行推理。

Hierholzer 算法

问题简述:现给出一个有向图,且为欧拉图。求欧拉回路。

Hierholzer 算法过程

  1. 选择任一顶点为起点,遍历所有相邻边。
  2. 深度搜索,访问相邻顶点。将经过的边都删除。
  3. 如果当前顶点没有相邻边,则将顶点入栈。
  4. 栈中的顶点倒序输出,就是从起点出发的欧拉回路。

从一个可能的起点出发,进行深度优先搜索,但是每次沿着辅助边从某个顶点移动到另外一个顶点的时候,都需要删除这个辅助边。如果没有可移动的路径,则将所在结点加入到栈中,并返回。

dfs(node, trace){
	while(!node.adj.isEmpty()){
		Node next = node.adj.removeLast();
		dfs(next, trace);
	}
	trace.addLast(node);
}

最后得到的栈中保存的就是整个欧拉闭迹中的顶点。(要恢复我们需要不断出栈,因此如果你用列表来存欧垃迹的话需要反转一次)。

重新安排行程:给你一份航线列表 tickets ,其中 tickets[i] = [fromi, toi] 表示飞机出发和降落的机场地点。请你对该行程进行重新规划排序。所有这些机票都属于一个从 JFK(肯尼迪国际机场)出发的先生,所以该行程必须从 JFK 开始。如果存在多种有效的行程,请你按字典排序返回最小的行程组合。例如,行程 [“JFK”, “LGA”] 与 [“JFK”, “LGB”] 相比就更小,排序更靠前。
假定所有机票至少存在一种合理的行程。且所有的机票 必须都用一次且只能用一次。

在这里插入图片描述
class Solution:
    ''' Hierholzer 算法
    '''
    def findItinerary(self, tickets):
        def dfs(cur, graph, res):
            while graph[cur]: # 遍历所有边
                dfs(graph[cur].pop(), graph, res) # 访问并删除
            if not graph[cur]: # 将顶点入栈
                res.append(cur)
        import collections
        # 建图
        graph = collections.defaultdict(list)
        for start, end in sorted(tickets)[::-1]:
            graph[start].append(end)
        res = []
        dfs("JFK", graph, res)
        return res[::-1] # 倒序输出

AlphaFold2论文

https://www.nature.com/articles/s41586-021-03819-2

https://github.com/deepmind/alphafold

沐神AlphaFold讲解

论文补充材料:https://www.biorxiv.org/content/10.1101/2021.10.04.463034v1

摘要

蛋白质对生命至关重要,了解它们的结构可以促进对其功能的系统理解。通过大量的实验,已经确定了大约 100,000 种独特的蛋白质结构,但这仅代表了数十亿已知蛋白质序列中的一小部分。

蛋白质结构是指蛋白质分子的空间结构。作为一类重要的生物大分子,蛋白质主要由化学元素组成。所有蛋白质都是由20种不同的L型α氨基酸连接形成的多聚体,在形成蛋白质后,这些氨基酸又被称为残基。
蛋白质一级结构:组成蛋白质多肽链的线性氨基酸序列。
蛋白质二级结构:依靠不同氨基酸之间的C=O和N-H基团间的氢键形成的稳定结构,主要为α螺旋和β折叠。
蛋白质三级结构:通过多个二级结构元素在三维空间的排列所形成的一个蛋白质分子的三维结构。
蛋白质四级结构:用于描述由不同多肽链(亚基)间相互作用形成具有功能的蛋白质复合物分子。

仅根据其氨基酸序列预测蛋白质的三级结构,这是“蛋白质折叠问题”, 50 多年来一直是一个重要的开放性研究问题。尽管最近取得了进展,但现有方法仍远未达到原子级准确度,尤其是当没有可用的同源结构时。

同源结构,那些不同物种因来自共同祖先而具有的相似性结构。 例如现代马经较长时间的修饰成为具有一个趾,鼹鼠及其他洞穴动物成为瘤状肢体,大象的肢体成为柱状,这样功能各异的前肢有一个共同来源,他们都来自原始陆生脊椎动物五趾型的肢体。

在这里,我们提供了第一种计算方法,即使在不知道相似结构的情况下,它也可以以原子精度定期预测蛋白质结构。我们在具有挑战性的第 14 次蛋白质结构预测关键评估 (CASP14) 中验证了我们基于神经网络模型的完全重新设计的AlphaFold,在大多数情况下表现出与实验相媲美的准确性,并且大大优于其他方法。支持最新版本的 AlphaFold 是一种新颖的机器学习方法,它将关于蛋白质结构的物理和生物学知识,利用多序列比对,融入深度学习算法的设计中。

序列比对指将两个或多个序列排列在一起,标明其相似之处。序列中可以插入间隔(通常用短横线“-”表示)。对应的相同或相似的符号(在核酸中是A, T(或U), C, G,在蛋白质中是氨基酸残基的单字母表示)排列在同一列上。
tcctctgcctctgccatcat—caaccccaaagt
|||| ||| ||||| ||||| ||||||||||||
tcctgtgcatctgcaatcatgggcaaccccaaagt
多序列比对是成对比对的延伸,是为了在一次比对里面处理多于两条的的序列。多序列比对方法试图比对一个指定序列集合里面的所有序列,这可以帮助确定这些序列的共同区段。进行多序列比对有几种方法,最常用的一种是Clustal程序集,它使用渐进多序列比对算法。Clustal在cladistics中被用来建立进化树,在PSI-BLAST和Hidden Markov model (HMM)中用来建立序列档案以在序列数据库中搜索更远的同源序列。

从蛋白质序列预测蛋白质3D结构的计算方法的发展沿着两条互补的路径前进,分别关注物理相互作用或进化历史。 物理相互作用方案将我们认知的分子驱动力(molecular driving forces)整合到物理热力学或动力学模拟或统计模型逼近中。 虽然理论上非常吸引人,但由于分子模拟的计算难度、蛋白质稳定性的上下文依赖性以及难以产生足够准确的蛋白质物理学模型,这种方法已被证明对即使是中等大小的蛋白质也极具挑战性。 近年来,进化方案提供了一种替代方案,其中蛋白质结构约束来自蛋白质进化历史的生物信息学分析、与已解决结构的同源性和成对进化相关性。

这种生物信息学方法极大地受益于蛋白质数据库 (PDB) 中存储的实验蛋白质结构的稳定增长、基因组测序的爆炸式增长以及用于解释这些相关性的深度学习技术的快速发展。 尽管取得了这些进展,基于物理和进化历史的方法产生的预测远低于实验准确性。

我们(DeepMind团队)开发了第一个能够在大多数情况下预测蛋白质结构接近实验准确性的计算方法。 我们开发的神经网络 AlphaFold 已进入 CASP14 评估(2020 年 5 月至 7 月)。 CASP 评估每两年进行一次,使用未在 PDB 中存放或公开披露的最近解决的结构,因此它是对参与方法的盲测,长期以来一直作为结构预测准确性的金标准评估。

蛋白质结构预测 (CASP) 实验的批判性评估旨在建立蛋白质结构预测的当前技术水平,确定已取得的进展,并突出未来可能最有成效的工作重点。

… 省略了部分内容

在图 2a 中证明,CASP14 中展示AlphaFold的高准确率扩展到最近大量PDB结构样本,其中所有结构在我们的训练数据截止后都存储在PDB中,并作为完整链进行分析。此外,当主链预测准确时,观察到高侧链准确性(图 2b),并且表明我们预测的局部距离差异测试(pLDDT)置信度可靠地预测了 Ca 局部距离差异测试(lDDT-Cα)准确度相应的预测(图 2c)。我们还发现可以准确估计全局叠加度量模板建模分数 (TM-score)(图 2d)。总体而言,这些验证了 AlphaFold 在 CASP14 蛋白上的高精度和可靠性也可以迁移到最近 PDB 提交的未经整理的数据集中。

Fig. 1 | AlphaFold produces highly accurate structures.

PDB蛋白质结构数据库(Protein Data Bank,简称PDB)是美国Brookhaven国家实验室于1971年创建的,由结构生物信息学研究合作组织(Research Collaboratory for Structural Bioinformatics,简称RCSB)维护。和核酸序列数据库一样,可以通过网络直接向PDB数据库提交数据。

AlphaFold网络

AlphaFold 通过结合基于蛋白质结构的进化、物理和几何约束的新型神经网络架构和训练程序,大大提高了结构预测的准确性。特别是,我们展示了一种联合嵌入多序列比对 (MSA) 和成对特征的新架构、一种新的输出表示和相关损失函数,可实现准确的端到端结构预测、新的等变注意架构、使用中间损失函数来实现预测的迭代改进,屏蔽 MSA 损失与结构联合训练,使用自蒸馏从未标记的蛋白质序列中学习,以及自我估计准确率。

蒸馏,就是知识蒸馏,将教师网络(teacher network)的知识迁移到学生网络(student network)上,使得学生网络的性能表现如教师网络一般;或者大型模型迁移到小型模型中,小型模型的参数规模小,运行速度快,但性能与大型模型参不多。

AlphaFold 网络使用一级氨基酸序列和同源物的比对序列作为输入,直接预测给定蛋白质的所有重原子的 3-D 坐标(图 1e,请参阅方法了解输入的详细信息,包括数据库、MSA 构建和使用的模板)。 最重要的方法和组件的描述如下在附件中提供了完整的网络架构和训练过程在附件的方法章节中。

该网络包括两个主要阶段。 首先,网络的主干通过我们称为 Evoformer 的新型神经网络块的重复层处理输入,以生成 Nseq × Nres 数组(Nseq:序列数,Nres:残基数),表示已处理的 MSA 和 Nres × Nres 数组,表示残基对。 MSA 表示是用原始 MSA 初始化的,但请参阅 附件-方法 1.2.7 了解处理非常深的 MSA 的详细信息。 Evoformer 块包含许多新颖的基于注意力和非基于注意力的组件。 我们在“可解释性”部分展示了证据,表明在 Evoformer 块中早期出现了具体的结构假设并不断完善。 Evoformer 模块的关键创新是在 MSA 内交换信息的新机制和允许直接推理空间和进化关系的配对表示。

网络的主干之后是结构模块,该模块以蛋白质的每个残基(全局刚体框架)的旋转和平移的形式引入了明确的 3-D 结构。这些表示在简单的状态下初始化,所有旋转设置为一致,所有位置设置为原点,但快速发展和完善,具有精确原子细节的高度准确的蛋白质结构。网络这一部分的关键创新包括打破链原子结构以允许同时对结构的所有部分进行局部细化,一种新颖的等变变换器允许网络隐式推理未表示的侧链原子,以及一个损失项代替残基的方向正确性的重要权重。在结构模块和整个网络中,我们通过反复将最终损失函数应用于输出,然后将输出递归地提供给相同的模块来强化迭代细化的概念。使用整个网络的迭代细化(我们称之为“循环”) 对准确性有显着贡献,而额外的训练时间很少(有关详细信息,请参阅附件-方法-1.8)。

Evoformer模块

名为 Evoformer(图 1e 和 3a)的网络构建块的关键原理是将蛋白质结构预测视为 3-D 空间中的图推理问题,其中图的边缘由邻近的残基定义。 配对表示的元素编码有关残基之间关系的信息(图 3b)。 MSA 表示的列编码输入序列的各个残基,而行表示这些残基出现的序列。 在这个框架内,我们定义了许多更新操作,这些更新操作应用于每个块中,其中不同的更新操作被串联应用。

Fig. 3 | Architectural details.

MSA 表示通过在 MSA 序列维度上求和的逐元素外积更新配对表示。 与之前的工作不同,此操作应用于每个块中,而不是在网络中应用一次,这使得从不断发展的 MSA 表示到配对表示的连续通信成为可能。

在配对表示中,有两种不同的更新模式。两者都受到配对表示一致性必要性的启发——为了将氨基酸的配对描述表示为单个 3-D 结构,必须满足许多约束,包括距离上的三角不等式。基于这种直觉,我们根据涉及三个不同节点的边三角形来安排对表示的更新操作(图 3c)。特别是,我们向轴向注意力添加了一个额外的 logit 偏差,以包括三角形的“缺失边”,并且我们定义了一个非注意力更新操作“三角形乘法更新”,它使用两条边来更新缺失的第三条边(参见 附件-方法-1.6.5 了解详情)。三角形乘法更新最初是作为注意力更对称且更便宜的替代品而开发的,仅使用注意力或乘法更新的网络都能够产生高精度结构。然而,两个更新的组合更准确。

我们还在 MSA 表示中使用了一种轴向注意力的变体。 在 MSA 中的 per-sequence attention 期间,我们从 pair stack 中投射额外的 logits 以偏置 MSA attention。 这通过提供从配对表示返回到 MSA 表示的信息流来关闭循环,确保整个 Evoformer 模块能够完全混合对和 MSA 表示之间的信息,并为结构模块中的结构生成做好准备。

端到端结构预测

结构模块(图 3d)使用对表示和来自主干的 MSA 表示的原始序列行(“single representation”,“单一表示”)在具体的 3-D 主干结构上运行。 3-D 主干结构表示为 Nres 独立的旋转和平移,每个旋转和平移相对于全局框架(residue gas,“残余气”?,图 3e)。这些旋转和平移,代表 N-Cα-C 原子的几何形状,优先考虑蛋白质骨架的方向,以便每个残基的侧链位置在该框架内受到高度限制。相反,肽键几何形状完全不受约束,并且在应用结构模块期间观察到网络经常违反链约束,因为打破此约束允许对链的所有部分进行局部细化,而无需解决复杂的闭环问题。在微调期间通过违规损失项鼓励满足肽键几何形状。只有在 Amber力场中的梯度下降结构的预测后松弛中才能实现肽键几何形状的精确执行。根据经验,这种最终松弛不会提高模型的准确性,如通过全局距离测试 (GDT) 或 IDDT-Cα34 测量的,但确实消除了分散注意力的立体化学违规而不会损失准确性。

AMBER力场是在生物大分子的模拟计算领域有着广泛应用的一个分子力场。AMBER力场的优势在于对生物大分子的计算,其对小分子体系的计算结果常常不能令人满意。

residue gas表示分两个阶段迭代更新(图 3d)。首先,我们称为不变点注意力 ( Point Attention) 的新型几何感知注意操作用于更新 Nres 神经激活集(single representation,“单一表示”)而不改变 3-D 位置,然后对residue gas使用更新的激活。不变点注意力通过在每个残基的局部框架中产生的 3-D 点来增强每个通常的注意力查询、键和值,这样最终值对全局旋转和平移是不变的(参见方法“不变点注意力(IPA)”了解详情)。 3-D 查询和键也对注意力施加了强烈的空间/局部性偏差,这非常适合蛋白质结构的迭代细化。在每个注意力操作和逐元素转换块之后,该模块计算每个主干帧的旋转和平移的更新。这些更新在每个残差的局部框架内的应用使得整体注意力和更新块成为对residue gas的等变操作。

侧链 chi 角的预测以及结构的最终每个残基精度 (pLDDT) 是在网络末端的最终激活上使用小的每个残基网络计算的。 TM 分数 (pTM) 的估计值是从成对错误预测中获得的,该预测被计算为最终对表示的线性投影。最后的损失(我们称之为帧对齐点误差(FAPE)(图 3f))将预测的原子位置与许多不同对齐下的真实位置进行比较。对于每个对齐,通过将预测帧 (Rk,tk) 对齐到相应的真实帧来定义,我们计算所有预测原子位置 xi 与真实原子位置的距离。由此产生的 Nframes × Natoms 距离受到限制的 L1 损失的惩罚。这对原子相对于每个残基的局部框架是正确的产生了强烈的偏见,因此在其侧链相互作用方面是正确的,并为 AlphaFold 提供了手性的主要来源(增刊-方法 1.9.3 和增刊-图 9)

使用标记和未标记数据进行训练

AlphaFold 架构能够仅使用对 PDB 数据的监督学习来训练到高精度,但我们能够使用类似于noisy-student自我蒸馏的方法来提高准确性(见图 4a)。 在这个过程中,我们使用一个训练有素的网络来预测来自 Uniclust30的约 350,000 个不同序列的结构,并将预测结构的新数据集过滤为高置信度子集。 然后,我们使用 PDB 和这个新的预测结构数据集的混合作为训练数据从头开始训练相同的架构,其中各种训练数据增强(例如裁剪和 MSA 子采样)使网络难以重现先前预测的结构。 这种自蒸馏过程有效地利用了未标记的序列数据,并显着提高了所得网络的准确性。

Self-training是最简单的半监督方法之一,其主要思想是找到一种方法,用未标记的数据集来扩充已标记的数据集。算法流程如下:
(1)首先,利用已标记的数据来训练一个好的模型,然后使用这个模型对未标记的数据进行标记。
(2)然后,进行伪标签的生成,因为我们知道,已训练好的模型对未标记数据的所有预测都不可能都是好的,因此对于经典的Self-training,通常是使用分数阈值过滤部分预测,以选择出未标记数据的预测标签的一个子集。
(3)其次,将生成的伪标签与原始的标记数据相结合,并在合并后数据上进行联合训练。
(4)整个过程可以重复n次,直到达到收敛。

此外,我们随机屏蔽或突变 MSA 中的单个残基,并具有来自 Transformers (BERT) 式目标的双向编码器表示来预测 MSA 序列的屏蔽元素。 这个目标鼓励网络学习解释系统发育和协变关系,而无需将特定的相关统计量硬编码到特征中。 与最近的独立工作相比,BERT 目标是在相同的训练示例上与正常的 PDB 结构损失联合训练的,并且没有进行预训练。

解释神经网络

为了了解 AlphaFold 如何预测蛋白质结构,我们为网络中的 48 个 Evoformer 块中的每一个训练了一个单独的结构模块,同时保持主网络的所有参数保持不变(补充方法 1.14)。包括我们的回收阶段,这提供了 192 个中间结构的轨迹,每个完整的 Evoformer 块一个,其中每个中间体代表网络对该块最可能结构的信念。在前几个块之后产生的轨迹出奇地平滑,表明 AlphaFold 对结构进行了不断的增量改进,直到它不能再改进为止(见图 4b 的准确度轨迹)。这些轨迹也说明了网络深度的作用。对于非常具有挑战性的蛋白质,如 SARS-CoV-2 Orf8 (T1064),网络搜索并重新排列多层的二级结构元素,然后再确定一个好的结构。对于 LmrP (T1024) 等其他蛋白质,网络会在前几层内找到最终结构。请参阅增刊-视频 1-4 的 CASP14 目标 T1024、T1044、T1064 和 T1091 的结构轨迹显示了一系列蛋白质大小和难度的清晰迭代构建过程。在补充-方法 1.16 和增刊-图12-13,我们解释了 AlphaFold 层产生的注意力图。

Fig. 4 | Interpreting the neural network

图 4a 包含 AlphaFold 组件的详细消融,表明各种不同的机制有助于 AlphaFold 的准确性。 请参阅增刊-方法 1.13 详细描述了每个消融模型、它们的训练细节、消融结果的扩展讨论以及 MSA 深度对每次消融的影响(补充-图 10)。

MSA 深度和跨链接触

虽然 AlphaFold 在绝大多数沉积的 PDB 结构中具有很高的准确性,但我们注意到仍然存在影响准确性或限制模型适用性的因素。当平均比对深度小于约 30 个序列时,该模型使用多个序列比对,准确度大幅下降(详见图 5a)。我们观察到阈值效应,其中 MSA 深度超过约 100 个序列的改进导致小增益。我们假设需要 MSA 信息在网络的早期阶段粗略地找到正确的结构,但是将该预测细化为高精度模型并不关键取决于 MSA 信息。我们观察到的另一个实质性限制是,与异型接触的数量相比,AlphaFold 对于链内或同型接触很少的蛋白质要弱得多。这通常发生在较大复合物中的桥接结构域,其中蛋白质的形状几乎完全由与复合物中其他链的相互作用产生。相反,AlphaFold 通常能够为同聚体提供高精度预测,即使链基本上交织在一起(例如图 5b)。我们希望 AlphaFold 的想法很容易适用于预测未来系统中的完整异质复合物,并且这将消除具有大量异质接触的蛋白质链的困难。

Fig. 5 | Effect of MSA depth and cross-chain contacts.

相关工作

蛋白质结构预测有一个漫长而多样的发展,在许多优秀的评论中得到了广泛的介绍。尽管将神经网络应用于结构预测的历史悠久,但它们最近才开始改进结构预测。这些方法通过将蛋白质结构预测问题处理为将进化耦合的“图像” 转换为蛋白质距离矩阵的“图像”,然后将距离预测集成到启发式系统中,从而有效地利用了计算机视觉系统的快速改进。产生最终的 3-D 坐标预测。最近开发了一些工作来直接预测 3-D 坐标 ,但这些方法的准确性与传统的手工结构预测管道不匹配 。同时,基于注意力的语言处理网络的成功 和最近的计算机视觉激发了对解释蛋白质序列的基于注意力的方法的探索。

讨论

我们在设计 AlphaFold 时采用的方法是生物信息学和物理方法的结合:我们使用物理和几何归纳偏差来构建组件,这些组件可以从 PDB 数据中学习,并最大限度地减少手工制作的特征(例如,AlphaFold 在没有氢的情况下有效地构建氢键债券评分函数)。这使得网络从 PDB 中的有限数据中更有效地学习,但能够应对结构数据的复杂性和多样性。特别是,AlphaFold 能够处理缺失的物理环境,并在具有挑战性的情况下生成准确的模型,例如交织的同源异构体或仅在未知血红素组存在时才会折叠的蛋白质。处理未指定结构条件的能力对于从 PDB 结构中学习至关重要,因为 PDB 代表了结构已被求解的所有条件。一般来说,AlphaFold 被训练产生最有可能作为 PDB 结构的一部分出现的蛋白质结构。在特定化学计量或配体/离子可单独从序列预测的情况下,AlphaFold 可能会产生一种隐式遵守这些约束的结构。

AlphaFold 已经向实验界展示了它的实用性,包括分子置换和解释低温电子显微镜 (cryo-EM) 图 。 此外,由于 AlphaFold 直接输出蛋白质坐标,因此 AlphaFold 会根据蛋白质序列的长度以图形处理单元 (GPU) 分钟到 GPU 小时生成预测(例如,对于 384 个残基,每个模型大约 1 个 GPU 分钟,请参阅方法了解详细信息 )。 这开辟了在蛋白质组规模及其他范围内预测结构的令人兴奋的可能性。

可用基因组测序技术和数据的爆炸式增长彻底改变了生物信息学,但实验结构测定的内在挑战阻止了我们结构知识的类似扩展。 通过开发准确的蛋白质结构预测算法,再加上由实验社区组装的现有大型且精心准备的结构和序列数据库,我们希望加速结构生物信息学的进步,以跟上基因组学革命的步伐。 我们希望 AlphaFold 以及将其技术应用于其他生物物理问题的计算方法将成为现代生物学的重要工具。

在线内容

任何方法、附加参考资料、自然研究报告摘要、源数据、扩展数据、补充信息、致谢、同行评审信息; 作者贡献和竞争利益的详细信息; 数据和代码可用性声明可在 Highly accurate protein structure prediction with AlphaFold – Nature 上获得。

 Highly accurate protein structure prediction with AlphaFold – Nature

dijkstra算法求最短路径

基本思想

dijkstra的算法思想是从以上最短距离数组中每次选择一个最近的点,将其作为下一个点,然后重新计算从起始点经过该点到其他所有点的距离,更新最短距离数据。已经选取过的点就是确定了最短路径的点,不再参与下一次计算。

  1. 通过Dijkstra计算图G中的最短路径时,需要指定起点s(即从顶点s开始计算)。
  2. 此外,引进两个集合S和U。S的作用是记录已求出最短路径的顶点(以及相应的最短路径长度),而U则是记录还未求出最短路径的顶点(以及该顶点到起点s的距离)。
  3. 初始时,S中只有起点s;U中是除s之外的顶点,并且U中顶点的路径是”起点s到该顶点的路径”。然后,从U中找出路径最短的顶点,并将其加入到S中;接着,更新U中的顶点和顶点对应的路径。 然后,再从U中找出路径最短的顶点,并将其加入到S中;接着,更新U中的顶点和顶点对应的路径。 … 重复该操作,直到遍历完所有顶点。

操作步骤

  1. 初始时,S只包含起点s;U包含除s外的其他顶点,且U中顶点的距离为”起点s到该顶点的距离”[例如,U中顶点v的距离为(s,v)的长度,然后s和v不相邻,则v的距离为∞]。
  2. 从U中选出”距离最短的顶点k”,并将顶点k加入到S中;同时,从U中移除顶点k。
  3. 更新U中各个顶点到起点s的距离。之所以更新U中顶点的距离,是由于上一步中确定了k是求出最短路径的顶点,从而可以利用k来更新其它顶点的距离;例如,(s,v)的距离可能大于(s,k)+(k,v)的距离。
  4. 重复步骤(2)和(3),直到遍历完所有顶点。

第1步:初始化距离,其实指与D直接连接的点的距离。dis[c]代表D到C点的最短距离,因而初始dis[C]=3,dis[E]=4,dis[D]=0,其余为无穷大。设置集合S用来表示已经找到的最短路径。此时,S={D}。现在得到D到各点距离{D(0),C(3),E(4),F(*),G(*),B(*),A(*)},其中*代表未知数也可以说是无穷大,括号里面的数值代表D点到该点的最短距离。

第2步:不考虑集合S中的值,因为dis[C]=3,是当中距离最短的,所以此时更新S,S={D,C}。接着我们看与C连接的点,分别有B,E,F,已经在集合S中的不看,dis[C-B]=10,因而dis[B]=dis[C]+10=13,dis[F]=dis[C]+dis[C-F]=9,dis[E]=dis[C]+dis[C-E]=3+5=8>4(初始化时的dis[E]=4)不更新。此时{D(0),C(3),E(4),F(9),G(*),B(13),A(*)}。

第3步:在第2步中,E点的值4最小,更新S={D,C,E},此时看与E点直接连接的点,分别有F,G。dis[F]=dis[E]+dis[E-F]=4+2=6(比原来的值小,得到更新),dis[G]=dis[E]+dis[E-G]=4+8=12(更新)。此时{D(0),C(3),E(4),F(6),G(12),B(13),A(*)}

第4步:在第3步中,F点的值6最小,更新S={D,C,E,F},此时看与F点直接连接的点,分别有B,A,G。dis[B]=dis[F]+dis[F-B]=6+7=13,dis[A]=dis[F]+dis[F-A]=6+16=22,dis[G]=dis[F]+dis[F-G]=6+9=15>12(不更新)。此时{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第5步:在第4步中,G点的值12最小,更新S={D,C,E,F,G},此时看与G点直接连接的点,只有A。dis[A]=dis[G]+dis[G-A]=12+14=26>22(不更新)。{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第6步:在第5步中,B点的值13最小,更新S={D,C,E,F,G,B},此时看与B点直接连接的点,只有A。dis[A]=dis[B]+dis[B-A]=13+12=25>22(不更新)。{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第6步:最后只剩下A值,直接进入集合S={D,C,E,F,G,B,A},此时所有的点都已经遍历结束,得到最终结果{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

python实现:

将以上的过程使用python来实现。
首先总结一个Dijkstra算法的核心思想,分成两步走:

  1. 构造一个最短路径数组,每次找到数组中未访问的节点里最小的点
  2. 以上一步的节点为最新节点,更新起始点到所有点的距离

使用python就是实现这两步即可



MAX= float('inf')

matrix = [
    [0,10,MAX,4,MAX,MAX],
    [10,0,8,2,6,MAX],
    [MAX,8,10,15,1,5],
    [4,2,15,0,6,MAX],
    [MAX,6,1,6,0,12],
    [MAX,MAX,5,MAX,12,0]
    ]


def dijkstra(matrix, start_node):
    
    #矩阵一维数组的长度,即节点的个数
    matrix_length = len(matrix)

    #访问过的节点数组
    used_node = [False] * matrix_length

    #最短路径距离数组
    distance = [MAX] * matrix_length

    #初始化,将起始节点的最短路径修改成0
    distance[start_node] = 0
    
    #将访问节点中未访问的个数作为循环值,其实也可以用个点长度代替。
    while used_node.count(False):
        min_value = float('inf')
        min_value_index = 999
        
        #在最短路径节点中找到最小值,已经访问过的不在参与循环。
        #得到最小值下标,每循环一次肯定有一个最小值
        for index in range(matrix_length):
            if not used_node[index] and distance[index] < min_value:
                min_value = distance[index]
                min_value_index = index
        
        #将访问节点数组对应的值修改成True,标志其已经访问过了
        used_node[min_value_index] = True

        #更新distance数组。
        #以B点为例:distance[x] 起始点达到B点的距离,
        #distance[min_value_index] + matrix[min_value_index][index] 是起始点经过某点达到B点的距离,比较两个值,取较小的那个。
        for index in range(matrix_length):
            distance[index] = min(distance[index], distance[min_value_index] + matrix[min_value_index][index])

    return distance




start_node = int(input('请输入起始节点:'))
result = dijkstra(matrix,start_node)
print('起始节点到其他点距离:%s' % result)

数据结构 — 单调队列

顾名思义,单调队列的重点分为 “单调” 和 “队列”

“单调” 指的是元素的的 “规律”——递增(或递减)

“队列” 指的是元素只能从队头和队尾进行操作

  • 单调递增队列:保证队列头元素一定是当前队列的最小值,用于维护区间的最小值。
  • 单调递减队列:保证队列头元素一定是当前队列的最大值,用于维护区间的最大值。

1、对于单调递增队列,对于一个元素如果它大于等于队尾元素, 那么直接把元素进队, 如果这个元素小于队尾元素,则将队尾元素出队列,直到满足这个元素大于等于队尾元素。

2、对于单调递减队列,对于一个元素如果它小于等于队尾元素, 那么直接把元素进队, 如果这个元素大于队尾元素,则将队尾元素出队列,直到满足这个元素小于等于队尾元素。

3、判断队头元素是否在待求解的区间之内,如果不在,就将其删除。(这个很好理解呀,因为单调队列的队头元素就是待求解区间的极值)

4、经过上面两个操作,取出 队列的头元素 ,就是 当前区间的极值

可以发现队列中的元素都是单调递减的(不然也就不叫单调递减队列啦),同时既有入队列的操作、也有出队列的操作。

实现单调队列,主要分为三个部分:

  • 去尾操作队尾元素出队列。当队列有新元素待入队,需要从队尾开始,删除影响队列单调性的元素,维护队列的单调性。(删除一个队尾元素后,就重新判断新的队尾元素)

去尾操作结束后,将该新元素入队列。

  • 删头操作队头元素出队列。判断队头元素是否在待求解的区间之内,如果不在,就将其删除。(这个很好理解呀,因为单调队列的队头元素就是待求解区间的极值)
  • 取解操作 :经过上面两个操作,取出 队列的头元素 ,就是 当前区间的极值

参考代码

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#define maxn 1000100
using namespace std;
int q[maxn], a[maxn];
int n, k;

void getmin() {  // 得到这个队列里的最小值,直接找到最后的就行了
  int head = 0, tail = 0;
  for (int i = 1; i < k; i++) {
    while (head <= tail && a[q[tail]] >= a[i]) tail--;
    q[++tail] = i;
  }
  for (int i = k; i <= n; i++) {
    while (head <= tail && a[q[tail]] >= a[i]) tail--;
    q[++tail] = i;
    while (q[head] <= i - k) head++;
    printf("%d ", a[q[head]]);
  }
}

void getmax() {  // 和上面同理
  int head = 0, tail = 0;
  for (int i = 1; i < k; i++) {
    while (head <= tail && a[q[tail]] <= a[i]) tail--;
    q[++tail] = i;
  }
  for (int i = k; i <= n; i++) {
    while (head <= tail && a[q[tail]] <= a[i]) tail--;
    q[++tail] = i;
    while (q[head] <= i - k) head++;
    printf("%d ", a[q[head]]);
  }
}

int main() {
  scanf("%d%d", &n, &k);
  for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
  getmin();
  printf("\n");
  getmax();
  printf("\n");
  return 0;
}

torch.nn.utils.rnn中的pack_padded_sequence和 pad_packed_sequence and pad_sequence

这两个是互逆的操作。

先来看为什么需要pad和pack操作:

先看一个例子,这个batch中有5个sample

如果不用pack和pad操作会有一个问题,什么问题呢?比如上图,句子“Yes”只有一个单词,但是padding了多余的pad符号,这样会导致LSTM对它的表示通过了非常多无用的字符,这样得到的句子表示就会有误差,更直观的如下图:

那么我们正确的做法应该是怎么样呢?

在上面这个例子,我们想要得到的表示仅仅是LSTM过完单词”Yes”之后的表示,而不是通过了多个无用的“Pad”得到的表示:如下图:

torch.nn.utils.rnn.pack_padded_sequence()

这里的pack,理解成压紧比较好。 将一个 填充过的变长序列 压紧。(填充时候,会有冗余,所以压紧一下)

其中pack的过程为:(注意pack的形式,不是按行压,而是按列压)

pack之后,原来填充的 PAD(一般初始化为0)占位符被删掉了。

输入的形状可以是(T×B×* )。T是最长序列长度,Bbatch size*代表任意维度(可以是0)。如果batch_first=True的话,那么相应的 input size 就是 (B×T×*)

Variable中保存的序列,应该按序列长度的长短排序,长的在前,短的在后。即input[:,0]代表的是最长的序列,input[:, B-1]保存的是最短的序列。

NOTE: 只要是维度大于等于2的input都可以作为这个函数的参数。你可以用它来打包labels,然后用RNN的输出和打包后的labels来计算loss。通过PackedSequence对象的.data属性可以获取 Variable

参数说明:

  • input (Variable) – 变长序列 被填充后的 batch
  • lengths (list[int]) – Variable 中 每个序列的长度。
  • batch_first (bool, optional) – 如果是True,input的形状应该是B*T*size
  • 参数 enforce_sorted ,如果是 True ,则输入应该是按长度降序排序的序列。如果是 False ,会在函数内部进行排序。默认值为 True 。也就是说在输入 pack_padded_sequence 前,我们也可以不对数据进行排序。

返回值:

一个PackedSequence 对象。

torch.nn.utils.rnn.pad_packed_sequence()

填充packed_sequence

上面提到的函数的功能是将一个填充后的变长序列压紧。 这个操作和pack_padded_sequence()是相反的。把压紧的序列再填充回来。填充时会初始化为0。

返回的Varaible的值的size是 T×B×*T 是最长序列的长度,B 是 batch_size,如果 batch_first=True,那么返回值是B×T×*

Batch中的元素将会以它们长度的逆序排列。

参数说明:

  • sequence (PackedSequence) – 将要被填充的 batch
  • batch_first (bool, optional) – 如果为True,返回的数据的格式为 B×T×*

返回值: 一个tuple,包含被填充后的序列,和batch中序列的长度列表

实例 代码:



>>> from torch.nn.utils.rnn import pack_padded_sequence, pad_packed_sequence
>>> seq = torch.tensor([[1,2,0], [3,0,0], [4,5,6]])
#seq的维度是3*3,第一个 维度表示T,is the length of the longest sequence 
#第二维表示B,批次大小,也就是说有3个长度为3的向量,其中每列表示一个序列,序列都是等长的,因为短的序列已经用0补齐了
>>> lens = [2, 1, 3]
>>> packed = pack_padded_sequence(seq, lens, batch_first=True, enforce_sorted=False)
#将补齐的序列压紧成一个序列,将0去掉
>>> packed
PackedSequence(data=tensor([4, 1, 3, 5, 2, 6]), batch_sizes=tensor([3, 2, 1]),
               sorted_indices=tensor([2, 0, 1]), unsorted_indices=tensor([1, 2, 0]))
>>> seq_unpacked, lens_unpacked = pad_packed_sequence(packed, batch_first=True)
>>> seq_unpacked
tensor([[1, 2, 0],
        [3, 0, 0],
        [4, 5, 6]])
>>> lens_unpacked
tensor([2, 1, 3])

pad_sequence

参数

sequences:表示输入样本序列,为 list 类型,list 中的元素为 tensor 类型。 tensor 的 size 为 L * F 。其中,L 为单个序列的长度,F 为序列中每个时间步(time step)特征的个数,根据任务的不同 F 的维度会有所不同。

batch_first:为 True 对应 [batch_size, seq_len, feature];False 对应[seq_len, batch_size, feature],从习惯上来讲一般设置为 True 比较符合我们的认知。

padding_value:填充值,默认值为 0 。

说明

主要用来对样本进行填充,填充值一般为 0 。我们在训练网络时,一般会采用一个一个 mini-batch 的方式,将训练样本数据喂给网络。在 PyTorch 里面数据都是以 tensor 的形式存在,一个 mini-batch 实际上就是一个高维的 tensor ,每个序列数据的长度必须相同才能组成一个 tensor 。为了使网络可以处理 mini-batch 形式的数据,就必须对序列样本进行填充,保证一个 mini-batch 里面的数据长度是相同的。

在 PyTorch 里面一般是使用 DataLoader 进行数据加载,返回 mini-batch 形式的数据,再将此数据喂给网络进行训练。我们一般会自定义一个 collate_fn 函数,完成对数据的填充。

import torch
from torch.utils.data import Dataset, DataLoader
from torch.nn.utils.rnn import pad_sequence,pack_padded_sequence,pack_sequence,pad_packed_sequence

class MyData(Dataset):
    def __init__(self, data):
        self.data = data

    def __len__(self):
        return len(self.data)

    def __getitem__(self, idx):
        return self.data[idx]

def collate_fn(data):
    data.sort(key=lambda x: len(x), reverse=True)
    data = pad_sequence(data, batch_first=True, padding_value=0)
    return data

a = torch.tensor([1,2,3,4])
b = torch.tensor([5,6,7])
c = torch.tensor([7,8])
d = torch.tensor([9])
train_x = [a, b, c, d]

data = MyData(train_x)
data_loader = DataLoader(data, batch_size=2, shuffle=True, collate_fn=collate_fn)
# 采用默认的 collate_fn 会报错
#data_loader = DataLoader(data, batch_size=2, shuffle=True) 
batch_x = iter(data_loader).next()

运行程序,得到 batch_x 的值:

# batch_x
tensor([[1, 2, 3, 4],
        [9, 0, 0, 0]])

从 batch_x 的值可以看出,第二行填充了三个 0 ,使其长度和第一行保持一致。

需要说明的是,对于长度不同的序列,使用默认的 collate_fn 函数,不自定义 collate_fn 函数完成对序列的填充,上面的程序就会报错。

Python3 迭代器与生成器

迭代器

迭代是Python最强大的功能之一,是访问集合元素的一种方式。

迭代器是一个可以记住遍历的位置的对象。

迭代器对象从集合的第一个元素开始访问,直到所有的元素被访问完结束。迭代器只能往前不会后退。

迭代器有两个基本的方法:iter() 用于生成迭代器 和 next() 求迭代器的值

字符串,列表或元组对象都可用于创建迭代器:

>>> list=[1,2,3,4]
>>> it = iter(list)    # 创建迭代器对象
>>> print (next(it))   # 输出迭代器的下一个元素
1
>>> print (next(it))
2
>>>

迭代器对象可以使用常规for语句进行遍历:

list=[1,2,3,4]
it = iter(list)    # 创建迭代器对象
for x in it:
    print (x, end=" ")

执行以上程序,输出结果如下:

1 2 3 4

也可以使用 next() 函数:

import sys         # 引入 sys 模块
 
list=[1,2,3,4]
it = iter(list)    # 创建迭代器对象
 
while True:
    try:
        print (next(it))
    except StopIteration:
        sys.exit()

执行以上程序,输出结果如下:

1
2
3
4

创建一个迭代器

把一个类作为一个迭代器使用需要在类中实现两个方法 __iter__() 与 __next__() 。

如果你已经了解的面向对象编程,就知道类都有一个构造函数,Python 的构造函数为 __init__(), 它会在对象初始化的时候执行。

更多内容查阅:Python3 面向对象

__iter__() 方法返回一个特殊的迭代器对象, 这个迭代器对象实现了 __next__() 方法并通过 StopIteration 异常标识迭代的完成。

__next__() 方法(Python 2 里是 next())会返回下一个迭代器对象。

创建一个返回数字的迭代器,初始值为 1,逐步递增 1:

class MyNumbers:
  def __iter__(self):
    self.a = 1
    return self
 
  def __next__(self):
    x = self.a
    self.a += 1
    return x
 
myclass = MyNumbers()
myiter = iter(myclass)
 
print(next(myiter))
print(next(myiter))
print(next(myiter))
print(next(myiter))
print(next(myiter))

执行输出结果为:

1
2
3
4
5

StopIteration

StopIteration 异常用于标识迭代的完成,防止出现无限循环的情况,在 __next__() 方法中我们可以设置在完成指定循环次数后触发 StopIteration 异常来结束迭代。

在 20 次迭代后停止执行:

class MyNumbers:
  def __iter__(self):
    self.a = 1
    return self
 
  def __next__(self):
    if self.a <= 20:
      x = self.a
      self.a += 1
      return x
    else:
      raise StopIteration
 
myclass = MyNumbers()
myiter = iter(myclass)
 
for x in myiter:
  print(x)

执行输出结果为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

生成器

yield是python的一个关键字,本质上是一个生成器generator。

生成器是一种特殊的函数,它会返回一个迭代器。定义一个生成器函数同定义一个普通函数没有什么区别,特殊之处在于生成器函数内部会包含yield表达式,专门用于生成一个序列。当一个生成器函数被调用时,它会返回一个迭代器。之后就由这个迭代器来控制生成器函数的执行。当生成器函数被调用后,首先会执行到第一个yield表达式处,然后会将生成器函数挂起,将yield生成的表达式的值返回给生成器函数的调用者。当生成器函数被挂起时,它的所有局部状态都会被保存起来,包括当前绑定的局部变量、指令指针、内部栈和异常处理的状态。当调用迭代器的方法时,一般都是调用next()方法,将会恢复生成器函数的执行,并且是从上次被挂起的地方继续执行,直到遇到另外一次yield调用,生成器函数将再次被挂起。

在一个生成器函数中,如果没有 return,则默认执行至函数完毕,如果在执行过程中 return,则直接抛出 StopIteration 终止迭代。

可以中断一个函数的执行,跳转到调用的地方

在 Python 中,使用了 yield 的函数被称为生成器(generator)。

跟普通函数不同的是,生成器是一个返回迭代器的函数,只能用于迭代操作,更简单点理解生成器就是一个迭代器。

在调用生成器运行的过程中,每次遇到 yield 时函数会暂停并保存当前所有的运行信息,返回 yield 的值, 并在下一次执行 next() 方法时从当前位置继续运行。

调用一个生成器函数,返回的是一个迭代器对象。

以下实例使用 yield 实现斐波那契数列:

import sys
 
def fibonacci(n): # 生成器函数 - 斐波那契
    a, b, counter = 0, 1, 0
    while True:
        if (counter > n): 
            return
        yield a
        a, b = b, a + b
        counter += 1
f = fibonacci(10) # f 是一个迭代器,由生成器返回生成
 
while True:
    try:
        print (next(f), end=" ")
    except StopIteration:
        sys.exit()

执行以上程序,输出结果如下:

0 1 1 2 3 5 8 13 21 34 55

数据结构 — 堆

python中直接使用heapq建立 小根堆

http://139.9.1.231/index.php/2022/02/22/python3heapq/

堆简介

堆是一棵树,其每个节点都有一个键值,且每个节点的键值都大于等于/小于等于其父亲的键值。

每个节点的键值都大于等于其父亲键值的堆叫做小根堆,否则叫做大根堆。STL 中的 priority_queue 其实就是一个大根堆。

(小根)堆主要支持的操作有:插入一个数、查询最小值、删除最小值、合并两个堆、减小一个元素的值。

一些功能强大的堆(可并堆)还能(高效地)支持 merge 等操作。

一些功能更强大的堆还支持可持久化,也就是对任意历史版本进行查询或者操作,产生新的版本。

二叉堆

结构

从二叉堆的结构说起,它是一棵二叉树,并且是完全二叉树,每个结点中存有一个元素(或者说,有个权值)。

堆性质:父亲的权值不小于儿子的权值(大根堆)。同样的,我们可以定义小根堆。本文以大根堆为例。

由堆性质,树根存的是最大值(getmax 操作就解决了)。

插入操作

插入操作是指向二叉堆中插入一个元素,要保证插入后也是一棵完全二叉树。

最简单的方法就是,最下一层最右边的叶子之后插入。

如果最下一层已满,就新增一层。

插入之后可能会不满足堆性质?

向上调整:如果这个结点的权值大于它父亲的权值,就交换,重复此过程直到不满足或者到根。

可以证明,插入之后向上调整后,没有其他结点会不满足堆性质。

向上调整的时间复杂度是  的。

二叉堆的插入操作

删除操作

删除操作指删除堆中最大的元素,即删除根结点。

但是如果直接删除,则变成了两个堆,难以处理。

所以不妨考虑插入操作的逆过程,设法将根结点移到最后一个结点,然后直接删掉。

然而实际上不好做,我们通常采用的方法是,把根结点和最后一个结点直接交换。

于是直接删掉(在最后一个结点处的)根结点,但是新的根结点可能不满足堆性质……

向下调整:在该结点的儿子中,找一个最大的,与该结点交换,重复此过程直到底层。

可以证明,删除并向下调整后,没有其他结点不满足堆性质。

减小某个点的权值

很显然,直接修改后,向上调整一次即可

参考代码:


void up(int x) {
  while (x > 1 && h[x] > h[x / 2]) {
    swap(h[x], h[x / 2]);
    x /= 2;
  }
}

void down(int x) {
  while (x * 2 <= n) {
    t = x * 2;
    if (t + 1 <= n && h[t + 1] > h[t]) t++;
    if (h[t] <= h[x]) break;
    std::swap(h[x], h[t]);
    x = t;
  }
}

建堆

方法一:使用 decreasekey(即,向上调整)¶
从根开始,按 BFS 序进行。


void build_heap_1() {
  for (i = 1; i <= n; i++) up(i);
}
方法二:使用向下调整¶
这时换一种思路,从叶子开始,逐个向下调整


void build_heap_2() {
  for (i = n; i >= 1; i--) down(i);
}

并查集

一、概念及其介绍

并查集是一种树型的数据结构,用于处理一些不相交集合的合并及查询问题。

并查集的思想是用一个数组表示了整片森林(parent),树的根节点唯一标识了一个集合,我们只要找到了某个元素的的树根,就能确定它在哪个集合里。

二、适用说明

并查集用在一些有 N 个元素的集合应用问题中,我们通常是在开始时让每个元素构成一个单元素的集合,然后按一定顺序将属于同一组的元素所在的集合合并,其间要反复查找一个元素在哪个集合中。这个过程看似并不复杂,但数据量极大,若用其他的数据结构来描述的话,往往在空间上过大,计算机无法承受,也无法在短时间内计算出结果,所以只能用并查集来处理。

三、并查集的基本数据表示

如上图 0-4 下面都是 05-9 下面都是 1,表示 0、1、2、3、4 这五个元素是相连接的,5、6、7、8、9 这五个元素是相连的。

再如上图 0、2、4、6、8 下面都是 0 这个集合,表示 0、2、4、6、8 这五个元素是相连接的,1、3、5、7、9 下面都是 1 这个集合,表示 0,1、3、5、7、9 这五个元素是相连的。

构造一个类 UnionFind,初始化, 每一个id[i]指向自己, 没有合并的元素:

public UnionFind1(int n) {
        count = n;
        id = new int[n];
        // 初始化, 每一个id[i]指向自己, 没有合并的元素
        for (int i = 0; i < n; i++)
            id[i] = i;
    }

查找

通俗地讲一个故事:几个家族进行宴会,但是家族普遍长寿,所以人数众多。由于长时间的分离以及年龄的增长,这些人逐渐忘掉了自己的亲人,只记得自己的爸爸是谁了,而最长者(称为「祖先」)的父亲已经去世,他只知道自己是祖先。为了确定自己是哪个家族,他们想出了一个办法,只要问自己的爸爸是不是祖先,一层一层的向上问,直到问到祖先。如果要判断两人是否在同一家族,只要看两人的祖先是不是同一人就可以了。

在这样的思想下,并查集的查找算法诞生了。

# Python Version
fa = [0] * MAXN # 记录某个人的爸爸是谁,特别规定,祖先的爸爸是他自己

# 递归
def find(x):
    # 寻找x的祖先
    if fa[x] == x:
        return x # 如果x是祖先则返回
    else:
        return find(fa[x]) # 如果不是则 x 的爸爸问 x 的爷爷

# 非递归
def find(x):
    while x != fa[x]: # 如果 x 不是祖先,就一直往上一辈找
        x = fa[x]
    return x # 如果 x 是祖先则返回

这样的确可以达成目的,但是显然效率实在太低。为什么呢?因为我们使用了太多没用的信息,我的祖先是谁与我父亲是谁没什么关系,这样一层一层找太浪费时间,不如我直接当祖先的儿子,问一次就可以出结果了。甚至祖先是谁都无所谓,只要这个人可以代表我们家族就能得到想要的效果。把在路径上的每个节点都直接连接到根上,这就是路径压缩。

# Python Version
def find(x):
    if x != fa[x]: # x 不是自身的父亲,即 x 不是该集合的代表
        fa[x] = find(fa[x]) # 查找 x 的祖先直到找到代表,于是顺手路径压缩
    return fa[x]

合并

宴会上,一个家族的祖先突然对另一个家族说:我们两个家族交情这么好,不如合成一家好了。另一个家族也欣然接受了。
我们之前说过,并不在意祖先究竟是谁,所以只要其中一个祖先变成另一个祖先的儿子就可以了。

# Python Version
def unionSet(x, y):
    # x 与 y 所在家族合并
    x = find(x)
    y = find(y)
    fa[x] = y # 把 x 的祖先变成 y 的祖先的儿子