{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
module Test.Orthogonal (testsVar) where

import qualified Test.Divide as Divide
import qualified Test.Generator as Gen
import qualified Test.Utility as Util
import Test.Generator ((<#*#>), (<#\#>), (<-*#>), (<#*|>), (<|=|>))
import Test.Utility
         (approx, approxReal, approxArrayTol, approxMatrix, approxVectorTol,
          Tagged, isIdentity, isUnitary, maybeConjugate)

import qualified Numeric.LAPACK.Orthogonal.Householder as HH
import qualified Numeric.LAPACK.Orthogonal as Ortho
import qualified Numeric.LAPACK.Matrix.Hermitian as Herm
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Square (Square)
import Numeric.LAPACK.Matrix
         (General, ShapeInt, (#*#), (##*#), (#*##), (#\##), (#*|))
import Numeric.LAPACK.Vector (Vector, (|+|), (|-|))
import Numeric.LAPACK.Scalar (RealOf, absolute, selectReal)

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Shape ((:+:))

import Control.Applicative (liftA2, (<$>))
import Data.Semigroup ((<>))

import qualified Test.QuickCheck as QC


pseudoInverseProjection ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
pseudoInverseProjection a =
   let ainv = snd $ Ortho.pseudoInverseRCond 1e-5 a
       tol = selectReal 1e-1 1e-5
   in approxArrayTol tol a (a <> ainv <> a) &&
      approxArrayTol tol ainv (ainv <> a <> ainv)

pseudoInverseHermitian ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
pseudoInverseHermitian a =
   let ainv = snd $ Ortho.pseudoInverseRCond 1e-5 a
       tol = selectReal 1e-2 1e-5
       aainv = a <> ainv
       ainva = ainv <> a
   in approxMatrix tol aainv (Matrix.adjoint aainv) &&
      approxMatrix tol ainva (Matrix.adjoint ainva)

pseudoInverseFactored ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a,
    Matrix.Wide ShapeInt ShapeInt a) -> Bool
pseudoInverseFactored (a,b) =
   let pinv x = snd $ Ortho.pseudoInverseRCond 1e-5 x
   in approxMatrix (selectReal 1e-1 1e-5)
         (pinv (a #*# b)) (pinv b #*# pinv a)

pseudoInverseInverse ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Square ShapeInt a -> Bool
pseudoInverseInverse a =
   approxMatrix (selectReal 1e-1 1e-5)
      (Matrix.inverse a)
      (snd $ Ortho.pseudoInverseRCond 1e-5 a)


determinant ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Square ShapeInt a -> Bool
determinant a =
   let detSquare = Square.determinant a
       detOrtho = Ortho.determinant a
   in approx
         (1e-3 * max 1 (max (absolute detSquare) (absolute detOrtho)))
         detSquare detOrtho

determinantAbsolute ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Square ShapeInt a -> Bool
determinantAbsolute a =
   let det = absolute $ Ortho.determinant a
       detAbs = Ortho.determinantAbsolute a
   in approxReal (1e-5 * max 1 (max det detAbs)) det detAbs

gramianDeterminant ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   General ShapeInt ShapeInt a -> Bool
gramianDeterminant a =
   let gram = Herm.gramian a
       Shape.ZeroBased n = Matrix.width a
       estimate = (Vector.sum (Herm.takeDiagonal gram) / fromIntegral n) ^ n
   in approxReal (1e-5 * max 1 estimate)
         (Herm.determinant gram)
         (Ortho.determinantAbsolute a ^ (2::Int))


multiplyDeterminantRight ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (General ShapeInt ShapeInt a, Square ShapeInt a) -> Bool
multiplyDeterminantRight (a,b) =
   let detA = Ortho.determinantAbsolute a
       detB = absolute $ Ortho.determinant b
   in approxReal
         (selectReal 1e-1 1e-5 * max 1 detA * max 1 detB)
         (Ortho.determinantAbsolute (a##*#b))
         (detA * detB)

multiplyDeterminantLeft ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Square ShapeInt a, General ShapeInt ShapeInt a) -> Bool
multiplyDeterminantLeft (a,b) =
   let detA = absolute $ Ortho.determinant a
       detB = Ortho.determinantAbsolute b
   in approxReal
         (selectReal 1e-1 1e-5 * max 1 detA * max 1 detB)
         (Ortho.determinantAbsolute (a#*##b))
         (detA * detB)


genFullRankTallRHS ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Gen.MatrixInt a
      (Matrix.Tall ShapeInt ShapeInt a,
       Matrix.General ShapeInt ShapeInt a)
genFullRankTallRHS = (,) <$> Gen.fullRankTall <#\#> Gen.matrix


normalEquationLeastSquares ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
normalEquationLeastSquares (a, b) =
   approxArrayTol
      (selectReal 10 1e-3)
      (Ortho.leastSquares a b)
      (Herm.solve (Herm.gramian $ Matrix.fromFull a) $
       Matrix.adjoint a #*# b)

specializedLeastSquares ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
specializedLeastSquares (a, b) =
   approxArrayTol
      (selectReal 1e-1 1e-5)
      (Ortho.leastSquares a b)
      (snd $ Ortho.leastSquaresMinimumNormRCond 1e-5 (Matrix.fromFull a) b)

householderLeastSquares ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
householderLeastSquares (a, b) =
   approxArrayTol
      (selectReal 1e-1 1e-5)
      (Ortho.leastSquares a b)
      (HH.leastSquares (HH.fromMatrix a) b)

triangularLeastSquares ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
triangularLeastSquares (a, b) =
   approxArrayTol
      (selectReal 1e-1 1e-5)
      (Ortho.leastSquares a b)
      (let (q,r) = Ortho.householderTall a
       in r #\## (Matrix.adjoint q #*# b))



genFullRankWideRHS ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Gen.MatrixInt a
      (Matrix.Wide ShapeInt ShapeInt a,
       Matrix.General ShapeInt ShapeInt a)
genFullRankWideRHS = (,) <$> Gen.fullRankWide <#\#> Gen.matrix


normalEquationMinimumNorm ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
normalEquationMinimumNorm (a, b) =
   approxArrayTol
      (selectReal 10 1e-3)
      (Ortho.minimumNorm a b)
      (Matrix.adjoint a #*#
       Herm.solve (Herm.gramian $ Matrix.fromFull $ Matrix.adjoint a) b)

specializedMinimumNorm ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
specializedMinimumNorm (a, b) =
   approxArrayTol
      (selectReal 1e-1 1e-5)
      (Ortho.minimumNorm a b)
      (snd $ Ortho.leastSquaresMinimumNormRCond 1e-5 (Matrix.fromFull a) b)

householderMinimumNorm ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
householderMinimumNorm (a, b) =
   approxArrayTol
      (selectReal 1e-1 1e-5)
      (Ortho.minimumNorm a b)
      (HH.minimumNorm (HH.fromMatrix $ Matrix.adjoint a) b)

triangularMinimumNorm ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
triangularMinimumNorm (a, b) =
   approxArrayTol
      (selectReal 1e-1 1e-5)
      (Ortho.minimumNorm a b)
      (let (q,r) = Ortho.householderTall $ Matrix.adjoint a
       in q #*# (Triangular.adjoint r #\## b))


complementDimension ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Tall ShapeInt ShapeInt a -> Bool
complementDimension a =
   let b = Matrix.fromFull a Matrix.||| Matrix.fromFull (Ortho.complement a)
   in Shape.size (Matrix.height b) == Shape.size (Matrix.width b)

complementBiorthogonal ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Tall ShapeInt ShapeInt a -> Bool
complementBiorthogonal a =
   all (approx 1e-3 0) $ Array.toList $ ArrMatrix.toVector $
   Matrix.adjoint a #*# Ortho.complement a

complementOrthogonal ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Tall ShapeInt ShapeInt a -> Bool
complementOrthogonal = isUnitary (selectReal 1e-3 1e-7) . Ortho.complement


affineSpanFromKernel ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
affineSpanFromKernel (a, by) =
   let b = Vector.take (Shape.size $ Matrix.height a) by
       y = Vector.drop (Shape.size $ Matrix.height a) by
       (c,d) = Ortho.affineSpanFromKernel a b
   in approxVectorTol
         (selectReal 1e-3 1e-7)
         b
         (a#*|(c#*|y|+|d))

affineKernelFromSpan ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Vector ShapeInt a, Vector ShapeInt a) ->
   Bool
affineKernelFromSpan (c,y,d) =
   let (a,b) = Ortho.affineKernelFromSpan c d
   in approxVectorTol
         (selectReal 1e-3 1e-7)
         b
         (a#*|(c#*|y|+|d))


projectHit ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt ShapeInt a,
    Vector ShapeInt a,
    Vector ShapeInt a) ->
   Bool
projectHit (b,x,d) =
   approxVectorTol
      (selectReal 1e-3 1e-9)
      d
      (b #*| Ortho.project b d x)


leastSquaresNoConstraint ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
leastSquaresNoConstraint (a, b) =
   approxVectorTol
      (selectReal 0.1 1e-7)
      (ArrMatrix.unliftColumn MatrixShape.ColumnMajor (Ortho.leastSquares a) b)
      (Ortho.leastSquaresConstraint
         (Matrix.fromFull a) b
         (Matrix.zero $
          MatrixShape.wide MatrixShape.ColumnMajor Shape.Zero (Matrix.width a))
         (Vector.zero Shape.Zero))

leastSquaresConstraintUnique ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Vector ShapeInt a,
    Matrix.General ShapeInt ShapeInt a, Matrix.Square ShapeInt a,
    Vector ShapeInt a) ->
   Bool
leastSquaresConstraintUnique (c,a,b,d) =
   approxVectorTol
      (selectReal 0.1 1e-7)
      d
      (b #*| Ortho.leastSquaresConstraint a c (Matrix.generalizeTall b) d)

splitLSCStack ::
   (Shape.C height, Shape.C constraints, Shape.C width, Class.Floating a) =>
   Matrix.Tall (constraints:+:height) width a ->
   Vector (constraints:+:height) a ->
   ((Matrix.General height width a, Vector height a),
    (Matrix.Wide constraints width a, Vector constraints a))
splitLSCStack baTall dc =
   let ba = Matrix.fromFull baTall
       b = Matrix.wideFromGeneral $ Matrix.takeTop ba
       a = Matrix.takeBottom ba
       (d,c) = (Vector.takeLeft dc, Vector.takeRight dc)
   in ((a,c),(b,d))

leastSquaresConstraintAdmissible ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a,
    Vector (ShapeInt:+:ShapeInt) a) ->
   Bool
leastSquaresConstraintAdmissible (ba,dc) =
   let ((a,c),(b,d)) = splitLSCStack ba dc
   in approxVectorTol
         (selectReal 0.1 1e-7)
         d
         (b #*| Ortho.leastSquaresConstraint a c b d)

leastSquaresConstraintMinimal ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a,
    Vector ShapeInt a,
    Vector (ShapeInt:+:ShapeInt) a) ->
   Bool
leastSquaresConstraintMinimal (ba,x,dc) =
   let ((a,c),(b,d)) = splitLSCStack ba dc
   in Vector.norm2 (c |-| a #*| Ortho.leastSquaresConstraint a c b d)
      <=
      Vector.norm2 (c |-| a #*| Ortho.project b d x) + selectReal 1e-1 1e-10


gaussMarkovLinearModelMinimumNorm ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar, Eq a) =>
   (Matrix.Wide ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
gaussMarkovLinearModelMinimumNorm (b, d) =
   let (x,y) =
         Ortho.gaussMarkovLinearModel
            (Matrix.zero $
             MatrixShape.tall MatrixShape.ColumnMajor
               (Matrix.height b) Shape.Zero)
            (Matrix.fromFull b) d
   in x == Vector.zero Shape.Zero
      &&
      approxVectorTol
         (selectReal 0.1 1e-7)
         y
         (ArrMatrix.unliftColumn MatrixShape.ColumnMajor
            (Ortho.minimumNorm b) d)

gaussMarkovLinearModelUnique ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar, Eq a) =>
   (Matrix.Square ShapeInt a, Vector ShapeInt a) -> Bool
gaussMarkovLinearModelUnique (a,d) =
   let (x,y) =
         Ortho.gaussMarkovLinearModel
            (Matrix.generalizeWide a)
            (Matrix.zero $
             MatrixShape.general MatrixShape.ColumnMajor
               (Matrix.height a) Shape.Zero)
            d
   in y == Vector.zero Shape.Zero
      &&
      approxVectorTol (selectReal 0.1 1e-7) d (a #*| x)

gaussMarkovLinearModelAdmissible ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt (ShapeInt:+:ShapeInt) a, Vector ShapeInt a) -> Bool
gaussMarkovLinearModelAdmissible (abWide,d) =
   let ab = Matrix.fromFull abWide
       a = Matrix.tallFromGeneral $ Matrix.takeLeft ab
       b = Matrix.takeRight ab
   in approxVectorTol
         (selectReal 0.1 1e-7)
         d
         (ab #*| uncurry Vector.append (Ortho.gaussMarkovLinearModel a b d))

gaussMarkovLinearModelMinimal ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Wide ShapeInt (ShapeInt:+:ShapeInt) a,
    Vector ShapeInt a,
    Vector (ShapeInt:+:ShapeInt) a) -> Bool
gaussMarkovLinearModelMinimal (abWide,d,xy) =
   let ab = Matrix.fromFull abWide
       a = Matrix.tallFromGeneral $ Matrix.takeLeft ab
       b = Matrix.takeRight ab
   in Vector.norm2 (snd $ Ortho.gaussMarkovLinearModel a b d)
      <=
      Vector.norm2 (Vector.takeRight $ Ortho.project abWide d xy)
         + selectReal 1e-3 1e-10


householderReconstruction ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.General ShapeInt ShapeInt a -> Bool
householderReconstruction a =
   approxArrayTol (selectReal 1e-3 1e-7)
      a (uncurry (#*##) (Ortho.householder a))

householderDeterminant ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Square ShapeInt a -> Bool
householderDeterminant a =
   let detOrtho = Ortho.determinant a
       detHH = HH.determinant $ HH.fromMatrix a
   in approx 1e-5 detOrtho detHH


maybeTriTranspose ::
   (Shape.C size, Class.Floating a, MatrixShape.TriDiag diag,
    MatrixShape.Content lo, MatrixShape.Content up) =>
   HH.Transposition ->
   Triangular.Triangular up diag lo size a -> Square size a
maybeTriTranspose HH.NonTransposed = Triangular.toSquare
maybeTriTranspose HH.Transposed = Triangular.toSquare . Triangular.transpose

householderSolveRR ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (HH.Transposition, HH.Conjugation) ->
   Matrix.Tall ShapeInt ShapeInt a -> Bool
householderSolveRR (trans,conj) a =
   let qr = HH.fromMatrix a
   in  isIdentity (selectReal 1e-3 1e-7) $
         HH.tallSolveR trans conj qr $
         maybeTriTranspose trans $ maybeConjugate conj $ HH.tallExtractR qr


householderMultiplyR ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   HH.Transposition ->
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) ->
   Bool
householderMultiplyR trans (a,b) =
   let qr = HH.fromMatrix a
       r = HH.tallExtractR qr
   in approxArrayTol
         (selectReal 1e-3 1e-7)
         (HH.tallMultiplyR trans qr b)
         (Matrix.multiplySquare trans r b)


householderQOrthogonal ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.General ShapeInt ShapeInt a -> Bool
householderQOrthogonal =
   isUnitary (selectReal 1e-3 1e-7) . HH.extractQ . HH.fromMatrix


householderMultiplyQ ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (HH.Transposition, HH.Conjugation) ->
   (Matrix.General ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) ->
   Bool
householderMultiplyQ (trans,conj) (a,b) =
   let qr = HH.fromMatrix a
   in approxArrayTol
         (selectReal 1e-3 1e-7)
         (Matrix.multiplySquare trans (maybeConjugate conj $ HH.extractQ qr) b)
         (HH.multiplyQ trans conj qr b)


householderTallQOrthogonal ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix.Tall ShapeInt ShapeInt a -> Bool
householderTallQOrthogonal =
   isUnitary (selectReal 1e-3 1e-7) . HH.tallExtractQ . HH.fromMatrix

householderTallMultiplyQ ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
householderTallMultiplyQ (a,b) =
   let qr = HH.fromMatrix a
   in approxArrayTol
         (selectReal 1e-3 1e-7)
         (HH.tallExtractQ qr #*# b)
         (HH.tallMultiplyQ qr b)

householderTallMultiplyQAdjoint ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Matrix.Tall ShapeInt ShapeInt a, Matrix.General ShapeInt ShapeInt a) -> Bool
householderTallMultiplyQAdjoint (a,b) =
   let qr = HH.fromMatrix a
   in approxArrayTol
         (selectReal 1e-3 1e-7)
         (Matrix.adjoint (HH.tallExtractQ qr) #*# b)
         (HH.tallMultiplyQAdjoint qr b)



checkForAll ::
   (Show a, QC.Testable test) =>
   Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property
checkForAll gen = Util.checkForAll (Gen.run gen 3 5)

checkForAllExtra ::
   (Show a, Show b, QC.Testable test) =>
   QC.Gen a -> Gen.T dim tag b ->
   (a -> b -> test) -> Tagged tag QC.Property
checkForAllExtra = Gen.withExtra checkForAll


testsVar ::
   (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
   [(String, Tagged a QC.Property)]
testsVar =
   ("pseudoInverseProjection",
      checkForAll Gen.matrix pseudoInverseProjection) :
   ("pseudoInverseHermitian",
      checkForAll Gen.matrix pseudoInverseHermitian) :
   ("pseudoInverseFactored",
      checkForAll
         ((,) <$> Gen.fullRankTall <#*#> Gen.fullRankWide)
         pseudoInverseFactored) :
   ("pseudoInverseInverse",
      checkForAll Gen.invertible pseudoInverseInverse) :

   ("determinant",
      checkForAll Gen.square determinant) :
   ("determinantAbsolute",
      checkForAll Gen.square determinantAbsolute) :
   ("gramianDeterminant",
      checkForAll Gen.matrix gramianDeterminant) :
   ("multiplyDeterminantRight",
      checkForAll
         ((,) <$> Gen.matrix <#*#> Gen.square) multiplyDeterminantRight) :
   ("multiplyDeterminantLeft",
      checkForAll
         ((,) <$> (fst . Ortho.householder <$> Gen.square) <#*#> Gen.matrix)
         multiplyDeterminantLeft) :
   ("normalEquationLeastSquares",
      checkForAll genFullRankTallRHS normalEquationLeastSquares) :
   ("normalEquationMinimumNorm",
      checkForAll genFullRankWideRHS normalEquationMinimumNorm) :
   ("specializedLeastSquares",
      checkForAll genFullRankTallRHS specializedLeastSquares) :
   ("specializedMinimumNorm",
      checkForAll genFullRankWideRHS specializedMinimumNorm) :

   ("complementDimension",
      checkForAll Gen.tall complementDimension) :
   ("complementBiorthogonal",
      checkForAll Gen.tall complementBiorthogonal) :
   ("complementOrthogonal",
      checkForAll Gen.tall complementOrthogonal) :

   ("affineSpanFromKernel",
      checkForAll
         ((,) <$> Gen.fullRankWide <#*|> Gen.vector) affineSpanFromKernel) :
   ("affineKernelFromSpan",
      checkForAll
         ((,,) <$> Gen.tall <#*|> Gen.vector <|=|> Gen.vector)
         affineKernelFromSpan) :

   ("projectHit",
      checkForAll
         ((,,) <$> Gen.fullRankWide <#*|> Gen.vector <|=|> Gen.vector)
         projectHit) :
   ("leastSquaresNoConstraint",
      checkForAll
         ((,) <$> Gen.transpose Gen.fullRankTall <#*|> Gen.vector)
         leastSquaresNoConstraint) :
   ("leastSquaresConstraintUnique",
      checkForAll
         ((,,,) <$>
            Gen.vector <-*#> Gen.matrix <-*#> Gen.invertible <|=|> Gen.vector)
         leastSquaresConstraintUnique) :
   ("leastSquaresConstraintAdmissible",
      checkForAll
         ((,) <$> Gen.transpose Gen.lscStack <#*|> Gen.vector)
         leastSquaresConstraintAdmissible) :
   ("leastSquaresConstraintMinimal",
      checkForAll
         ((,,) <$> Gen.lscStack <#*|> Gen.vector <|=|> Gen.vector)
         leastSquaresConstraintMinimal) :
   ("gaussMarkovLinearModelMinimumNorm",
      checkForAll
         ((,) <$> Gen.transpose Gen.fullRankWide <#*|> Gen.vector)
         gaussMarkovLinearModelMinimumNorm) :
   ("gaussMarkovLinearModelUnique",
      checkForAll
         ((,) <$> Gen.invertible <#*|> Gen.vector)
         gaussMarkovLinearModelUnique) :
   ("gaussMarkovLinearModelAdmissible",
      checkForAll
         ((,) <$> fmap Matrix.transpose Gen.lscStack <#*|> Gen.vector)
         gaussMarkovLinearModelAdmissible) :
   ("gaussMarkovLinearModelMinimal",
      checkForAll
         ((,,) <$> fmap Matrix.transpose Gen.lscStack
               <#*|> Gen.vector <|=|> Gen.vector)
         gaussMarkovLinearModelMinimal) :

   ("triangularLeastSquares",
      checkForAll genFullRankTallRHS triangularLeastSquares) :
   ("triangularMinimumNorm",
      checkForAll genFullRankWideRHS triangularMinimumNorm) :

   ("householderReconstruction",
      checkForAll Gen.matrix householderReconstruction) :
   ("householderDeterminant",
      checkForAll Gen.square householderDeterminant) :
   ("householderLeastSquares",
      checkForAll genFullRankTallRHS householderLeastSquares) :
   ("householderMinimumNorm",
      checkForAll genFullRankWideRHS householderMinimumNorm) :
   ("householderSolveRR",
      checkForAllExtra
         (liftA2 (,) QC.arbitraryBoundedEnum QC.arbitraryBoundedEnum)
         Gen.fullRankTall householderSolveRR) :
   ("householderMultiplyR",
      checkForAllExtra QC.arbitraryBoundedEnum
         ((,) <$> Gen.tall <#*#> Gen.matrix)
         householderMultiplyR) :
   ("householderQOrthogonal",
      checkForAll Gen.matrix householderQOrthogonal) :
   ("householderMultiplyQ",
      checkForAllExtra
         (liftA2 (,) QC.arbitraryBoundedEnum QC.arbitraryBoundedEnum)
         ((,) <$> Gen.matrix <#\#> Gen.matrix)
         householderMultiplyQ) :
   ("householderTallQOrthogonal",
      checkForAll Gen.tall householderTallQOrthogonal) :
   ("householderTallMultiplyQ",
      checkForAll ((,) <$> Gen.tall <#*#> Gen.matrix) householderTallMultiplyQ) :
   ("householderTallMultiplyQAdjoint",
      checkForAll
         ((,) <$> Gen.tall <#\#> Gen.matrix) householderTallMultiplyQAdjoint) :
   Divide.testsVar (HH.fromMatrix <$> Gen.invertible) ++
   []