5

只是在搞乱 F#，我试图基于这个 C# 版本创建一个基本的拉格朗日插值函数（从 C++ wiki 条目复制）：

    double Lagrange(double[] pos, double[] val, double desiredPos)
    {
        double retVal = 0;

        for (int i = 0; i < val.Length; ++i)
        {
            double weight = 1;

            for (int j = 0; j < val.Length; ++j)
            {
                // The i-th term has to be skipped
                if (j != i)
                {
                    weight *= (desiredPos - pos[j]) / (pos[i] - pos[j]);
                }
            }

            retVal += weight * val[i];
        }

        return retVal;
    }

使用我有限的 F# 和函数式编程知识，我能想到的最好的方法是：

let rec GetWeight desiredPos i j (pos : float[]) weight = 
   match i with
   | i when j = pos.Length -> weight
   | i when i = j -> GetWeight desiredPos i (j+1) pos weight 
   | i -> GetWeight desiredPos i (j+1) pos (weight * (desiredPos - pos.[j])/(pos.[i] - pos.[j]) ) 

let rec Lagrange (pos : float[]) (vals : float[]) desiredPos result counter = 
   match counter with
   | counter when counter = pos.Length -> result
   | counter -> Lagrange pos vals desiredPos (result + (GetWeight desiredPos counter 0 pos 1.0)* vals.[counter]) (counter+1)

有人可以基于相同的C# 代码提供更好/更整洁的 F# 版本吗？

4

7 回答 7

4

折叠序列是用累加器替换循环的常用方法。

let Lagrange(pos:_[], v:_[], desiredPos) =
  seq {0 .. v.Length-1} 
  |> Seq.fold (fun retVal i -> 
      seq {for j in 0 .. pos.Length-1 do if i <> j then yield j} 
      |> Seq.fold (fun w j -> w * (desiredPos - pos.[j]) / (pos.[i] - pos.[j])) 1.0
      |> (fun weight -> weight * v.[i] + retVal)) 0.0
于 2009-10-21T18:53:48.537 回答
4

使您的功能解决方案变得丑陋的部分是跳过第 i 个元素，这意味着索引。将其提取到可重用函数中，以便隔离所有丑陋的索引处理。我称我为 RoundRobin。

let RoundRobin l = seq {
  for i in {0..Seq.length l - 1} do
    yield (Seq.nth i l, Seq.take i l |> Seq.append <| Seq.skip (i+1) l)
}

不过，如果你想制作一个高效的版本，它可能会更丑陋。

我在 Seq 模块中找不到product，所以我自己写了。

let prod (l : seq<float>) = Seq.reduce (*) l

现在生成代码相当简单：

