{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}

-- SPDX-License-Identifier: MPL-2.0

{- | Immutable dynamically sized and type-indexed vectors. All indices are
zero based. Checked operations reject mismatched dimensions and zero norms.
-}
module HMath.Vector (
    Vector,
    vector,
    vectorValues,
    vectorDimension,
    FixedVector,
    SparseVector,
    fixedVector,
    fixedValues,
    sparseVector,
    sparseVectorValues,
    vectorSlice,
    vectorConcat,
    vectorAdd,
    vectorSubtract,
    vectorScale,
    vectorDot,
    cross,
    outer,
    triple,
    norm0,
    norm1,
    norm2,
    normP,
    normInfinity,
    normalizeVector,
    project,
    reflect,
    refract,
    distance,
    angle,
    cosineSimilarity,
    gramSchmidt,
) where

import qualified Data.Foldable as Foldable (foldl')
import qualified Data.IntMap.Strict as IM
import Data.Proxy (Proxy (..))
import GHC.TypeLits (KnownNat, Nat, natVal)
import HMath.Error
import HMath.Matrix (Matrix, fromListChecked)
import HMath.Matrix.Internal (finiteValues, stableNorm)

{- | Immutable dynamically sized coordinate list. Construct with vector; checked binary
operations reject dimension mismatch.
-}
newtype Vector a = Vector [a] deriving (Vector a -> Vector a -> Bool
(Vector a -> Vector a -> Bool)
-> (Vector a -> Vector a -> Bool) -> Eq (Vector a)
forall a. Eq a => Vector a -> Vector a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Vector a -> Vector a -> Bool
== :: Vector a -> Vector a -> Bool
$c/= :: forall a. Eq a => Vector a -> Vector a -> Bool
/= :: Vector a -> Vector a -> Bool
Eq, Int -> Vector a -> ShowS
[Vector a] -> ShowS
Vector a -> String
(Int -> Vector a -> ShowS)
-> (Vector a -> String) -> ([Vector a] -> ShowS) -> Show (Vector a)
forall a. Show a => Int -> Vector a -> ShowS
forall a. Show a => [Vector a] -> ShowS
forall a. Show a => Vector a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Vector a -> ShowS
showsPrec :: Int -> Vector a -> ShowS
$cshow :: forall a. Show a => Vector a -> String
show :: Vector a -> String
$cshowList :: forall a. Show a => [Vector a] -> ShowS
showList :: [Vector a] -> ShowS
Show)

{- | Opaque vector whose length matches its KnownNat type index. Construct with fixedVector
for a checked dynamic/static boundary.
-}
newtype FixedVector (n :: Nat) a = FixedVector [a] deriving (FixedVector n a -> FixedVector n a -> Bool
(FixedVector n a -> FixedVector n a -> Bool)
-> (FixedVector n a -> FixedVector n a -> Bool)
-> Eq (FixedVector n a)
forall (n :: Nat) a.
Eq a =>
FixedVector n a -> FixedVector n a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall (n :: Nat) a.
Eq a =>
FixedVector n a -> FixedVector n a -> Bool
== :: FixedVector n a -> FixedVector n a -> Bool
$c/= :: forall (n :: Nat) a.
Eq a =>
FixedVector n a -> FixedVector n a -> Bool
/= :: FixedVector n a -> FixedVector n a -> Bool
Eq, Int -> FixedVector n a -> ShowS
[FixedVector n a] -> ShowS
FixedVector n a -> String
(Int -> FixedVector n a -> ShowS)
-> (FixedVector n a -> String)
-> ([FixedVector n a] -> ShowS)
-> Show (FixedVector n a)
forall (n :: Nat) a. Show a => Int -> FixedVector n a -> ShowS
forall (n :: Nat) a. Show a => [FixedVector n a] -> ShowS
forall (n :: Nat) a. Show a => FixedVector n a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall (n :: Nat) a. Show a => Int -> FixedVector n a -> ShowS
showsPrec :: Int -> FixedVector n a -> ShowS
$cshow :: forall (n :: Nat) a. Show a => FixedVector n a -> String
show :: FixedVector n a -> String
$cshowList :: forall (n :: Nat) a. Show a => [FixedVector n a] -> ShowS
showList :: [FixedVector n a] -> ShowS
Show)

{- | Opaque dynamic vector with canonical nonzero coordinate entries. Construct with
sparseVector; duplicate indices add.
-}
data SparseVector a = SparseVector Int (IM.IntMap a) deriving (SparseVector a -> SparseVector a -> Bool
(SparseVector a -> SparseVector a -> Bool)
-> (SparseVector a -> SparseVector a -> Bool)
-> Eq (SparseVector a)
forall a. Eq a => SparseVector a -> SparseVector a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => SparseVector a -> SparseVector a -> Bool
== :: SparseVector a -> SparseVector a -> Bool
$c/= :: forall a. Eq a => SparseVector a -> SparseVector a -> Bool
/= :: SparseVector a -> SparseVector a -> Bool
Eq, Int -> SparseVector a -> ShowS
[SparseVector a] -> ShowS
SparseVector a -> String
(Int -> SparseVector a -> ShowS)
-> (SparseVector a -> String)
-> ([SparseVector a] -> ShowS)
-> Show (SparseVector a)
forall a. Show a => Int -> SparseVector a -> ShowS
forall a. Show a => [SparseVector a] -> ShowS
forall a. Show a => SparseVector a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> SparseVector a -> ShowS
showsPrec :: Int -> SparseVector a -> ShowS
$cshow :: forall a. Show a => SparseVector a -> String
show :: SparseVector a -> String
$cshowList :: forall a. Show a => [SparseVector a] -> ShowS
showList :: [SparseVector a] -> ShowS
Show)

{- | Wrap a list as an immutable dynamic vector without numeric validation. The list order
is preserved. Example: vectorValues (vector [1,2]) == [1,2].
-}
vector :: [a] -> Vector a
vector :: forall a. [a] -> Vector a
vector = [a] -> Vector a
forall a. [a] -> Vector a
Vector

{- | Extract the immutable list of vector values in order; no shape conversion or numeric
validation occurs.
-}
vectorValues :: Vector a -> [a]
vectorValues :: forall a. Vector a -> [a]
vectorValues (Vector [a]
xs) = [a]
xs

{- | Return the number of coordinates in a dynamic vector. Empty vectors have dimension
zero.
-}
vectorDimension :: Vector a -> Int
vectorDimension :: forall a. Vector a -> Int
vectorDimension = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> Int) -> (Vector a -> [a]) -> Vector a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
vectorValues

