Module.hs

{-# LANGUAGE AllowAmbiguousTypes  #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE TemplateHaskell      #-}
{-# LANGUAGE UndecidableInstances #-}

{-# OPTIONS_GHC -Wno-orphans #-}

module Math.Module
  ( Module(..), lerp
  , Inner(..), Cross(..)
  , Interpolatable
  , norm, squaredNorm, normalise
  , quadrance, distance
  , proj, projC, closestPointOnSegment
  , strictlyParallel, convexCombination

  , ViaModule(..)
  )
  where

-- base
import Control.Monad
  ( guard )
import Data.Coerce
  ( coerce )
import Data.Kind
  ( Type, Constraint )
import Data.Monoid
  ( Ap(..), Sum(..) )

-- acts
import Data.Act
  ( Act
    ( (•) )
  , Torsor
    ( (-->) )
  )

-- groups
import Data.Group
  ( Group(..) )

-- brush-strokes
import Math.Epsilon
  ( epsilon )
import Math.Linear
import Math.Module.Internal
  ( Module(..), Inner(..) )
import Math.Ring
  ( Ring )
import qualified Math.Ring as Ring

--------------------------------------------------------------------------------

instance ( Applicative f, Module r m ) => Module r ( Ap f m ) where
  origin = pure origin
  (^+^)  = liftA2 (^+^)
  (^-^)  = liftA2 (^-^)
  (*^) r = fmap ( r *^ )

lerp :: forall v r p. ( Module r v, Torsor v p ) => r -> p -> p -> p
lerp t p0 p1 = ( t *^ ( p0 --> p1 :: v ) ) • p0

class Module r m => Cross r m where
  (×) :: m -> m -> r

-- | Norm of a vector, computed using the inner product.
norm :: forall m r. ( Floating r, Inner r m ) => m -> r
norm = sqrt . squaredNorm

-- | Squared norm of a vector, computed using the inner product.
squaredNorm :: forall m r. Inner r m => m -> r
squaredNorm v = v ^.^ v

normalise :: ( Floating r, Inner r m ) => m -> m
normalise v = recip ( norm v ) *^ v

-- | Squared distance between two points.
quadrance :: forall v r p. ( Inner r v, Torsor v p ) => p -> p -> r
quadrance p1 p2 = squaredNorm ( p1 --> p2 :: v )

-- | Distance between two points.
distance :: forall v r p. ( Floating r, Inner r v, Torsor v p ) => p -> p -> r
distance p1 p2 = sqrt ( quadrance @v p1 p2 )

-- | Projects the first argument onto the second.
proj :: forall m r. ( Inner r m, Fractional r ) => m -> m -> m
proj x y = projC x y *^ y

-- | Projection constant: how far along the projection of the first vector lands along the second vector.
projC :: forall m r. ( Inner r m, Fractional r ) => m -> m -> r
projC x y = x ^.^ y / squaredNorm y

closestPointOnSegment
  :: forall v r p
  . ( Inner r v, Torsor v p, Fractional r, Ord r )
  => p -> Segment p -> ( r, p )
closestPointOnSegment c ( Segment p0 p1 )
  | t <= 0
  = ( 0, p0 )
  | t >= 1
  = ( 1, p1 )
  | otherwise
  = ( t, ( t *^ v01 ) • p0 )
  where
    v01 :: v
    v01 = p0 --> p1
    t :: r
    t = projC ( p0 --> c ) v01

--------------------------------------------------------------------------------

-- | A convenient constraint synonym for types that support interpolation.
type Interpolatable :: Type -> Type -> Constraint
class    ( Torsor ( T u ) u, Module r ( T u ) ) => Interpolatable r u
instance ( Torsor ( T u ) u, Module r ( T u ) ) => Interpolatable r u

--------------------------------------------------------------------------------

instance Ring a => Inner a ( Sum a ) where
  Sum a ^.^ Sum b = a Ring.* b

instance Inner Double ( T ( ℝ 2 ) ) where
  V2 x1 y1 ^.^ V2 x2 y2 = x1 Ring.* x2 + y1 Ring.* y2

instance Cross Double ( T ( ℝ 2 ) ) where
  V2 x1 y1 × V2 x2 y2 = x1 Ring.* y2 Ring.- x2 Ring.* y1

-- | Compute whether two vectors point in the same direction,
-- that is, whether each vector is a (strictly) positive multiple of the other.
--
-- Returns @False@ if either of the vectors is zero (or very close to zero).
strictlyParallel :: T ( ℝ 2 ) -> T ( ℝ 2 ) -> Bool
strictlyParallel u v
  =  abs ( u  ×  v ) < tol -- vectors are collinear
  &&       u ^.^ v   > tol -- vectors point in the same direction (parallel and not anti-parallel)
  where
    tol = norm u * norm v * epsilon

-- | Finds whether the query vector @ u @ is a convex combination of the two provided vectors @ v0 @, @ v1 @.
--
-- If so, returns @ t @ in @ [ 0, 1 ] @ such that @ ( 1 - t ) v0 + t v1 @ is a positive multiple of @ u @.
convexCombination
  :: T ( ℝ 2 ) -- ^ first vector
  -> T ( ℝ 2 ) -- ^ second vector
  -> T ( ℝ 2 ) -- ^ query vector
  -> Maybe Double
convexCombination v0 v1 u
  | abs c10 < epsilon
  = if strictlyParallel u v0
    then Just 0
    else if strictlyParallel u v1
      then Just 1
      else Nothing
  | otherwise
  = do
      let
        t :: Double
        t = c0 / c10
      guard ( t > -epsilon && t < 1 + epsilon )
      guard ( epsilon < u ^.^ ( lerp @( T ( ℝ 2 ) ) t v0 v1 ) )
      Just $ min 1 ( max 0 t )

  where
    c0, c10 :: Double
    c0  = v0 × u
    c10 = ( v0 ^-^ v1 ) × u

--------------------------------------------------------------------------------
-- Not sure how to set things up to automate the following...

instance Module Double ( T ( ℝ 0 ) ) where
  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )

instance Module Double ( T ( ℝ 1 ) ) where
  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )

instance Module Double ( T ( ℝ 2 ) ) where
  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )

