{-# LANGUAGE TemplateHaskell       #-}

module Bio.Utils.Geometry
    ( R
    , V3R
    , Ray (..)
    , AffineTransformable(..)
    , Epsilon (..)
    , zoRay
    , cross, dot
    , norm , normalize
    , distance, angle, dihedral
    , svd3
    ) where

import           Control.Lens
import           Linear.V3                      ( V3
                                                , cross
                                                )
import           Linear.Vector                  ( zero )
import           Linear.Epsilon                 ( Epsilon (..) )
import           Linear.Matrix                  ( M33 )
import           Linear.Metric                  ( dot
                                                , norm
                                                , normalize
                                                , distance
                                                )
import qualified Linear.Quaternion             as Q
                                                ( rotate
                                                , axisAngle
                                                )

-- | Default floating point type, switch here to move to Doubles
--
type R = Float

-- | Defalut type of 3D vectors
--
type V3R = V3 R

-- | Ray has an origin and a direction
--
data Ray a = Ray { forall a. Ray a -> a
_origin    :: a
                 , forall a. Ray a -> a
_direction :: a
                 }

makeLenses ''Ray

-- | Zero-origin ray
zoRay :: V3R -> Ray V3R
zoRay :: V3R -> Ray V3R
zoRay = forall a. a -> a -> Ray a
Ray forall (f :: * -> *) a. (Additive f, Num a) => f a
zero forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a (f :: * -> *).
(Floating a, Metric f, Epsilon a) =>
f a -> f a
normalize

-- | Affine transformations for vectors and sets of vectors
--
class AffineTransformable a where
    -- | Rotate an object around the vector by some angle
    --
    rotate    :: V3R -> R -> a -> a

    -- | Rotate an object around the ray by some angle
    --
    rotateR   :: Ray V3R -> R -> a -> a

    -- | Translocate an object by some vectors
    --
    translate :: V3R -> a -> a

-- | We can apply affine transformations to vectors
--
instance AffineTransformable V3R where
    rotate :: V3R -> R -> V3R -> V3R
rotate V3R
v R
a = forall a.
(Conjugate a, RealFloat a) =>
Quaternion a -> V3 a -> V3 a
Q.rotate (forall a. (Epsilon a, Floating a) => V3 a -> a -> Quaternion a
Q.axisAngle V3R
v R
a)
    rotateR :: Ray V3R -> R -> V3R -> V3R
rotateR Ray V3R
r R
a V3R
x = forall a. AffineTransformable a => V3R -> R -> a -> a
rotate (Ray V3R
r forall s a. s -> Getting a s a -> a
^. forall a. Lens' (Ray a) a
direction) R
a (V3R
x forall a. Num a => a -> a -> a
- Ray V3R
r forall s a. s -> Getting a s a -> a
^. forall a. Lens' (Ray a) a
origin) forall a. Num a => a -> a -> a
+ Ray V3R
r forall s a. s -> Getting a s a -> a
^. forall a. Lens' (Ray a) a
origin
    translate :: V3R -> V3R -> V3R
translate V3R
v = (V3R
v forall a. Num a => a -> a -> a
+)

-- | If we have any collection of vectors, than we can transform it too
--
instance Functor f => AffineTransformable (f V3R) where
    rotate :: V3R -> R -> f V3R -> f V3R
rotate V3R
v R
a = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. AffineTransformable a => V3R -> R -> a -> a
rotate V3R
v R
a)
    rotateR :: Ray V3R -> R -> f V3R -> f V3R
rotateR Ray V3R
r R
a = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. AffineTransformable a => Ray V3R -> R -> a -> a
rotateR Ray V3R
r R
a)
    translate :: V3R -> f V3R -> f V3R
translate V3R
v = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. AffineTransformable a => V3R -> a -> a
translate V3R
v)

-- | Measure angle between vectors
--
angle :: V3R -> V3R -> R
angle :: V3R -> V3R -> R
angle V3R
a V3R
b = forall a. RealFloat a => a -> a -> a
atan2 (forall (f :: * -> *) a. (Metric f, Floating a) => f a -> a
norm (V3R
a forall a. Num a => V3 a -> V3 a -> V3 a
`cross` V3R
b)) (V3R
a forall (f :: * -> *) a. (Metric f, Num a) => f a -> f a -> a
`dot` V3R
b)

-- | Measure dihedral between four points
-- by https://math.stackexchange.com/a/47084
--
dihedral :: V3R -> V3R -> V3R -> V3R -> R
dihedral :: V3R -> V3R -> V3R -> V3R -> R
dihedral V3R
x V3R
y V3R
z V3R
w = let b1 :: V3R
b1 = V3R
y forall a. Num a => a -> a -> a
- V3R
x
                       b2 :: V3R
b2 = V3R
z forall a. Num a => a -> a -> a
- V3R
y
                       b3 :: V3R
b3 = V3R
w forall a. Num a => a -> a -> a
- V3R
z
                       n1 :: V3R
n1 = forall a (f :: * -> *).
(Floating a, Metric f, Epsilon a) =>
f a -> f a
normalize forall a b. (a -> b) -> a -> b
$ V3R
b1 forall a. Num a => V3 a -> V3 a -> V3 a
`cross` V3R
b2
                       n2 :: V3R
n2 = forall a (f :: * -> *).
(Floating a, Metric f, Epsilon a) =>
f a -> f a
normalize forall a b. (a -> b) -> a -> b
$ V3R
b2 forall a. Num a => V3 a -> V3 a -> V3 a
`cross` V3R
b3
                       m1 :: V3R
m1 = V3R
n1 forall a. Num a => V3 a -> V3 a -> V3 a
`cross` forall a (f :: * -> *).
(Floating a, Metric f, Epsilon a) =>
f a -> f a
normalize V3R
b2
                   in  forall a. RealFloat a => a -> a -> a
atan2 (V3R
m1 forall (f :: * -> *) a. (Metric f, Num a) => f a -> f a -> a
`dot` V3R
n2) (V3R
n1 forall (f :: * -> *) a. (Metric f, Num a) => f a -> f a -> a
`dot` V3R
n2)

data SVD a = SVD { forall a. SVD a -> a
svdU :: a
                 , forall a. SVD a -> a
svdS :: a
                 , forall a. SVD a -> a
svdV :: a
                 }
  deriving (Int -> SVD a -> ShowS
forall a. Show a => Int -> SVD a -> ShowS
forall a. Show a => [SVD a] -> ShowS
forall a. Show a => SVD a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [SVD a] -> ShowS
$cshowList :: forall a. Show a => [SVD a] -> ShowS
show :: SVD a -> String
$cshow :: forall a. Show a => SVD a -> String
showsPrec :: Int -> SVD a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> SVD a -> ShowS
Show, SVD a -> SVD a -> Bool
forall a. Eq a => SVD a -> SVD a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: SVD a -> SVD a -> Bool
$c/= :: forall a. Eq a => SVD a -> SVD a -> Bool
== :: SVD a -> SVD a -> Bool
$c== :: forall a. Eq a => SVD a -> SVD a -> Bool
Eq)

-- | Singular value decomposition
-- for 3x3 matricies
svd3 :: M33 R -> SVD (M33 R)
svd3 :: M33 R -> SVD (M33 R)
svd3 = forall a. HasCallStack => a
undefined