{- | Construct a vector whose list length equals the type-level KnownNat size. A mismatch
returns Left InvalidArgument. For example, fixedVector @3 requires exactly three
values.
-}
fixedVector :: forall n a. (KnownNat n) => [a] -> Either MathError (FixedVector n a)
fixedVector :: forall (n :: Nat) a.
KnownNat n =>
[a] -> Either MathError (FixedVector n a)
fixedVector [a]
xs
    | Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Integer
forall a. Integral a => a -> Integer
toInteger (Proxy n -> Integer
forall (n :: Nat) (proxy :: Nat -> *).
KnownNat n =>
proxy n -> Integer
natVal (Proxy n
forall {k} (t :: k). Proxy t
Proxy :: Proxy n)) = FixedVector n a -> Either MathError (FixedVector n a)
forall a b. b -> Either a b
Right ([a] -> FixedVector n a
forall (n :: Nat) a. [a] -> FixedVector n a
FixedVector [a]
xs)
    | Bool
otherwise = MathError -> Either MathError (FixedVector n a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"fixedVector: dimension differs from type-level size")

{- | Extract the ordered coordinates of a type-indexed vector; its length is the type-level
size.
-}
fixedValues :: FixedVector n a -> [a]
fixedValues :: forall (n :: Nat) a. FixedVector n a -> [a]
fixedValues (FixedVector [a]
xs) = [a]
xs

