-- SPDX-License-Identifier: MPL-2.0

{- | Matrix-free sparse power algorithms. Neither a dense matrix nor A^T A is
formed. Sparse multiplication, explicit vector workspaces and orthogonal
deflation determine storage. Power convergence requires a starting component
in the desired singular/eigen direction and a separated dominant magnitude;
an exact residual certifies an eigen/singular pair, not spectral dominance.
Every returned pair has residuals checked against the original sparse input.
See the power-method SVD derivation in
<https://ee263.stanford.edu/archive/ee263_course_reader.pdf Stanford EE263, exercise 15.10>.
-}
module HMath.Sparse.Spectral (
    SparseSpectralOptions (..),
    defaultSparseSpectralOptions,
    SparseEigenResult (..),
    sparseDominantEigen,
    SparseSVD (..),
    sparseTruncatedSVD,
) where

import Control.Monad (foldM)
import qualified Data.Map.Strict as Map
import HMath.Error (MathError (..))
import HMath.Matrix.Internal (dotV, finiteValues, stableNorm)
import HMath.Sparse (Sparse, sparseEntries, sparseShape, sparseTranspose, spmv)

{- | Residual tolerances, per-pair iteration budget and an explicit workspace
limit measured in scalar elements. Tolerances are finite and nonnegative.
-}
data SparseSpectralOptions a = SparseSpectralOptions
    { forall a. SparseSpectralOptions a -> a
sparseAbsoluteTolerance :: a
    -- ^ Absolute norm tolerance, finite and nonnegative.
    , forall a. SparseSpectralOptions a -> a
sparseRelativeTolerance :: a
    -- ^ Relative tolerance against max(abs eigen/singular value,maxAbs input).
    , forall a. SparseSpectralOptions a -> Int
sparseIterationLimit :: Int
    -- ^ Nonnegative number of power updates per pair.
    , forall a. SparseSpectralOptions a -> Integer
sparseWorkspaceLimit :: Integer
    -- ^ Maximum estimated simultaneous vector elements, nonnegative.
    }
    deriving (SparseSpectralOptions a -> SparseSpectralOptions a -> Bool
(SparseSpectralOptions a -> SparseSpectralOptions a -> Bool)
-> (SparseSpectralOptions a -> SparseSpectralOptions a -> Bool)
-> Eq (SparseSpectralOptions a)
forall a.
Eq a =>
SparseSpectralOptions a -> SparseSpectralOptions a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a.
Eq a =>
SparseSpectralOptions a -> SparseSpectralOptions a -> Bool
== :: SparseSpectralOptions a -> SparseSpectralOptions a -> Bool
$c/= :: forall a.
Eq a =>
SparseSpectralOptions a -> SparseSpectralOptions a -> Bool
/= :: SparseSpectralOptions a -> SparseSpectralOptions a -> Bool
Eq, Int -> SparseSpectralOptions a -> ShowS
[SparseSpectralOptions a] -> ShowS
SparseSpectralOptions a -> String
(Int -> SparseSpectralOptions a -> ShowS)
-> (SparseSpectralOptions a -> String)
-> ([SparseSpectralOptions a] -> ShowS)
-> Show (SparseSpectralOptions a)
forall a. Show a => Int -> SparseSpectralOptions a -> ShowS
forall a. Show a => [SparseSpectralOptions a] -> ShowS
forall a. Show a => SparseSpectralOptions a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> SparseSpectralOptions a -> ShowS
showsPrec :: Int -> SparseSpectralOptions a -> ShowS
$cshow :: forall a. Show a => SparseSpectralOptions a -> String
show :: SparseSpectralOptions a -> String
$cshowList :: forall a. Show a => [SparseSpectralOptions a] -> ShowS
showList :: [SparseSpectralOptions a] -> ShowS
Show)

{- | Absolute tolerance 1e-12, relative tolerance 1e-10, 1000 iterations per
pair and a workspace limit of ten million scalar elements.
-}
defaultSparseSpectralOptions :: (Fractional a) => SparseSpectralOptions a
defaultSparseSpectralOptions :: forall a. Fractional a => SparseSpectralOptions a
defaultSparseSpectralOptions = a -> a -> Int -> Integer -> SparseSpectralOptions a
forall a. a -> a -> Int -> Integer -> SparseSpectralOptions a
SparseSpectralOptions a
1e-12 a
1e-10 Int
1000 Integer
10000000

