Reescribiendo el código C # en F #
-
05-07-2019 - |
Pregunta
Solo estaba jugando con F # e intentaba crear una función básica de interpolación de Lagrange basada en esta versión de C # (copiada de una entrada wiki de C ++):
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 pude encontrar usando 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 versión F # mejor / más ordenada basada en el mismo código C #?
Solución
Plegar secuencias es una forma común de reemplazar 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
Otros consejos
La parte que hace que su solución funcional sea fea es omitir el elemento i, lo que significa índices. Saque eso a una función reutilizable para que todo el manejo del índice feo 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 desea 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 se incluye con el resto de pos en el bucle interno. Para incluir la matriz val
, la comprimí con la matriz redondeada pos
.
La lección aquí es que la indexación es muy fea en un estilo funcional.
También descubrí un truco genial: |> Seq.append <|
te da la sintaxis infija para agregar secuencias. Sin embargo, no es tan agradable como ^
.
Creo que esto funciona bien como código 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 necesitas llevar 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 usé algunos valores aleatorios para probar (con suerte) para asegurarme de que no arruiné nada:
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
Aquí hay una solución no recursiva. Es un poco extraño porque el algoritmo requiere índices, pero con suerte muestra cómo se pueden componer las funciones de 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)
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 las construcciones simples de if / elif / else se ven aquí mucho mejor sin gastos generales como
match i with
|i when i=...
Si solo está jugando, aquí hay una versión similar a la de Brian que utiliza la función de curry y el operador de tubería de tupla.
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))
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
Refactorizar haciendo que el bucle interno sea una función. También podemos hacer que el código sea más sencillo y & Quot; comprensible & Quot; definiendo 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:
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 div-by-zero. Esto lleva a una solución genérica e indizada:
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