Skip to content

Conversione non banale di un algoritmo da imperativo a funzionale

La guida o il codice passo-passo che vedrai in questo post è la soluzione più semplice ed efficace che troviamo alla tua domanda o problema.

Soluzione:

Sebbene sia certamente possibile tradurre algoritmi numerici in Haskell partendo dallo stile Fortran/Matlab/C "tutto è un array" (e, usando vettori mutabili non inscatolati, le prestazioni non saranno in genere molto peggiori), questo non è il punto di partenza per l'uso di un linguaggio funzionale. Il matematica sottostante è in realtà molto più vicina alla programmazione funzionale che a quella imperativa, quindi la cosa migliore è iniziare proprio da lì. In particolare, la formula di ricorrenza

Equazione 2.5 dal libro Nurbs

può essere tradotta quasi letteralmente in Haskell, molto meglio che in un linguaggio imperativo:

baseFuncs :: [Double]  -- ^ Knots, ({u_i}_i)
          -> Int       -- ^ Index (i) at which to evaluate
          -> Int       -- ^ Spline degree (p)
          -> Double    -- ^ Position (u) at which to evaluate
          -> Double
baseFuncs us i 0 u
  | u >= us!!i, u < us!!(i+1)  = 1
  | otherwise                  = 0
baseFuncs us i p u
      = (u - us!!i)/(us!!(i+p) - us!!i) * baseFuncs us i (p-1) u
       + (us!!(i+p+1) - u)/(us!!(i+p+1) - us!!(i+1)) * baseFuncs us (i+1) (p-1) u

Sfortunatamente, però, questa soluzione non sarà efficiente per diversi motivi.

Primo, le liste fanno schifo per l'accesso casuale. Una soluzione semplice è quella di passare a vettori unboxed (ma puri). Già che ci siamo, avvolgiamoli in un newtype, perché l'opzione ui si suppone che siano strettamente crescenti. Parlando di tipi: gli accessi diretti sono unsafe; si potrebbe risolvere questo problema portando p e il numero di segmenti al livello del tipo e permettendo solo gli indici i < n-p, ma non lo approfondirò in questa sede.

Inoltre, è scomodo passare il parametro us e u per tutto il percorso della ricorsione, meglio legarlo una volta e poi usare una funzione ausiliaria per scendere:

import Data.Vector.Unboxed (Vector, (!))
import qualified Data.Vector.Unboxed as VU

newtype Knots = Knots {getIncreasingKnotsSeq :: Vector Double}

baseFuncs :: Knots     -- ^ ({u_i}_i)
          -> Int       -- ^ Index (i) at which to evaluate
          -> Int       -- ^ Spline degree (p)
          -> Double    -- ^ Position (u) at which to evaluate
          -> Double
baseFuncs (Knots us) i₀ p₀ u = go i₀ p₀
 where go i 0
        | u >= us!i
        , i>=VU.length us-1 || u < us!(i+1)  = 1
        | otherwise                          = 0
       go i p
           = (u - us!i)/(us!(i+p) - us!i) * go i (p-1)
            + (us!(i+p+1) - u)/(us!(i+p+1) - us!(i+1)) * go (i+1) (p-1)

L'altra cosa che non è ottimale è che non condividiamo le valutazioni di livello inferiore tra chiamate ricorsive vicine. (La valutazione si estende effettivamente su un grafo diretto con p22 ma lo valutiamo come un albero con 2 nodi p nodi). Questo è un massiccio inefficienza per grandi p ma in realtà è abbastanza innocua per le tipiche spline a basso grado.

Il modo per evitare questa inefficienza è quello di memorizzare. La versione C lo fa esplicitamente con il parametro N ma, trattandosi di Haskell, possiamo essere pigri e risparmiare lo sforzo di allocare la dimensione corretta utilizzando una libreria di memoizzazione generica, ad esempio memo-trie:

import Data.MemoTrie (memo2)

baseFuncs (Knots us) i₀ p₀ u = go' i₀ p₀
 where go i 0
        | u >= us!i
        , i>=VU.length us || u < us!(i+1)  = 1
        | otherwise                        = 0
       go i p
           = (u - us!i)/(us!(i+p) - us!i) * go' i (p-1)
            + (us!(i+p+1) - u)/(us!(i+p+1) - us!(i+1)) * go' (i+1) (p-1)
       go' = memo2 go

Questa era la versione senza cervello ("basta memorizzare l'intero dominio di go"). Come osserva dfeuer, è abbastanza facile memorizzare esplicitamente solo la regione che viene effettivamente valutata, e quindi si può usare di nuovo un efficiente vettore unboxed:

