使用不可变数据结构变换数据

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根本不提供突变是一个常见的神话.实际上,它提供了一种非常特殊的突变:一个值可以改变一次,从未评估到评估.利用这种特殊类型的突变的艺术被称为打结.我们将从像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.对于现时,我们将要求computeV辜负你所说,即非正式的合同,它可以看看position和neighbors任何领域Point,而不是v或w领域.(同样地,computeW可以看看除了w任何Point它的领域之外的任何东西.)实际上可以在没有太多体操的情况下在类型级别强制执行此操作,但是现在让我们跳过它.

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

编辑:...在开头解释"特殊类型的变异"评论:在评估过程中,你可以期待当你要求w一个Point事物将按此顺序发生时:computeW将迫使该v领域; 然后computeV将迫使该neighbors领域; 然后该neighbors领域将从未评估变为评估; 然后该v领域将从未评估变为评估; 然后该w字段将从未评估变为评估.最后三个步骤看起来非常类似于C++算法的三个变异步骤!

double edit:我觉得我想看到这个东西运行,所以我用虚拟实现实例化了上面提到的所有东西.我也想看到它只评估一次,因为我甚至不确定我做得对!所以我打了一些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

如您所见,每个字段最多只计算一次.