module Math.Spline.Knots
( Knots
, empty, isEmpty
, knot, multipleKnot
, mkKnots, fromList
, numKnots, lookupKnot
, toList, numDistinctKnots, lookupDistinctKnot
, knots, knotsVector
, distinctKnots, distinctKnotsVector
, toMap
, fromMap
, toVector
, fromVector
, splitLookup
, takeKnots, dropKnots, splitKnotsAt
, takeDistinctKnots, dropDistinctKnots, splitDistinctKnotsAt
, maxMultiplicity
, knotMultiplicity, setKnotMultiplicity
, fromAscList, fromDistinctAscList
, valid
, knotSpan
, knotsInSpan
, knotSpans
, knotDomain
, uniform
) where
import Prelude hiding (sum, maximum)
import Control.Monad (guard)
import Data.Foldable (Foldable(foldMap), sum, maximum)
import qualified Data.Map as M
import Data.Monoid (Monoid(..))
import Data.Maybe (fromMaybe)
import qualified Data.Set as S (Set)
import qualified Data.Vector as V
import Data.VectorSpace
data Knots a = Knots !Int (M.Map a Int) deriving (Eq, Ord)
instance Show a => Show (Knots a) where
showsPrec p ks@(Knots 0 _) = showString "empty"
showsPrec p ks@(Knots 1 _) = showParen (p > 10)
( showString "knot "
. showsPrec 11 (head $ knots ks)
)
showsPrec p ks = showParen (p > 10)
( showString "mkKnots "
. showsPrec 11 (knots ks)
)
instance (Ord a) => Monoid (Knots a) where
mempty = empty
mappend (Knots n1 v1) (Knots n2 v2) =
Knots (n1 + n2) (M.filter (/=0) (M.unionWith (+) v1 v2))
instance Foldable Knots where
foldMap f = foldMap f . knots
empty :: Knots a
empty = Knots 0 M.empty
isEmpty :: Knots a -> Bool
isEmpty (Knots 0 _) = True
isEmpty _ = False
knot :: Ord a => a -> Knots a
knot x = multipleKnot x 1
multipleKnot :: Ord a => a -> Int -> Knots a
multipleKnot k n
| n <= 0 = Knots 0 (M.empty)
| otherwise = Knots n (M.singleton k n)
mkKnots :: (Ord a) => [a] -> Knots a
mkKnots ks = fromList (map (\k -> (k,1)) ks)
fromList :: (Ord k) => [(k, Int)] -> Knots k
fromList ks = Knots (sum kMap) kMap
where kMap = M.fromListWith (+) (filter ((>0).snd) ks)
fromAscList :: Eq k => [(k, Int)] -> Knots k
fromAscList ks = Knots (sum kMap) kMap
where kMap = M.fromAscListWith (+) (filter ((>0).snd) ks)
fromDistinctAscList :: [(k, Int)] -> Knots k
fromDistinctAscList ks = Knots (sum kMap) kMap
where kMap = M.fromDistinctAscList (filter ((>0).snd) ks)
fromMap :: M.Map k Int -> Knots k
fromMap ks = Knots (sum kMap) kMap
where
kMap = mFilter (>0) ks
mFilter p = M.fromDistinctAscList . filter (p.snd) . M.toAscList
fromVector :: Ord k => V.Vector (k,Int) -> Knots k
fromVector = fromList . V.toList
toList :: Knots k -> [(k, Int)]
toList = M.toList . toMap
toVector :: Knots k -> V.Vector (k, Int)
toVector = V.fromList . toList
toMap :: Knots k -> M.Map k Int
toMap (Knots _ ks) = ks
numKnots :: Knots t -> Int
numKnots (Knots n _) = n
numDistinctKnots :: Knots t -> Int
numDistinctKnots (Knots _ ks) = M.size ks
maxMultiplicity :: Knots t -> Int
maxMultiplicity (Knots 0 _) = 0
maxMultiplicity (Knots _ ks) = maximum ks
lookupKnot :: Int -> Knots a -> Maybe a
lookupKnot k kts
| k < 0 = Nothing
| k < numKnots kts = fmap fst mbKt
| otherwise = Nothing
where (_, mbKt, _) = splitLookup k kts
lookupDistinctKnot :: Int -> Knots a -> Maybe a
lookupDistinctKnot k (Knots _ ks)
| k < 0 = Nothing
| k < M.size ks = Just (fst (M.elemAt k ks))
| otherwise = Nothing
splitLookup :: Int -> Knots a -> (Knots a, Maybe (a, Int), Knots a)
splitLookup k (Knots n ks) = scan 0 M.empty n ks
where
scan nPre pre nPost post
| nPost <= 0 = (Knots nPre pre, Nothing, Knots nPost post)
| nPre + m > k = (Knots nPre pre, Just kt, Knots nNewPost newPost)
| otherwise = scan (nPre + m) (pre `ascSnoc` kt) nNewPost newPost
where
Just (kt@(x,m), newPost) = M.minViewWithKey post
nNewPost = nPost m
done x = (Knots nPre pre, x, Knots nPost post)
ascCons x m = M.fromDistinctAscList (x : M.toAscList m)
ascSnoc m x = M.fromDistinctAscList (M.toAscList m ++ [x])
ascConsKnot (_,0) kts = kts
ascConsKnot kt@(k,m) (Knots n ks) = Knots (n+m) (kt `ascCons` ks)
ascSnocKnot kts (_,0) = kts
ascSnocKnot (Knots n ks) kt@(k,m) = Knots (n+m) (ks `ascSnoc` kt)
clamp lo hi = max lo . min hi
dropKnots :: Int -> Knots a -> Knots a
dropKnots k kts = fromMaybe post $ do
(x,xAvail) <- mbKt
let xWanted = numKnots kts (numKnots post + k)
return ((x, clamp 0 xAvail xWanted) `ascConsKnot` post)
where
(pre, mbKt, post) = splitLookup k kts
takeKnots :: Int -> Knots a -> Knots a
takeKnots k kts = fromMaybe pre $ do
(x,xAvail) <- mbKt
let xWanted = k numKnots pre
return (pre `ascSnocKnot` (x, clamp 0 xAvail xWanted))
where
(pre, mbKt, post) = splitLookup k kts
splitKnotsAt :: Int -> Knots a -> (Knots a, Knots a)
splitKnotsAt k kts = fromMaybe (pre, post) $ do
(x,xAvail) <- mbKt
let xWanted = k numKnots pre
xTaken = clamp 0 xAvail xWanted
return ( pre `ascSnocKnot` (x,xTaken)
, (x, xAvail xTaken) `ascConsKnot` post
)
where
(pre, mbKt, post) = splitLookup k kts
takeDistinctKnots :: Int -> Knots a -> Knots a
takeDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap
where
kMap = M.fromDistinctAscList (take k (M.toAscList ks))
dropDistinctKnots :: Int -> Knots a -> Knots a
dropDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap
where
kMap = M.fromDistinctAscList (drop k (M.toAscList ks))
splitDistinctKnotsAt :: Int -> Knots a -> (Knots a, Knots a)
splitDistinctKnotsAt k (Knots n ks) = (Knots sz1 kMap1, Knots (n sz1) kMap2)
where
(ks1, ks2) = splitAt k (M.toAscList ks)
kMap1 = M.fromDistinctAscList ks1
kMap2 = M.fromDistinctAscList ks2
sz1 = sum kMap1
knots :: Knots t -> [t]
knots (Knots _ ks) = concat [replicate n k | (k,n) <- M.toAscList ks]
knotsVector :: Knots t -> V.Vector t
knotsVector (Knots _ ks) = V.concat [V.replicate n k | (k,n) <- M.toAscList ks]
distinctKnots :: Knots t -> [t]
distinctKnots (Knots _ ks) = M.keys ks
distinctKnotsVector :: Knots t -> V.Vector t
distinctKnotsVector = V.fromList . distinctKnots
distinctKnotsSet :: Knots k -> S.Set k
distinctKnotsSet (Knots _ ks) = M.keysSet ks
knotMultiplicity :: (Ord k) => k -> Knots k -> Int
knotMultiplicity k (Knots _ ks) = fromMaybe 0 (M.lookup k ks)
setKnotMultiplicity :: Ord k => k -> Int -> Knots k -> Knots k
setKnotMultiplicity k n (Knots m ks)
| n <= 0 = Knots (m n') (M.delete k ks)
| otherwise = Knots (m + n n') (M.insert k n ks)
where
n' = knotMultiplicity k (Knots m ks)
valid :: Ord k => Knots k -> Bool
valid (Knots n ks) = and
[ M.valid ks
, n == sum ks
, all (>0) (M.elems ks)
]
knotSpan :: Knots a -> Int -> Int -> Maybe (a, a)
knotSpan ks i j = do
guard (i <= j)
lo <- lookupKnot i ks
hi <- lookupKnot j ks
return (lo,hi)
knotsInSpan :: Knots a -> Int -> Int -> Knots a
knotsInSpan kts i j = takeKnots (j i) (dropKnots i kts)
knotSpans :: Knots a -> Int -> [(a,a)]
knotSpans ks w
| w <= 0 = error "knotSpans: width must be positive"
| otherwise = zip kts (drop w kts)
where kts = knots ks
knotDomain :: Knots a -> Int -> Maybe (a,a)
knotDomain ks@(Knots n _) p = knotSpan ks p (np1)
uniform :: (Ord s, Fractional s) => Int -> Int -> (s,s) -> Knots s
uniform deg nPts (lo,hi) = ends `mappend` internal
where
ends = fromList [(lo,deg), (hi,deg)]
n = nPts + deg numKnots ends
f i = (fromIntegral i * lo + fromIntegral (n i) * hi) / fromIntegral n
internal = mkKnots [f i | i <- [0..n]]