-- | A normalized real symmetric eigenpair and its independently checked residual.
data SparseEigenResult a = SparseEigenResult
    { forall a. SparseEigenResult a -> a
sparseEigenvalue :: a
    -- ^ Rayleigh quotient; negative dominant eigenvalues retain their sign.
    , forall a. SparseEigenResult a -> [a]
sparseEigenvector :: [a]
    -- ^ Unit Euclidean norm eigenvector.
    , forall a. SparseEigenResult a -> a
sparseEigenResidual :: a
    -- ^ Original norm ||A v-lambda v||_2.
    , forall a. SparseEigenResult a -> Int
sparseEigenIterations :: Int
    -- ^ Number of completed power updates.
    }
    deriving (SparseEigenResult a -> SparseEigenResult a -> Bool
(SparseEigenResult a -> SparseEigenResult a -> Bool)
-> (SparseEigenResult a -> SparseEigenResult a -> Bool)
-> Eq (SparseEigenResult a)
forall a.
Eq a =>
SparseEigenResult a -> SparseEigenResult a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a.
Eq a =>
SparseEigenResult a -> SparseEigenResult a -> Bool
== :: SparseEigenResult a -> SparseEigenResult a -> Bool
$c/= :: forall a.
Eq a =>
SparseEigenResult a -> SparseEigenResult a -> Bool
/= :: SparseEigenResult a -> SparseEigenResult a -> Bool
Eq, Int -> SparseEigenResult a -> ShowS
[SparseEigenResult a] -> ShowS
SparseEigenResult a -> String
(Int -> SparseEigenResult a -> ShowS)
-> (SparseEigenResult a -> String)
-> ([SparseEigenResult a] -> ShowS)
-> Show (SparseEigenResult a)
forall a. Show a => Int -> SparseEigenResult a -> ShowS
forall a. Show a => [SparseEigenResult a] -> ShowS
forall a. Show a => SparseEigenResult a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> SparseEigenResult a -> ShowS
showsPrec :: Int -> SparseEigenResult a -> ShowS
$cshow :: forall a. Show a => SparseEigenResult a -> String
show :: SparseEigenResult a -> String
$cshowList :: forall a. Show a => [SparseEigenResult a] -> ShowS
showList :: [SparseEigenResult a] -> ShowS
Show)

{- | Symmetric power iteration targeting the eigenvalue of largest magnitude.
Exact symmetry and finite input are required. The caller supplies a nonzero
finite initial vector; its projection onto the desired eigenspace must be
nonzero. An initial eigenvector can converge with a zero budget. Equal opposite
dominant eigenvalues may cause non-convergence. Returned residuals do not prove
that the converged eigenvalue is the largest when the initial projection is zero.
-}
sparseDominantEigen ::
    (RealFloat a) =>
    SparseSpectralOptions a ->
    Sparse a ->
    [a] ->
    Either MathError (SparseEigenResult a)
sparseDominantEigen :: forall a.
RealFloat a =>
SparseSpectralOptions a
-> Sparse a -> [a] -> Either MathError (SparseEigenResult a)
sparseDominantEigen SparseSpectralOptions a
options Sparse a
matrix [a]
initial = do
    SparseSpectralOptions a -> Either MathError ()
forall a.
RealFloat a =>
SparseSpectralOptions a -> Either MathError ()
validateOptions SparseSpectralOptions a
options
    let (Int
m, Int
n) = Sparse a -> (Int, Int)
forall a. Sparse a -> (Int, Int)
sparseShape Sparse a
matrix
    if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError ()
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n)) else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"empty matrix has no eigenvector") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    SparseSpectralOptions a -> Integer -> Either MathError ()
forall a. SparseSpectralOptions a -> Integer -> Either MathError ()
checkWorkspace SparseSpectralOptions a
options (Integer
6 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n)
    Sparse a -> Either MathError ()
