分享哈斯克尔的力量计算

Question

我正在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,然后更新各自Particle的ps加速度为我做到这一点.

我想过要做点什么foldr f ps combs.起始累加器将是当前的ps并且f将是一些函数,它接收每个comb并更新相关Particle的内容ps,并返回该累加器.对于这样一个简单的过程来说,这似乎非常耗费内存并且非常复

有任何想法吗？

Answer 1

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

\n\n

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

\n\n

并修改它，以便在矩阵（好吧，列表列表......）中计算出加速度，与更新每个粒子分开，我们得到......

\n\n

updateParticle :: Model -> (Particle, [Acc]) -> Particle\nupdateParticle ps (p@(Particle m pos vel acc), accs) =\n    let acc\' = foldr (\\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs\n        vel\' = (fst vel + fst acc, snd vel + snd acc)\n        pos\' = (fst pos + fst vel, snd pos + snd vel)\n    in  Particle m pos\' vel\' acc\'\n\nstep :: ViewPort -> Float -> Model -> Model\nstep _ _ ps = map (updateParticle ps) $ zip ps accsMatrix where\n    accsMatrix = [map (gravitate p) ps | p <- ps]\n

\n\n

...所以问题本质上是如何利用=的事实来减少gravitate锻炼时调用的次数。accsMatrixgravitate a b-1 * gravitate b a

\n\n

如果我们打印出来accsMatrix，它看起来像......

\n\n

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

\n\n

...所以我们看到了accsMatrix !! i !! j == -1 * accsMatrix !! j !! i。

\n\n

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

\n\n

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

\n\n

...并用列表理解替换内部列表...

\n\n

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

\n\n

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

\n\n

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

\n\n

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

\n\n

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]\n

\n\n

我们也可以将其拆分一下，如下所示......

\n\n

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

\n\n

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

\n\n

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

\n\n

...或者通过使用列表 monad...

\n\n

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

\n\n

...尽管（对我来说）更难看出其中发生了什么。

\n\n

这种列表列表方法可能存在一些可怕的性能问题，尤其是使用，并且由于遍历而!!您仍在运行O(n\xc2\xb2)操作这一事实，以及 O (n \xc2 \xb7 (n \xc2\xad\xe2\x80\x93 1)) \xe2\x89\xa1 O (n\xc2\xb2) 正如@leftaroundabout提到的，但每次迭代都应该调用gravitate n * (n-1) / 2times.

\n