{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TypeFamilies  #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.TwoD.Segment.Bernstein
-- Copyright   :  (c) 2014-2015 diagrams-lib team (see LICENSE)
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  diagrams-discuss@googlegroups.com
--
-- Bernstein polynomials, used internally by code to find
-- intersections of paths.  This module is probably not of any
-- relevance to most users of diagrams.
-----------------------------------------------------------------------------
module Diagrams.TwoD.Segment.Bernstein
  ( BernsteinPoly (..)
  , listToBernstein
  , evaluateBernstein

  , degreeElevate
  , bernsteinDeriv
  , evaluateBernsteinDerivs
  ) where

import           Data.List           (tails)
import           Diagrams.Core.V
import           Diagrams.Parametric
import           Linear.V1

-- | Compute the binomial coefficients of degree n.
binomials :: Num n => Int -> [n]
binomials :: forall n. Num n => Int -> [n]
binomials Int
n = forall a b. (a -> b) -> [a] -> [b]
map forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\Int
x Int
m -> Int
x forall a. Num a => a -> a -> a
* (Int
n forall a. Num a => a -> a -> a
- Int
mforall a. Num a => a -> a -> a
+Int
1) forall a. Integral a => a -> a -> a
`quot` Int
m) Int
1 [Int
1..Int
n]

data BernsteinPoly n = BernsteinPoly
  { forall n. BernsteinPoly n -> Int
bernsteinDegree :: Int
  , forall n. BernsteinPoly n -> [n]
bernsteinCoeffs :: [n]
  } deriving (Int -> BernsteinPoly n -> ShowS
forall n. Show n => Int -> BernsteinPoly n -> ShowS
forall n. Show n => [BernsteinPoly n] -> ShowS
forall n. Show n => BernsteinPoly n -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [BernsteinPoly n] -> ShowS
$cshowList :: forall n. Show n => [BernsteinPoly n] -> ShowS
show :: BernsteinPoly n -> String
$cshow :: forall n. Show n => BernsteinPoly n -> String
showsPrec :: Int -> BernsteinPoly n -> ShowS
$cshowsPrec :: forall n. Show n => Int -> BernsteinPoly n -> ShowS
Show, forall a b. a -> BernsteinPoly b -> BernsteinPoly a
forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
<$ :: forall a b. a -> BernsteinPoly b -> BernsteinPoly a
$c<$ :: forall a b. a -> BernsteinPoly b -> BernsteinPoly a
fmap :: forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
$cfmap :: forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
Functor)

type instance V        (BernsteinPoly n) = V1
type instance N        (BernsteinPoly n) = n
type instance Codomain (BernsteinPoly n) = V1

-- | Create a bernstein polynomial from a list of coëfficients.
listToBernstein :: Fractional n => [n] -> BernsteinPoly n
listToBernstein :: forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein [] = BernsteinPoly n
0
listToBernstein [n]
l  = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly (forall (t :: * -> *) a. Foldable t => t a -> Int
length [n]
l forall a. Num a => a -> a -> a
- Int
1) [n]
l

-- | Degree elevate a bernstein polynomial a number of times.
degreeElevate :: Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate :: forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
b                    Int
0     = BernsteinPoly n
b
degreeElevate (BernsteinPoly Int
lp [n]
p) Int
times =
  forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate (forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly (Int
lpforall a. Num a => a -> a -> a
+Int
1) (forall a. [a] -> a
head [n]
pforall a. a -> [a] -> [a]
:[n] -> n -> [n]
inner [n]
p n
1)) (Int
timesforall a. Num a => a -> a -> a
-Int
1)
  where
    n :: n
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
lp

    inner :: [n] -> n -> [n]
inner []         n
_ = [n
0]
    inner [n
a]        n
_ = [n
a]
    inner (n
a:n
b:[n]
rest) n
i = (n
iforall a. Num a => a -> a -> a
*n
aforall a. Fractional a => a -> a -> a
/(n
nforall a. Num a => a -> a -> a
+n
1) forall a. Num a => a -> a -> a
+ n
bforall a. Num a => a -> a -> a
*(n
1 forall a. Num a => a -> a -> a
- n
iforall a. Fractional a => a -> a -> a
/(n
nforall a. Num a => a -> a -> a
+n
1))) forall a. a -> [a] -> [a]
: [n] -> n -> [n]
inner (n
bforall a. a -> [a] -> [a]
:[n]
rest) (n
iforall a. Num a => a -> a -> a
+n
1)

-- | Evaluate the bernstein polynomial.
evaluateBernstein :: Fractional n => BernsteinPoly n -> n -> n
evaluateBernstein :: forall n. Fractional n => BernsteinPoly n -> n -> n
evaluateBernstein (BernsteinPoly Int
_ [])       n
_ = n
0
evaluateBernstein (BernsteinPoly Int
_ [n
b])      n
_ = n
b
evaluateBernstein (BernsteinPoly Int
lp (n
b':[n]
bs)) n
t = n -> n -> n -> n -> [n] -> n
go n
t n
n (n
b'forall a. Num a => a -> a -> a
*n
u) n
2 [n]
bs
  where
    u :: n
