hs-memoised-splines.hs

{-# LANGUAGE OverloadedLists #-}

import qualified Data.Vector.Unboxed as VU
import Data.Vector.Unboxed (Vector, (!))
import Data.MemoTrie (memo2)
import Graphics.Dynamic.Plot.R2
import Text.Printf

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

main :: IO ()
main = do
  plotWindow
   [ plotLatest
    $ plotDelay 4 mempty :
    [ legendName (printf "𝑖=%i, 𝑝=%i" i p)
      . plotDelay 1
      . continFnPlot $ baseFuncs (Knots us) i p
    | p <- [0..4]
    , i <- [0..n-p-2]
    ]
   , dynamicAxes
   ]
  return ()
 where us = [0, 0.2, 0.3, 0.6, 0.9, 1]
       n = VU.length us