forall a. RealFloat a => Sparse a -> Either MathError ()
validateSparse Sparse a
matrix
    if [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch Int
n ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial)) else [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
initial
    let values :: Map (Int, Int) a
values = [((Int, Int), a)] -> Map (Int, Int) a
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [((Int
i, Int
j), a
v) | (Int
i, Int
j, a
v) <- Sparse a -> [(Int, Int, a)]
forall a. Sparse a -> [(Int, Int, a)]
sparseEntries Sparse a
matrix]
    if (((Int, Int), a) -> Bool) -> [((Int, Int), a)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\((Int
i, Int
j), a
v) -> a -> (Int, Int) -> Map (Int, Int) a -> a
forall k a. Ord k => a -> k -> Map k a -> a
Map.findWithDefault a
0 (Int
j, Int
i) Map (Int, Int) a
values a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
v) (Map (Int, Int) a -> [((Int, Int), a)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (Int, Int) a
values)
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"sparseDominantEigen requires exact symmetry")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    q <- [a] -> Either MathError [a]
forall a. RealFloat a => [a] -> Either MathError [a]
unitVector [a]
initial
    go 0 q
  where
    scale :: a
scale = Sparse a -> a
forall a. (Num a, Ord a) => Sparse a -> a
inputScale Sparse a
matrix
    go :: Int -> [a] -> Either MathError (SparseEigenResult a)
go Int
k [a]
q = do
        aq <- Sparse a -> [a] -> Either MathError [a]
forall a. RealFloat a => Sparse a -> [a] -> Either MathError [a]
checkedProduct Sparse a
matrix [a]
q
        let lambda = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
q [a]
aq
            r = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
a a
v -> a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
lambda a -> a -> a
forall a. Num a => a -> a -> a
* a
v) [a]
aq [a]
q
        finiteValues (lambda : r)
        rn <- checkedNorm r
        if within options (max scale (abs lambda)) rn
            then Right (SparseEigenResult lambda q rn k)
            else
                if k == sparseIterationLimit options
                    then Left (NonConvergence "sparse symmetric power iteration" k)
                    else do
                        next <- unitVector aq
                        go (k + 1) next

{- | Compact singular triplets. Left and right vectors are stored as lists of
columns, requiring O(k*(rows+columns)) storage instead of a dense matrix.
-}
data SparseSVD a = SparseSVD
    { forall a. SparseSVD a -> [a]
sparseSingularValues :: [a]
    -- ^ Nonnegative singular estimates in extraction order.
    , forall a. SparseSVD a -> [[a]]
sparseLeftVectors :: [[a]]
    -- ^ Left unit vectors, one list per singular triplet.
    , forall a. SparseSVD a -> [[a]]
sparseRightVectors :: [[a]]
    -- ^ Right unit vectors (columns of V, rather than rows of V^T).
    , forall a. SparseSVD a -> [(a, a)]
sparseSingularResiduals :: [(a, a)]
    -- ^ Original (||A v-s u||_2,||A^T u-s v||_2) residuals.
    , forall a. SparseSVD a -> [Int]
sparseSingularIterations :: [Int]
    -- ^ Power updates used by each triplet.
    }
    deriving (SparseSVD a -> SparseSVD a -> Bool
(SparseSVD a -> SparseSVD a -> Bool)
-> (SparseSVD a -> SparseSVD a -> Bool) -> Eq (SparseSVD a)
forall a. Eq a => SparseSVD a -> SparseSVD a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => SparseSVD a -> SparseSVD a -> Bool
== :: SparseSVD a -> SparseSVD a -> Bool
$c/= :: forall a. Eq a => SparseSVD a -> SparseSVD a -> Bool
/= :: SparseSVD a -> SparseSVD a -> Bool
Eq, Int -> SparseSVD a -> ShowS
[SparseSVD a] -> ShowS
SparseSVD a -> String
(Int -> SparseSVD a -> ShowS)
-> (SparseSVD a -> String)
-> ([SparseSVD a] -> ShowS)
-> Show (SparseSVD a)
forall a. Show a => Int -> SparseSVD a -> ShowS
forall a. Show a => [SparseSVD a] -> ShowS
forall a. Show a => SparseSVD a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> SparseSVD a -> ShowS
showsPrec :: Int -> SparseSVD a -> ShowS
$cshow :: forall a. Show a => SparseSVD a -> String
show :: SparseSVD a -> String
$cshowList :: forall a. Show a => [SparseSVD a] -> ShowS
showList :: [SparseSVD a] -> ShowS
Show)

{- | Truncated SVD by two-sided power iteration and two-pass orthogonal deflation.
The seed initializes a local deterministic stream, with no global RNG state.
Both original singular residuals must meet the configured tolerance. The
algorithm targets descending singular values if starting vectors have nonzero
projection and the successive singular values are separated; it does not claim
a certified rank or ordering when these power-method assumptions fail.
Exact null images report RankDeficient, rather than inventing zero singular
vectors. Small positive singular values may be returned within tolerance;
requesting more vectors than the numerical rank can also exhaust the budget.
k=0 returns empty factors, including for empty/zero inputs.
-}
sparseTruncatedSVD ::
    (RealFloat a) =>
    SparseSpectralOptions a ->
    Int ->
    Integer ->
    Sparse a ->
    Either MathError (SparseSVD a)