u = n
1forall a. Num a => a -> a -> a
-n
t
    n :: n
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
lp

    go :: n -> n -> n -> n -> [n] -> n
go n
tn n
bc n
tmp n
_ [n
b]      = n
tmp forall a. Num a => a -> a -> a
+ n
tnforall a. Num a => a -> a -> a
*n
bcforall a. Num a => a -> a -> a
*n
b
    go n
tn n
bc n
tmp n
i (n
b:[n]
rest) =
      n -> n -> n -> n -> [n] -> n
go (n
tnforall a. Num a => a -> a -> a
*n
t)              -- tn
         (n
bcforall a. Num a => a -> a -> a
*(n
n forall a. Num a => a -> a -> a
- n
iforall a. Num a => a -> a -> a
+n
1)forall a. Fractional a => a -> a -> a
/n
i)    -- bc
         ((n
tmp forall a. Num a => a -> a -> a
+ n
tnforall a. Num a => a -> a -> a
*n
bcforall a. Num a => a -> a -> a
*n
b)forall a. Num a => a -> a -> a
*n
u) -- tmp
         (n
iforall a. Num a => a -> a -> a
+n
1)               -- i
         [n]
rest
    go n
_ n
_ n
_ n
_ []           = forall a. HasCallStack => String -> a
error String
"evaluateBernstein: impossible"

-- | Evaluate the bernstein polynomial and its derivatives.
evaluateBernsteinDerivs :: Fractional n => BernsteinPoly n -> n -> [n]
evaluateBernsteinDerivs :: forall n. Fractional n => BernsteinPoly n -> n -> [n]
evaluateBernsteinDerivs BernsteinPoly n
b n
t
  | forall n. BernsteinPoly n -> Int
bernsteinDegree BernsteinPoly n
b forall a. Eq a => a -> a -> Bool
== Int
0 = [forall n. Fractional n => BernsteinPoly n -> n -> n
evaluateBernstein BernsteinPoly n
b n
t]
  | Bool
otherwise              = forall n. Fractional n => BernsteinPoly n -> n -> n
evaluateBernstein BernsteinPoly n
b n
t forall a. a -> [a] -> [a]
: forall n. Fractional n => BernsteinPoly n -> n -> [n]
evaluateBernsteinDerivs (forall n. Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv BernsteinPoly n
b) n
t

-- | Find the derivative of a bernstein polynomial.
bernsteinDeriv :: Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv :: forall n. Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv (BernsteinPoly Int
0 [n]
_)  = BernsteinPoly n
0
bernsteinDeriv (BernsteinPoly Int
lp [n]
p) =
  -- BernsteinPoly (lp-1) $ map (* fromIntegral lp) $ zipWith (-) (tail p) p
  forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly (Int
lpforall a. Num a => a -> a -> a
-Int
1) forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\n
a n
b -> (n
a forall a. Num a => a -> a -> a
- n
b) forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
lp) (forall a. [a] -> [a]
tail [n]
p) [n]
p

instance Fractional n => Parametric (BernsteinPoly n) where
  atParam :: BernsteinPoly n
-> N (BernsteinPoly n)
-> Codomain (BernsteinPoly n) (N (BernsteinPoly n))
atParam BernsteinPoly n
b = forall a. a -> V1 a
V1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall n. Fractional n => BernsteinPoly n -> n -> n
evaluateBernstein BernsteinPoly n
b
instance Num n        => DomainBounds (BernsteinPoly n)
instance Fractional n => EndValues    (BernsteinPoly n)
instance Fractional n => Sectionable  (BernsteinPoly n) where
  splitAtParam :: BernsteinPoly n
-> N (BernsteinPoly n) -> (BernsteinPoly n, BernsteinPoly n)
splitAtParam  = forall n.
Num n =>
BernsteinPoly n -> n -> (BernsteinPoly n, BernsteinPoly n)
bernsteinSplit
  reverseDomain :: BernsteinPoly n -> BernsteinPoly n
reverseDomain (BernsteinPoly Int
i [n]
xs) = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
i (forall a. [a] -> [a]
reverse [n]
xs)

-- | Split a bernstein polynomial.
bernsteinSplit :: Num n => BernsteinPoly n -> n -> (BernsteinPoly n, BernsteinPoly n)
bernsteinSplit :: forall n.
Num n =>
BernsteinPoly n -> n -> (BernsteinPoly n, BernsteinPoly n)
bernsteinSplit (BernsteinPoly Int
lp [n]
p) n
t =
  (forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
lp forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a. [a] -> a
head [[n]]
controls,
   forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
lp forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a. [a] -> a
last [[n]]
controls)
  where
    interp :: n -> n -> n
