splines 0.3 → 0.5.0.1

raw patch · 14 files changed

+595/−236 lines, 14 filesdep +QuickCheckdep +criteriondep +splinesdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, criterion, splines, test-framework, test-framework-quickcheck2

Dependency ranges changed: vector

API changes (from Hackage documentation)

Files

• benchmark/DeBoor.hs +63/−0
• splines.cabal +30/−2
• src/Math/NURBS.hs +9/−7
• src/Math/Spline.hs +2/−0
• src/Math/Spline/BSpline.hs +40/−37
• src/Math/Spline/BSpline/Internal.hs +168/−19
• src/Math/Spline/BSpline/Reference.hs +45/−13
• src/Math/Spline/BezierCurve.hs +1/−1
• src/Math/Spline/Class.hs +17/−4
• src/Math/Spline/Hermite.hs +89/−0
• src/Math/Spline/ISpline.hs +2/−0
• src/Math/Spline/Knots.hs +112/−153
• src/Math/Spline/MSpline.hs +3/−0
• test/Main.hs +14/−0

@@ -0,0 +1,63 @@+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE FlexibleContexts #-}++import Criterion.Main++import Data.List (find)+import qualified Data.Vector.Generic.Safe as V+import qualified Data.Vector.Safe as BV+import qualified Data.Vector.Unboxed.Safe as UV+import Math.Spline.BSpline+import Math.Spline.BSpline.Reference+import Math.Spline.Knots+import Math.Polynomial++import Debug.Trace+import Control.Monad++kts = mkKnots $ [0,0,0] ++ [0..10] ++ [11,11,11]+ctPts = map sin [0..12]++unboxedSpline :: BSpline UV.Vector Double+unboxedSpline = bSpline kts (V.fromList ctPts)++boxedSpline :: BSpline BV.Vector Double+boxedSpline = bSpline kts (V.fromList ctPts)+++intervalPoly :: [(Double, Double, Poly Double)]+intervalPoly = map f $ zip3 dkts (tail dkts) (basisPolynomials kts)+  where+    dkts = distinctKnots kts+    f (begin, end, k) = (begin, end, ) . sumPolys $ zipWith scalePoly ctPts (k !! 3)++applyDeBoor :: V.Vector v Double => BSpline v Double -> Double -> Double+applyDeBoor s = evalBSpline s++{-# SPECIALIZE applyNaturalDeBoor :: BSpline BV.Vector Double -> Double -> Double #-}+{-# SPECIALIZE applyNaturalDeBoor :: BSpline UV.Vector Double -> Double -> Double #-}+applyNaturalDeBoor :: V.Vector v Double =>  BSpline v Double -> Double -> Double+applyNaturalDeBoor s = evalNaturalBSpline s++applyPoly :: Double -> Double+applyPoly x = maybe 0 (\(_,_,p) -> evalPoly p x) $ find (\(b,e,_) -> x >= b && x < e) intervalPoly++applyAndSum :: (Double -> Double) -> [Double] -> Double+applyAndSum f = sum . map f++main = defaultMain+       [ bgroup "Boxed"+         [ bench "deBoor 1000" $ whnf (applyAndSum (applyDeBoor boxedSpline)) [0,0.01..10]+         , bench "deBoor 10000" $ whnf (applyAndSum (applyDeBoor boxedSpline)) [0,0.001..10]+         , bench "natural 1000" $ whnf (applyAndSum (applyNaturalDeBoor boxedSpline)) [0,0.01..10]+         , bench "natural 10000" $ whnf (applyAndSum (applyNaturalDeBoor boxedSpline)) [0,0.001..10]+         ]+       , bgroup "Unboxed"+         [ bench "deBoor 1000" $ whnf (applyAndSum (applyDeBoor unboxedSpline)) [0,0.01..10]+         , bench "deBoor 10000" $ whnf (applyAndSum (applyDeBoor unboxedSpline)) [0,0.001..10]+         , bench "natural 1000" $ whnf (applyAndSum (applyNaturalDeBoor unboxedSpline)) [0,0.01..10]+         , bench "natural 10000" $ whnf (applyAndSum (applyNaturalDeBoor unboxedSpline)) [0,0.001..10]+         ]+       , bench "poly 1000" $ whnf (applyAndSum applyPoly) [0,0.01..10]+       , bench "poly 10000" $ whnf (applyAndSum applyPoly) [0,0.001..10]+       ]

@@ -1,8 +1,8 @@ name:                   splines-version:                0.3+version:                0.5.0.1 stability:              provisional -cabal-version:          >= 1.6+cabal-version:          >= 1.9.2 build-type:             Simple  author:                 James Cook <mokus@deepbondi.net>@@ -24,11 +24,13 @@  Library   hs-source-dirs:       src+  ghc-options:          -Wall   exposed-modules:      Math.Spline                         Math.Spline.BezierCurve                         Math.Spline.BSpline                         Math.Spline.BSpline.Reference                         Math.Spline.Class+                        Math.Spline.Hermite                         Math.Spline.ISpline                         Math.Spline.Knots                         Math.Spline.MSpline@@ -37,5 +39,31 @@   build-depends:        base >= 3 && < 5,                         containers,                         polynomial,+                        vector >= 0.8,+                        vector-space++Test-Suite splines-test+  type:                 exitcode-stdio-1.0+  hs-source-dirs:       test+  main-is:              Main.hs+  +  build-depends:        base >= 3 && <5,+                        containers,+                        polynomial,+                        splines,+                        test-framework,+                        test-framework-quickcheck2,+                        QuickCheck >= 2,                         vector,                         vector-space++Benchmark splines-bench+  type:                 exitcode-stdio-1.0+  hs-source-dirs:       benchmark+  main-is:              DeBoor.hs+  +  build-depends:        base >= 3 && < 5,+                        criterion,+                        polynomial,+                        splines,+                        vector

