我想实现一个特定的算法,但我找不到适合这项工作的良好数据结构。该算法的一个更简单版本的工作原理如下:
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 步,使其更加清晰。
最佳答案
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
如您所见,每个字段最多为每个位置计算一次。
关于haskell - 使用不可变数据结构改变数据,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/19176849/