interp n
a n
b = (n
1forall a. Num a => a -> a -> a
-n
t)forall a. Num a => a -> a -> a
*n
a forall a. Num a => a -> a -> a
+ n
tforall a. Num a => a -> a -> a
*n
b

    terp :: [n] -> [[n]]
terp [n
_] = []
    terp [n]
l   = let ctrs :: [n]
ctrs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith n -> n -> n
interp [n]
l (forall a. [a] -> [a]
tail [n]
l)
               in  [n]
ctrs forall a. a -> [a] -> [a]
: [n] -> [[n]]
terp [n]
ctrs
    controls :: [[n]]
controls = [n]
p forall a. a -> [a] -> [a]
: [n] -> [[n]]
terp [n]
p

instance Fractional n => Num (BernsteinPoly n) where
  ba :: BernsteinPoly n
ba@(BernsteinPoly Int
la [n]
a) + :: BernsteinPoly n -> BernsteinPoly n -> BernsteinPoly n
+ bb :: BernsteinPoly n
bb@(BernsteinPoly Int
lb [n]
b)
    | Int
la forall a. Ord a => a -> a -> Bool
< Int
lb   = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
lb forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) (forall n. BernsteinPoly n -> [n]
bernsteinCoeffs forall a b. (a -> b) -> a -> b
$ forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
ba forall a b. (a -> b) -> a -> b
$ Int
lb forall a. Num a => a -> a -> a
- Int
la) [n]
b
    | Int
la forall a. Ord a => a -> a -> Bool
> Int
lb   = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
la forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) [n]
a (forall n. BernsteinPoly n -> [n]
bernsteinCoeffs forall a b. (a -> b) -> a -> b
$ forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
bb forall a b. (a -> b) -> a -> b
$ Int
la forall a. Num a => a -> a -> a
- Int
lb)
    | Bool
otherwise = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
la forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) [n]
a [n]
b

  ba :: BernsteinPoly n
ba@(BernsteinPoly Int
la [n]
a) - :: BernsteinPoly n -> BernsteinPoly n -> BernsteinPoly n
- bb :: BernsteinPoly n
bb@(BernsteinPoly Int
lb [n]
b)
    | Int
la forall a. Ord a => a -> a -> Bool
< Int
lb   = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
lb forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (forall n. BernsteinPoly n -> [n]
bernsteinCoeffs forall a b. (a -> b) -> a -> b
$ forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
ba (Int
lb forall a. Num a => a -> a -> a
- Int
la)) [n]
b
    | Int
la forall a. Ord a => a -> a -> Bool
> Int
lb   = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
la forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [n]
a (forall n. BernsteinPoly n -> [n]
bernsteinCoeffs forall a b. (a -> b) -> a -> b
$ forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
bb (Int
la forall a. Num a => a -> a -> a
- Int
lb))
    | Bool
otherwise = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
la forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [n]
a [n]
b

  (BernsteinPoly Int
la [n]
a) * :: BernsteinPoly n -> BernsteinPoly n -> BernsteinPoly n
* (BernsteinPoly Int
lb [n]
b) =
    forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly (Int
laforall a. Num a => a -> a -> a
+Int
lb) forall a b. (a -> b) -> a -> b
$
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Fractional a => a -> a -> a
(/)) (forall n. Num n => Int -> [n]
binomials (Int
la forall a. Num a => a -> a -> a
+ Int
lb)) forall a b. (a -> b) -> a -> b
$
                   forall a. [a] -> [a]
init forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$
                   forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [n]
a') (forall {a}. [a] -> [[a]]
down [n]
b') forall a. [a] -> [a] -> [a]
++
                   forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) (forall a. [a] -> [a]
reverse [n]
b')) (forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall {a}. [a] -> [[a]]
tails [n]
a')
                   -- zipWith (zipWith (*)) (tail $ tails a') (repeat $ reverse b')
    where down :: [a] -> [[a]]
down [a]
l = forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (forall a b c. (a -> b -> c) -> b -> a -> c
flip (:)) [] [a]
l -- [[1], [2, 1], [3, 2, 1], ...
          a' :: [n]
a' = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [n]
a (forall n. Num n => Int -> [n]
binomials Int
la)
          b' :: [n]
b' = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [n]
b (forall n. Num n => Int -> [n]
binomials Int
lb)

  fromInteger :: Integer -> BernsteinPoly n
fromInteger Integer
a = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
0 [forall a. Num a => Integer -> a
fromInteger Integer
a]

  signum :: BernsteinPoly n -> BernsteinPoly n
signum (BernsteinPoly Int
_ [])    = BernsteinPoly n
0
  signum (BernsteinPoly Int
_ (n
a:[n]
_)) = forall n. Int -> [n] -> BernsteinPoly n
BernsteinPoly Int
0 [forall a. Num a => a -> a
signum n
a]

  abs :: BernsteinPoly n -> BernsteinPoly n
abs = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Num a => a -> a
abs