2009-10-21 10 views
5

Sólo perder el tiempo con F # y yo estaba tratando de crear una función básica interpolación de Lagrange sobre la base de esta versión C# (copiado de una entrada de wiki C++):Reescritura de código C# en Fa #

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

Lo mejor que podía venir con el uso de mi conocimiento limitado de F # y la programación funcional fue:

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) 

Puede alguien proporcionar una mejor versión más ordenado #/F basado en el mismo código C# ?

+3

Creo que este es un buen ejemplo de cuándo el código imperativo es más fácil de leer y mantener que funcional. – gradbot

Respuesta

4

plegable sobre secuencias es una forma común de sustituir los bucles con un acumulador.

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 
+0

Gracias por la edición de Brian. – gradbot

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) 

personalmente creo que simple si/elif/else construcciones mira aquí mucho mejor sin los gastos generales tales como

match i with 
|i when i=... 
2

Aquí hay una solución no recursiva. Es un poco raro porque el algoritmo requiere índices, pero esperemos que muestra cómo las funciones F # 's pueden estar compuestos:

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

Creo que este código funciona muy bien como imperativo:

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 

pero si quieres funcional, algunos pliegues (junto con MAPI ya que a menudo tienen que llevar a lo largo de los índices) funcionan bien:

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 

no sé las matemáticas aquí, así que utilizan algunos valores aleatorios para probar al correo (con suerte) segurar lo arruiné nada hasta:

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 
+0

También podría eliminar el mapi haciendo que su acumulador de pliegues sea una tupla que contenga el índice. v |> fold (diversión (retVal, i) posi -> newRetValuefunction, i + 1) (0.0, 0) – gradbot

+0

Esto produce las mismas respuestas que el código C# original para mis datos –

4

La parte que hace que su solución funcional sea fea omite el elemento i, lo que significa índices. Ponlo en una función reutilizable para que todo el feo manejo del índice esté aislado. Yo llamo al mío 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) 
} 

Sin embargo, podría ser mucho más feo si quiere producir una versión eficiente.

No pude encontrar product en el módulo Seq, así que escribí el mío.

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

Ahora producir el código es bastante 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 asegura que pos [i] no está incluido con el resto de pos en el bucle interior. Para incluir el conjunto de val, lo comprimí con la matriz pos de round-robinned.

La lección aquí es que la indexación es muy fea en un estilo funcional. También descubrí un truco genial: |> Seq.append <| le da sintaxis infijo para anexar secuencias. No es tan bueno como ^ sin embargo.

1

Si solo te estás metiendo en el camino, entonces esta es una versión similar a la de Brian que usa el currying de funciones y el operador de tuple pipe.

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

Mi intento:

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 

Refactor haciendo que el bucle interno de una función. También podemos hacer que el código sea más directo y "comprensible" al definir algunas funciones significativas.

Además, esta reescritura resalta un error en su código original (y todas las demás variantes). La función de distancia debería ser en realidad:

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

para evitar el error general de división por cero. Esto lleva a una solución genérica e indistinta:

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