从命令式到函数式的非平凡算法转换

Question

为了加强我的（弱）函数式编程技能，我正在学习Piegl和Tiller的 NURBS 书，将所有算法转换为 Haskell。这是一个非常好的和有指导意义的过程，但我被困在算法 2.2 上，这是（一个类似 C 的由我重新设计的版本）伪代码：

double[] BasisFuns( int i, double u, int p, double U[]) {
  double N[p+1];
  double left[p+1];
  double ritght[p+1];
  N[0]=1.0;
  for (j=1; j<=p; j++) {
    left[j]  = u - U[i+1-j];
    right[j] = U[i+j] - u;
    double saved = 0.0;
    for (r=O; r<j; r++) {
      double temp= N[r]/(right[r+1]+left[j-r]);
      N[r] = saved+right[r+1]*temp;
      saved = left[j-r]*temp;
    }
    N[j] = saved;
  }
  return N;
}

外循环看起来很简单，但内循环，对元素的所有那些必须有序的修改N让我头疼。

我开始这样设置：

baseFunc :: RealFrac a => Int -> Int -> a -> [a] -> [a]
baseFunc i p u knots
  [ ??? | j <- [1..p], r <- [0..j] ]
  where
  left = [ u - knots !! (i+1-j) | j <- [ 1 .. p ]
  right= [ knots !! (i+j) - u   | j <- [ 1 .. p ]

但我觉得我可能完全偏离了道路。

我已经根据方程编写了一个完全不同且效率低下的函数版本。2.5，因此在这里我希望保持命令式版本的性能。

Answer 1

虽然从 Fortran/Matlab/C“一切都是数组”风格开始，当然可以将数值算法转换成 Haskell（并且，使用未装箱的可变向量，性能通常不会更差），但这确实缺少使用函数式语言。在基础数学实际上是更接近的功能，而不是命令式编程，所以最好的办法是开始就在这里。具体来说，递推公式

几乎可以按字面意思翻译成 Haskell，比翻译成命令式语言要好得多：

baseFuncs :: [Double]  -- ^ Knots, \(\{u_i\}_i\)
          -> Int       -- ^ Index \(i\) at which to evaluate
          -> Int       -- ^ Spline degree \(p\)
          -> Double    -- ^ Position \(u\) at which to evaluate
          -> Double
baseFuncs us i 0 u
  | u >= us!!i, u < us!!(i+1)  = 1
  | otherwise                  = 0
baseFuncs us i p u
      = (u - us!!i)/(us!!(i+p) - us!!i) * baseFuncs us i (p-1) u
       + (us!!(i+p+1) - u)/(us!!(i+p+1) - us!!(i+1)) * baseFuncs us (i+1) (p-1) u

不幸的是，由于多种原因，这实际上效率不高。

首先，列表在随机访问方面很糟糕。一个简单的解决方法是切换到未装箱（但纯）的向量。虽然我们现在是让我们将它们包装在一个 newtype 中，因为u _i应该严格增加。谈类型：直接访问是不安全的；我们可以通过将p段数带到类型级别并只允许索引来解决这个问题i < n-p，但我不会在这里讨论。

此外，在递归过程中一直传递us和u循环很尴尬，最好只绑定一次，然后使用辅助函数向下传递：

import Data.Vector.Unboxed (Vector, (!))
import qualified Data.Vector.Unboxed as VU

newtype Knots = Knots {getIncreasingKnotsSeq :: Vector Double}

baseFuncs :: Knots     -- ^ \(\{u_i\}_i\)
          -> Int       -- ^ Index \(i\) at which to evaluate
          -> Int       -- ^ Spline degree \(p\)
          -> Double    -- ^ Position \(u\) at which to evaluate
          -> Double
baseFuncs (Knots us) i? p? u = go i? p?
 where go i 0
        | u >= us!i
        , i>=VU.length us-1 || u < us!(i+1)  = 1
        | otherwise                          = 0
       go i p
           = (u - us!i)/(us!(i+p) - us!i) * go i (p-1)
            + (us!(i+p+1) - u)/(us!(i+p+1) - us!(i+1)) * go (i+1) (p-1)

另一件不是最优的事情是我们不共享相邻递归调用之间的较低级别的评估。（评估有效地跨越了具有^{p 2} ? _{2 个}节点的有向图，但我们将其评估为具有 2 ^{p 个}节点的树。）对于大p 来说这是一个巨大的低效率，但实际上对于典型的低阶样条是无害的。

避免这种低效率的方法是记忆。C 版本使用N数组明确地做到了这一点，但是——这是 Haskell——我们可以懒惰地通过使用通用的记忆库来节省分配正确大小的工作，例如memo-trie：

import Data.MemoTrie (memo2)

baseFuncs (Knots us) i? p? u = go' i? p?
 where go i 0
        | u >= us!i
        , i>=VU.length us || u < us!(i+1)  = 1
        | otherwise                        = 0
       go i p
           = (u - us!i)/(us!(i+p) - us!i) * go' i (p-1)
            + (us!(i+p+1) - u)/(us!(i+p+1) - us!(i+1)) * go' (i+1) (p-1)
       go' = memo2 go

那是无脑版本（“只需记住 ”的整个域go）。正如 dfeuer 所说，显式地只记住实际被评估的 regian 是很容易的，然后我们可以再次使用一个有效的未装箱向量：

baseFuncs (Knots us) i? p? u = VU.unsafeHead $ gol i? p?
 where gol i 0 = VU.generate (p?+1) $ \j ->
        if u >= us!(i+j)
            && (i+j>=VU.length us || u < us!(i+j+1))
         then 1 else 0
       gol i p = case gol i (p-1) of
        res' -> VU.izipWith
         (\j l r -> let i' = i+j
            in (u - us!i')/(us!(i'+p) - us!i') * l
            + (us!(i'+p+1) - u)/(us!(i'+p+1) - us!(i'+1)) * r)
         res' (VU.unsafeTail res')

（我可以安全地使用unsafeHeadandunsafeTail在这里，因为在每个递归级别，压缩将长度减少 1，所以在顶层我仍然有p? - (p?-1) = 1元素。）

我认为，这个版本应该与 C 版本具有相同的渐近线。通过一些更小的改进，例如预先计算间隔长度和预先检查参数是否在允许的范围内，以便可以进行所有访问unsafe，它的性能可能与 C 版本非常接近。

正如 dfeuer 所说的那样，甚至可能没有必要在那里使用向量，因为我只是将结果压缩在一起。对于这类东西，即使使用普通列表，GHC 也可以非常擅长优化代码。但是，我不会在这里进一步调查性能。

我用来确认它确实有效的测试：

https://gist.github.com/leftaroundabout/4fd6ef8642029607e1b222783b9d1c1e