haskell - 使用不可变数据结构改变数据

Question

我想实现一个特定的算法，但我找不到适合这项工作的良好数据结构。该算法的一个更简单版本的工作原理如下:

Input: A set of points.
Output: A new set of points.
Step 1: For each point, calculate the closest points in a radius.
Step 2: For each point, calculate a value "v" from the closest points subset.
Step 3: For each point, calculate a new value "w" from the closest points and
        the values "v" from the previous step, i.e, "w" depends on the neighbors
        and "v" of each neighbor.
Step 4: Update points.

在 C++ 中，我可以这样解决:

struct Point {
    Vector position;
    double v, w;
    std::vector<Point *> neighbors;
};

std::vector<Point> points = initializePoints();
calculateNeighbors(points);
calculateV(points); // points[0].v = value; for example.
calculateW(points);

对于诸如点列表之类的朴素结构，我无法将值“v”更新为原始点集，并且需要计算邻居两次。我怎样才能避免这种情况并保持函数的纯净，因为计算邻居是算法中最昂贵的部分(超过 30% 的时间)？

PS.:对于那些在数值方法和 CFD 方面有经验的人来说，这是平滑粒子流体动力学方法的简化版本。

更新:更改了第 3 步，使其更加清晰。

Answer 1

Haskell 根本不提供突变是一个普遍的误解。实际上，它提供了一种非常特殊的变异:一个值只能变异一次，从未求值到求值。利用这种特殊突变的艺术被称为 tying the knot .我们将从一个数据结构开始，就像您的 C++ 数据结构一样:

data Vector -- held abstract

data Point = Point
    { position  :: Vector
    , v, w      :: Double
    , neighbors :: [Point]
    }

现在，我们要做的是构建一个Array Point，其neighbors 包含指向同一数组中其他元素的指针。以下代码中 Array 的主要特点是它是惰性的(它不会过早地强制其元素)并且具有快速的随机访问；如果愿意，您可以用这些属性替换您最喜欢的备用数据结构。

找邻居功能的界面有多种选择。为了具体和简化我自己的工作，我假设你有一个函数，它接受一个 Vector 和一个 Vectors 列表，并给出邻居的索引。

findNeighbors :: Vector -> [Vector] -> [Int]
findNeighbors = undefined

我们还为 computeV 和 computeW 设置一些类型。对于 nonce，我们将要求 computeV 遵守您声明的非正式契约(Contract)，即它可以查看 position 和 neighbors任何 Point 的字段，但不是 v 或 w 字段。 (类似地，computeW 可以查看任何东西，但可以查看任何 Point 的 w 字段。)实际上可以强制执行这在类型级别没有太多体操，但现在让我们跳过它。

computeV, computeW :: Point -> Double
(computeV, computeW) = undefined

现在我们已准备好构建我们的(标记的)内存中图。

buildGraph :: [Vector] -> Array Int Point
buildGraph vs = answer where
    answer = listArray (0, length vs-1) [point pos | pos <- vs]
    point pos = this where
        this = Point
            { position = pos
            , v = computeV this
            , w = computeW this
            , neighbors = map (answer!) (findNeighbors pos vs)
            }

就是这样，真的。现在你可以写你的

newPositions :: Point -> [Vector]
newPositions = undefined

newPositions 可以完全自由地检查它所传递的 Point 的任何字段，并将所有函数放在一起:

update :: [Vector] -> [Vector]
update = newPositions <=< elems . buildGraph

编辑:...解释开头的“特殊类型的突变”评论:在评估期间，您可以期望何时需要 Point 的 w 字段> 事情将按以下顺序发生:computeW 将强制 v 字段；然后 computeV 将强制 neighbors 字段；然后 neighbors 字段将从未评估变为已评估；然后 v 字段将从未评估变为已评估；然后 w 字段将从未评估变为已评估。这最后三个步骤看起来与您的 C++ 算法的三个变异步骤非常相似!

双重编辑:我决定让这个东西运行起来，所以我用虚拟实现实例化了上面所有抽象的东西。我还希望看到它只对事物进行一次评估，因为我什至不确定我是否做对了!所以我加入了一些 trace 调用。这是一个完整的文件:

import Control.Monad
import Data.Array
import Debug.Trace

announce s (Vector pos) = trace $ "computing " ++ s ++ " for position " ++ show pos

data Vector = Vector Double deriving Show

data Point = Point
    { position  :: Vector
    , v, w      :: Double
    , neighbors :: [Point]
    }

findNeighbors :: Vector -> [Vector] -> [Int]
findNeighbors (Vector n) vs = [i | (i, Vector n') <- zip [0..] vs, abs (n - n') < 1]

computeV, computeW :: Point -> Double
computeV (Point pos _ _ neighbors) = sum [n | Point { position = Vector n } <- neighbors]
computeW (Point pos v _ neighbors) = sum [v | Point { v = v } <- neighbors]

buildGraph :: [Vector] -> Array Int Point
buildGraph vs = answer where
    answer = listArray (0, length vs-1) [point pos | pos <- vs]
    point pos = this where { this = Point
        { position  = announce "position" pos $ pos
        , v         = announce "v" pos $ computeV this
        , w         = announce "w" pos $ computeW this
        , neighbors = announce "neighbors" pos $ map (answer!) (findNeighbors pos vs)
        } }

newPositions :: Point -> [Vector]
newPositions (Point { position = Vector n, v = v, w = w }) = [Vector (n*v), Vector w]

update :: [Vector] -> [Vector]
update = newPositions <=< elems . buildGraph

并在 ghci 中运行:

*Main> length . show . update . map Vector $ [0, 0.25, 0.75, 1.25, 35]
computing position for position 0.0
computing v for position 0.0
computing neighbors for position 0.0
computing position for position 0.25
computing position for position 0.75
computing w for position 0.0
computing v for position 0.25
computing neighbors for position 0.25
computing v for position 0.75
computing neighbors for position 0.75
computing position for position 1.25
computing w for position 0.25
computing w for position 0.75
computing v for position 1.25
computing neighbors for position 1.25
computing w for position 1.25
computing position for position 35.0
computing v for position 35.0
computing neighbors for position 35.0
computing w for position 35.0
123

如您所见，每个字段最多为每个位置计算一次。