{- | Given nonnegative dimension and (zero-based index,value) entries, sum duplicate indices
and remove explicit zeros. Reject out-of-range indices or negative dimension.
-}
sparseVector :: (Eq a, Num a) => Int -> [(Int, a)] -> Either MathError (SparseVector a)
sparseVector :: forall a.
(Eq a, Num a) =>
Int -> [(Int, a)] -> Either MathError (SparseVector a)
sparseVector Int
n [(Int, a)]
xs
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| ((Int, a) -> Bool) -> [(Int, a)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\(Int
i, a
_) -> Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n) [(Int, a)]
xs =
        MathError -> Either MathError (SparseVector a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"sparseVector: invalid dimension/index")
    | Bool
otherwise = SparseVector a -> Either MathError (SparseVector a)
forall a b. b -> Either a b
Right (Int -> IntMap a -> SparseVector a
forall a. Int -> IntMap a -> SparseVector a
SparseVector Int
n ((a -> Bool) -> IntMap a -> IntMap a
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0) ((a -> a -> a) -> [(Int, a)] -> IntMap a
forall a. (a -> a -> a) -> [(Int, a)] -> IntMap a
IM.fromListWith a -> a -> a
forall a. Num a => a -> a -> a
(+) [(Int, a)]
xs)))

{- | Expand a sparse vector to a dense list of its declared dimension, filling absent
entries with zero.
-}
sparseVectorValues :: (Num a) => SparseVector a -> [a]
sparseVectorValues :: forall a. Num a => SparseVector a -> [a]
sparseVectorValues (SparseVector Int
n IntMap a
xs) = [a -> Int -> IntMap a -> a
forall a. a -> Int -> IntMap a -> a
IM.findWithDefault a
0 Int
i IntMap a
xs | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]

{- | Given start, count, then vector, return the half-open slice [start,start+count). Reject
negative or out-of-range bounds; an empty end slice is valid.
-}
vectorSlice :: Int -> Int -> Vector a -> Either MathError (Vector a)
vectorSlice :: forall a. Int -> Int -> Vector a -> Either MathError (Vector a)
vectorSlice Int
start Int
count (Vector [a]
xs)
    | Int
start Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
start Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n Bool -> Bool -> Bool
|| Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
start =
        MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"vectorSlice: invalid bounds")
    | Bool
otherwise = Vector a -> Either MathError (Vector a)
forall a b. b -> Either a b
Right ([a] -> Vector a
forall a. [a] -> Vector a
Vector (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
count (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
start [a]
xs)))
  where
    n :: Int
n = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs

{- | Concatenate a list of vectors in order. Return Left InvalidArgument if the combined
dimension would overflow Int; an empty list gives an empty vector.
-}
vectorConcat :: [Vector a] -> Either MathError (Vector a)
vectorConcat :: forall a. [Vector a] -> Either MathError (Vector a)
vectorConcat [Vector a]
vs
    | [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Vector a -> Integer) -> [Vector a] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Integer) -> (Vector a -> Int) -> Vector a -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> Int
forall a. Vector a -> Int
vectorDimension) [Vector a]
vs) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int) =
        MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"vectorConcat: size overflow")
    | Bool
otherwise = Vector a -> Either MathError (Vector a)
forall a b. b -> Either a b
Right ([a] -> Vector a
forall a. [a] -> Vector a
Vector ((Vector a -> [a]) -> [Vector a] -> [a]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Vector a -> [a]
forall a. Vector a -> [a]
vectorValues [Vector a]
vs))

zipVector :: (a -> b -> c) -> Vector a -> Vector b -> Either MathError (Vector c)
zipVector :: forall a b c.
(a -> b -> c)
-> Vector a -> Vector b -> Either MathError (Vector c)
zipVector a -> b -> c
f (Vector [a]
a) (Vector [b]
b)
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [b] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
b = MathError -> Either MathError (Vector c)
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a) ([b] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
b))
    | Bool
otherwise = Vector c -> Either MathError (Vector c)
forall a b. b -> Either a b
Right ([c] -> Vector c
forall a. [a] -> Vector a
Vector ((a -> b -> c) -> [a] -> [b] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> b -> c
f [a]
a [b]
b))