sparseTruncatedSVD :: forall a.
RealFloat a =>
SparseSpectralOptions a
-> Int -> Integer -> Sparse a -> Either MathError (SparseSVD a)
sparseTruncatedSVD SparseSpectralOptions a
options Int
requested Integer
seed Sparse a
matrix = do
    SparseSpectralOptions a -> Either MathError ()
forall a.
RealFloat a =>
SparseSpectralOptions a -> Either MathError ()
validateOptions SparseSpectralOptions a
options
    let (Int
m, Int
n) = Sparse a -> (Int, Int)
forall a. Sparse a -> (Int, Int)
sparseShape Sparse a
matrix
    if Int
requested Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
requested Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"invalid sparse SVD rank")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    SparseSpectralOptions a -> Integer -> Either MathError ()
forall a. SparseSpectralOptions a -> Integer -> Either MathError ()
checkWorkspace
        SparseSpectralOptions a
options
        (if Int
requested Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Integer
0 else (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
requested Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
10) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n))
    Sparse a -> Either MathError ()
forall a. RealFloat a => Sparse a -> Either MathError ()
validateSparse Sparse a
matrix
    (triplets, _) <- (([([a], a, [a], (a, a), Int)], Integer)
 -> Int -> Either MathError ([([a], a, [a], (a, a), Int)], Integer))
-> ([([a], a, [a], (a, a), Int)], Integer)
-> [Int]
-> Either MathError ([([a], a, [a], (a, a), Int)], Integer)
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM ([([a], a, [a], (a, a), Int)], Integer)
-> Int -> Either MathError ([([a], a, [a], (a, a), Int)], Integer)
extract ([], Integer
seed Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
4294967296) [Int
0 .. Int
requested Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    pure
        ( SparseSVD
            [s | (_, s, _, _, _) <- triplets]
            [u | (u, _, _, _, _) <- triplets]
            [v | (_, _, v, _, _) <- triplets]
            [r | (_, _, _, r, _) <- triplets]
            [k | (_, _, _, _, k) <- triplets]
        )
  where
    transposed :: Sparse a
transposed = Sparse a -> Sparse a
forall a. Sparse a -> Sparse a
sparseTranspose Sparse a
matrix
    (Int
_, Int
columns) = Sparse a -> (Int, Int)
forall a. Sparse a -> (Int, Int)
sparseShape Sparse a
matrix
    scale :: a
scale = Sparse a -> a
forall a. (Num a, Ord a) => Sparse a -> a
inputScale Sparse a
matrix
    extract :: ([([a], a, [a], (a, a), Int)], Integer)
-> Int -> Either MathError ([([a], a, [a], (a, a), Int)], Integer)
extract ([([a], a, [a], (a, a), Int)]
prior, Integer
state) Int
index = do
        let ([a]
candidate, Integer
nextSeed) = Integer -> Int -> ([a], Integer)
forall a. Fractional a => Integer -> Int -> ([a], Integer)
randomVector Integer
state Int
columns
            lefts :: [[a]]
lefts = [[a]
u | ([a]
u, a
_, [a]
_, (a, a)
_, Int
_) <- [([a], a, [a], (a, a), Int)]
prior]
            rights :: [[a]]
rights = [[a]
v | ([a]
_, a
_, [a]
v, (a, a)
_, Int
_) <- [([a], a, [a], (a, a), Int)]
prior]
        start <- [a] -> Either MathError [a]
forall a. RealFloat a => [a] -> Either MathError [a]
unitVector ([[a]] -> [a] -> [a]
forall a. RealFloat a => [[a]] -> [a] -> [a]
projectOut [[a]]
rights [a]
candidate)
        triplet <- iteratePair index lefts rights 0 start
        pure (prior ++ [triplet], nextSeed)
    iteratePair :: Int
-> [[a]]
-> [[a]]
-> Int
-> [a]
-> Either MathError ([a], a, [a], (a, a), Int)
iteratePair Int
index [[a]]
lefts [[a]]
rights Int
k [a]
v
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== SparseSpectralOptions a -> Int
forall a. SparseSpectralOptions a -> Int
sparseIterationLimit SparseSpectralOptions a
options = MathError -> Either MathError ([a], a, [a], (a, a), Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"sparse two-sided SVD" Int
k)
        | Bool
otherwise = do
            av <- Sparse a -> [a] -> Either MathError [a]
forall a. RealFloat a => Sparse a -> [a] -> Either MathError [a]
checkedProduct Sparse a
matrix [a]
v
            let left = [[a]] -> [a] -> [a]
forall a. RealFloat a => [[a]] -> [a] -> [a]
projectOut [[a]]
lefts [a]
av
            if all (== 0) left then Left (RankDeficient index) else Right ()
            u <- unitVector left
            atu <- checkedProduct transposed u
            let right = [[a]] -> [a] -> [a]
forall a. RealFloat a => [[a]] -> [a] -> [a]
projectOut [[a]]
rights [a]
atu
            if all (== 0) right then Left (RankDeficient index) else Right ()
            vn <- unitVector right
            avn <- checkedProduct matrix vn
            let sigma = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
u [a]
avn
                leftResidual = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
a a
z -> a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
sigma a -> a -> a
forall a. Num a => a -> a -> a
* a
z) [a]
avn [a]
u
                rightResidual = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