let Lagrange pos value desiredPos = Seq.sum (seq {
  for (v,(p,rest)) in Seq.zip value (RoundRobin pos) do
    yield v * prod (seq { for p' in rest do yield (desiredPos - p') / (p - p') })
})

RoundRobin 确保 pos[i] 不包含在内循环中的其余 pos 中。为了包含val数组，我用循环pos数组压缩了它。

这里的教训是，索引在函数式风格中非常难看。我还发现了一个很酷的技巧：|> Seq.append <|为您提供附加序列的中缀语法。没有那么好^。

于 2009-10-21T19:19:58.490 回答
3

我认为这作为命令式代码可以正常工作：

let LagrangeI(pos:_[], v:_[], desiredPos) =
    let mutable retVal = 0.0
    for i in 0..v.Length-1 do
        let mutable weight = 1.0
        for j in 0..pos.Length-1 do
            // The i-th term has to be skipped
            if j <> i then
                weight <- weight * (desiredPos - pos.[j]) / (pos.[i] - pos.[j])
        retVal <- retVal + weight * v.[i]
    retVal

但是如果你想要功能，一些折叠（以及 mapi 因为你经常需要携带索引）效果很好：

let LagrangeF(pos:_[], v:_[], desiredPos) =
    v |> Seq.mapi (fun i x -> i, x)
      |> Seq.fold (fun retVal (i,vi) ->
        let weight = 
            pos |> Seq.mapi (fun j x -> j<>i, x) 
                |> Seq.fold (fun weight (ok, posj) ->
                    if ok then
                        weight * (desiredPos - posj) / (pos.[i] - posj)
                    else
                        weight) 1.0
        retVal + weight * vi) 0.0

我不知道这里的数学，所以我使用了一些随机值来测试（希望）确保我没有搞砸：

let pos = [| 1.0; 2.0; 3.0 |]
let v = [|8.0; 4.0; 9.0 |]

printfn "%f" (LagrangeI(pos, v, 2.5))  // 5.375
printfn "%f" (LagrangeF(pos, v, 2.5))  // 5.375
于 2009-10-21T17:53:29.470 回答
2

这是一个非递归的解决方案。这有点古怪，因为该算法需要索引，但希望它显示了 F# 的函数是如何组合的：

let Lagrange (pos : float[]) (vals : float[]) desiredPos = 
    let weight pos desiredPos (i,v) =
        let w = pos |> Array.mapi (fun j p -> j,p)
                    |> Array.filter (fun (j,p) -> i <> j)
                    |> Array.fold (fun acc (j,p) -> acc * (desiredPos - p)/(pos.[i] - p)) 1.
        w * v
    vals |> Array.mapi (fun i v -> i,v)
         |> Array.sumBy (weight pos desiredPos)
于 2009-10-21T14:57:10.330 回答
1 
            let rec GetWeight desiredPos i j (pos : float[]) weight = 
               if j = pos.Length then weight
               elif i = j then GetWeight desiredPos i (j+1) pos weight 
               else GetWeight desiredPos i (j+1) pos (weight * (desiredPos - pos.[j])/(pos.[i] - pos.[j]) ) 

            let rec Lagrange (pos : float[]) (vals : float[]) desiredPos result counter = 
               if counter = pos.Length then result
               else Lagrange pos vals desiredPos (result + (GetWeight desiredPos counter 0 pos 1.0)* vals.[counter]) (counter+1)

就我个人而言，我认为简单的 if/elif/else 构造在这里看起来要好得多，没有这样的开销

match i with   
|i when i=...
于 2009-10-21T14:26:05.683 回答
1

如果您只是在搞乱，那么这里有一个类似于 Brian 的版本，它使用函数柯里化和元组管道运算符。

let Lagrange(pos:_[], v:_[], desiredPos) =
    let foldi f state = Seq.mapi (fun i x -> i, x) >> Seq.fold f state
    (0.0, v) ||> foldi (fun retVal (i, posi) -> 
        (1.0, pos) ||> foldi (fun weight (j, posj) -> 
            if j <> i then
                (desiredPos - posj) / (posi - posj)
            else
                1.0)
        |> (fun weight -> weight * posi + retVal))
于 2009-10-28T05:43:22.663 回答
0

我的尝试：

let Lagrange(p:_[], v, desiredPos) =
    let Seq_multiply = Seq.fold (*) 1.0
    let distance i j = if (i=j) then 1.0 else (desiredPos-p.[j])/(p.[i]-p.[j])
    let weight i = p |> Seq.mapi (fun j _ -> distance i j) |> Seq_multiply
    v |> Seq.mapi (fun i vi -> (weight i)*vi) |> Seq.sum

通过使内部循环成为函数来重构。我们还可以通过定义一些有意义的函数来使代码更直接和“易于理解”。

此外，此重写突出显示了您的原始代码（以及所有其他变体）中的一个错误。距离函数实际上应该是：

let distance i j = if (p.[i]=p.[j]) then 1.0 else (desiredPos-p.[j])/(p.[i]-p.[j])

以避免一般的 div-by-zero 错误。这导致了一个通用且无索引的解决方案：

let Lagrange(p, v, desiredPos) =
    let distance pi pj = if (pi=pj) then 1.0 else (desiredPos-pj)/(pi-pj)
    let weight pi vi = p |> Seq.map (distance pi) |> Seq.fold (*) vi
    Seq.map2 weight p v |> Seq.sum
于 2012-03-14T22:07:57.763 回答