{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
--------------------------------------------------------------------------------
-- |
-- Module      :  HGeometry.Vector
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- D-dimensional Vectors
--
--------------------------------------------------------------------------------
module HGeometry.Vector
  ( Vector(..)
  , HasComponents(..)
  , module  HGeometry.Vector.Class
  , isScalarMultipleOf
  , scalarMultiple
  , sameDirection
  -- * Special vector operations
  , cross
  ) where

import           GHC.TypeNats(KnownNat, natVal)
import           Control.Lens
import           Data.Proxy
import           Data.Semigroup
import qualified Data.Vector.Generic as GV
import           Data.Vector.Generic.Mutable (MVector(basicInitialize))
import qualified Data.Vector.Generic.Mutable as GMV
import qualified Data.Vector.Unboxed as UV
import qualified Data.Vector.Unboxed.Mutable as UMV
import           HGeometry.Vector.Class
import           HGeometry.Vector.Type

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

-- | Given two colinar vectors, u and v, test if they point in the same direction, i.e.
-- iff scalarMultiple' u v == Just lambda, with lambda > 0
--
-- pre: u and v are colinear, u and v are non-zero
sameDirection     :: ( Additive_ vector d r
                     , Num r, Eq r
                     ) => vector -> vector -> Bool
sameDirection :: forall vector (d :: Nat) r.
(Additive_ vector d r, Num r, Eq r) =>
vector -> vector -> Bool
sameDirection vector
u vector
v = All -> Bool
getAll (All -> Bool) -> All -> Bool
forall a b. (a -> b) -> a -> b
$ (r -> r -> All) -> vector -> vector -> All
forall vector (d :: Nat) r m.
(Additive_ vector d r, Semigroup m) =>
(r -> r -> m) -> vector -> vector -> m
foldMapZip (\r
ux r
vx -> Bool -> All
All (Bool -> All) -> Bool -> All
forall a b. (a -> b) -> a -> b
$ r -> r
forall a. Num a => a -> a
signum r
ux r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r -> r
forall a. Num a => a -> a
signum r
vx) vector
u vector
v
{-# INLINE sameDirection #-}

-- | 'isScalarmultipleof u v' test if v is a scalar multiple of u.
--
-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 10
-- True
-- >>> Vector3 1 1 2 `isScalarMultipleOf` Vector3 10 10 20
-- True
-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 1
-- False
-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 (-1) (-1)
-- True
-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 11.1 11.1
-- True
-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 11.1 11.2
-- False
-- >>> Vector2 2 1 `isScalarMultipleOf` Vector2 11.1 11.2
-- False
-- >>> Vector2 2 1 `isScalarMultipleOf` Vector2 4 2
-- True
-- >>> Vector2 2 1 `isScalarMultipleOf` Vector2 4 0
-- False
-- >>> Vector3 2 1 0 `isScalarMultipleOf` Vector3 4 0 5
-- False
-- >>> Vector3 0 0 0 `isScalarMultipleOf` Vector3 4 0 5
-- True
isScalarMultipleOf       :: ( Eq r, Num r
                            , Metric_ vector d r
                            )
                         => vector -> vector  -> Bool
vector
u isScalarMultipleOf :: forall r vector (d :: Nat).
(Eq r, Num r, Metric_ vector d r) =>
vector -> vector -> Bool
`isScalarMultipleOf` vector
v = let d :: r
d = vector
u vector -> vector -> r
forall vector (d :: Nat) r.
(Metric_ vector d r, Num r) =>
vector -> vector -> r
`dot` vector
v
                               num :: r
num = vector -> r
forall vector (d :: Nat) r.
(Metric_ vector d r, Num r) =>
vector -> r
quadrance vector
u r -> r -> r
forall a. Num a => a -> a -> a
* vector -> r
forall vector (d :: Nat) r.
(Metric_ vector d r, Num r) =>
vector -> r
quadrance vector
v
                           in r
num r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0 Bool -> Bool -> Bool
|| r
num r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
dr -> r -> r
forall a. Num a => a -> a -> a
*r
d
-- u `isScalarMultipleOf` v = isJust $ scalarMultiple u v
{-# INLINE isScalarMultipleOf #-}

-- | scalarMultiple u v computes the scalar labmda s.t. v = lambda * u (if it exists)
scalarMultiple     :: (Eq r, Fractional r, Additive_ vector d r)
                   => vector -> vector -> Maybe r
scalarMultiple :: forall r vector (d :: Nat).
(Eq r, Fractional r, Additive_ vector d r) =>
vector -> vector -> Maybe r
scalarMultiple vector
u vector
v
      | vector -> Bool
allZero vector
u Bool -> Bool -> Bool
|| vector -> Bool
allZero vector
v = r -> Maybe r
forall a. a -> Maybe a
Just r
0
      | Bool
otherwise              = vector -> vector -> Maybe r
forall r vector (d :: Nat).
(Eq r, Fractional r, Additive_ vector d r) =>
vector -> vector -> Maybe r
scalarMultiple' vector
u vector
v
  where
    -- allZero :: (Eq r, Num r, Vector_ vector d r) => vector -> Bool
    allZero :: vector -> Bool
allZero = Getting All vector (IxValue vector)
-> (IxValue vector -> Bool) -> vector -> Bool
forall s a. Getting All s a -> (a -> Bool) -> s -> Bool
allOf Getting All vector (IxValue vector)
forall vector vector'.
HasComponents vector vector' =>
IndexedTraversal1
  Int vector vector' (IxValue vector) (IxValue vector')
IndexedTraversal1
  Int vector vector (IxValue vector) (IxValue vector)
components (IxValue vector -> IxValue vector -> Bool
forall a. Eq a => a -> a -> Bool
== IxValue vector
0)
{-# INLINE scalarMultiple #-}

-- | Helper type for implementing scalarMultiple
data ScalarMultiple r = No | Maybe | Yes r deriving (ScalarMultiple r -> ScalarMultiple r -> Bool
(ScalarMultiple r -> ScalarMultiple r -> Bool)
-> (ScalarMultiple r -> ScalarMultiple r -> Bool)
-> Eq (ScalarMultiple r)
forall r. Eq r => ScalarMultiple r -> ScalarMultiple r -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall r. Eq r => ScalarMultiple r -> ScalarMultiple r -> Bool
== :: ScalarMultiple r -> ScalarMultiple r -> Bool
$c/= :: forall r. Eq r => ScalarMultiple r -> ScalarMultiple r -> Bool
/= :: ScalarMultiple r -> ScalarMultiple r -> Bool
Eq,Int -> ScalarMultiple r -> ShowS
[ScalarMultiple r] -> ShowS
ScalarMultiple r -> String
(Int -> ScalarMultiple r -> ShowS)
-> (ScalarMultiple r -> String)
-> ([ScalarMultiple r] -> ShowS)
-> Show (ScalarMultiple r)
forall r. Show r => Int -> ScalarMultiple r -> ShowS
forall r. Show r => [ScalarMultiple r] -> ShowS
forall r. Show r => ScalarMultiple r -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall r. Show r => Int -> ScalarMultiple r -> ShowS
showsPrec :: Int -> ScalarMultiple r -> ShowS
$cshow :: forall r. Show r => ScalarMultiple r -> String
show :: ScalarMultiple r -> String
$cshowList :: forall r. Show r => [ScalarMultiple r] -> ShowS
showList :: [ScalarMultiple r] -> ShowS
Show)

instance Eq r => Semigroup (ScalarMultiple r) where
  ScalarMultiple r
No      <> :: ScalarMultiple r -> ScalarMultiple r -> ScalarMultiple r
<> ScalarMultiple r
_       = ScalarMultiple r
forall r. ScalarMultiple r
No
  ScalarMultiple r
_       <> ScalarMultiple r
No      = ScalarMultiple r
forall r. ScalarMultiple r
No
  ScalarMultiple r
Maybe   <> ScalarMultiple r
x       = ScalarMultiple r
x
  ScalarMultiple r
x       <> ScalarMultiple r
Maybe   = ScalarMultiple r
x
  (Yes r
x) <> (Yes r
y)
     | r
x r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
y               = r -> ScalarMultiple r
forall r. r -> ScalarMultiple r
Yes r
x
     | Bool
otherwise            = ScalarMultiple r
forall r. ScalarMultiple r
No


instance Eq r => Monoid (ScalarMultiple r) where
  mempty :: ScalarMultiple r
mempty = ScalarMultiple r
forall r. ScalarMultiple r
Maybe
  mappend :: ScalarMultiple r -> ScalarMultiple r -> ScalarMultiple r
mappend = ScalarMultiple r -> ScalarMultiple r -> ScalarMultiple r
forall a. Semigroup a => a -> a -> a
(<>)

-- | Actual implementation of scalarMultiple
scalarMultiple'      :: (Eq r, Fractional r, Additive_ vector d r)
                     => vector -> vector -> Maybe r
scalarMultiple' :: forall r vector (d :: Nat).
(Eq r, Fractional r, Additive_ vector d r) =>
vector -> vector -> Maybe r
scalarMultiple' vector
u vector
v = ScalarMultiple r -> Maybe r
forall {a}. ScalarMultiple a -> Maybe a
g (ScalarMultiple r -> Maybe r) -> ScalarMultiple r -> Maybe r
forall a b. (a -> b) -> a -> b
$ (r -> r -> ScalarMultiple r)
-> vector -> vector -> ScalarMultiple r
forall vector (d :: Nat) r m.
(Additive_ vector d r, Semigroup m) =>
(r -> r -> m) -> vector -> vector -> m
foldMapZip r -> r -> ScalarMultiple r
forall {r}. (Eq r, Fractional r) => r -> r -> ScalarMultiple r
f vector
u vector
v
  where
    f :: r -> r -> ScalarMultiple r
f r
0  r
0  = ScalarMultiple r
forall r. ScalarMultiple r
Maybe -- we don't know lambda yet, but it may still be a scalar mult.
    f r
_  r
0  = ScalarMultiple r
forall r. ScalarMultiple r
No      -- Not a scalar multiple
    f r
ui r
vi = r -> ScalarMultiple r
forall r. r -> ScalarMultiple r
Yes (r -> ScalarMultiple r) -> r -> ScalarMultiple r
forall a b. (a -> b) -> a -> b
$ r
ui r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
vi -- can still be a scalar multiple

    g :: ScalarMultiple a -> Maybe a
g ScalarMultiple a
No      = Maybe a
forall a. Maybe a
Nothing
    g ScalarMultiple a
Maybe   = String -> Maybe a
forall a. HasCallStack => String -> a
error String
"scalarMultiple': found a Maybe, which means the vectors either have length zero, or one of them is all Zero!"
    g (Yes a
x) = a -> Maybe a
forall a. a -> Maybe a
Just a
x


--------------------------------------------------------------------------------
-- * unboxed vectors

-- | elements of the vector are stored consecutively
newtype instance UMV.MVector s (Vector d r) = MV_VectorD (UMV.MVector s r)
newtype instance UV.Vector     (Vector d r) = V_VectorD  (UV.Vector     r)

-- | get the dimension as an Int
natVal' :: forall d. KnownNat d => Int
natVal' :: forall (d :: Nat). KnownNat d => Int
natVal' = Nat -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Nat -> Int) -> Nat -> Int
forall a b. (a -> b) -> a -> b
$ Proxy d -> Nat
forall (n :: Nat) (proxy :: Nat -> *). KnownNat n => proxy n -> Nat
natVal (forall (t :: Nat). Proxy t
forall {k} (t :: k). Proxy t
Proxy @d)

instance ( GMV.MVector UMV.MVector r
         , Vector_ (Vector d r) d r
         ) => GMV.MVector UMV.MVector (Vector d r) where

  basicLength :: forall s. MVector s (Vector d r) -> Int
basicLength (MV_VectorD MVector s r
v) = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                               in MVector s r -> Int
forall s. MVector s r -> Int
forall (v :: * -> * -> *) a s. MVector v a => v s a -> Int
GMV.basicLength MVector s r
v Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
d
  {-# INLINE basicLength #-}
  basicUnsafeSlice :: forall s.
Int -> Int -> MVector s (Vector d r) -> MVector s (Vector d r)
basicUnsafeSlice Int
s Int
l (MV_VectorD MVector s r
v) = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                                        in MVector s r -> MVector s (Vector d r)
forall s (d :: Nat) r. MVector s r -> MVector s (Vector d r)
MV_VectorD (MVector s r -> MVector s (Vector d r))
-> MVector s r -> MVector s (Vector d r)
forall a b. (a -> b) -> a -> b
$ Int -> Int -> MVector s r -> MVector s r
forall s. Int -> Int -> MVector s r -> MVector s r
forall (v :: * -> * -> *) a s.
MVector v a =>
Int -> Int -> v s a -> v s a
GMV.basicUnsafeSlice (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
s) (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
l) MVector s r
v
  {-# INLINE basicUnsafeSlice #-}
  basicOverlaps :: forall s. MVector s (Vector d r) -> MVector s (Vector d r) -> Bool
basicOverlaps  (MV_VectorD MVector s r
v) (MV_VectorD MVector s r
v') = MVector s r -> MVector s r -> Bool
forall s. MVector s r -> MVector s r -> Bool
forall (v :: * -> * -> *) a s.
MVector v a =>
v s a -> v s a -> Bool
GMV.basicOverlaps MVector s r
v MVector s r
v'
  {-# INLINE basicOverlaps #-}
  basicUnsafeNew :: forall s. Int -> ST s (MVector s (Vector d r))
basicUnsafeNew Int
n = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                     in MVector s r -> MVector s (Vector d r)
forall s (d :: Nat) r. MVector s r -> MVector s (Vector d r)
MV_VectorD (MVector s r -> MVector s (Vector d r))
-> ST s (MVector s r) -> ST s (MVector s (Vector d r))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> ST s (MVector s r)
forall s. Int -> ST s (MVector s r)
forall (v :: * -> * -> *) a s. MVector v a => Int -> ST s (v s a)
GMV.basicUnsafeNew (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
  {-# INLINE basicUnsafeNew #-}
  basicInitialize :: forall s. MVector s (Vector d r) -> ST s ()
basicInitialize (MV_VectorD MVector s r
v) = MVector s r -> ST s ()
forall s. MVector s r -> ST s ()
forall (v :: * -> * -> *) a s. MVector v a => v s a -> ST s ()
GMV.basicInitialize MVector s r
v
  {-# INLINE basicInitialize#-}
  basicUnsafeRead :: forall s. MVector s (Vector d r) -> Int -> ST s (Vector d r)
basicUnsafeRead (MV_VectorD MVector s r
v) Int
i = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                                     in (Int -> ST s r) -> ST s (Vector d r)
forall vector (d :: Nat) r (f :: * -> *).
(Vector_ vector d r, Applicative f) =>
(Int -> f r) -> f vector
forall (f :: * -> *).
Applicative f =>
(Int -> f r) -> f (Vector d r)
generateA ((Int -> ST s r) -> ST s (Vector d r))
-> (Int -> ST s r) -> ST s (Vector d r)
forall a b. (a -> b) -> a -> b
$ \Int
j -> MVector s r -> Int -> ST s r
forall s. MVector s r -> Int -> ST s r
forall (v :: * -> * -> *) a s.
MVector v a =>
v s a -> Int -> ST s a
GMV.basicUnsafeRead MVector s r
v (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j)
  {-# INLINE basicUnsafeRead #-}
  basicUnsafeWrite :: forall s. MVector s (Vector d r) -> Int -> Vector d r -> ST s ()
basicUnsafeWrite (MV_VectorD MVector s r
v) Int
i Vector d r
w = IndexedGetting Int (Sequenced () (ST s)) (Vector d r) r
-> (Int -> r -> ST s ()) -> Vector d r -> ST s ()
forall (m :: * -> *) i r s a.
Monad m =>
IndexedGetting i (Sequenced r m) s a
-> (i -> a -> m r) -> s -> m ()
imapMOf_ IndexedGetting Int (Sequenced () (ST s)) (Vector d r) r
Indexed
  Int
  (IxValue (Vector d r))
  (Const (Sequenced () (ST s)) (IxValue (Vector d r)))
-> Vector d r -> Const (Sequenced () (ST s)) (Vector d r)
forall vector vector'.
HasComponents vector vector' =>
IndexedTraversal1
  Int vector vector' (IxValue vector) (IxValue vector')
IndexedTraversal1
  Int
  (Vector d r)
  (Vector d r)
  (IxValue (Vector d r))
  (IxValue (Vector d r))
components Int -> r -> ST s ()
f Vector d r
w
    where
      f :: Int -> r -> ST s ()
f Int
j r
x = MVector s r -> Int -> r -> ST s ()
forall s. MVector s r -> Int -> r -> ST s ()
forall (v :: * -> * -> *) a s.
MVector v a =>
v s a -> Int -> a -> ST s ()
GMV.basicUnsafeWrite MVector s r
v (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j) r
x
      d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
  {-# INLINE basicUnsafeWrite #-}


instance ( GV.Vector UV.Vector r
         , Vector_ (Vector d r) d r
         ) => GV.Vector UV.Vector (Vector d r) where

  basicUnsafeFreeze :: forall s.
Mutable Vector s (Vector d r) -> ST s (Vector (Vector d r))
basicUnsafeFreeze (MV_VectorD MVector s r
mv) = Vector r -> Vector (Vector d r)
forall (d :: Nat) r. Vector r -> Vector (Vector d r)
V_VectorD (Vector r -> Vector (Vector d r))
-> ST s (Vector r) -> ST s (Vector (Vector d r))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Mutable Vector s r -> ST s (Vector r)
forall s. Mutable Vector s r -> ST s (Vector r)
forall (v :: * -> *) a s. Vector v a => Mutable v s a -> ST s (v a)
GV.basicUnsafeFreeze MVector s r
Mutable Vector s r
mv
  {-# INLINE basicUnsafeFreeze #-}
  basicUnsafeThaw :: forall s.
Vector (Vector d r) -> ST s (Mutable Vector s (Vector d r))
basicUnsafeThaw (V_VectorD Vector r
v) = MVector s r -> MVector s (Vector d r)
forall s (d :: Nat) r. MVector s r -> MVector s (Vector d r)
MV_VectorD (MVector s r -> MVector s (Vector d r))
-> ST s (MVector s r) -> ST s (MVector s (Vector d r))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Vector r -> ST s (Mutable Vector s r)
forall s. Vector r -> ST s (Mutable Vector s r)
forall (v :: * -> *) a s. Vector v a => v a -> ST s (Mutable v s a)
GV.basicUnsafeThaw Vector r
v
  {-# INLINE basicUnsafeThaw #-}
  basicLength :: Vector (Vector d r) -> Int
basicLength (V_VectorD Vector r
v) = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                              in Vector r -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
GV.basicLength Vector r
v Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
d
  {-# INLINE basicLength #-}
  basicUnsafeSlice :: Int -> Int -> Vector (Vector d r) -> Vector (Vector d r)
basicUnsafeSlice Int
s Int
l (V_VectorD Vector r
v) = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                                       in Vector r -> Vector (Vector d r)
forall (d :: Nat) r. Vector r -> Vector (Vector d r)
V_VectorD (Vector r -> Vector (Vector d r))
-> Vector r -> Vector (Vector d r)
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Vector r -> Vector r
forall (v :: * -> *) a. Vector v a => Int -> Int -> v a -> v a
GV.basicUnsafeSlice (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
s) (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
l) Vector r
v
  {-# INLINE basicUnsafeSlice #-}
  basicUnsafeIndexM :: Vector (Vector d r) -> Int -> Box (Vector d r)
basicUnsafeIndexM (V_VectorD Vector r
v) Int
i = let d :: Int
d = forall (d :: Nat). KnownNat d => Int
natVal' @d
                                      in (Int -> Box r) -> Box (Vector d r)
forall vector (d :: Nat) r (f :: * -> *).
(Vector_ vector d r, Applicative f) =>
(Int -> f r) -> f vector
forall (f :: * -> *).
Applicative f =>
(Int -> f r) -> f (Vector d r)
generateA ((Int -> Box r) -> Box (Vector d r))
-> (Int -> Box r) -> Box (Vector d r)
forall a b. (a -> b) -> a -> b
$ \Int
j -> Vector r -> Int -> Box r
forall (v :: * -> *) a. Vector v a => v a -> Int -> Box a
GV.basicUnsafeIndexM Vector r
v (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j)
  {-# INLINE basicUnsafeIndexM #-}

instance ( UV.Unbox r, Vector_ (Vector d r) d r
         ) => UV.Unbox (Vector d r) where


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