a a
z -> a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
sigma a -> a -> a
forall a. Num a => a -> a -> a
* a
z) [a]
atu [a]
vn
            finiteValues (sigma : leftResidual ++ rightResidual)
            if sigma < 0
                then Left (DomainError "sparse SVD lost nonnegative singular orientation")
                else if sigma == 0 then Left (RankDeficient index) else Right ()
            lr <- checkedNorm leftResidual
            rr <- checkedNorm rightResidual
            if within options (max scale sigma) lr && within options (max scale sigma) rr
                then Right (u, sigma, vn, (lr, rr), k + 1)
                else iteratePair index lefts rights (k + 1) vn

validateOptions :: (RealFloat a) => SparseSpectralOptions a -> Either MathError ()
validateOptions :: forall a.
RealFloat a =>
SparseSpectralOptions a -> Either MathError ()
validateOptions SparseSpectralOptions a
options
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any
        (\a
x -> Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
finite a
x) Bool -> Bool -> Bool
|| a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0)
        [SparseSpectralOptions a -> a
forall a. SparseSpectralOptions a -> a
sparseAbsoluteTolerance SparseSpectralOptions a
options, SparseSpectralOptions a -> a
forall a. SparseSpectralOptions a -> a
sparseRelativeTolerance SparseSpectralOptions a
options] =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | SparseSpectralOptions a -> Int
forall a. SparseSpectralOptions a -> Int
sparseIterationLimit SparseSpectralOptions a
options Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative sparse spectral iteration limit")
    | SparseSpectralOptions a -> Integer
forall a. SparseSpectralOptions a -> Integer
sparseWorkspaceLimit SparseSpectralOptions a
options Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative sparse spectral workspace limit")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()

checkWorkspace :: SparseSpectralOptions a -> Integer -> Either MathError ()
checkWorkspace :: forall a. SparseSpectralOptions a -> Integer -> Either MathError ()
checkWorkspace SparseSpectralOptions a
options Integer
count
    | Integer
count Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> SparseSpectralOptions a -> Integer
forall a. SparseSpectralOptions a -> Integer
sparseWorkspaceLimit SparseSpectralOptions a
options Bool -> Bool -> Bool
|| Integer
count 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 ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"sparse spectral vector workspace limit exceeded")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()

validateSparse :: (RealFloat a) => Sparse a -> Either MathError ()
validateSparse :: forall a. RealFloat a => Sparse a -> Either MathError ()
validateSparse = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a] -> Either MathError ())
-> (Sparse a -> [a]) -> Sparse a -> Either MathError ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int, a) -> a) -> [(Int, Int, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
_, Int
_, a
v) -> a
v) ([(Int, Int, a)] -> [a])
-> (Sparse a -> [(Int, Int, a)]) -> Sparse a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Sparse a -> [(Int, Int, a)]
forall a. Sparse a -> [(Int, Int, a)]
sparseEntries

inputScale :: (Num a, Ord a) => Sparse a -> a
inputScale :: forall a. (Num a, Ord a) => Sparse a -> a
inputScale = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([a] -> a) -> (Sparse a -> [a]) -> Sparse a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [a]) -> (Sparse a -> [a]) -> Sparse a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int, a) -> a) -> [(Int, Int, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
_, Int
_, a
v) -> a -> a
forall a. Num a => a -> a
abs a
v) ([(Int, Int, a)] -> [a])
-> (Sparse a -> [(Int, Int, a)]) -> Sparse a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Sparse a -> [(Int, Int, a)]
forall a. Sparse a -> [(Int, Int, a)]
sparseEntries