{- | Add corresponding values of equally sized vectors, returning a new vector. Unequal
dimensions return VectorLengthMismatch. Generic Num arithmetic is unchecked.
-}
vectorAdd :: (Num a) => Vector a -> Vector a -> Either MathError (Vector a)
vectorAdd :: forall a.
Num a =>
Vector a -> Vector a -> Either MathError (Vector a)
vectorAdd = (a -> a -> a)
-> Vector a -> Vector a -> Either MathError (Vector a)
forall a b c.
(a -> b -> c)
-> Vector a -> Vector b -> Either MathError (Vector c)
zipVector a -> a -> a
forall a. Num a => a -> a -> a
(+)

{- | Subtract corresponding right-hand coordinates from left-hand coordinates. Unequal
dimensions return VectorLengthMismatch; Num arithmetic is unchecked.
-}
vectorSubtract :: (Num a) => Vector a -> Vector a -> Either MathError (Vector a)
vectorSubtract :: forall a.
Num a =>
Vector a -> Vector a -> Either MathError (Vector a)
vectorSubtract = (a -> a -> a)
-> Vector a -> Vector a -> Either MathError (Vector a)
forall a b c.
(a -> b -> c)
-> Vector a -> Vector b -> Either MathError (Vector c)
zipVector (-)

{- | Given scalar then vector, multiply every coordinate by the scalar. Return a new vector,
retaining the original dimension and Num arithmetic behavior.
-}
vectorScale :: (Num a) => a -> Vector a -> Vector a
vectorScale :: forall a. Num a => a -> Vector a -> Vector a
vectorScale a
s (Vector [a]
xs) = [a] -> Vector a
forall a. [a] -> Vector a
Vector ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
s a -> a -> a
forall a. Num a => a -> a -> a
*) [a]
xs)

{- | Sum coordinatewise products of equally sized vectors. Empty vectors give zero; unequal
dimensions return VectorLengthMismatch. No floating-point finiteness check is added.
-}
vectorDot :: (Num a) => Vector a -> Vector a -> Either MathError a
vectorDot :: forall a. Num a => Vector a -> Vector a -> Either MathError a
vectorDot Vector a
a Vector a
b = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> (Vector a -> [a]) -> Vector a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
vectorValues (Vector a -> a)
-> Either MathError (Vector a) -> Either MathError a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (a -> a -> a)
-> Vector a -> Vector a -> Either MathError (Vector a)
forall a b c.
(a -> b -> c)
-> Vector a -> Vector b -> Either MathError (Vector c)
zipVector a -> a -> a
forall a. Num a => a -> a -> a
(*) Vector a
a Vector a
b

{- | Compute the right-handed cross product of two three-dimensional vectors. Any other
dimension returns Left InvalidArgument.
-}
cross :: (Num a) => Vector a -> Vector a -> Either MathError (Vector a)
cross :: forall a.
Num a =>
Vector a -> Vector a -> Either MathError (Vector a)
cross (Vector [a
a, a
b, a
c]) (Vector [a
d, a
e, a
f]) = Vector a -> Either MathError (Vector a)
forall a b. b -> Either a b
Right ([a] -> Vector a
forall a. [a] -> Vector a
Vector [a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
f a -> a -> a
forall a. Num a => a -> a -> a
- a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
e, a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
f, a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
e a -> a -> a
forall a. Num a => a -> a -> a
- a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
d])
cross Vector a
_ Vector a
_ = MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"cross: vectors must have dimension three")

{- | Return the matrix of pairwise products a[i]*b[j], with shape (dimension a,dimension b).
Matrix construction can return shape/count errors.
-}
outer :: (Num a) => Vector a -> Vector a -> Either MathError (Matrix a)
outer :: forall a.
Num a =>
Vector a -> Vector a -> Either MathError (Matrix a)
outer (Vector [a]
a) (Vector [a]
b) = Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a) ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
b) [a
x a -> a -> a
forall a. Num a => a -> a -> a
* a
y | a
x <- [a]
a, a
y <- [a]
b]