@@ -8,13 +8,13 @@     ) where  import qualified Data.Vector as V-import Data.VectorSpace+import Data.VectorSpace hiding (project) import Math.Spline.Class (Spline, toBSpline) import Math.Spline.BSpline.Internal import Math.Spline.BSpline import Math.Spline.Knots -newtype NURBS v = NURBS (BSpline (Scalar v, v))+newtype NURBS v = NURBS (BSpline V.Vector (Scalar v, v))  deriving instance (Eq   v, Eq   (Scalar v), Eq   (Scalar (Scalar v))) => Eq   (NURBS v) deriving instance (Ord  v, Ord  (Scalar v), Ord  (Scalar (Scalar v))) => Ord  (NURBS v)@@ -34,16 +34,18 @@  -- |Constructs the homogeneous-coordinates B-spline that corresponds to this -- NURBS curve+nurbsAsSpline :: VectorSpace v => NURBS v -> BSpline V.Vector (Scalar v, v) nurbsAsSpline (NURBS spline) = spline      { controlPoints = V.map homogenize (controlPoints spline) }     where-        homogenize (w,v) = (w, w *^ v)+        homogenize (w,v) = (w, v ^* w)  -- |Constructs the NURBS curve corresponding to a homogeneous-coordinates B-spline+splineAsNURBS :: (VectorSpace v, Fractional (Scalar v)) => BSpline V.Vector (Scalar v, v) -> NURBS v splineAsNURBS spline = NURBS spline      { controlPoints = V.map unHomogenize (controlPoints spline) }     where-        unHomogenize (w,v) = (w, recip w *^ v)+        unHomogenize (w,v) = (w, v ^/ w)   evalNURBS@@ -51,7 +53,7 @@       VectorSpace w, Scalar w ~ w,       Fractional w, Ord w) =>      NURBS v -> w -> v-evalNURBS nurbs = project . evalBSpline (nurbsAsSpline nurbs)+evalNURBS f = project . evalBSpline (nurbsAsSpline f)     where         project (w,v) = recip w *^ v @@ -75,6 +77,6 @@                VectorSpace w, Scalar w ~ w,                Ord w, Fractional w)     => NURBS v -> Scalar v -> Maybe (NURBS v, NURBS v)-splitNURBS nurbs t = do-    (s0, s1) <- splitBSpline (nurbsAsSpline nurbs) t+splitNURBS f t = do+    (s0, s1) <- splitBSpline (nurbsAsSpline f) t     return (splineAsNURBS s0, splineAsNURBS s1)

@@ -7,6 +7,7 @@     , BSpline, bSpline     , MSpline, mSpline, toMSpline     , ISpline, iSpline, toISpline+    , CSpline, cSpline     ) where  import Math.Spline.Class@@ -15,3 +16,4 @@ import Math.Spline.BSpline import Math.Spline.MSpline import Math.Spline.ISpline+import Math.Spline.Hermite