checkedProduct :: (RealFloat a) => Sparse a -> [a] -> Either MathError [a]
checkedProduct :: forall a. RealFloat a => Sparse a -> [a] -> Either MathError [a]
checkedProduct Sparse a
a [a]
v = do
    productVector <- Sparse a -> [a] -> Either MathError [a]
forall a. Num a => Sparse a -> [a] -> Either MathError [a]
spmv Sparse a
a [a]
v
    finiteValues productVector
    pure productVector

-- Normalize by max magnitude first. The true norm can overflow for entirely
-- finite input; scaling before normalization avoids converting it to a zero vector.
unitVector :: (RealFloat a) => [a] -> Either MathError [a]
unitVector :: forall a. RealFloat a => [a] -> Either MathError [a]
unitVector [a]
xs = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs
    let scale :: a
scale = [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 [a]
xs)
    if a
scale a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
        then MathError -> Either MathError [a]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero sparse spectral starting or power vector")
        else do
            let scaled :: [a]
scaled = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale) [a]
xs
                norm :: a
norm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
scaled
                q :: [a]
q = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
norm) [a]
scaled
            [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
q
            [a] -> Either MathError [a]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure [a]
q

checkedNorm :: (RealFloat a) => [a] -> Either MathError a
checkedNorm :: forall a. RealFloat a => [a] -> Either MathError a
checkedNorm [a]
xs = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs
    let norm :: a
norm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
xs
    if a -> Bool
forall a. RealFloat a => a -> Bool
finite a
norm then a -> Either MathError a
forall a b. b -> Either a b
Right a
norm else MathError -> Either MathError a
forall a b. a -> Either a b
Left MathError
NonFiniteResult

within :: (RealFloat a) => SparseSpectralOptions a -> a -> a -> Bool
within :: forall a. RealFloat a => SparseSpectralOptions a -> a -> a -> Bool
within SparseSpectralOptions a
options a
scale a
residual
    | a
residual a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= SparseSpectralOptions a -> a
forall a. SparseSpectralOptions a -> a
sparseAbsoluteTolerance SparseSpectralOptions a
options = Bool
True
    | a
scale a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = Bool
False
    | Bool
otherwise =
        a
residual a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= SparseSpectralOptions a -> a
forall a. SparseSpectralOptions a -> a
sparseRelativeTolerance SparseSpectralOptions a
options a -> a -> a
forall a. Num a => a -> a -> a
+ SparseSpectralOptions a -> a
forall a. SparseSpectralOptions a -> a
sparseAbsoluteTolerance SparseSpectralOptions a
options a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale

projectOut :: (RealFloat a) => [[a]] -> [a] -> [a]
projectOut :: forall a. RealFloat a => [[a]] -> [a] -> [a]
projectOut [[a]]
basis = [a] -> [a]
pass ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
pass
  where
    pass :: [a] -> [a]
pass [a]
xs =
        ([a] -> [a] -> [a]) -> [a] -> [[a]] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
            ( \[a]
v [a]
q ->
                let coefficient :: a
coefficient = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
q [a]
v
                 in (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
a a
b -> a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
coefficient a -> a -> a
forall a. Num a => a -> a -> a
* a
b) [a]
v [a]
q
            )
            [a]
xs
            [[a]]
basis

randomVector :: (Fractional a) => Integer -> Int -> ([a], Integer)
randomVector :: forall a. Fractional a => Integer -> Int -> ([a], Integer)
randomVector Integer
seed Int
n =
    let states :: [Integer]
states = Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
take Int
n (Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
drop Int
1 ((Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (\Integer
x -> (Integer
1664525 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1013904223) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
4294967296) Integer
seed))
        values :: [a]
values = [(Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.5) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
4294967296 a -> a -> a
forall a. Num a => a -> a -> a
- a
0.5 | Integer
x <- [Integer]
states]
     in ([a]
values, case [Integer] -> [Integer]
forall a. [a] -> [a]
reverse [Integer]
states of Integer
x : [Integer]
_ -> Integer
x; [] -> Integer
seed)

finite :: (RealFloat a) => a -> Bool
finite :: forall a. RealFloat a => a -> Bool
finite a
x = Bool -> Bool
not (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)