我的NW算法的实现,这个东西蛮有意思的。@jit/@njit来做加速函数的部分,发现整体上来看反而不快,比较奇怪。
在代码层面,numba装饰的函数在循环体中进行运转反而不快,要优化还是得看并行化编程,或者我直接用C写一遍,抑或着是信息办法改掉数据结构,让循环体也能被@jit修饰。在算法层面,我应该尝试优化算法结构和数据结构,我使用了python原生的list对象,这导致我在进行numba加速的时候被拖累了。
实验室新人来做NW算法了,我来尝试在上次的基础上用C语言重写一遍。上次的代码就其功能而言是完备的,也能处理复杂序列的情况,在不考虑时间消耗的情况下。 但是,当我们必须考虑时间及产生的序列的可用性的时候,之前的代码就不合适了。
我首先想到的就是把python代码转换为C的代码,这部分的问题依旧出在性能上。我使用的用来做测试的序列如下:
// 两个需要比对的DNA序列
char A[] = "ATCGGGCTACATCGGGCTACGGATCGGGCTACGAAAAAAAAAAAAAAAAAAAAAA";
char B[] = "ATCGGATCGGGCTACGATCGGGCTACGATCGCGTTTTTTTAAAAAGAAAAAAAAGGGGGGGGTGTATTGTA";
在实践中我发现,如果我继续用遍历所有可能路径的方式来做这个序列的比对,其可能的路径的数量级在“亿”的级别。 即使用C语言来编写搜索算法也很慢,我开始寻求并行化的方法。现在的做法是我们依次对每个可能的路径依次搜索,作为并行化的更改我们在这一步进行拆分,并行的进行每个路径的遍历。
在并行化这一步遇到了新的问题,并行化的时候遇到了Segmentation fault (core dumped)
的问题。看起来是在内存分配上遇到了问题。在这一步,我选择了“剪枝”的思路来处理这个问题。即考虑到比对结果的可用性,在比对的前期对可能的路径进行选择,去除一部分结果。
然而,这样的思路导致了另外一个问题,一些在前期表现较好的比对结果更容易出现在最终的结果中,而那些在稍微靠后的位置上才表现优异的比对结果就无法体现了。
我本来就是实现贪婪算法,再“贪婪”一点也不是问题
这个问题的本质是: - 早期剪枝:可能会提前剪掉一些后续整体得分较高的路径,从而导致最终结果不佳。 - 局部最优:在路径生成过程中,可能会优先保留局部最优的路径,而忽略了全局最优的路径。
为了避免这个问题,可以考虑以下几种策略:
-
延迟剪枝: 先生成所有可能的路径,计算它们的得分,然后在所有路径中选择得分最高的前5000个。这种方法避免了提前剪枝带来的问题,但代价是需要更多的内存和计算时间。
-
分段剪枝: 在每个部分(即从矩阵的某一行到某一行,或某一列到某一列)进行剪枝,而不是在每个步骤进行剪枝。这样可以让算法有机会探索更多的路径组合,从而增加全局最优解出现的可能性。
-
比如优先队列剪枝(Best-First Search): 使用优先队列存储路径,路径的优先级基于当前的总得分(已计算部分的得分 + 对未计算部分的一个估计得分)。每次扩展时,从队列中取出优先级最高的路径进行扩展,这样可以尽量保留潜在得分高的路径。 但是这样做也有显著的缺点,那就是这样的结果并不一定是最优的
说到这里,就不得不提路径搜索的三大 baseline 算法 Dijkstra(荷兰语 迪斯科特拉)、Best-first 以及 A* search 算法。
今天不进行展开,只先尝试用best-first先完成这个算法作业。可能对于人类而言,只需要前10个比对结果就OK👌了。
代码附在文件NW实现3.c
中。
教授让新人重新搞一下NW的算法实现,新人的程序和讲解都是问题。于是我也来重写一下。主要是填上之前的几个坑。
- 逻辑优化:我看别人搞得都非常简洁,我开始怀疑自己的实现的思路。结果我一跑他们的程序就知道了,根本不是那回事,果然复杂性不会消失,只会转移。那些简洁的代码是有代价的,他们的空间复杂度要比我高,我的实现方式相较之下非常节约内存,我以后也应该首先关注算法的时间和空间复杂度,再考虑实现方式。
- 逻辑优化2:坦白讲,我写的时候根本没什么空间优化的想法在,我只是模拟了我进行手动比对的过程,即我是首先用伪代码模拟了人类计算的过程。人类这种生物先天就内存小,所以模拟人类的话内存也小?
- 实现几个baseline算法,我查阅了一下他们的实现方法,如果我真的那样做了,是有一些背离本意的。另外一个比较巧的是,我过去猜测的,对矩阵分块进行路径回溯:即局部的最优的和等于全局最优的思路,竟然就是A*搜索的主要思想,那就不言者这个方向继续了。
- 对于并行化的处理:之前技术路线歪到剪枝的主要原因就是C语言的并行化计算时候的内存操作我没整明白,导致超限。对这个问题的最直接的解决办法是换用Rust语言来实现。Rust在与C语言性能相当的同时实现了自动化的内存管理。
- 对于并行化的处理2:并行化的主要推动力是在python实现NW算法的时候产生的,对于长序列,python太慢了,同时python不太好实现完整的并行,因为python有GIL机制的存在。我按照一般的思路,对计算方向的部分进行并行化处理遇到了很多问题,什么数据竞争,什么数据封闭之类的,也还产生了字符串拼接的问题。不得不放弃。
- 对于并行化的处理3:我放弃的是那种处理并行化的方式,并不是放弃并行化处理本身。我引入了子任务调度的概念,进行分层并行化。即当且仅当需要计算的任务超过一定限度的时候,再去独立子进程,且子进程的数量不是无限的,超限的任务排队。所有计算完成之后再合并结果。 这样处理结果出人意料的好,我莫不是个天才🤣。
- 对于并行化的处理4:并行化处理大成功!!!! Rust好强!
- 对于并行化的处理5:计算完成后把写过写在文件里也是问题,这里如果结果很多就会像是卡住一下下。因为采取的方式是先写在小文件里再合并起来。那么这里的优化思路是明显的,人类可读的结果当且仅当的是十几个结果。所以会有两个版本
- 这次没有引入兼并碱基之类的东西,也没有考虑到进化上的关系。下次再把这个东西实现,做成主程序、配置文件的形式。
最近始终在想怎么实现并行化,一个思路就是对可能经过的矩阵的区域进行拆分,这样可以同时计算所有矩阵。这样我们减轻了每次计算的粒度,并且增加了并行计算的宽度。可以借助硬件进行进一步的大规模并行计算。现在借助python实现了这个算法,但是性能有待进一步的优化。这里也可以导向rust来重写。当然也有借助cuda进行改写的计划,主要的思路还是和之前一致,借助GPU的并行计算的能力进行加速,尽可能降低沟通的损耗。 只不过这样带来了新的问题,即如何把产生的那些新的路径(短又多)拼接起来,这个问题非常类似短reads拼接的问题。解决的思路也很一致。 本阶段的算法练习让我一步步学到了相当多的算法和技术知识。 在nw实现7.8中,我们完成了主要的算法实现,可以利用新的算法计算出所有可能的路径。但在路径拼接和转换上还有欠缺,会在接下里的时间中进一步完善。但也可能先去搞CUDA函数的编写,GPU的异构调用非常符合我的算法的直觉。
BTW,接下来在算法上的革新主要会做“矩阵的启发式稀疏计算”,尝试通过启发式的方式只计算需要计算的小矩阵。即通过什么方式预测路径在矩阵中的边界,这样来减少计算量。
矩阵拆分图如下
进一步修订了并行化的思路。 之前的思路需要拆分小矩阵,原计划是按照一定的方式对矩阵进行拆分,但这里有两个问题,一个是如果按照最小因数拆分,存在质数无法进一步拆分的问题;此外,这样的拆分方式增加了处理矩阵接缝区域的复杂度。我进一步深化拆分的粒度,已知对于NW算法这样的矩阵,反对角线的元素相互之间没有依赖关系,因此从这个角度进行并行化的拆分。
这里更新了两类算法实现,
- 第一是在之前的NW实现7.8的基础上进行并行化方式的修改,改为放对角线上元素的并行计算,且只保留每个元素的指针(就是方向)。回溯的时候只对着最高分的方向回溯(潜在问题是分数相同怎么办,这里是记录了多种路线)。形成的文件在NW算法实现8.0
- 第二是在NW算法实现9.0的基础上进行修改, 让其能够把CUDA的优点利用起来,把主要的计算放在GPU上执行,速度非常快,上百倍的加速。
在GPU的帮助下,我的代码终于和现有的成熟工具可以在同一个baseline,这波是大力出奇迹。在这里发现了,当我们要求非常多的路径的时候,主要的限速步骤是在回溯的地方,回溯也是接下来的加速的一个着力点。形成的文件在NW算法实现9.0_CUDA- 调用方式也进行了修正,
python nw9.py -i input.fasta -o alignment_output.txt -g 10.0 -e 0.5 -p 5
- 调用方式也进行了修正,
下面是使用的截图
接下来的重点
- 把已有的成熟的打分矩阵规则引入NW算法实现9.0_CUDA中
- cpu版本进行修正,写一个NW8.0_rust的版本
48GB显存也就是75kb长度的序列,目标是1Mb,这个还有些距离。 内存和内核的优化可以继续尝试