2009-10-21 3 views
5

Просто возиться с F #, и я пытался создать основную функцию Лагранжа Интерполяция на основе этой C# версии (скопированного из записи вики C++):Переписывая C# код в F #

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

Лучшее, что я мог бы прийти с помощью моего ограниченного знания F # и функционального программирования было:

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) 

Может кто-то обеспечить лучшее/опрятнее F # версии на основе же C# код?

+3

Я думаю, что это хороший пример того, когда императивный код легче читать и поддерживать, чем функциональный. – gradbot

ответ

4

Складные последовательности - это обычный способ замены петель аккумулятором.

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

Спасибо за редактирование Брайана. – 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) 

Лично я считаю, что просто, если/Элиф/остальное конструкции смотрите здесь намного лучше без таких накладных расходов, как

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

Вот нерекурсивна решение. Это немного напуганный, потому что алгоритм требует индексов, но, надеюсь, она показывает, как функции F # 'S может состоять:

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

Я думаю, что это работает отлично, как императивный код:

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 

, но если вы хотите функциональный, некоторые складки (наряду с MAPI, так как часто необходимо проводить индексы вместе) работают хорошо:

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 

я не знаю математику, поэтому я использовал некоторые случайные значения для проверки на (надеюсь) е nsure, я ничего не испортил:

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

Вы также можете покончить с mapi, если ваш складной аккумулятор будет кортежем, содержащим индекс. v |> fold (fun (retVal, i) posi -> newRetValuefunction, i + 1) (0.0, 0) – gradbot

+0

Это дает те же ответы, что и исходный код C# для моих данных. –

4

Часть, которая делает ваше функциональное решение уродливым, пропускает i-й элемент, что означает индексы. Вытяните это в функцию многократного использования, чтобы исключить всю уродливую обработку индекса. Я называю мой 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) 
} 

Это может быть намного уродливее, если вы хотите создать эффективную версию.

Не удалось найти product в модуле Seq, поэтому я написал свой собственный.

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

Теперь производить код достаточно прост:

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 гарантирует, что позиция [я] не входит в остальных поз во внутреннем цикле. Чтобы включить массив val, я закрепил его массивным размером pos.

Урок здесь заключается в том, что индексирование очень уродливое в функциональном стиле. Также я обнаружил отличный трюк: |> Seq.append <| дает вам синтаксис infix для добавления последовательностей. Не совсем так хорошо, как ^.

1

Если вы просто возитесь, то вот версия, подобная Брайану, которая использует функцию currying и оператор 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

Моя попытка:

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 путем внутреннего контура функции. Также мы можем сделать код более понятным и понятным, определив некоторые значимые функции.

Кроме того, эта переписка подчеркивает ошибку в исходном коде (и всех других вариантах). Функция расстояния должна быть фактически:

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

, чтобы избежать общей ошибки div-by-zero. Это приводит к общему и безразмерному решению:

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 
Смежные вопросы