工程中有句谚语说“在至少有 3 个实现之前不要一概而论”。这是有道理的——尤其是在在步行几次之前寻找新的功能迭代概念时。
“步行”在这里的意思是,如果没有你知道的友好辅助函数,你应该使用递归。递归地编写您的“特殊情况”。最好采用尾递归形式。然后,如果您开始看到循环模式,您可能会想出一种方法来重构某些循环迭代方案及其“内核”。
为了澄清上述内容,假设您从未听说过foldl,并且您想从列表上的迭代中累积结果...然后,您可以编写如下内容:
myAvg values =
total / (length values)
where
mySum acc [] = acc
mySum acc (x:xs) = mySum (acc + x) xs
total = mySum 0 values
在这样做几次之后,模式可能会显示,那些 where 子句中的递归总是看起来非常相似。然后,您可能会为该内部递归 sn-p 想出一个类似“fold”或“reduce”的名称,并最终得到:
myAvg values = (foldl (+) 0.0 values) / fromIntegral (length values) :: Float
所以,如果您正在寻找对您的用例有帮助的辅助函数,我的建议是您首先编写一些实例作为递归函数,然后寻找模式。
所以,说了这么多,让我们动动手指,看看 Jacobi 算法如何 转化为 Haskell。就这样我们有话要说。现在 - 通常我不会将 Haskell 用于任何需要数组(具有 O(1) 元素访问权限的容器)的东西,因为我知道至少有 5 个数组包,我必须阅读 2 天来决定哪个适合我的应用程序。 TL;博士;)。所以我坚持使用列表,并且在下面的代码中除了序曲之外没有包依赖项。但那是 - 鉴于我们试图解决的示例方程的大小很小 - 这根本不是一件坏事。此外,代码表明,lazy Haskell 中的列表推导允许对单元集(例如在矩阵中)进行非强制性但高性能的操作,而无需显式循环。
type Matrix = [[Double]]
-- sorry - my mind went blank while looking for a better name for this...
-- but it is useful nonetheless
idefix nr nc =
[ [(r,c) | c <- [0..nc-1]] | r <- [0..nr-1]]
matElem m (r,c) = (m !! r) !! c
transpose (r,c) = (c,r)
matrixDim m = (length m, length . head $ m)
-- constructs a Matrix by enumerating the indices and querying
-- 'unfolder' for a value.
-- try "unfoldMatrix 3 3 id" and you see how indices relate to
-- cells in the matrix.
unfoldMatrix nr nc unfolder =
fmap (\row -> fmap (\cell -> unfolder cell) row) $ idefix nr nc
-- Not really needed for Jacobi problem but good
-- training to get our fingers wet with unfoldMatrix.
transposeMatrix m =
let (nr,nc) = matrixDim m in
unfoldMatrix nc nr (matElem m . transpose)
addMatrix m1 m2
| (matrixDim m1) == (matrixDim m2) =
let (nr,nc) = matrixDim m1 in
unfoldMatrix nr nc (\idx -> matElem m1 idx + matElem m2 idx)
subMatrix m1 m2
| (matrixDim m1) == (matrixDim m2) =
let (nr,nc) = matrixDim m1 in
unfoldMatrix nr nc (\idx -> matElem m1 idx - matElem m2 idx)
dluMatrix :: Matrix -> (Matrix,Matrix,Matrix)
dluMatrix m
| (fst . matrixDim $ m) == (snd . matrixDim $ m) =
let n = fst . matrixDim $ m in
(unfoldMatrix n n (\(r,c) -> if r == c then matElem m (r,c) else 0.0)
,unfoldMatrix n n (\(r,c) -> if r > c then matElem m (r,c) else 0.0)
,unfoldMatrix n n (\(r,c) -> if c > r then matElem m (r,c) else 0.0)
)
mulMatrix m1 m2
| (snd . matrixDim $ m1) == (fst . matrixDim $ m2) =
let (nr, nc) = ((fst . matrixDim $ m1),(snd . matrixDim $ m2)) in
unfoldMatrix nr nc
(\(ro,co) ->
sum [ matElem m1 (ro,i) * matElem m2 (i,co) | i <- [0..nr-1]]
)
isSquareMatrix m = let (nr,nc) = matrixDim m in nr == nc
jacobi :: Double -> Matrix -> Matrix -> Matrix -> Matrix
jacobi errMax a b x0
| isSquareMatrix a && (snd . matrixDim $ a) == (fst . matrixDim $ b) =
approximate x0
-- We could possibly avoid our hand rolled recursion
-- with the help of 'loop' from Control.Monad.Extra
-- according to hoogle. But it would not look better at all.
-- loop (\x -> let x' = jacobiStep x in if converged x' then Right x' else Left x') x0
where
(nra, nca) = matrixDim a
(d,l,u) = dluMatrix a
dinv = unfoldMatrix nra nca (\(r,c) ->
if r == c
then 1.0 / matElem d (r,c)
else 0.0)
lu = addMatrix l u
converged x =
let delta = (subMatrix (mulMatrix a x) b) in
let (nrd,ncd) = matrixDim delta in
let err = sum (fmap (\idx -> let v = matElem delta idx in v * v)
(concat (idefix nrd ncd))) in
err < errMax
jacobiStep x =
(mulMatrix dinv (subMatrix b (mulMatrix lu x)))
approximate x =
let x' = jacobiStep x in
if converged x' then x' else approximate x'
wikiExample errMax =
let a = [[ 2.0, 1.0],[5.0,7.0]] in
let b = [[11], [13]] in
jacobi errMax a b [[1.0],[1.0]]
函数idefix,尽管它的名字很傻,但恕我直言,对于来自非懒惰语言的人来说,这是一个大开眼界。他们的第一反应是害怕:“什么——他用索引创建了一个列表而不是写循环?真是浪费!”但是浪费,它不是在懒惰的语言中。你在这个函数中看到的(列表推导)会产生一个 lazy 列表。它不是真正创建的。幕后发生的事情在本质上类似于 LINQ 在 C# 中所做的事情 - IEnumerator<T> juggling。
当我们想要对delta 中的所有元素求和时,我们第二次使用idefix。在那里,我们不关心矩阵的具体结构。因此,我们使用标准的前奏函数concat 将矩阵展平为线性列表。当然,也很懒惰。这就是美。
命令式维基百科伪代码的下一个显着区别是,与嵌套循环和在单个单元格上操作相比,使用矩阵表示法要简单得多。幸运的是,维基百科文章显示了两者。因此,我们只需要一个等效于最外层的 while 循环,而不是具有 2 个嵌套循环的 while 循环。我们的 2 线性递归函数 approximate 涵盖了这一点。
经验教训:
- 列表和列表推导可以帮助简化代码,否则需要嵌套循环。 (在惰性语言中)。
- Ocaml 和 Common Lisp 具有可变性并内置数组和循环。当从命令式语言或命令式伪代码翻译算法时,这使得一个包非常方便。
- Haskell 具有不变性,没有内置数组,也没有循环,但它有一套同样强大的工具,即惰性、尾调用优化和简洁的语法。这种组合需要更多的计划(并编写一些通常很短的辅助函数),而不是“让我们全部在 main() 中编写”的经典 C 方法。
- 有时,编写一个 2 行长的递归函数比考虑如何抽象它更容易。