在 python 中实现基本元胞自动机时如何利用并行性？

Question

这更像是一个了解实现并行性的适当方法的练习。我不一定需要一个最佳程序来运行基本的元胞自动机。

假设我们在 1d 中布置了 N=10000 单元格。 ECA 通过在每个单元格 i 上应用本地规则 f 来工作，给定它的邻居状态和它自己的状态：

c[t][i] = f(c[t-1][i-1], c[t-1][i], c[t-1][i+1])

在 9999 个时间步长上进化它的蛮力方法是这样的：

N = 10000
c[0] = np.random.randint(2, size=N)
for t in range(1, 9999):
    for i in range(N):
        temp[i] = f(c[t-1][i-1], c[t-1][i], c[t-1][i+1])
    c[t] = temp.copy()

给予或接受一些变量初始化。

有没有办法在每个时间步并行化 N=10000 应用程序？这应该使用线程来完成吗？还是那些地方规则太快太简单而无法真正利用任何东西？我对利用代码的并行性知之甚少，但这似乎是它的自然设置。

Answer 1

手头的问题

最近我一直在围绕这个问题进行一段时间的研究，以了解 ECA 的多线程应用程序的可行性，并使用一些关于多线程的基本直觉和基本元胞自动机的性质，多线程不太可能将完全有助于程序，并且在大多数情况下只会减慢模拟速度，尤其是在像 python 这样相对未优化的语言案例的情况下。

这个问题的关键问题是多线程需要一个能够分布式工作负载的问题，理论上通用元胞自动机是，但是，一维元胞自动机依赖于上一代的数据（i-1）为了推导出当前一代 (i)。结果，如果我们将 n 个线程专用于长度为 n*256 的 1 和 0 的数组，每个线程可以计算规则 258 次（+2，因为每个分区的边界）。起初，这似乎是一个完全合理的解决方案，但是，如果目标是速度/效率，那么根据手头的语言，它会在不同程度上表现得更差。

为什么不用多线程？

因为在现代处理器上计算 258 次比较的速度非常快（无论您是使用逐位比较还是以编程方式进行），而且即使使用固定大小的蜂窝，启动线程的计算成本也非常高您可以在您选择使用的任何数据结构上显式定义每个线程的边界的自动机，通过并行化这个简单的计算可​​以获得的性能只是不存在。

除了启动线程的成本之外，除非您使用为并发而构建的语言（即 go-lang）进行编程，其中内置了用于加速线程初始化和处理内存访问冲突的工具（您可以必须处理给定每个新单元取决于其重叠的最近邻居），您必须在循环中写入一些警告以处理冲突（主要在内存中）和潜在的系统故障（例如线程被某些人锁定）系统进程弄乱了每个线程必须比较多少数组）。你的程序中的所有这些边缘情况都会显着减慢它，尤其是在 python 中，除非你的元胞自动机的初始条件或即时进化空间在一次进化中是 >= n*500000 的顺序（根据我的经验），其中计算对于每 500000 个单元，单个线程的时间大约在 10 秒以上，因此不会提高性能，并且很可能会因使用多线程而降低性能。

计算不可约性

换句话说，每个线程一次只能对一段进化进行计算，我们不能使用线程来垂直而不是水平地计算元胞自动机。这是 wolfram 的 computational irreducibility 理念的核心（即，必须按顺序计算整个自动机，以便程序可以访问所有之前的数据）。尽管在某些情况下，规则可以推广到布尔表达式（参考），这可能是可能的，因为结果是可预测的，但在一般情况下，这不是由于规则 30、110 和 54 等显示独特的计算不可约行为的规则。

我们能做什么？

话虽如此，有可能使用或不使用布尔表达式（如下所述）对元胞自动机进行伪多线程处理。 Wolfram 的新科学 (reference) 中描述的其中一种方法是利用 CPU 的内部寄存器大小来并行化比较，书中描述的方法称为“位切片”。

您可能经常听到人们提到计算机是 64 位与 32 位系统，这取决于内部寄存器的大小。我说我们有一台 64 位计算机，而不是使用 1/64 的寄存器大小顺序比较每个 1 和 0，我们可以比较 64 位或 8 字节的单元，从而在单核上获得理论上的 8 倍速度收益。我们可以通过在 8 位计算机上将一维元胞自动机的初始条件垂直排列为 3 个 8 位数字来做到这一点。

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

1  1  1  0  1  1  0  1  1  0  1  1  0  1  1  0  0  1  1  0  1  1  1  0
     R1|      R2|      R3|      R4|      R5|     R6|      R7|      R8|

->

   C1 C2 C3 
 
R1| 1  1  1
R2| 0  1  1
R3| 0  1  1
R4| 0  1  1
R5| 0  1  1
R6| 0  0  1
R7| 1  0  1
R8| 1  1  0    

C1 = 131
C2 = 249
C3 = 254

在此配置中，每个位的左右索引对应于下一个十进制数的同一行上的位。使用这种配置，我们可以使用规则的布尔表达式和这些数字上的按位运算符（同时移动第一个和最后一个数字）来评估规则。 比如关于Wolfram Alpha可以描述为...

(p, q, r) -> (q AND (NOT p)) OR (q XOR r)  

如果我们使用在二维模拟中预先分配所有前导零的位切片数组，在 go-lang 中看起来像这样。我不太熟悉 Python 访问裸机指令，但我可以肯定地说这可以在 python 中复制，但在 go-lang 等静态语言上非常快。

func historicallyAware(evolutions int) {
     for i := 1; i < evolutions+1; i++ {
         sim[i][0] = ((^(sim[i-1][len(sim[i-1])-1]) << 1) & sim[i-1][0]) | (sim[i-1][0] ^ sim[i-1][1])
         for j := 1; j < len(sim[i])-1; j++ {
            sim[i][j] = ((^(sim[i-1][j-1])) & sim[i-1][j]) | (sim[i-1][j] ^ sim[i-1][j+1])
         }
         sim[i][len(sim[i])-1] = ((^(sim[i-1][len(sim[i])-2])) & sim[i-1][len(sim[i])-1]) | (sim[i-1][len(sim[i])-1] ^ sim[i-1][0]>>1)
    }
}