在Haskell中共享力计算

Question

在Haskell中共享力计算

8

我正在使用Haskell实现一个N-Body模拟。 https://github.com/thorlucas/N-Body-Simulation

目前，每个粒子都会计算自己的力量，然后针对其他每个粒子计算加速度。换句话说，需要进行O(n²)次力量计算。如果我只计算每种组合一次，可以将其减少为O(n选择2)。

let combs = [(a, b) | (a:bs) <- tails ps, b <- bs ]
    force = map (\comb -> gravitate (fst comb) (snd comb)) combs

但是我不知道如何在没有使用状态的情况下应用它们到粒子上。在上面的例子中，ps 是一个 [Particle] 数组，其中

data Particle = Particle Mass Pos Vel Acc deriving (Eq, Show)

从理论上讲，在一个有状态的语言中，我可以很容易地遍历组合，计算每个“a”和“b”的力所产生的加速度，然后在这样做时更新“ps”中每个“Particle”的加速度。

我曾经考虑过类似于“foldr f ps combs”的东西。起始累加器将是当前的“ps”，而“f”将是一些函数，它接受每个“comb”，并更新“ps”中相关的“Particle”，并返回该累加器。这似乎对于如此简单的过程来说非常占用内存且相当复杂。

有什么想法吗？

- Thor Correia

4

请注意，O (n · (n – 1)) ≡ O (n²)。仅计算每个力一次并不能使您在真正称之为“多体”模拟时获得可接受的性能。要在这种计算中真正节省渐近成本，您需要使用逼近方法；通常可以使用类似于八叉树的东西来实现（参见：https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation）。 - leftaroundabout

1

@ThorCorreia：O(n!/(2(n-2)!)) = O(n!/(n-2)!) = O(n · (n - 1)) = O(n²) - Ry-

我不太理解 @Ryan 的观点。以选择 200 中的 2 个数为例，结果是 19,900，而使用 200(199) 则是 39,800。计算量几乎增加了一倍。随着问题规模变大，这种差异会更加明显。例如选择 2,000 中的 2 个数仅约为 2e6，而使用 2000(1999) 却是 4e6。 - Thor Correia

1

@ThorCorreia的O（发音为“大O”）符号实际上并不是某个输入所需计算次数。它表示随着n趋近于无穷大时计算次数的增长率。另一种思考方式是“当n变得很大时，哪个因素最终占主导地位”。虽然在评论中深入讨论可能有些困难！ - Michal Charemza

1

在规模上，我们不关心常数乘数。如果“两倍”始终是“两倍”，那么它就相当无害。例如，O(n)和O(2^n)之间的区别在于，如果您将输入大小加倍，则前者的时间也会加倍，而如果您仅将输入大小增加1，则后者的时间会加倍。这是一个规模问题，而不是实际数量问题。 - Silvio Mayolo

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Michal Charemza · Accepted Answer

从https://github.com/thorlucas/N-Body-Simulation获取代码。

updateParticle :: Model -> Particle -> Particle
updateParticle ps p@(Particle m pos vel acc) =
    let accs = map (gravitate p) ps
        acc' = foldr (\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
        vel' = (fst vel + fst acc, snd vel + snd acc)
        pos' = (fst pos + fst vel, snd pos + snd vel)
    in  Particle m pos' vel' acc'

step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) ps

通过将加速度计算出来的矩阵（实际上是列表的列表）与每个粒子的更新分开处理并进行修改，我们得到...

updateParticle :: Model -> (Particle, [Acc]) -> Particle
updateParticle ps (p@(Particle m pos vel acc), accs) =
    let acc' = foldr (\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
        vel' = (fst vel + fst acc, snd vel + snd acc)
        pos' = (fst pos + fst vel, snd pos + snd vel)
    in  Particle m pos' vel' acc'

step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) $ zip ps accsMatrix where
    accsMatrix = [map (gravitate p) ps | p <- ps]

所以问题本质上是如何通过使用gravitate a b = -1 * gravitate b a的事实，在计算accsMatrix时减少对gravitate的调用次数。

如果我们打印出accsMatrix，它会像这样...

[[( 0.0,  0.0), ( 1.0,  2.3), (-1.0,  0.0), ...
[[(-1.0, -2.3), ( 0.0,  0.0), (-1.2,  5.3), ...
[[( 1.0,  0.0), ( 1.2, -5.3), ( 0.0,  0.0), ...
...

...因此我们可以看到accsMatrix !! i !! j == -1 * accsMatrix !! j !! i。

因此，要使用上述事实，我们需要访问一些索引。首先，我们对外部列表进行索引...

accsMatrix = [map (gravitate p) ps | (i,p) <- zip [0..] ps]

...并使用列表推导式替换内部列表...

accsMatrix = [[ gravitate p p' | p' <- ps] | (i,p) <- zip [0..] ps]

"

...通过zip获取更多可用的索引...

"

accsMatrix = [[ gravitate p p' | (j, p') <- zip [0..] ps] | (i,p) <- zip [0..] ps]

然后，关键是让accsMatrix在矩阵的一半中依赖于自身...

accsMatrix = [[ if i == j then 0 else if i < j then gravitate p p' else -1 * accsMatrix !! j !! i | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]

我们也可以将其分成几部分，如下所示...

accsMatrix = [[ accs (j, p') (i, p) | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]
accs (j, p') (i, p) 
    | i == j    = 0
    | i < j     = gravitate p p'
    | otherwise = -1 * accsMatrix !! j !! i

... 或者通过使用map避免列表理解。

accsMatrix = map (flip map indexedPs) $ map accs indexedPs  
indexedPs = zip [0..] ps
accs (i, p) (j, p')
    | i == j    = 0
    | i < j     = gravitate p p'
    | otherwise = -1 * accsMatrix !! j !! i

...或者通过使用列表单子（list monad）...

accsMatrix = map accs indexedPs >>= (:[]) . flip map indexedPs

虽然对我来说，这些东西更难理解。

这种列表嵌套的方法可能会存在一些性能问题，特别是使用!!，以及由于遍历而仍然运行O(n²)操作的事实，正如@leftaroundabout所提到的那样，O(n·(n–1)) ≡ O(n²)，但每次迭代应该调用gravitaten*(n-1)/2次。