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 番目の要素をスキップすることです。それを再利用可能な関数に引き出して、醜いインデックス処理がすべて分離されるようにします。私はラウンドロビンと呼んでいます。

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

あなたがただいじっているなら、関数カリー化とタプルパイプ演算子を使用する 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 に答える
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 に答える
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 に答える