Domanda

Stavo solo scherzando con F # e stavo cercando di creare una funzione base di interpolazione Lagrange basata su questa versione C # (copiata da una voce wiki 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;
    }

Il meglio che ho potuto inventare usando la mia conoscenza limitata di F # e la programmazione funzionale è stato:

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)

Qualcuno può fornire una versione F # migliore / più ordinata basata sullo stesso codice C #?

È stato utile?

Soluzione

Ripiegare le sequenze è un modo comune per sostituire i loop con un accumulatore.

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

Altri suggerimenti

La parte che rende brutta la tua soluzione funzionale sta saltando l'ennesimo elemento, il che significa indici. Tiralo fuori in una funzione riutilizzabile in modo che tutta la brutta gestione dell'indice sia isolata. Chiamo il mio 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)
}

Potrebbe essere molto più brutto se si desidera produrre una versione efficiente.

Non sono riuscito a trovare product nel modulo Seq, quindi ho scritto il mio.

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

Ora produrre il codice è abbastanza semplice:

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 assicura che pos [i] non sia incluso con il resto di pos nel loop interno. Per includere l'array val, l'ho compresso con l'array pos a struttura circolare

La lezione qui è che l'indicizzazione è molto brutta in uno stile funzionale. Inoltre ho scoperto un bel trucco: |> Seq.append <| ti dà la sintassi di infisso per le sequenze aggiunte. Non abbastanza carino come ^ però.

Penso che funzioni perfettamente come codice 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

ma se vuoi funzionale, alcune pieghe (insieme a mapi poiché spesso devi portare gli indici) funzionano bene:

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

Non conosco la matematica qui, quindi ho usato alcuni valori casuali per testare (si spera) per assicurarsi di aver rovinato nulla:

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

Ecco una soluzione non ricorsiva. È un po 'strano perché l'algoritmo richiede indici, ma si spera che mostri come possono essere composte le funzioni di 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 penso che i semplici costrutti if / elif / else abbiano un aspetto molto migliore senza spese generali come

match i with   
|i when i=...

Se stai solo scherzando, ecco una versione simile a quella di Brian che utilizza la funzione curry e l'operatore 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))

Il mio tentativo:

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 rendendo il loop interno una funzione. Inoltre possiamo rendere il codice più semplice e & Quot; comprensibile & Quot; definendo alcune funzioni significative.

Inoltre, questa riscrittura evidenzia un bug nel codice originale (e in tutte le altre varianti). La funzione di distanza dovrebbe effettivamente essere:

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

per evitare l'errore generale div-by-zero. Questo porta a una soluzione generica e senza indice:

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
Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top