1

我编写了一些代码，使用 repa计算距离矩阵：

distance :: Int -> Int -> Mat -> Double
distance aindx bindx arr = let a = slice arr (Any :. aindx :. All)
                               b = slice arr (Any :. bindx :. All)-
                               sqdiff = R.map (\x -> x*x) $ R.zipWith (-) a b 
                            in sqrt $ sumAllS sqdiff

buildDistanceMatrix :: Mat -> Mat 
buildDistanceMatrix m = let (Z :. height :. width) = R.extent m
                            cords = fromListUnboxed (Z :. (height * height) )  [ (x,y) | x <- [0..height-1], y <- [0..height-1]]
                            dist = R.smap (\(a,b) -> distance a b m) cords
                            dmat = R.reshape (Z :. height  :. height ) dist
                         in R.computeS dmat

它似乎工作。但后来我添加了一个快速检查：

prop_distmat :: Double -> Bool
prop_distmat d = let dvec = [d,0,0,0]
                     dmat = R.fromListUnboxed (Z :. (2::Int) :. (2::Int)) dvec
                     dist = buildDistanceMatrix dmat
                  in (R.toList dist) == [0.0, d, d, 0.0 ]

换句话说，被距离 D 分开的两个点应该产生一个看起来像 [0,D,D,0] 的距离矩阵。在我的临时手动测试中，确实如此。但 QuickCheck 很快发现 5.0e-324 的距离会产生 [0,0,0,0] 的距离矩阵

distance matrix  *** Failed! Falsifiable (after 2 tests and 1074 shrinks):    
5.0e-324

这仅仅是因为双打的精度吗？我是否需要限制 QuickCheck 将发送的可能值？或者这是一个真正的错误？

4

2 回答 2

5

您正在测试浮点数是否相等，这通常应该避免（在任何语言中，这不是 Haskell 特定的）。大双打也会溢出到无穷大。即使对于那些没有上溢或下溢的双打，sqrt (x*x) == x一般也不成立。因此，您需要同时==检查差异最多是一些合理的 epsilon 并限制可能的值（或检查属性中的溢出）。

于 2014-08-20T05:34:51.073 回答
1

天真地计算向量的 L2 范数可以在应用平方根函数之前很久就给出欠流或溢出。我引用一位知道的人的话：“Fortran 中两个范数的最稳健计算有超过 200 行代码，但那是 25 年前的事了”。我建议搜索 Fortran 实现，然后使用可能出错的知识将其应用于您的 Haskell 实现。数字很​​棘手；好消息是，大约 50 年前，大多数问题可能都在 Fortran 中得到了解决。

于 2014-08-21T09:54:55.170 回答