Question

Je me contentais de jouer avec F # et j'essayais de créer une fonction d'interpolation Lagrange de base basée sur cette version C # (copiée à partir d'une entrée de 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;
    }

Le mieux que je pouvais trouver en utilisant mes connaissances limitées en F # et en programmation fonctionnelle était:

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)

Quelqu'un peut-il fournir une version F # meilleure / plus ordonnée basée sur le même code C #?

Était-ce utile?

La solution

Le repliement de séquences est un moyen courant de remplacer les boucles par un accumulateur.

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

Autres conseils

La partie qui rend votre solution fonctionnelle laide est de sauter le ième élément, ce qui signifie les indices. Tirez cela dans une fonction réutilisable afin que toute la gestion des index laids soit isolée. J'appelle le mien 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)
}

Cela pourrait être beaucoup plus laid si vous voulez produire une version efficace, cependant.

Je ne pouvais pas trouver product dans le module Seq, j'ai donc écrit le mien.

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

Produire le code est assez 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 garantit que pos [i] n'est pas inclus avec le reste de pos dans la boucle interne. Pour inclure le tableau val, je l'ai compressé avec le tableau pos circulaire.

La leçon à tirer est que l’indexation est très laide dans un style fonctionnel. J'ai aussi découvert un truc intéressant: |> Seq.append <| vous donne la syntaxe infixe pour ajouter des séquences. Pas tout à fait aussi bien que ^ bien.

Je pense que cela fonctionne bien en tant que code impératif:

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

mais si vous voulez être fonctionnel, certains plis (ainsi que mapi, car vous avez souvent besoin de transporter les index) fonctionnent 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

Je ne connais pas les calculs ici, alors j’ai utilisé des valeurs aléatoires pour tester (espérons-le) s’assurer de ne rien foirer:

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

Voici une solution non récursive. C'est un peu génial parce que l'algorithme nécessite des index, mais j'espère qu'il montre comment les fonctions de F # peuvent être composées:

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)

Personnellement, je pense que les constructions simples si / elif / else regardent ici beaucoup mieux sans frais généraux tels que

match i with   
|i when i=...

Si vous vous contentez de déconner, voici une version similaire à celle de Brian qui utilise la fonction currying et l'opérateur de tubes de tuples.

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

Ma tentative:

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

Refactoriser en faisant de la boucle interne une fonction. Nous pouvons également rendre le code plus simple et & "Compréhensible &"; en définissant des fonctions significatives.

De plus, cette réécriture met en évidence un bogue dans votre code original (et toutes les autres variantes). La fonction de distance devrait en réalité être:

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

pour éviter l'erreur générale div par zéro. Cela conduit à une solution générique et sans index:

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
Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top