{- | Given a,b,c, compute scalar triple product a dot (b cross c). All three vectors must
have dimension three; dimension errors propagate.
-}
triple :: (Num a) => Vector a -> Vector a -> Vector a -> Either MathError a
triple :: forall a.
Num a =>
Vector a -> Vector a -> Vector a -> Either MathError a
triple Vector a
a Vector a
b Vector a
c = Vector a -> Vector a -> Either MathError (Vector a)
forall a.
Num a =>
Vector a -> Vector a -> Either MathError (Vector a)
cross Vector a
b Vector a
c Either MathError (Vector a)
-> (Vector a -> Either MathError a) -> Either MathError a
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Vector a -> Vector a -> Either MathError a
forall a. Num a => Vector a -> Vector a -> Either MathError a
vectorDot Vector a
a

{- | Count nonzero coordinates, returning Int. This is the usual sparsity count rather than
a norm satisfying homogeneity.
-}
norm0 :: (Eq a, Num a) => Vector a -> Int
norm0 :: forall a. (Eq a, Num a) => Vector a -> Int
norm0 = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> Int) -> (Vector a -> [a]) -> Vector a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0) ([a] -> [a]) -> (Vector a -> [a]) -> Vector a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
vectorValues

{- | Sum absolute coordinate values. Empty vector returns zero; arithmetic follows the
chosen Num type.
-}
norm1 :: (Num a) => Vector a -> a
norm1 :: forall a. Num a => Vector a -> a
norm1 = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> (Vector a -> [a]) -> Vector a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs ([a] -> [a]) -> (Vector a -> [a]) -> Vector a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
vectorValues

{- | Compute Euclidean length using a scaled stable norm to reduce intermediate
overflow/underflow. Empty vector returns zero; the result itself is not wrapped in
Either.
-}
norm2 :: (RealFloat a) => Vector a -> a
norm2 :: forall a. RealFloat a => Vector a -> a
norm2 = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm ([a] -> a) -> (Vector a -> [a]) -> Vector a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> [a]
forall a. Vector a -> [a]
vectorValues

{- | Return maximum absolute coordinate, or zero for an empty vector. NaN coordinates
propagate rather than being silently discarded.
-}
normInfinity :: (Ord a, Num a) => Vector a -> a
normInfinity :: forall a. (Ord a, Num a) => Vector a -> a
normInfinity Vector a
v = case [a
x | a
x <- Vector a -> [a]
forall a. Vector a -> [a]
vectorValues Vector a
v, a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
x] of
    a
x : [a]
_ -> a
x
    [] -> [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs (Vector a -> [a]
forall a. Vector a -> [a]
vectorValues Vector a
v))

{- | Given positive exponent p then finite vector, compute the scaled p-norm. Positive
infinity selects normInfinity. NaN/nonpositive p or non-finite coordinates fail; p\<1 is
accepted as a quasi-norm.
-}
normP :: (RealFloat a) => a -> Vector a -> Either MathError a
normP :: forall a. RealFloat a => a -> Vector a -> Either MathError a
normP a
p v :: Vector a
v@(Vector [a]
xs)
    | a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
p Bool -> Bool -> Bool
|| a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 = MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"normP: p must be positive")
    | a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
p = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> a -> Either MathError a
forall a b. b -> Either a b
Right (Vector a -> a
forall a. (Ord a, Num a) => Vector a -> a
normInfinity Vector a
v)
    | Bool
otherwise = do
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs
        let s :: a
s = Vector a -> a
forall a. (Ord a, Num a) => Vector a -> a
normInfinity Vector a
v
        a -> Either MathError a
forall a b. b -> Either a b
Right (if a
s a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
0 else a
s a -> a -> a
forall a. Num a => a -> a -> a
* ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a -> a
forall a. Num a => a -> a
abs (a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) a -> a -> a
forall a. Floating a => a -> a -> a
** a
p | a
x <- [a]
xs]) a -> a -> a
forall a. Floating a => a -> a -> a
** (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
p))