@@ -1,11 +1,13 @@-{-# LANGUAGE MultiParamTypeClasses, StandaloneDeriving, FlexibleContexts, UndecidableInstances, TypeFamilies, ParallelListComp #-}+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, TypeFamilies, ParallelListComp #-} module Math.Spline.BSpline     ( BSpline     , bSpline     , evalBSpline+    , evalNaturalBSpline     , insertKnot     , splitBSpline     , differentiateBSpline, integrateBSpline+    , deBoor     ) where  import Math.Spline.Knots@@ -13,20 +15,19 @@  import Data.Maybe (fromMaybe) import Data.VectorSpace-import qualified Data.Vector as V+import qualified Data.Vector.Generic as V  -- |@bSpline kts cps@ creates a B-spline with the given knot vector and control  -- points.  The degree is automatically inferred as the difference between the  -- number of spans in the knot vector (@numKnots kts - 1@) and the number of  -- control points (@length cps@).-bSpline :: Knots (Scalar a) -> V.Vector a -> BSpline a-bSpline kts cps-    | V.null cps    = error "bSpline: no control points"-    | otherwise     = fromMaybe-        (error "bSpline: too few knots")-        (maybeSpline kts cps)+bSpline :: V.Vector v a => Knots (Scalar a) -> v a -> BSpline v a+bSpline kts cps = fromMaybe+    (error "bSpline: too few knots")+    (maybeSpline kts cps) -maybeSpline :: Knots (Scalar a) -> V.Vector a -> Maybe (BSpline a)+-- not exported; precondition: n > 0+maybeSpline :: V.Vector v a => Knots (Scalar a) -> v a -> Maybe (BSpline v a) maybeSpline kts cps      | n > m     = Nothing     | otherwise = Just (Spline (m - n) kts cps)@@ -34,23 +35,18 @@         n = V.length cps         m = numKnots kts - 1 -deriving instance (Eq   (Scalar v), Eq   v) => Eq   (BSpline v)-deriving instance (Ord  (Scalar v), Ord  v) => Ord  (BSpline v)-instance (Show (Scalar v), Show v) => Show (BSpline v) where-    showsPrec p (Spline _ kts cps) = showParen (p>10) -        ( showString "bSpline "-        . showsPrec 11 kts-        . showChar ' '-        . showsPrec 11 cps-        )- differentiateBSpline-  :: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v+    :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a), V.Vector v a+       , V.Vector v (Scalar a)) => BSpline v a -> BSpline v a differentiateBSpline spline-    | numKnots ks  < 2  = spline-    | numKnots ks == 2  = bSpline ks (V.singleton zeroV)-    | otherwise         = bSpline ks' ds'+    | V.null ds = error "differentiateBSpline: no control points"+    | m  < 1    = spline+    | p == 0    = bSpline ks (V.replicate n zeroV)+    | otherwise = bSpline ks' ds'     where+        n = V.length ds+        m = numKnots ks - 1+                 ks' = mkKnots . init . tail $ ts         ds' = V.zipWith (*^) (V.tail cs) (V.zipWith (^-^) (V.tail ds) ds)         @@ -58,10 +54,12 @@         ds = controlPoints spline                  p  = degree spline-        cs = V.fromList [fromIntegral p / (t1 - t0) | (t0,t1) <- spans p ts]+        cs = V.fromList [ if t1 /= t0 then fromIntegral p / (t1 - t0) else 0 | (t0,t1) <- spans p ts]  integrateBSpline-  :: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v+  :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a), V.Vector v a+     , V.Vector v (Scalar a)) =>+     BSpline v a -> BSpline v a integrateBSpline spline = bSpline (mkKnots ts') (V.scanl (^+^) zeroV ds')     where         ds' = V.zipWith (*^) cs (controlPoints spline)@@ -70,26 +68,31 @@         p = degree spline + 1         cs = V.fromList [(t1 - t0) / fromIntegral p | (t0,t1) <- spans p ts] +spans :: Int -> [a] -> [(a,a)] spans n xs = zip xs (drop n xs) +-- |Split a B-spline at the specified point (which must be inside the spline's domain),+-- returning two disjoint splines, the sum of which is equal to the original.  The domain+-- of the first will be below the split point and the domain of the second will be above. splitBSpline-  :: (VectorSpace v, Ord (Scalar v), Fractional (Scalar v)) =>-     BSpline v -> Scalar v -> Maybe (BSpline v, BSpline v)-splitBSpline spline@(Spline p kv ds) t -    | inDomain  = Just (Spline p (mkKnots us0) ds0, Spline p (mkKnots us1) ds1)+  :: ( VectorSpace a, Ord (Scalar a), Fractional (Scalar a), V.Vector v a+     , V.Vector v (Scalar a)) =>+     BSpline v a -> Scalar a -> Maybe (BSpline v a, BSpline v a)+splitBSpline spline@(Spline p kv _) t +    | inDomain  = Just (Spline p us0 ds0, Spline p us1 ds1)     | otherwise = Nothing     where         inDomain = case knotDomain kv p of             Nothing         -> False-            Just (t0, t1)   -> t >= t0 || t <= t1+            Just (t0, t1)   -> t >= t0 && t <= t1         -        us = knots kv+        (lt, _, gt) = splitFind t kv         dss = deBoor spline t         -        us0 = takeWhile (<t) us ++ replicate (p+1) t-        ds0 = V.fromList (trimTo (drop (p+1) us0) (map V.head dss))+        us0 = setKnotMultiplicity t (p+1) lt+        ds0 = trimTo us0 (map V.head dss)         -        us1 = replicate (p+1) t ++ dropWhile (<=t) us-        ds1 = V.reverse (V.fromList (trimTo (drop (p+1) us1) (map V.last dss)))--        trimTo list  xs = zipWith const xs list+        us1 = setKnotMultiplicity t (p+1) gt+        ds1 = V.reverse (trimTo us1 (map V.last dss))+        +        trimTo kts = V.fromList . take (numKnots kts - p - 1)

@@ -1,61 +1,156 @@+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-} module Math.Spline.BSpline.Internal-    (BSpline(..), mapControlPoints, evalBSpline, insertKnot, deBoor) where+    ( BSpline(..)+    , mapControlPoints+    , evalBSpline+    , evalNaturalBSpline+    , insertKnot+    , deBoor+    , slice+    ) where  import Math.Spline.Knots +import Control.Monad.ST import Data.Monoid-import Data.Vector as V+import qualified Data.Vector.Generic as V+import qualified Data.Vector.Generic.Mutable as MV+import qualified Data.Vector as BV (Vector)+import qualified Data.Vector.Unboxed as UV (Vector, Unbox) import Data.VectorSpace import Prelude as P -data BSpline t = Spline+-- | A B-spline, defined by a knot vector (see 'Knots') and a sequence of control points.+data BSpline v t = Spline     { degree        :: !Int     , knotVector    :: Knots (Scalar t)-    , controlPoints :: Vector t+    , controlPoints :: v t     } +deriving instance (Eq (Scalar a), Eq (v a)) => Eq   (BSpline v a)+deriving instance (Ord (Scalar a), Ord (v a)) => Ord  (BSpline v a)+instance (Show (Scalar a), Show a, Show (v a)) => Show (BSpline v a) where+    showsPrec p (Spline _ kts cps) = showParen (p>10) +        ( showString "bSpline "+        . showsPrec 11 kts+        . showChar ' '+        . showsPrec 11 cps+        )++mapControlPoints :: (Scalar a ~ Scalar b, V.Vector v a, V.Vector v b) => (a -> b) -> BSpline v a -> BSpline v b+{-# SPECIALIZE mapControlPoints :: (Scalar a ~ Scalar b) => (a -> b)+ -> BSpline BV.Vector a -> BSpline BV.Vector b #-} mapControlPoints f spline = spline     { controlPoints = V.map f (controlPoints spline)     , knotVector = knotVector spline     } -evalBSpline spline = V.head . P.last . deBoor spline+-- |Evaluate a B-spline at the given point.  This uses a slightly modified version of +-- de Boor's algorithm which is only strictly correct inside the domain of the spline.+-- Unlike the standard algorithm, the basis functions always sum to 1, even outside the+-- domain of the spline.  This is mainly useful for \"clamped\" splines - the values at+-- or outside the endpoints will always be the value of the nearest control point.+--+-- For a standard implementation of de Boor's algorithm, see 'evalNaturalBSpline'.+-- For a (much slower) strictly mathematically correct evaluation, see 'evalReferenceBSpline'.+evalBSpline :: ( VectorSpace a, Fractional (Scalar a), Ord (Scalar a)+               , V.Vector v a, V.Vector v (Scalar a))+     => BSpline v a -> Scalar a -> a+{-# SPECIALIZE evalBSpline :: ( VectorSpace a, Fractional (Scalar a)+  , Ord (Scalar a)) => BSpline BV.Vector a -> Scalar a -> a#-}+evalBSpline spline+     | V.null (controlPoints spline) = zeroV+     | otherwise = V.head . P.last . deBoor spline --- |Insert one knot into a 'BSpline'+-- | Evaluate a B-spline at the given point.  This uses de Boor's algorithm, which is +-- only strictly correct inside the domain of the spline.+-- +-- For a (much slower) strictly mathematically correct evaluation, see 'evalReferenceBSpline'.+{-# INLINE evalNaturalBSpline #-}+{-# SPECIALIZE evalNaturalBSpline :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a) +    => BSpline BV.Vector a -> Scalar a -> a #-}+{-# SPECIALIZE evalNaturalBSpline :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, UV.Unbox a) +    => BSpline UV.Vector a -> Scalar a -> a #-}+evalNaturalBSpline :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, V.Vector v a)+    => BSpline v a -> Scalar a -> a+evalNaturalBSpline spline x = runST $ do+    let Spline p kvec cps = slice spline x+        kts = V.tail (knotsVector kvec)+        s = p - 1 - V.length (V.takeWhile (x==) kts)+    +    ds <- V.thaw cps+    +    sequence_+        [ do+            -- No need to check whether u0 < x < u1:+            --  x > u0 is guaranteed by choice of 's'+            --  x < u1 is guaranteed by 'slice'+            let !u0 = kts V.! (j + i)+                !u1 = kts V.! (j + p)+                !du = u1 - u0+                !a  = if du <= 0 then 1 else (x - u0) / du+            +            d0 <- MV.read ds  j+            d1 <- MV.read ds (j + 1)+            MV.write ds j (lerp d0 d1 a)+        | i <- [0 .. s]+        , j <- [0 .. s - i]+        ]+    +    MV.read ds 0++-- |Insert one knot into a 'BSpline' without changing the spline's shape. insertKnot-  :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a)) =>-     BSpline a -> Scalar a -> BSpline a+  :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a), V.Vector v a, V.Vector v (Scalar a)) =>+     BSpline v a -> Scalar a -> BSpline v a+{-# SPECIALIZE insertKnot :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a)) =>+  BSpline BV.Vector a -> Scalar a -> BSpline BV.Vector a #-} insertKnot spline x = spline     { knotVector    = knotVector spline `mappend` knot x     , controlPoints = V.zipWith4 (interp x) us (V.drop p us) ds (V.tail ds)     }     where-        us = knotsVector (knotVector spline)+        us = V.convert $ knotsVector (knotVector spline)         p  = degree spline         ds = extend (controlPoints spline) - -- duplicate the endpoints of a list; for example, -- extend [1..5] -> [1,1,2,3,4,5,5]+{-# INLINE extend #-}+extend :: V.Vector v t => v t -> v t extend vec     | V.null vec    = V.empty-    | otherwise     = V.cons (V.head vec) (V.snoc vec (V.last vec)) +    | otherwise     = V.cons (V.head vec) (V.snoc vec (V.last vec)) -deBoor spline x = go us (controlPoints spline)+{-# SPECIALIZE deBoor :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a)+  => BSpline BV.Vector a -> Scalar a -> [BV.Vector a] #-}+{-# SPECIALIZE deBoor :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a,+                          UV.Unbox a, UV.Unbox (Scalar a))+  => BSpline UV.Vector a -> Scalar a -> [UV.Vector a] #-}++-- | The table from de Boor's algorithm, calculated for the entire spline.  If that is not necessary+-- (for example, if you are only evaluating the spline), then use 'slice' on the spline first.+-- 'splitBSpline' currently uses the whole table.  It is probably not necessary there, but it +-- greatly simplifies the definition and makes the similarity to splitting Bezier curves very obvious.+deBoor :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, V.Vector v a, V.Vector v (Scalar a))+    => BSpline v a -> Scalar a -> [v a]+deBoor spline x = go us0 (controlPoints spline)     where-        us = knotsVector (knotVector spline)-        +        us0 = V.convert $ knotsVector (knotVector spline)         -- Upper endpoints of the intervals are the same for         -- each row in the table (they just line up differently         -- with the lower endpoints):-        uHi = V.drop (degree spline + 1) us-        -        -- On each pass, the lower endpoints of the -        -- interpolation intervals advance and the new +        uHi = V.drop (degree spline + 1) us0++        -- On each pass, the lower endpoints of the+        -- interpolation intervals advance and the new         -- coefficients are given by linear interpolation         -- on the current intervals:-        go us ds +        go us ds             | V.null ds = []             | otherwise = ds : go uLo ds'             where@@ -63,10 +158,64 @@                 ds' = V.zipWith4 (interp x) uLo uHi                                             ds (V.tail ds) +interp :: (Fractional (Scalar v), Ord (Scalar v), VectorSpace v)+    => Scalar v -> Scalar v -> Scalar v -> v -> v -> v interp x x0 x1 y0 y1     |  x <  x0  = y0     |  x >= x1  = y1     | otherwise = lerp y0 y1 a     where         a = (x - x0) / (x1 - x0)++-- "slice" a spline to contain only those knots and control points that +-- actually influence the value at 'x'.+--+-- It should be true for any valid BSpline that:+-- degree (slice f x) == degree f+-- slice (slice f x) x == slice f x+-- {x in domain of f} => {x in domain of slice f x}+-- {x in domain of f} => evalBSpline (slice f x) x == evalBSpline f x+{-# INLINE slice #-}+slice :: (Num (Scalar a), Ord (Scalar a), AdditiveGroup a, V.Vector v a)+     => BSpline v a -> Scalar a -> BSpline v a+slice spline x = spline+    { knotVector    = stakeKnots (n + n) . sdropKnots (l - n) $ knotVector spline+    , controlPoints = vtake       n      . vdrop      (l - n) $ controlPoints spline+    }+    where+        l = maybe 0 id $ V.findIndex (> x) us+        n = degree spline + 1+        +        us = knotsVector (knotVector spline)++-- Try to take n, but if there's not enough, pad the rest with 0s+vtake :: (V.Vector v t, AdditiveGroup t) => Int -> v t -> v t+{-# SPECIALIZE vtake :: AdditiveGroup t => Int -> BV.Vector t -> BV.Vector t #-}+vtake n xs+    | n <= V.length xs = V.take n xs+    | otherwise = xs V.++ V.replicate (n - V.length xs) zeroV++-- Try to drop n, but if n is negative, pad the beginning with 0s+vdrop :: (V.Vector v t, AdditiveGroup t) => Int -> v t -> v t+{-# SPECIALIZE vdrop :: AdditiveGroup t => Int -> BV.Vector t -> BV.Vector t #-}+vdrop n xs+    | n >= 0 = V.drop n xs+    | otherwise = V.replicate (-n) zeroV V.++ xs++-- Try to take n knots, but if there aren't enough, increase the multiplicity of the last knot+stakeKnots :: (Num k, Ord k) => Int -> Knots k -> Knots k+stakeKnots n kts+    | n <= nKts = takeKnots n kts+    | otherwise = case maxKnot kts of+        Nothing     -> multipleKnot 0 (n - nKts)+        Just (k, m) -> setKnotMultiplicity k (m + n - nKts) kts+    where nKts = numKnots kts++-- Try to drop n knots, but if n is negative, increase the multiplicity of the first knot by @abs n@+sdropKnots :: (Num k, Ord k) => Int -> Knots k -> Knots k+sdropKnots n kts+    | n >= 0    = dropKnots n kts+    | otherwise = case minKnot kts of+        Nothing     -> multipleKnot 0 (-n)+        Just (k, m) -> setKnotMultiplicity k (m - n) kts 

@@ -1,4 +1,6 @@+{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ParallelListComp #-}+{-# LANGUAGE FlexibleContexts #-} -- |Reference implementation of B-Splines; very inefficient but \"obviously\" -- correct. module Math.Spline.BSpline.Reference@@ -6,15 +8,40 @@     , basisFunctions     , basisPolynomials     , basisPolynomialsAt+    , evalReferenceBSpline+    , fitPolyToBSplineAt     ) where +import qualified Data.Vector.Generic as V+import Data.VectorSpace (VectorSpace, Scalar, (^*), sumV) import Math.Spline.Knots+import Math.Spline.BSpline.Internal import Math.Polynomial (Poly) import qualified Math.Polynomial as Poly +-- | This is a fairly slow function which computes the value of a B-spline at a given point,+-- using the mathematical definition of B-splines.  It is mainly for testing purposes, as a+-- reference against which the other evaluation functions are checked.+evalReferenceBSpline :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a), V.Vector v a) +    => BSpline v a -> Scalar a -> a+evalReferenceBSpline (Spline deg kts cps) x =+    sumV (zipWith (^*) (V.toList cps) (bases kts x !! deg))++-- | This is a fairly slow function which computes one polynomial segment of a B-spline (the +-- one containing the given point), using the mathematical definition of B-splines.  It is +-- mainly for testing purposes, as a reference against which the other evaluation functions+-- are checked.+fitPolyToBSplineAt :: (Fractional a, Ord a, Scalar a ~ a, V.Vector v a)+    => BSpline v a -> a -> Poly a+fitPolyToBSplineAt (Spline deg kts cps) x = +    Poly.sumPolys (zipWith Poly.scalePoly (V.toList cps) (basisPolynomialsAt kts x !! deg))++ind :: Num a => Bool -> a ind True  = 1 ind False = 0 +-- | The values of all the B-spline basis functions for the given knot vector at the given+-- point, ordered by degree; \"b_{i,j}(x)\" is @bases kts x !! i !! j@. bases :: (Fractional a, Ord a) => Knots a -> a -> [[a]] bases kts x = coxDeBoor interp initial kts     where@@ -26,8 +53,8 @@             = (if d0 == 0 then 0 else (x       - t_j) / d0) * b_nm1_j             + (if d1 == 0 then 0 else (t_jpnp1 -   x) / d1) * b_nm1_jp1 --- Alternate version constructing table of functions rather than computing--- table of values+-- | All the B-spline basis functions for the given knot vector at the given+-- point, ordered by degree; \"b_{i,j}\" is @basisFunctions kts x !! i !! j@. basisFunctions :: (Fractional a, Ord a) => Knots a -> [[a -> a]] basisFunctions kts = coxDeBoor interp initial kts     where@@ -39,13 +66,15 @@             = (if d0 == 0 then 0 else (x       - t_j) / d0) * b_nm1_j   x             + (if d1 == 0 then 0 else (t_jpnp1 -   x) / d1) * b_nm1_jp1 x --- compute all the basis polynomials for a knot vector, ordered by knot span.+-- | All the B-spline basis polynomials for the given knot vector, ordered first +-- by knot span and then by degree. basisPolynomials :: (Fractional a, Ord a) => Knots a -> [[[Poly a]]] basisPolynomials kts     | isEmpty kts   = []     | otherwise     = [basisPolynomialsAt kts kt | kt <- init (distinctKnots kts)] --- compute all the basis polynomials for the knot span containing a given location.+-- | All the B-spline basis polynomials for the given knot vector at the given+-- point, ordered by degree; \"b_{i,j}\" is @basisPolynomialsAt kts x !! i !! j@. basisPolynomialsAt :: (Fractional a, Ord a) => Knots a -> a -> [[Poly a]] basisPolynomialsAt kts x = coxDeBoor interp initial kts     where@@ -57,22 +86,23 @@             | (t_j, t_jp1) <- knotSpans kts 1             ]         interp t_j d0 b_nm1_j t_jpnp1 d1 b_nm1_jp1-            = (if d0 == 0 then Poly.zero else (Poly.x                 - Poly.constPoly t_j) / d0) * b_nm1_j-            + (if d1 == 0 then Poly.zero else (Poly.constPoly t_jpnp1 -             Poly.x) / d1) * b_nm1_jp1+            =   (if d0 == 0 then Poly.zero else (Poly.x                 ^-^ Poly.constPoly t_j) ^/ d0) ^*^ b_nm1_j+            ^+^ (if d1 == 0 then Poly.zero else (Poly.constPoly t_jpnp1 ^-^             Poly.x) ^/ d1) ^*^ b_nm1_jp1             where-                infixl 6 +, --                p + q   = Poly.addPoly p q-                p - q   = p + (Poly.negatePoly q)+                infixl 6 ^+^, ^-^+                p ^+^ q   = Poly.addPoly p q+                p ^-^ q   = p ^+^ (Poly.negatePoly q)                 -                infixl 7 *, /-                p * q   = Poly.multPoly p q-                p / s   = Poly.scalePoly (recip s) p+                infixl 7 ^*^, ^/+                p ^*^ q   = Poly.multPoly p q+                p ^/  s   = Poly.scalePoly (recip s) p --- This is a straightforward implementation of the Cox-De Boor recursion scheme+-- | This is a straightforward implementation of the Cox-De Boor recursion scheme -- generalized in a slightly strange way; the initial vector is a parameter  -- and the actual computation of the recursion step is a function parameter. -- The purpose is to allow the same recursion to be applied when computing basis -- function values and  basis polynomials.+coxDeBoor :: Num a => (a -> a -> b -> a -> a -> b -> b) -> [b] -> Knots a -> [[b]] coxDeBoor interp initial kts = table     where         ts = knots kts@@ -90,4 +120,6 @@ spans     = spansWith (,) spanDiffs :: Num a => Int -> [a] -> [a] spanDiffs = spansWith subtract++spansWith :: (a -> a -> b) -> Int -> [a] -> [b] spansWith f n ts = zipWith f ts (drop n ts)

@@ -12,7 +12,7 @@ import qualified Data.Vector as V import Data.VectorSpace --- |A BezierCurve curve on @0 <= x <= 1@.+-- |A Bezier curve on @0 <= x <= 1@. data BezierCurve t = BezierCurve !Int !(V.Vector t) deriving (Eq, Ord)  -- |Construct a Bezier curve from a list of control points.  The degree

@@ -1,4 +1,5 @@ {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts, UndecidableInstances #-}+{-# LANGUAGE OverlappingInstances, IncoherentInstances #-} module Math.Spline.Class where  import Control.Applicative@@ -6,6 +7,7 @@ import qualified Math.Spline.BSpline.Internal as BSpline  import qualified Data.Vector as V+import qualified Data.Vector.Generic as G import Data.VectorSpace  -- |A spline is a piecewise polynomial vector-valued function.  The necessary@@ -30,16 +32,27 @@     knotVector :: s v -> Knots (Scalar v)     knotVector = knotVector . toBSpline     -    toBSpline :: s v -> BSpline.BSpline v+    toBSpline :: s v -> BSpline.BSpline V.Vector v +-- TODO: this class should probably go away.  all it really does is overload something that doesn't really have any implementation-independent semantics (or does it?). class Spline s v => ControlPoints s v where     controlPoints :: s v -> V.Vector v -instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline BSpline.BSpline v where-    evalSpline spline = V.head . last . BSpline.deBoor spline+instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline (BSpline.BSpline V.Vector) v where+    evalSpline = BSpline.evalBSpline     splineDegree = BSpline.degree     knotVector = BSpline.knotVector     toBSpline = id -instance Spline BSpline.BSpline v => ControlPoints BSpline.BSpline v where+instance ( VectorSpace a, Fractional (Scalar a), Ord (Scalar a), G.Vector v a+         , G.Vector v (Scalar a)) => Spline (BSpline.BSpline v) a where+    evalSpline = BSpline.evalBSpline+    splineDegree = BSpline.degree+    knotVector = BSpline.knotVector+    toBSpline (BSpline.Spline deg ks ctp) = BSpline.Spline deg ks (G.convert $ ctp)++instance Spline (BSpline.BSpline V.Vector) a => ControlPoints (BSpline.BSpline V.Vector) a where     controlPoints = BSpline.controlPoints++instance (Spline (BSpline.BSpline v) a, G.Vector v a) => ControlPoints (BSpline.BSpline v) a where+    controlPoints = V.convert . BSpline.controlPoints

@@ -0,0 +1,89 @@+{-# LANGUAGE ParallelListComp #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+module Math.Spline.Hermite+    ( CSpline, cSpline+    , evalSpline+    ) where++import Data.List+import Data.Ord+import Math.Polynomial hiding (x)+import Math.Spline.BSpline+import Math.Spline.Class+import Math.Spline.Knots+import qualified Data.Vector as V+import Data.VectorSpace++-- | Cubic Hermite splines.  These are cubic splines defined by a +-- sequence of control points and derivatives at those points.+newtype CSpline a = CSpline [(Scalar a,a,a)]++-- | Cubic splines specified by a list of control points, +-- where each control point is given by a triple of parameter value, +-- position of the spline at that parameter value,+-- and derivative of the spline at that parameter value.+cSpline :: Ord (Scalar a) => [(Scalar a,a,a)] -> CSpline a+cSpline = CSpline . sortBy (comparing fst3)+    where fst3 (a,_,_) = a++h00, h10, h01, h11 :: (Num a, Eq a) => Poly a+h00 = poly BE [ 2,-3, 0, 1]+h10 = poly BE [ 1,-2, 1, 0]+h01 = poly BE [-2, 3, 0, 0]+h11 = poly BE [ 1,-1, 0, 0]++evalHermite+  :: (Eq (Scalar v), Num (Scalar v), VectorSpace v) =>+     v -> v -> v -> v -> Scalar v -> v+evalHermite y0 m0 y1 m1 x = sumV+    [ evalPoly h x *^ p+    | p <- [ y0,  m0,  y1,  m1]+    | h <- [h00, h10, h01, h11]+    ]++evalCSpline+  :: (VectorSpace v, Ord (Scalar v), Fractional (Scalar v)) =>+     CSpline v -> Scalar v -> v+evalCSpline (CSpline cps) = loop cps+    where+        loop []         _ = zeroV+        loop [(_,y,_)]  _ = y+        loop ((x0,y0,m0):rest@((x1,y1,m1):_)) x+            | x <= x0   = y0+            | x <  x1   = let dx = x1 - x0+                           in evalHermite y0 (m0 ^* dx) y1 (m1 ^* dx) ((x - x0) / dx)+            | otherwise = loop rest x++instance (VectorSpace a, Fractional (Scalar a), Ord (Scalar a)) => Spline CSpline a where+    splineDegree _ = 3+    +    splineDomain (CSpline  []) = Nothing+    splineDomain (CSpline cps) = Just (head xs, last xs) where (xs, _, _) = unzip3 cps+    +    evalSpline = evalCSpline+    +    -- TODO: check.  Also work out a more compact translation, taking advantage of the+    -- known continuity on the interior knots.  It should be possible to work out an +    -- equivalent b-spline with only 'n+4' knots.  If that translation isn't ill-+    -- conditioned, it might be a good thing to implement.+    toBSpline (CSpline cSpl) = bSpline kts (V.fromList cps)+        where +            kts = fromList [(x,4) | x <- xs]+            +            (xs, _, _) = unzip3 cSpl+            cps = concat +                [ [ y0+                  , y0 ^+^ dy0 ^* dx3+                  , y1 ^-^ dy1 ^* dx3+                  , y1+                  ]+                | ((x0, y0, dy0), (x1, y1, dy1)) <- spans 1 cSpl+                , let dx3 = (x1 - x0) / 3+                ]++spans :: Int -> [a] -> [(a,a)]+spans n xs = zip xs (drop n xs)+

@@ -62,6 +62,8 @@ toISpline :: (Spline s v, Eq v) => s v -> ISpline v toISpline = fromBSpline . toBSpline +fromBSpline :: (Eq v, VectorSpace v, Fractional (Scalar v), Ord (Scalar v))+    => BSpline V.Vector v -> ISpline v fromBSpline spline     | V.head ds == zeroV      && numKnots ks >= 2 = iSpline (mkKnots (init (tail ts))) (V.tail ds')

@@ -10,7 +10,9 @@     , toList, numDistinctKnots, lookupDistinctKnot          , knots, knotsVector-    , distinctKnots, distinctKnotsVector+    , distinctKnots, multiplicities+    , distinctKnotsVector, multiplicitiesVector+    , distinctKnotsSet          , toMap     , fromMap@@ -25,6 +27,8 @@     , maxMultiplicity     , knotMultiplicity, setKnotMultiplicity     +    , splitFind+         , fromAscList, fromDistinctAscList     , valid     @@ -34,27 +38,26 @@     , knotDomain          , uniform+    +    , minKnot+    , maxKnot     ) where  import Prelude hiding (sum, maximum)+import Control.Arrow ((***)) import Control.Monad (guard)-import Data.Foldable (Foldable(foldMap), sum, maximum)+import Data.Foldable (Foldable(foldMap), maximum)+import Data.List (sortBy, sort, unfoldr) import qualified Data.Map as M import Data.Monoid (Monoid(..))-import Data.Maybe (fromMaybe)-import qualified Data.Set as S (Set)+import Data.Ord+import qualified Data.Set as S (Set, fromAscList) import qualified Data.Vector as V-import Data.VectorSpace  -- |Knot vectors - multisets of points in a 1-dimensional space.-data Knots a = Knots !Int (M.Map a Int) deriving (Eq, Ord)+newtype Knots a = Knots (V.Vector a) 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)@@ -62,20 +65,19 @@  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))+    mappend (Knots v1) (Knots v2) =+      Knots . V.fromList . sort . V.toList $ v1 V.++ v2  instance Foldable Knots where-    foldMap f = foldMap f . knots+    foldMap f = foldMap f . knotsVector   -- |An empty knot vector empty :: Knots a-empty = Knots 0 M.empty+empty = Knots V.empty  isEmpty :: Knots a -> Bool-isEmpty (Knots 0 _) = True-isEmpty  _          = False+isEmpty (Knots v) = V.null v  -- |Create a knot vector consisting of one knot. knot :: Ord a => a -> Knots a@@ -83,79 +85,74 @@  -- |Create a knot vector consisting of one knot with the specified multiplicity. multipleKnot :: Ord a => a -> Int -> Knots a-multipleKnot k n -    | n <= 0    = Knots 0 (M.empty)-    | otherwise = Knots n (M.singleton k n)+multipleKnot k n = Knots $ V.replicate n k  -- |Create a knot vector consisting of all the knots in a list. mkKnots :: (Ord a) => [a] -> Knots a-mkKnots ks = fromList (map (\k -> (k,1)) ks)+mkKnots = Knots . V.fromList . sort  -- |Create a knot vector consisting of all the knots and corresponding  -- multiplicities in a list. fromList :: (Ord k) => [(k, Int)] -> Knots k-fromList ks = Knots (sum kMap) kMap-    where kMap = M.fromListWith (+) (filter ((>0).snd) ks)+fromList ks = Knots v+    where v = V.concat . map (\(k, mult) -> V.replicate mult k) .+              sortBy (comparing fst) $ filter ((>0).snd) ks  -- |Create a knot vector consisting of all the knots and corresponding  -- multiplicities in a list ordered by the knots' 'Ord' instance.  The -- ordering precondition is not checked. fromAscList :: Eq k => [(k, Int)] -> Knots k-fromAscList ks = Knots (sum kMap) kMap-    where kMap = M.fromAscListWith (+) (filter ((>0).snd) ks)+fromAscList ks = Knots v+    where v = V.concat . map (\(k, mult) -> V.replicate mult k)+                 $ filter ((>0).snd) ks  -- |Create a knot vector consisting of all the knots and corresponding  -- multiplicities in a list ordered by the knots' 'Ord' instance with no -- duplicates.  The preconditions are not checked.-fromDistinctAscList :: [(k, Int)] -> Knots k-fromDistinctAscList ks = Knots (sum kMap) kMap-    where kMap = M.fromDistinctAscList (filter ((>0).snd) ks)+fromDistinctAscList :: Eq k => [(k, Int)] -> Knots k+fromDistinctAscList = fromAscList -fromMap :: M.Map k Int -> Knots k-fromMap ks = Knots (sum kMap) kMap-    where-        kMap = mFilter (>0) ks-        -- filter is monotonic, I have no idea why M.filter requires Ord on the key-        mFilter p = M.fromDistinctAscList . filter (p.snd) . M.toAscList+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  -- |Returns a list of all distinct knots in ascending order along with -- their multiplicities.-toList :: Knots k -> [(k, Int)]-toList = M.toList . toMap+toList :: Eq k => Knots k -> [(k, Int)]+toList = unfoldr $ \kts -> do+    kt <- minKnot kts+    return (kt, dropKnots (snd kt) kts) -toVector :: Knots k -> V.Vector (k, Int)-toVector = V.fromList . toList+toVector :: Eq k => Knots k -> V.Vector (k, Int)+toVector = V.unfoldr $ \kts -> do+    kt <- minKnot kts+    return (kt, dropKnots (snd kt) kts) -toMap :: Knots k -> M.Map k Int-toMap (Knots _ ks) = ks+toMap :: Ord k => Knots k -> M.Map k Int+toMap = M.fromAscListWith (+) . toList  -- |Returns the number of knots (not necessarily distinct) in a knot vector. numKnots :: Knots t -> Int-numKnots (Knots n _) = n+numKnots (Knots v) = V.length v  -- |Returns the number of distinct knots in a knot vector.-numDistinctKnots :: Knots t -> Int-numDistinctKnots (Knots _ ks) = M.size ks+numDistinctKnots :: Eq t => Knots t -> Int+numDistinctKnots = V.length . distinctKnotsVector -maxMultiplicity :: Knots t -> Int-maxMultiplicity (Knots 0  _) = 0-maxMultiplicity (Knots _ ks) = maximum ks+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 kts-    | k < 0             = Nothing-    | k < numKnots kts  = fmap fst mbKt-    | otherwise         = Nothing-    where (_, mbKt, _) = splitLookup k kts+lookupKnot k (Knots kts) = do+    guard (0 <= k && k < V.length kts)+    V.indexM kts k -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+lookupDistinctKnot :: Eq a => Int -> Knots a -> Maybe a+lookupDistinctKnot k kts = lookupKnot k . Knots $ distinctKnotsVector kts  -- |@splitLookup n kts@: Split a knot vector @kts@ into 3 parts @(pre, mbKt, post)@ -- such that:@@ -165,133 +162,84 @@ --  * Putting the 3 parts back together yields exactly the original knot vector --  * The @n@'th knot, if one exists, will be in @mbKt@ along with its multiplicity ---splitLookup :: Int -> Knots a -> (Knots a, Maybe (a, Int), Knots a)-splitLookup k (Knots n ks) = scan 0 M.empty n ks-    where-        -- The general plan: iteratively pull the smallest knot out of "post",-        -- either moving it to "pre" or terminating by returning it along with-        -- current values of "pre" and "post"-        -        -- invariants:-        --   nPre  = sum pre-        --   nPost = sum post-        --   M.union pre post = ks-        --   every key in pre < every key in post-        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)---- Prepend or append an element to a map, without checking the precondition--- that the new pair's key is less than (greater than, resp.) all keys in --- the map.-ascCons x m = M.fromDistinctAscList (x : M.toAscList m)-ascSnoc m x = M.fromDistinctAscList (M.toAscList m ++ [x])---- Prepend or append an knot to a knot vector, without checking the--- precondition that the new knot's location is less than (greater than,--- resp.) all knots in the vector.-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+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 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+dropKnots k (Knots v) = Knots $ V.drop k v  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+takeKnots k (Knots v) = Knots $ V.take k v  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+splitKnotsAt k (Knots v) = Knots *** Knots $ V.splitAt k v +-- |Count the number of knots less than the n'th distinct knot.+findDistinctKnot :: Eq a => Int -> Knots a -> Int+findDistinctKnot n = V.last . V.take (1 + max 0 n) . V.scanl (+) 0 . multiplicitiesVector -takeDistinctKnots :: Int -> Knots a -> Knots a-takeDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap-    where-        kMap = M.fromDistinctAscList (take k (M.toAscList ks))+takeDistinctKnots :: (Ord a) => Int -> Knots a -> Knots a+takeDistinctKnots k kts = takeKnots (findDistinctKnot k kts) kts -dropDistinctKnots :: Int -> Knots a -> Knots a-dropDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap-    where-        kMap = M.fromDistinctAscList (drop k (M.toAscList ks))+dropDistinctKnots :: (Ord a) => Int -> Knots a -> Knots a+dropDistinctKnots k kts = dropKnots (findDistinctKnot k kts) kts -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+splitDistinctKnotsAt :: (Ord a, Eq a) => Int -> Knots a -> (Knots a, Knots a)+splitDistinctKnotsAt k kts = splitKnotsAt (findDistinctKnot k kts) kts  -- |Returns a list of all knots (not necessarily distinct) of a knot vector in ascending order knots :: Knots t -> [t]-knots (Knots _ ks) = concat [replicate n k | (k,n) <- M.toAscList ks]+knots = V.toList . knotsVector  -- |Returns a vector of all knots (not necessarily distinct) of a knot vector in ascending order knotsVector :: Knots t -> V.Vector t-knotsVector (Knots _ ks) = V.concat [V.replicate n k | (k,n) <- M.toAscList ks]+knotsVector (Knots v) = v  -- |Returns a list of all distinct knots of a knot vector in ascending order-distinctKnots :: Knots t -> [t]-distinctKnots (Knots _ ks) = M.keys ks+distinctKnots :: Eq t => Knots t -> [t]+distinctKnots = map fst . toList +multiplicities :: Eq t => Knots t -> [Int]+multiplicities = map snd . toList+ -- |Returns a vector of all distinct knots of a knot vector in ascending order-distinctKnotsVector :: Knots t -> V.Vector t-distinctKnotsVector = V.fromList . distinctKnots+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+ -- |Returns a 'S.Set' of all distinct knots of a knot vector-distinctKnotsSet :: Knots k -> S.Set k-distinctKnotsSet (Knots _ ks) = M.keysSet ks+distinctKnotsSet :: Eq k => Knots k -> S.Set k+distinctKnotsSet (Knots k) = S.fromAscList $ V.toList k  -- |Looks up the multiplicity of a knot (which is 0 if the point is not a knot) knotMultiplicity :: (Ord k) => k -> Knots k -> Int-knotMultiplicity k (Knots _ ks) = fromMaybe 0 (M.lookup k ks)+knotMultiplicity k (Knots ks) = V.length $ V.elemIndices k ks  -- |Returns a new knot vector with the given knot set to the specified  -- multiplicity and all other knots unchanged. 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)+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+ -- |Check the internal consistency of a knot vector-valid :: Ord k => Knots k -> Bool-valid (Knots n ks) = and-    [ M.valid ks-    , n == sum ks-    , all (>0) (M.elems ks)-    ]+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 kts i j@ returns the knot span extending from the @i@'th knot -- to the @j@'th knot, if  @i <= j@ and both knots exist.@@ -324,7 +272,7 @@ -- the basis functions sum to 1, which is only true on this range, and so -- this is also precisely the domain on which de Boor's algorithm is valid. knotDomain :: Knots a -> Int -> Maybe (a,a)-knotDomain ks@(Knots n _) p = knotSpan ks p (n-p-1)+knotDomain ks@(Knots v) p = knotSpan ks p (V.length v-p-1)  -- |@uniform deg nPts (lo,hi)@ constructs a uniformly-spaced knot vector over -- the interval from @lo@ to @hi@ which, when used to construct a B-spline @@ -336,3 +284,14 @@         n = nPts + deg - numKnots ends         f i = (fromIntegral i * lo + fromIntegral (n - i) * hi) / fromIntegral n         internal = mkKnots [f i | i <- [0..n]]++{-# INLINE minKnot #-}+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++{-# INLINE maxKnot #-}+maxKnot :: Eq a => Knots a -> Maybe (a, Int)+maxKnot (Knots v) = minKnot (Knots (V.reverse v))

@@ -46,6 +46,7 @@         n = V.length cps         m = numKnots kts - 1 +spans :: Int -> V.Vector a -> V.Vector (a,a) spans n xs = V.zip xs (V.drop n xs)  instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline MSpline v where@@ -63,6 +64,8 @@ toMSpline :: Spline s v => s v -> MSpline v toMSpline = fromBSpline . toBSpline +fromBSpline :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a))+    => BSpline (V.Vector) a -> MSpline a fromBSpline spline = mSpline ks cs     where         n = splineDegree spline + 1; n' = fromIntegral n

@@ -0,0 +1,14 @@+#!/usr/bin/env runhaskell+module Main where++import Test.Framework (defaultMain, testGroup)++import Tests.BSpline.Reference (referenceBSplineTests)+import Tests.BSpline (bSplineTests)+import Tests.Knots (knotsTests)++main = defaultMain +    [ testGroup "Math.Spline.BSpline.Reference" referenceBSplineTests+    , testGroup "Math.Spline.BSpline"           bSplineTests+    , testGroup "Math.Spline.Knots"             knotsTests+    ]