instance Module Double ( T ( ℝ 3 ) ) where
  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )

instance Module Double ( T ( ℝ 4 ) ) where
  origin      = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
  T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
  T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
  k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )

deriving via ViaModule Double ( T ( ℝ n ) )
  instance Module Double ( T ( ℝ n ) ) => Semigroup ( T ( ℝ n ) )
deriving via ViaModule Double ( T ( ℝ n ) )
  instance Module Double ( T ( ℝ n ) ) => Monoid ( T ( ℝ n ) )
deriving via ViaModule Double ( T ( ℝ n ) )
  instance Module Double ( T ( ℝ n ) ) => Group ( T ( ℝ n ) )

deriving via ViaModule Double ( ℝ n )
  instance Module Double ( T ( ℝ n ) ) => Act ( T ( ℝ n ) ) ( ℝ n )
deriving via ( ViaModule Double ( ℝ n ) )
  instance Module Double ( T ( ℝ n ) ) => Torsor ( T ( ℝ n ) ) ( ℝ n )

--------------------------------------------------------------------------------

newtype ViaModule r m = ViaModule { unViaModule :: m }

instance Module r m => Semigroup ( ViaModule r m ) where
  (<>) = coerce $ (^+^) @r @m
instance Module r m => Monoid ( ViaModule r m ) where
  mempty = coerce $ origin @r @m
instance Module r m => Group ( ViaModule r m ) where
  invert = coerce $ ( \ x -> (^-^) @r @m ( origin @r @m ) x )

instance ( Semigroup ( T m ), Module r ( T m ) ) => Act    ( T m ) ( ViaModule r m ) where
  g • ViaModule m = ViaModule ( unT $ g ^+^ T m )
instance ( Group     ( T m ), Module r ( T m ) ) => Torsor ( T m ) ( ViaModule r m ) where
  ViaModule a --> ViaModule b = T ( unT $ T b ^-^ T a )