module Math.Spline.Knots
( Knots
, empty, isEmpty
, knot, multipleKnot
, mkKnots, fromList
, numKnots, lookupKnot
, toList, numDistinctKnots, lookupDistinctKnot
, knots, knotsVector
, distinctKnots, multiplicities
, distinctKnotsVector, multiplicitiesVector
, distinctKnotsSet
, toMap
, fromMap
, toVector
, fromVector
, splitLookup
, takeKnots, dropKnots, splitKnotsAt
, takeDistinctKnots, dropDistinctKnots, splitDistinctKnotsAt
, maxMultiplicity
, knotMultiplicity, setKnotMultiplicity
, splitFind
, fromAscList, fromDistinctAscList
, valid
, knotSpan
, knotsInSpan
, knotSpans
, knotDomain
, uniform
, minKnot
, maxKnot
) where
import Prelude hiding (sum, maximum)
import Control.Arrow ((***))
import Control.Monad (guard)
import Data.Foldable (Foldable(foldMap), maximum)
import Data.List (sortBy, sort, unfoldr)
import qualified Data.Map as M
import Data.Monoid (Monoid(..))
import Data.Ord
import qualified Data.Set as S (Set, fromAscList)
import qualified Data.Vector as V
newtype Knots a = Knots (V.Vector a) deriving (Eq, Ord)
instance Show a => Show (Knots a) where
showsPrec p ks = showParen (p > 10)
( showString "mkKnots "
. showsPrec 11 (knots ks)
)
instance (Ord a) => Monoid (Knots a) where
mempty = empty
mappend (Knots v1) (Knots v2) =
Knots . V.fromList . sort . V.toList $ v1 V.++ v2
instance Foldable Knots where
foldMap f = foldMap f . knotsVector
empty :: Knots a
empty = Knots V.empty
isEmpty :: Knots a -> Bool
isEmpty (Knots v) = V.null v
knot :: Ord a => a -> Knots a
knot x = multipleKnot x 1
multipleKnot :: Ord a => a -> Int -> Knots a
multipleKnot k n = Knots $ V.replicate n k
mkKnots :: (Ord a) => [a] -> Knots a
mkKnots = Knots . V.fromList . sort
fromList :: (Ord k) => [(k, Int)] -> Knots k
fromList ks = Knots v
where v = V.concat . map (\(k, mult) -> V.replicate mult k) .
sortBy (comparing fst) $ filter ((>0).snd) ks
fromAscList :: Eq k => [(k, Int)] -> Knots k
fromAscList ks = Knots v
where v = V.concat . map (\(k, mult) -> V.replicate mult k)
$ filter ((>0).snd) ks
fromDistinctAscList :: Eq k => [(k, Int)] -> Knots k
fromDistinctAscList = fromAscList
fromMap :: Eq k => M.Map k Int -> Knots k
fromMap = fromAscList . M.toAscList
fromVector :: Ord k => V.Vector (k,Int) -> Knots k
fromVector = fromList . V.toList
toList :: Eq k => Knots k -> [(k, Int)]
toList = unfoldr $ \kts -> do
kt <- minKnot kts
return (kt, dropKnots (snd kt) kts)
toVector :: Eq k => Knots k -> V.Vector (k, Int)
toVector = V.unfoldr $ \kts -> do
kt <- minKnot kts
return (kt, dropKnots (snd kt) kts)
toMap :: Ord k => Knots k -> M.Map k Int
toMap = M.fromAscListWith (+) . toList
numKnots :: Knots t -> Int
numKnots (Knots v) = V.length v
numDistinctKnots :: Eq t => Knots t -> Int
numDistinctKnots = V.length . distinctKnotsVector
maxMultiplicity :: Ord t => Knots t -> Int
maxMultiplicity kts@(Knots v)
| V.length v == 0 = 0
| otherwise = maximum $ toMap kts
lookupKnot :: Int -> Knots a -> Maybe a
lookupKnot k (Knots kts) = do
guard (0 <= k && k < V.length kts)
V.indexM kts k
lookupDistinctKnot :: Eq a => Int -> Knots a -> Maybe a
lookupDistinctKnot k kts = lookupKnot k . Knots $ distinctKnotsVector kts
splitLookup :: Int -> Knots a -> (Knots a, Maybe a, Knots a)
splitLookup k (Knots v)
| V.null gt = (Knots lt, Nothing, Knots V.empty)
| otherwise = (Knots lt, Just $ V.head gt, Knots $ V.tail gt)
where
(lt, gt) = V.splitAt k v
dropKnots :: Int -> Knots a -> Knots a
dropKnots k (Knots v) = Knots $ V.drop k v
takeKnots :: Int -> Knots a -> Knots a
takeKnots k (Knots v) = Knots $ V.take k v
splitKnotsAt :: Int -> Knots a -> (Knots a, Knots a)
splitKnotsAt k (Knots v) = Knots *** Knots $ V.splitAt k v
findDistinctKnot :: Eq a => Int -> Knots a -> Int
findDistinctKnot n = V.last . V.take (1 + max 0 n) . V.scanl (+) 0 . multiplicitiesVector
takeDistinctKnots :: (Ord a) => Int -> Knots a -> Knots a
takeDistinctKnots k kts = takeKnots (findDistinctKnot k kts) kts
dropDistinctKnots :: (Ord a) => Int -> Knots a -> Knots a
dropDistinctKnots k kts = dropKnots (findDistinctKnot k kts) kts
splitDistinctKnotsAt :: (Ord a, Eq a) => Int -> Knots a -> (Knots a, Knots a)
splitDistinctKnotsAt k kts = splitKnotsAt (findDistinctKnot k kts) kts
knots :: Knots t -> [t]
knots = V.toList . knotsVector
knotsVector :: Knots t -> V.Vector t
knotsVector (Knots v) = v
distinctKnots :: Eq t => Knots t -> [t]
distinctKnots = map fst . toList
multiplicities :: Eq t => Knots t -> [Int]
multiplicities = map snd . toList
distinctKnotsVector :: Eq t => Knots t -> V.Vector t
distinctKnotsVector = V.map fst . toVector
multiplicitiesVector :: Eq a => Knots a -> V.Vector Int
multiplicitiesVector = V.map snd . toVector
distinctKnotsSet :: Eq k => Knots k -> S.Set k
distinctKnotsSet (Knots k) = S.fromAscList $ V.toList k
knotMultiplicity :: (Ord k) => k -> Knots k -> Int
knotMultiplicity k (Knots ks) = V.length $ V.elemIndices k ks
setKnotMultiplicity :: Ord k => k -> Int -> Knots k -> Knots k
setKnotMultiplicity k n kts@(Knots v)
| n <= 0 = Knots (V.filter (/= k) v)
| otherwise = Knots $ V.concat [lt, V.replicate n k, gt]
where (Knots lt, _, Knots gt) = splitFind k kts
splitFind :: Ord k => k -> Knots k -> (Knots k, Knots k, Knots k)
splitFind k (Knots v) = (Knots lt, Knots eq, Knots gt)
where
(lt, xs) = V.span (<k) v
(eq, gt) = V.span (==k) xs
valid :: (Ord k, Num k) => Knots k -> Bool
valid (Knots v)
| V.null v = True
| otherwise = V.and $ V.zipWith (>=) (V.tail v) v
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 v) p = knotSpan ks p (V.length vp1)
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]]
minKnot :: (Eq a) => Knots a -> Maybe (a, Int)
minKnot (Knots v)
| V.null v = Nothing
| otherwise = Just (kt, V.length (V.takeWhile (kt ==) v))
where kt = V.head v
maxKnot :: Eq a => Knots a -> Maybe (a, Int)
maxKnot (Knots v) = minKnot (Knots (V.reverse v))