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 #?

¿Fue útil?

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
Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top