baseFuncs (Knots us) i₀ p₀ u = VU.unsafeHead $ gol i₀ p₀
 where gol i 0 = VU.generate (p₀+1) $ j ->
        if u >= us!(i+j)
            && (i+j>=VU.length us || u < us!(i+j+1))
         then 1 else 0
       gol i p = case gol i (p-1) of
        res' -> VU.izipWith
         (j l r -> let i' = i+j
            in (u - us!i')/(us!(i'+p) - us!i') * l
            + (us!(i'+p+1) - u)/(us!(i'+p+1) - us!(i'+1)) * r)
         res' (VU.unsafeTail res')

(posso tranquillamente usare unsafeHead e unsafeTail qui, perché ad ogni livello di ricorsione lo zipping riduce la lunghezza di 1, quindi al livello superiore ho ancora p₀ - (p₀-1) = 1 elementi).

Questa versione dovrebbe, credo, avere gli stessi asintoti della versione C. Con altri piccoli miglioramenti, come il precalcolo delle lunghezze degli intervalli e il precontrollo che gli argomenti siano nell'intervallo consentito, in modo che tutti gli accessi possano essere effettuati unsafe, probabilmente le prestazioni si avvicinano molto a quelle della versione C.

Come osserva ancora dfeuer, potrebbe anche non essere necessario usare i vettori, perché mi limito a comprimere il risultato. Per questo tipo di cose, GHC può essere molto bravo a ottimizzare il codice anche quando si usano liste semplici. Ma non indagherò ulteriormente sulle prestazioni in questa sede.


Il test che ho usato per confermare che funziona davvero:

Trama delle spline di base

https://gist.github.com/leftaroundabout/4fd6ef8642029607e1b222783b9d1c1e

(Disclaimer: non ho idea di cosa venga calcolato qui).

Lo schema di accesso all'array U sembra essere il seguente: dall'indice i verso l'esterno, si consumano valori a sinistra e a destra. Potremmo immaginare di avere due liste che consistono proprio in queste sequenze di elementi. In effetti, potremmo costruire tali liste a partire da un elenco di partenza come questo:

pryAt :: Int -> [a] -> ([a], [a])
pryAt i xs = go i ([], xs)
  where
    go 0 a = a
    go n (us, v : vs) = go (pred n) (v : us, vs)
-- pryAt 5 ['a'..'x']
-- ("edcba","fghijklmnopqrstuvwx")

Per i contenitori ad accesso casuale, potremmo avere versioni specializzate di pryAtperché attraversando l'intero elenco fino a raggiungere il punto isarà inefficiente.

Nel ciclo esterno, abbiamo gli array N, left e right che crescono a ogni iterazione (N sembra essere completamente ricostruito a ogni iterazione). Potremmo rappresentarli come liste. In ogni iterazione, consumiamo anche una coppia di elementi di Ua sinistra e a destra.

Il seguente tipo di dato rappresenta la situazione all'inizio di un'iterazione del ciclo esterno:

data State = State
  { ns :: [Double],
    left :: [Double],
    right :: [Double],
    us :: ([Double], [Double])
  }

Supponendo di avere già i dati outerStep :: State -> State potremmo semplicemente girare la manovella p volte:

basis :: Int -> Double -> Int -> [Double] -> [Double]
basis i u p us =
  ns $ iterate outerStep initial !! p
  where
    initial =
      State
        { ns = [1.0],
          left = [],
          right = [],
          us = pryAt i us
        }

Ciò che viene fatto a outerStep? Aggiungiamo nuovi elementi a left e righte poi ricreiamo l'intero N dall'inizio, portando con sé un elemento saved lungo il percorso. Questo è un mapAccumR. Abbiamo bisogno di qualche informazione in più: l'accumulatore right (nella stessa direzione di N) e i valori left (in direzione inversa), quindi è necessario zipparli in anticipo:

    outerStep (State {ns, left, right, us = (ul : uls, ur : urs)}) =
      let left' = u - ul : left
          right' = ur - u : right
          (saved', ns') = mapAccumR innerStep 0.0 $ zip3 ns right' (reverse left')
       in State
            { ns = saved' : ns',
              left = left',
              right = right',
              us = (uls, urs)
            }

Ed ecco i calcoli del passo interno:

    innerStep saved (n, r, l) = 
        let temp = n / (r - l)
            n' = saved + r
            saved' = l * temp
         in (saved', n')

Oltre a correggere i possibili bug, rimarrebbe altro lavoro, perché il basis nella sua forma attuale rischia di perdere memoria (in particolare, la funzione mapAccumR creerà molti thunk). Forse si potrebbe riscrivere la funzione in modo da usare funzioni come iterate' o foldl' che mantengono i loro accumulatori rigorosi.

Vi invitiamo ad aggiungere valore alle nostre informazioni collaborando con la vostra anzianità nelle recensioni.



Utilizzate il nostro motore di ricerca

Ricerca
Generic filters

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.