Rewriting C# code in F#

Question

Just messing about with F# and I was trying to create a basic Lagrange Interpolation function based on this C# version (copied from a C++ wiki entry):

    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;
    }

The best I could come up with using my limited knowledge of F# and functional programming was:

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)

Can someone provide a better/tidier F# version based on the same C# code?

Answer 1

Folding over sequences is a common way to replace loops with an accumulator.

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

Answer 2

The part that makes your functional solution ugly is skipping the i'th element, which means indices. Pull that out into a reusable function so that all the ugly index handling is isolated. I call mine 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)
}

It could be a lot uglier if you want to produce an efficient version, though.

I couldn't find product in the Seq module, so I wrote my own.

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

Now producing the code is fairly simple:

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 ensures that pos[i] is not included with the rest of pos in the inner loop. To include the val array, I zipped it with the round-robinned pos array.

The lesson here is that indexing is very ugly in a functional style. Also I discovered a cool trick: |> Seq.append <| gives you infix syntax for appending sequences. Not quite as nice as ^ though.

Answer 3

I think this works fine as imperative code:

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

but if you want functional, some folds (along with mapi since you often need to carry the indices along) work well:

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

I don't know the math here, so I used some random values to test to (hopefully) ensure I screwed nothing up:

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