-- | Not sure why I wanted these functions, but here they are.
module Why (
    fareyApproximate
  , divideOnIndex
  , bhattacharyyaDistance
  , bhattacharyyaCoefficient
  , hellingerDistance
) where



import           Data.Ratio
import           Data.Vector (Vector)
import qualified Data.Vector as V



-- | Approximate a rational number with one that has a maximum denominator.
--
-- This can be used to implement explitit rounding when building geometry with
-- rational numbers, which have the advantage of exact calculations, in contrast to
-- Double: the algorithm can be exact, and rounding can happen after the fact.
fareyApproximate
    :: Integer -- ^ Maximum denominator
    -> Ratio Integer -- ^ Number to approximate
    -> Ratio Integer
fareyApproximate :: Integer -> Ratio Integer -> Ratio Integer
fareyApproximate Integer
maxD Ratio Integer
exact = Ratio Integer -> Ratio Integer -> Ratio Integer
go Ratio Integer
0 (Ratio Integer -> Ratio Integer
forall {a}. Integral a => Ratio a -> Ratio a
ceilRatio Ratio Integer
exact)
  where
    ceilRatio :: Ratio a -> Ratio a
ceilRatio Ratio a
x = Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
x a -> a -> a
forall a. Num a => a -> a -> a
* (a -> a -> a
forall a. Integral a => a -> a -> a
quot (Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
x) (Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
x) a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) a -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
% a
1
    mediant :: Ratio a -> Ratio a -> Ratio a
mediant Ratio a
x Ratio a
y = (Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
x a -> a -> a
forall a. Num a => a -> a -> a
+ Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
y) a -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
% (Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
x a -> a -> a
forall a. Num a => a -> a -> a
+ Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
y)
    go :: Ratio Integer -> Ratio Integer -> Ratio Integer
go Ratio Integer
x Ratio Integer
y
      = let m :: Ratio Integer
m = Ratio Integer -> Ratio Integer -> Ratio Integer
forall {a}. Integral a => Ratio a -> Ratio a -> Ratio a
mediant Ratio Integer
x Ratio Integer
y
            dx :: Integer
dx = Ratio Integer -> Integer
forall a. Ratio a -> a
denominator Ratio Integer
x
            dy :: Integer
dy = Ratio Integer -> Integer
forall a. Ratio a -> a
denominator Ratio Integer
y
        in case Ratio Integer -> Ratio Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Ratio Integer
exact Ratio Integer
m of
            Ordering
_ | Integer
dx Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
maxD -> Ratio Integer
y
              | Integer
dy Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
maxD -> Ratio Integer
x
            Ordering
LT -> Ratio Integer -> Ratio Integer -> Ratio Integer
go Ratio Integer
m Ratio Integer
y
            Ordering
GT -> Ratio Integer -> Ratio Integer -> Ratio Integer
go Ratio Integer
x Ratio Integer
m
            Ordering
EQ | Integer
dx Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
dy Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Ratio Integer -> Integer
forall a. Ratio a -> a
denominator Ratio Integer
exact -> Ratio Integer
m
               | Integer
dy Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
dx -> Ratio Integer
y
               | Bool
otherwise -> Ratio Integer
x

-- | Split a 'V.Vector' at an index into the slice before, the value at the index, and the slice after.
--
-- >>> divideOnIndex 0 (V.fromList [0..10])
-- ([],0,[1,2,3,4,5,6,7,8,9,10])
--
-- >>> divideOnIndex 3 (V.fromList [0..10])
-- ([0,1,2],3,[4,5,6,7,8,9,10])
--
-- >>> divideOnIndex 10 (V.fromList [0..10])
-- ([0,1,2,3,4,5,6,7,8,9],10,[])
divideOnIndex :: Int -> V.Vector a -> (V.Vector a, a, V.Vector a)
divideOnIndex :: forall a. Int -> Vector a -> (Vector a, a, Vector a)
divideOnIndex Int
ix Vector a
vec
    | Int
ix Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
vec = [Char] -> (Vector a, a, Vector a)
forall a. HasCallStack => [Char] -> a
error ([Char]
"divideOnIndex: index out of bounds. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"i = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
ix [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
", length = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show (Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
vec))
divideOnIndex Int
ix Vector a
vec
  = let (Vector a
before, Vector a
after') = Int -> Vector a -> (Vector a, Vector a)
forall a. Int -> Vector a -> (Vector a, Vector a)
V.splitAt Int
ix Vector a
vec
        Just (a
pivot, Vector a
after) = Vector a -> Maybe (a, Vector a)
forall a. Vector a -> Maybe (a, Vector a)
V.uncons Vector a
after'
    in (Vector a
before, a
pivot, Vector a
after)

-- | The Bhattacharyya distance measures the similarity of two probability distributions.
--
-- I’ve added it here because I’m hoping to use it to implement a distance metric
-- between two pictures at some point.
--
-- \[
--   D_{B}(p,q)=-\ln \left(BC(p,q)\right) \\
--   BC(p,q)=\sum _{{x\in X}}{\sqrt {p(x)q(x)}}
-- \]
--
-- See https://en.wikipedia.org/wiki/Bhattacharyya_distance
bhattacharyyaDistance :: Vector Double -> Vector Double -> Double
bhattacharyyaDistance :: Vector Double -> Vector Double -> Double
bhattacharyyaDistance Vector Double
xs Vector Double
ys = - Double -> Double
forall a. Floating a => a -> a
log (Vector Double -> Vector Double -> Double
bhattacharyyaCoefficient Vector Double
xs Vector Double
ys)

-- | The Bhattacharyya coefficient is an approximate measurement of the amount of overlap between two statistical samples.
--
-- \[
--   BC(p,q)=\sum _{{x\in X}}{\sqrt {p(x)q(x)}}
-- \]
--
-- See https://en.wikipedia.org/wiki/Bhattacharyya_distance#Bhattacharyya_coefficient
bhattacharyyaCoefficient :: Vector Double -> Vector Double -> Double
bhattacharyyaCoefficient :: Vector Double -> Vector Double -> Double
bhattacharyyaCoefficient Vector Double
xs Vector Double
ys = (Double -> Double -> Double) -> Double -> Vector Double -> Double
forall a b. (a -> b -> a) -> a -> Vector b -> a
V.foldl' Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 ((Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c. (a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (\Double
x Double
y -> Double -> Double
forall a. Floating a => a -> a
sqrt (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y)) Vector Double
xs Vector Double
ys)

-- | The Hellinger distance measures the similarity of two probability distributions.
--
-- I’ve added it here because I’m hoping to use it to implement a distance metric
-- between two pictures at some point.
--
-- \[
--   H(p,q)={\sqrt {1-BC(p,q)}} \\
--   BC(p,q)=\sum _{{x\in X}}{\sqrt {p(x)q(x)}}
-- \]
--
-- See https://en.wikipedia.org/wiki/Hellinger_distance
hellingerDistance :: Vector Double -> Vector Double -> Double
hellingerDistance :: Vector Double -> Vector Double -> Double
hellingerDistance Vector Double
xs Vector Double
ys = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Vector Double -> Vector Double -> Double
bhattacharyyaCoefficient Vector Double
xs Vector Double
ys)