{- | Return a unit-length vector in the same direction. Reject non-finite coordinates, zero
norm, or non-finite resulting norm. Empty input is a zero vector and fails.
-}
normalizeVector :: (RealFloat a) => Vector a -> Either MathError (Vector a)
normalizeVector :: forall a. RealFloat a => Vector a -> Either MathError (Vector a)
normalizeVector v :: Vector a
v@(Vector [a]
xs) = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs
    let n :: a
n = Vector a -> a
forall a. RealFloat a => Vector a -> a
norm2 Vector a
v
    if a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
        then MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"normalizeVector: zero vector")
        else if a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
n then MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult else Vector a -> Either MathError (Vector a)
forall a b. b -> Either a b
Right (a -> Vector a -> Vector a
forall a. Num a => a -> Vector a -> Vector a
vectorScale (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
n) Vector a
v)

{- | Given a then b, return the projection of a onto b. Dimensions must agree and b must
have nonzero finite norm; non-finite input/result fails.
-}
project :: (RealFloat a) => Vector a -> Vector a -> Either MathError (Vector a)
project :: forall a.
RealFloat a =>
Vector a -> Vector a -> Either MathError (Vector a)
project Vector a
a Vector a
b = do
    u <- Vector a -> Either MathError (Vector a)
forall a. RealFloat a => Vector a -> Either MathError (Vector a)
normalizeVector Vector a
b
    s <- vectorDot a u
    finiteValues (vectorValues a)
    checkedVector (vectorScale s u)

checkedVector :: (RealFloat a) => Vector a -> Either MathError (Vector a)
checkedVector :: forall a. RealFloat a => Vector a -> Either MathError (Vector a)
checkedVector Vector a
v
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\a
x -> a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x) (Vector a -> [a]
forall a. Vector a -> [a]
vectorValues Vector a
v) = MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | Bool
otherwise = Vector a -> Either MathError (Vector a)
forall a b. b -> Either a b
Right Vector a
v

{- | Given incident vector then surface normal, return incident minus twice its projection
onto the normal. The normal is normalized internally; invalid dimensions, zero normal,
or non-finite values fail.
-}
reflect :: (RealFloat a) => Vector a -> Vector a -> Either MathError (Vector a)
reflect :: forall a.
RealFloat a =>
Vector a -> Vector a -> Either MathError (Vector a)
reflect Vector a
incident Vector a
normal = do
    p <- Vector a -> Vector a -> Either MathError (Vector a)
forall a.
RealFloat a =>
Vector a -> Vector a -> Either MathError (Vector a)
project Vector a
incident Vector a
normal
    vectorSubtract incident (vectorScale 2 p) >>= checkedVector

{- | Given positive finite refractive-index ratio, incident vector, then normal, return the
refracted unit-direction expression. Normalize both inputs; reject invalid dimensions,
zero vectors, non-finite values, or total internal reflection.
-}
refract :: (RealFloat a) => a -> Vector a -> Vector a -> Either MathError (Vector a)
refract :: forall a.
RealFloat a =>
a -> Vector a -> Vector a -> Either MathError (Vector a)
refract a
eta Vector a
incident Vector a
normal
    | a
eta a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
eta Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
eta =
        MathError -> Either MathError (Vector a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"refract: positive finite index ratio required")
    | Bool
otherwise = do
        i <- Vector a -> Either MathError (Vector a)
forall a. RealFloat a => Vector a -> Either MathError (Vector a)
normalizeVector Vector a
incident
        n <- normalizeVector normal
        d <- vectorDot i n
        let k = a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
eta a -> a -> a
forall a. Num a => a -> a -> a
* a
eta a -> a -> a
forall a. Num a => a -> a -> a
* (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
d a -> a -> a
forall a. Num a => a -> a -> a
* a
d)
        if k < 0
            then Left (DomainError "refract: total internal reflection")
            else vectorSubtract (vectorScale eta i) (vectorScale (eta * d + sqrt k) n) >>= checkedVector

