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

Question

为了加强我（薄弱的）函数式编程技能，我正在研究 Piegl 和 Tiller 将所有算法转换为 Haskell 的 NURBS 书 .这是一个非常好的和有启发性的过程，但我被算法 2.2 卡住了，这里是伪代码（一个类似 C 的 reworked-by-me 版本）：

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 ]

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

我已经基于 Eq. 2.5，因此在这里我希望保持命令式版本的性能。

Answer 1

我不确定它是否会起作用，但是...：

baseFunc :: RealFrac a => Int -> Int -> a -> [a] -> [a]
baseFunc i p u knots =
    foldl' helper [1.0] [1..p]

    where
    left = [ u - knots !! (i+1-j) | j <- [ 1 .. p ] ]
    right= [ knots !! (i+j) - u   | j <- [ 1 .. p ] ]

    helper  N  j = outer_loop j N left right

    inner_loop :: RealFrac a => Int -> Int -> [a] -> [a] -> [a] -> a -> (a, a)
    inner_loop  r  j  N  left  right  saved =
        let temp = N !! r / (right !! (r+1) + left !! (j-r))
        in  (saved + right !! (r+1) * temp, left !! (j-r) * temp)

    outer_loop :: RealFrac a => Int -> [a] -> [a] -> [a] -> [a]
    outer_loop  j  N  left  right =
        let (new_N, saved) = foldl' helper (N, 0.0) [0..j-1]
            helper  (prev_N, saved)  r = 
                let (N_r, new_saved) = inner_loop r j prev_N left right saved
                in  (insertAt r N_r prev_N, new_saved)
        in  new_N ++ [saved]

Answer 2

虽然从 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

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

首先，列表不适合随机访问。一个简单的解决方法是切换到未装箱（但纯）的向量。当我们在这里时，让我们将它们包装在一个新类型中，因为 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)

另一件不是最优的事情是我们不会在相邻的递归调用之间共享较低级别的评估。 （评估实际上跨越了一个具有 ^p2⁄₂ 节点的有向图，但我们将其评估为具有 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 所说，很容易只显式记住实际被评估的区域，然后我们可以再次使用有效的未装箱向量：

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')

（我可以在这里安全地使用unsafeHead 和unsafeTail，因为在每个递归级别，压缩都会将长度减少1，所以在顶层我还剩下p₀ - (p₀-1) = 1 元素。）

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

正如 - 再次 - dfeuer 所说，甚至可能不需要在此处使用向量，因为我只是将结果压缩在一起。对于这类东西，即使使用普通列表，GHC 也可以非常擅长优化代码。但是，我不会在这里进一步研究性能。

---

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

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

Answer 3

（免责声明：我对这里计算的内容一无所知。）

数组U 中的访问模式似乎如下：从索引i 向外，我们消费左边和右边的值。我们可以想象有两个列表精确地包含在这些元素序列中。事实上，我们可以像这样从源列表中构造出这样的列表：

pryAt :: Int -> [a] -> ([a], [a])
pryAt i xs = go i ([], xs)
  where
    go 0 a = a
    go n (us, v : vs) = go (pred n) (v : us, vs)
-- pryAt 5 ['a'..'x']
-- ("edcba","fghijklmnopqrstuvwx")

对于随机访问容器，我们可以有专门的 pryAt 版本，因为遍历整个列表直到到达 ith 元素将是低效的。

在外部循环中，我们有数组N、left 和right，它们随着每次迭代而增长（N 似乎在每次迭代时也完全重建）。我们可以将它们表示为列表。在每次迭代中，我们还消耗了左右一对U的元素。

以下数据类型表示外循环迭代开始时的情况：

data State = State
  { ns :: [Double],
    left :: [Double],
    right :: [Double],
    us :: ([Double], [Double])
  }

假设我们已经实现了outerStep :: State -> State，我们可以简单地转动曲柄p 次：

basis :: Int -> Double -> Int -> [Double] -> [Double]
basis i u p us =
  ns $ iterate outerStep initial !! p
  where
    initial =
      State
        { ns = [1.0],
          left = [],
          right = [],
          us = pryAt i us
        }

outerStep 做了什么？我们向left 和right 添加新元素，然后从头开始重新创建整个N 列表，同时携带saved 累加器。这是mapAccumR。我们需要一些额外的信息：right 值（与N 方向相同）和left 值（反向）所以我们需要预先压缩它们：

    outerStep (State {ns, left, right, us = (ul : uls, ur : urs)}) =
      let left' = u - ul : left
          right' = ur - u : right
          (saved', ns') = mapAccumR innerStep 0.0 $ zip3 ns right' (reverse left')
       in State
            { ns = saved' : ns',
              left = left',
              right = right',
              us = (uls, urs)
            }

下面是内部步骤的计算：

    innerStep saved (n, r, l) = 
        let temp = n / (r - l)
            n' = saved + r
            saved' = l * temp
         in (saved', n')

除了纠正可能的错误外，还需要做更多的工作，因为当前形式的basis 函数可能会泄漏内存（特别是mapAccumR 会产生大量的thunk）。也许可以重写它以使用像 iterate' 或 foldl' 这样的函数来保持其累加器严格。