{- | Return Euclidean distance between equal-sized finite vectors. Dimension mismatch and
non-finite input/intermediate/result return Left MathError.
-}
distance :: (RealFloat a) => Vector a -> Vector a -> Either MathError a
distance :: forall a. RealFloat a => Vector a -> Vector a -> Either MathError a
distance Vector a
a Vector a
b = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues (Vector a -> [a]
forall a. Vector a -> [a]
vectorValues Vector a
a [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ Vector a -> [a]
forall a. Vector a -> [a]
vectorValues Vector a
b)
    difference <- Vector a -> Vector a -> Either MathError (Vector a)
forall a.
Num a =>
Vector a -> Vector a -> Either MathError (Vector a)
vectorSubtract Vector a
a Vector a
b Either MathError (Vector a)
-> (Vector a -> Either MathError (Vector a))
-> Either MathError (Vector a)
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Vector a -> Either MathError (Vector a)
forall a. RealFloat a => Vector a -> Either MathError (Vector a)
checkedVector
    let result = Vector a -> a
forall a. RealFloat a => Vector a -> a
norm2 Vector a
difference
    if isNaN result || isInfinite result then Left NonFiniteResult else Right result

{- | Return normalized dot product clamped to [-1,1]. Equal dimensions and nonzero finite
vectors are required; validation errors propagate.
-}
cosineSimilarity :: (RealFloat a) => Vector a -> Vector a -> Either MathError a
cosineSimilarity :: forall a. RealFloat a => Vector a -> Vector a -> Either MathError a
cosineSimilarity Vector a
a Vector a
b = do
    u <- Vector a -> Either MathError (Vector a)
forall a. RealFloat a => Vector a -> Either MathError (Vector a)
normalizeVector Vector a
a
    v <- normalizeVector b
    max (-1) . min 1 <$> vectorDot u v

{- | Return the angle in radians, acos of cosineSimilarity, between two nonzero finite
equally sized vectors. Invalid dimensions or zero/non-finite vectors fail.
-}
angle :: (RealFloat a) => Vector a -> Vector a -> Either MathError a
angle :: forall a. RealFloat a => Vector a -> Vector a -> Either MathError a
angle Vector a
a Vector a
b = a -> a
forall a. Floating a => a -> a
acos (a -> a) -> Either MathError a -> Either MathError a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Vector a -> Vector a -> Either MathError a
forall a. RealFloat a => Vector a -> Vector a -> Either MathError a
cosineSimilarity Vector a
a Vector a
b

{- | Given finite nonnegative relative tolerance then vectors, perform two-pass modified
Gram-Schmidt, returning an orthonormal basis in input order. Reject unequal dimensions,
zero/dependent vectors (RankDeficient), invalid tolerance, or non-finite data.
-}
gramSchmidt :: (RealFloat a) => a -> [Vector a] -> Either MathError [Vector a]
gramSchmidt :: forall a.
RealFloat a =>
a -> [Vector a] -> Either MathError [Vector a]
gramSchmidt a
tolerance [Vector a]
vs
    | a
tolerance a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tolerance Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tolerance = MathError -> Either MathError [Vector a]
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Bool
otherwise = (Either MathError [Vector a]
 -> Vector a -> Either MathError [Vector a])
-> Either MathError [Vector a]
-> [Vector a]
-> Either MathError [Vector a]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' Either MathError [Vector a]
-> Vector a -> Either MathError [Vector a]
step ([Vector a] -> Either MathError [Vector a]
forall a b. b -> Either a b
Right []) [Vector a]
vs
  where
    step :: Either MathError [Vector a]
-> Vector a -> Either MathError [Vector a]
step Either MathError [Vector a]
acc Vector a
v = do
        basis <- Either MathError [Vector a]
acc
        -- Two modified Gram-Schmidt passes reduce loss of orthogonality.
        w <- Foldable.foldl' remove (Right v) (basis ++ basis)
        if norm2 w <= tolerance * max 1 (norm2 v)
            then Left (RankDeficient (length basis))
            else do u <- normalizeVector w; Right (basis ++ [u])
    remove :: Either MathError (Vector a)
-> Vector a -> Either MathError (Vector a)
remove Either MathError (Vector a)
acc Vector a
u = do w <- Either MathError (Vector a)
acc; p <- project w u; vectorSubtract w p