{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Orthogonal (
   leastSquares,
   minimumNorm,
   leastSquaresMinimumNormRCond,
   pseudoInverseRCond,

   project,
   leastSquaresConstraint,
   gaussMarkovLinearModel,

   determinant,
   determinantAbsolute,
   complement,
   affineFrameFromFiber,
   affineFiberFromFrame,

   householder,
   householderTall,
   ) where

import qualified Numeric.LAPACK.Orthogonal.Householder as HH
import qualified Numeric.LAPACK.Orthogonal.Basic as Plain

import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Layout as Layout
import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Type as Matrix
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Shape as ExtShape
import Numeric.LAPACK.Matrix.Array (Full, General, Tall, Wide, Square)
import Numeric.LAPACK.Matrix.Private (ShapeInt)
import Numeric.LAPACK.Vector (Vector, (|-|))
import Numeric.LAPACK.Scalar (RealOf)

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Shape as Shape

import Data.Tuple.HT (mapSnd)


{- |
If @x = leastSquares a b@
then @x@ minimizes @Vector.norm2 (multiply a x `sub` b)@.

Precondition: @a@ must have full rank and @height a >= width a@.
-}
leastSquares ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full meas horiz Extent.Small height width a ->
   Full meas vert horiz height nrhs a ->
   Full meas vert horiz width nrhs a
leastSquares :: Full meas horiz Small height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
leastSquares = (PlainArray
   Unpacked Arbitrary Filled Filled meas horiz Small height width a
 -> PlainArray
      Unpacked Arbitrary Filled Filled meas vert horiz height nrhs a
 -> PlainArray
      Unpacked Arbitrary Filled Filled meas vert horiz width nrhs a)
-> Full meas horiz Small height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
forall packA propA lowerA upperA measA vertA horizA heightA widthA
       packB propB lowerB upperB measB vertB horizB heightB widthB packC
       propC lowerC upperC measC vertC horizC heightC widthC a b c.
(ToPlain
   packA propA lowerA upperA measA vertA horizA heightA widthA,
 ToPlain
   packB propB lowerB upperB measB vertB horizB heightB widthB,
 FromPlain
   packC propC lowerC upperC measC vertC horizC heightC widthC) =>
(PlainArray
   packA propA lowerA upperA measA vertA horizA heightA widthA a
 -> PlainArray
      packB propB lowerB upperB measB vertB horizB heightB widthB b
 -> PlainArray
      packC propC lowerC upperC measC vertC horizC heightC widthC c)
-> ArrayMatrix
     packA propA lowerA upperA measA vertA horizA heightA widthA a
-> ArrayMatrix
     packB propB lowerB upperB measB vertB horizB heightB widthB b
-> ArrayMatrix
     packC propC lowerC upperC measC vertC horizC heightC widthC c
ArrMatrix.lift2 PlainArray
  Unpacked Arbitrary Filled Filled meas horiz Small height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled meas vert horiz height nrhs a
-> PlainArray
     Unpacked Arbitrary Filled Filled meas vert horiz width nrhs a
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
Full meas horiz Small height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
Plain.leastSquares

{- |
The vector @x@ with @x = minimumNorm a b@
is the vector with minimal @Vector.norm2 x@
that satisfies @multiply a x == b@.

Precondition: @a@ must have full rank and @height a <= width a@.
-}
minimumNorm ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full meas Extent.Small vert height width a ->
   Full meas vert horiz height nrhs a ->
   Full meas vert horiz width nrhs a
minimumNorm :: Full meas Small vert height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
minimumNorm = (PlainArray
   Unpacked Arbitrary Filled Filled meas Small vert height width a
 -> PlainArray
      Unpacked Arbitrary Filled Filled meas vert horiz height nrhs a
 -> PlainArray
      Unpacked Arbitrary Filled Filled meas vert horiz width nrhs a)
-> Full meas Small vert height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
forall packA propA lowerA upperA measA vertA horizA heightA widthA
       packB propB lowerB upperB measB vertB horizB heightB widthB packC
       propC lowerC upperC measC vertC horizC heightC widthC a b c.
(ToPlain
   packA propA lowerA upperA measA vertA horizA heightA widthA,
 ToPlain
   packB propB lowerB upperB measB vertB horizB heightB widthB,
 FromPlain
   packC propC lowerC upperC measC vertC horizC heightC widthC) =>
(PlainArray
   packA propA lowerA upperA measA vertA horizA heightA widthA a
 -> PlainArray
      packB propB lowerB upperB measB vertB horizB heightB widthB b
 -> PlainArray
      packC propC lowerC upperC measC vertC horizC heightC widthC c)
-> ArrayMatrix
     packA propA lowerA upperA measA vertA horizA heightA widthA a
-> ArrayMatrix
     packB propB lowerB upperB measB vertB horizB heightB widthB b
-> ArrayMatrix
     packC propC lowerC upperC measC vertC horizC heightC widthC c
ArrMatrix.lift2 PlainArray
  Unpacked Arbitrary Filled Filled meas Small vert height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled meas vert horiz height nrhs a
-> PlainArray
     Unpacked Arbitrary Filled Filled meas vert horiz width nrhs a
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
Full meas Small vert height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
Plain.minimumNorm


{- |
If @(rank,x) = leastSquaresMinimumNormRCond rcond a b@
then @x@ is the vector with minimum @Vector.norm2 x@
that minimizes @Vector.norm2 (a #*| x `sub` b)@.

Matrix @a@ can have any rank
but you must specify the reciprocal condition of the rank-truncated matrix.
-}
leastSquaresMinimumNormRCond ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   RealOf a ->
   Full meas horiz vert height width a ->
   Full meas vert horiz height nrhs a ->
   (Int, Full meas vert horiz width nrhs a)
leastSquaresMinimumNormRCond :: RealOf a
-> Full meas horiz vert height width a
-> Full meas vert horiz height nrhs a
-> (Int, Full meas vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond Full meas horiz vert height width a
a Full meas vert horiz height nrhs a
b =
   (Full meas vert horiz width nrhs a
 -> Full meas vert horiz width nrhs a)
-> (Int, Full meas vert horiz width nrhs a)
-> (Int, Full meas vert horiz width nrhs a)
forall b c a. (b -> c) -> (a, b) -> (a, c)
mapSnd Full meas vert horiz width nrhs a
-> Full meas vert horiz width nrhs a
forall pack prop lower upper meas vert horiz height width a.
FromPlain pack prop lower upper meas vert horiz height width =>
PlainArray pack prop lower upper meas vert horiz height width a
-> ArrayMatrix pack prop lower upper meas vert horiz height width a
ArrMatrix.lift0 ((Int, Full meas vert horiz width nrhs a)
 -> (Int, Full meas vert horiz width nrhs a))
-> (Int, Full meas vert horiz width nrhs a)
-> (Int, Full meas vert horiz width nrhs a)
forall a b. (a -> b) -> a -> b
$
   RealOf a
-> Full meas horiz vert height width a
-> Full meas vert horiz height nrhs a
-> (Int, Full meas vert horiz width nrhs a)
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
RealOf a
-> Full meas horiz vert height width a
-> Full meas vert horiz height nrhs a
-> (Int, Full meas vert horiz width nrhs a)
Plain.leastSquaresMinimumNormRCond
      RealOf a
rcond (Full meas horiz vert height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled meas horiz vert height width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector Full meas horiz vert height width a
a) (Full meas vert horiz height nrhs a
-> PlainArray
     Unpacked Arbitrary Filled Filled meas vert horiz height nrhs a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector Full meas vert horiz height nrhs a
b)


pseudoInverseRCond ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Class.Floating a) =>
   RealOf a ->
   Full meas vert horiz height width a ->
   (Int, Full meas horiz vert width height a)
pseudoInverseRCond :: RealOf a
-> Full meas vert horiz height width a
-> (Int, Full meas horiz vert width height a)
pseudoInverseRCond RealOf a
rcond =
   (Full meas horiz vert (Unchecked width) (Unchecked height) a
 -> Full meas horiz vert width height a)
-> (Int,
    Full meas horiz vert (Unchecked width) (Unchecked height) a)
-> (Int, Full meas horiz vert width height a)
forall b c a. (b -> c) -> (a, b) -> (a, c)
mapSnd (Array (Full meas horiz vert width height) a
-> Full meas horiz vert width height a
forall pack prop lower upper meas vert horiz height width a.
FromPlain pack prop lower upper meas vert horiz height width =>
PlainArray pack prop lower upper meas vert horiz height width a
-> ArrayMatrix pack prop lower upper meas vert horiz height width a
ArrMatrix.lift0 (Array (Full meas horiz vert width height) a
 -> Full meas horiz vert width height a)
-> (Full meas horiz vert (Unchecked width) (Unchecked height) a
    -> Array (Full meas horiz vert width height) a)
-> Full meas horiz vert (Unchecked width) (Unchecked height) a
-> Full meas horiz vert width height a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full meas horiz vert (Unchecked width) (Unchecked height) a
-> Array (Full meas horiz vert width height) a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz (Unchecked height) (Unchecked width) a
-> Full meas vert horiz height width a
Basic.recheck) ((Int, Full meas horiz vert (Unchecked width) (Unchecked height) a)
 -> (Int, Full meas horiz vert width height a))
-> (Full meas vert horiz height width a
    -> (Int,
        Full meas horiz vert (Unchecked width) (Unchecked height) a))
-> Full meas vert horiz height width a
-> (Int, Full meas horiz vert width height a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   RealOf a
-> Full meas vert horiz (Unchecked height) (Unchecked width) a
-> (Int,
    Full meas horiz vert (Unchecked width) (Unchecked height) a)
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 Eq width, Floating a) =>
RealOf a
-> Full meas vert horiz height width a
-> (Int, Full meas horiz vert width height a)
Plain.pseudoInverseRCond RealOf a
rcond (Full meas vert horiz (Unchecked height) (Unchecked width) a
 -> (Int,
     Full meas horiz vert (Unchecked width) (Unchecked height) a))
-> (Full meas vert horiz height width a
    -> Full meas vert horiz (Unchecked height) (Unchecked width) a)
-> Full meas vert horiz height width a
-> (Int,
    Full meas horiz vert (Unchecked width) (Unchecked height) a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Full meas vert horiz height width a
-> Full meas vert horiz (Unchecked height) (Unchecked width) a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width a
-> Full meas vert horiz (Unchecked height) (Unchecked width) a
Basic.uncheck (Full meas vert horiz height width a
 -> Full meas vert horiz (Unchecked height) (Unchecked width) a)
-> (Full meas vert horiz height width a
    -> Full meas vert horiz height width a)
-> Full meas vert horiz height width a
-> Full meas vert horiz (Unchecked height) (Unchecked width) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full meas vert horiz height width a
-> Full meas vert horiz height width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector


{- |
@project b d x@ projects @x@ on the plane described by @B*x = d@.

@b@ must have full rank.
-}
project ::
   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   Wide height width a -> Vector height a ->
   Vector width a -> Vector width a
project :: Wide height width a
-> Vector height a -> Vector width a -> Vector width a
project Wide height width a
b Vector height a
d Vector width a
x =
   Vector width a
x
   Vector width a -> Vector width a -> Vector width a
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
|-|
   Order
-> (General height () a -> General width () a)
-> Vector height a
-> Vector width a
forall height0 a height1 b.
Order
-> (General height0 () a -> General height1 () b)
-> Vector height0 a
-> Vector height1 b
ArrMatrix.unliftColumn Order
Layout.ColumnMajor
      (Wide height width a -> General height () a -> General width () a
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
Full meas Small vert height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
minimumNorm Wide height width a
b) (Wide height width a -> Vector width a -> Vector height a
forall typ xl xu lower upper meas vert horiz height width a.
(MultiplyVector typ xl xu, Strip lower, Strip upper, Measure meas,
 C vert, C horiz, C height, C width, Eq width, Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Vector width a -> Vector height a
Multiply.matrixVector Wide height width a
b Vector width a
x Vector height a -> Vector height a -> Vector height a
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
|-| Vector height a
d)


{- |
@leastSquaresConstraint a c b d@ computes @x@
with minimal @|| c - A*x ||_2@ and constraint @B*x = d@.

@b@ must be wide and @a===b@ must be tall
and both matrices must have full rank.
-}
leastSquaresConstraint ::
   (Shape.C height, Eq height,
    Shape.C width, Eq width,
    Shape.C constraints, Eq constraints, Class.Floating a) =>
   General height width a -> Vector height a ->
   Wide constraints width a -> Vector constraints a ->
   Vector width a
leastSquaresConstraint :: General height width a
-> Vector height a
-> Wide constraints width a
-> Vector constraints a
-> Vector width a
leastSquaresConstraint General height width a
a Vector height a
c Wide constraints width a
b Vector constraints a
d =
   General height width a
-> Vector height a
-> Wide constraints width a
-> Vector constraints a
-> Vector width a
forall height width constraints a.
(C height, Eq height, C width, Eq width, C constraints,
 Eq constraints, Floating a) =>
General height width a
-> Vector height a
-> Wide constraints width a
-> Vector constraints a
-> Vector width a
Plain.leastSquaresConstraint
      (General height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Big height width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector General height width a
a) Vector height a
c
      (Wide constraints width a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Small Big constraints width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector Wide constraints width a
b) Vector constraints a
d

{- |
@gaussMarkovLinearModel a b d@ computes @(x,y)@
with minimal @|| y ||_2@ and constraint @d = A*x + B*y@.

@a@ must be tall and @a|||b@ must be wide
and both matrices must have full rank.
-}
gaussMarkovLinearModel ::
   (Shape.C height, Eq height,
    Shape.C width, Eq width,
    Shape.C opt, Eq opt, Class.Floating a) =>
   Tall height width a -> General height opt a -> Vector height a ->
   (Vector width a, Vector opt a)
gaussMarkovLinearModel :: Tall height width a
-> General height opt a
-> Vector height a
-> (Vector width a, Vector opt a)
gaussMarkovLinearModel Tall height width a
a General height opt a
b Vector height a
d =
   {-
   Fortran-LAPACK and OpenBLAS would leave Y uninitalized
   instead of setting Y to zeros.
   -}
   if height -> Int
forall sh. C sh => sh -> Int
Shape.size (Tall height width a -> height
forall typ meas vert horiz xl xu lower upper height width a.
(Box typ, Measure meas, C vert, C horiz) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> height
Matrix.height Tall height width a
a) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
      then (width -> Vector width a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero (Tall height width a -> width
forall typ meas vert horiz xl xu lower upper height width a.
(Box typ, Measure meas, C vert, C horiz) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> width
Matrix.width Tall height width a
a), opt -> Vector opt a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero (General height opt a -> opt
forall typ meas vert horiz xl xu lower upper height width a.
(Box typ, Measure meas, C vert, C horiz) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> width
Matrix.width General height opt a
b))
      else Tall height width a
-> General height opt a
-> Vector height a
-> (Vector width a, Vector opt a)
forall height width opt a.
(C height, Eq height, C width, Eq width, C opt, Eq opt,
 Floating a) =>
Tall height width a
-> General height opt a
-> Vector height a
-> (Vector width a, Vector opt a)
Plain.gaussMarkovLinearModel
               (Tall height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Small height width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector Tall height width a
a) (General height opt a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Big height opt a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector General height opt a
b) Vector height a
d


{-
@(q,r) = householder a@
means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.
-}
householder ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    ExtShape.Permutable height, Shape.C width, Class.Floating a) =>
   Full meas vert horiz height width a ->
   (Square height a, Unpacked.UpperTrapezoid meas vert horiz height width a)
householder :: Full meas vert horiz height width a
-> (Square height a, UpperTrapezoid meas vert horiz height width a)
householder Full meas vert horiz height width a
a =
   let hh :: Householder meas vert horiz height width a
hh = Full meas vert horiz height width a
-> Householder meas vert horiz height width a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Householder meas vert horiz height width a
HH.fromMatrix Full meas vert horiz height width a
a
   in  (Householder meas vert horiz height width a -> Square height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, Permutable height, C width,
 Floating a) =>
Householder meas vert horiz height width a -> Square height a
HH.extractQ Householder meas vert horiz height width a
hh, Householder meas vert horiz height width a
-> UpperTrapezoid meas vert horiz height width a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, Permutable height, C width,
 Floating a) =>
Householder meas vert horiz height width a
-> UpperTrapezoid meas vert horiz height width a
HH.extractR Householder meas vert horiz height width a
hh)

householderTall ::
   (Extent.Measure meas, Extent.C vert,
    Shape.C height, ExtShape.Permutable width, Class.Floating a) =>
   Full meas vert Extent.Small height width a ->
   (Full meas vert Extent.Small height width a, Triangular.Upper width a)
householderTall :: Full meas vert Small height width a
-> (Full meas vert Small height width a, Upper width a)
householderTall Full meas vert Small height width a
a =
   let hh :: Householder meas vert Small height width a
hh = Full meas vert Small height width a
-> Householder meas vert Small height width a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Householder meas vert horiz height width a
HH.fromMatrix Full meas vert Small height width a
a
   in  (Householder meas vert Small height width a
-> Full meas vert Small height width a
forall meas vert height width a.
(Measure meas, C vert, C height, Permutable width, Floating a) =>
Householder meas vert Small height width a
-> Full meas vert Small height width a
HH.tallExtractQ Householder meas vert Small height width a
hh, Householder meas vert Small height width a -> Upper width a
forall meas vert height width a.
(Measure meas, C vert, C height, Permutable width, Floating a) =>
Householder meas vert Small height width a -> Upper width a
HH.tallExtractR Householder meas vert Small height width a
hh)


determinant :: (Shape.C sh, Class.Floating a) => Square sh a -> a
determinant :: Square sh a -> a
determinant = HouseholderFlex Filled Filled Shape Small Small sh sh a -> a
forall sh a lower upper.
(C sh, Floating a) =>
HouseholderFlex lower upper Shape Small Small sh sh a -> a
HH.determinant (HouseholderFlex Filled Filled Shape Small Small sh sh a -> a)
-> (Square sh a
    -> HouseholderFlex Filled Filled Shape Small Small sh sh a)
-> Square sh a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Square sh a
-> HouseholderFlex Filled Filled Shape Small Small sh sh a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Householder meas vert horiz height width a
HH.fromMatrix

{-|
Gramian determinant -
works also for non-square matrices, but is sensitive to transposition.

> determinantAbsolute a = sqrt (Herm.determinant (Herm.gramian a))
-}
determinantAbsolute ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Class.Floating a) =>
   Full meas vert horiz height width a -> RealOf a
determinantAbsolute :: Full meas vert horiz height width a -> RealOf a
determinantAbsolute = Full meas vert horiz height width a -> RealOf a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a -> RealOf a
Plain.determinantAbsolute (Full meas vert horiz height width a -> RealOf a)
-> (Full meas vert horiz height width a
    -> Full meas vert horiz height width a)
-> Full meas vert horiz height width a
-> RealOf a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full meas vert horiz height width a
-> Full meas vert horiz height width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector


{- |
For an m-by-n-matrix @a@ with m>=n
this function computes an m-by-(m-n)-matrix @b@
such that @Matrix.multiply (adjoint b) a@ is a zero matrix.
The function does not try to compensate a rank deficiency of @a@.
That is, @a|||b@ has full rank if and only if @a@ has full rank.

For full-rank matrices you might also call this @kernel@ or @nullspace@.
-}
complement ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   Tall height width a -> Tall height ShapeInt a
complement :: Tall height width a -> Tall height ShapeInt a
complement = (PlainArray
   Unpacked Arbitrary Filled Filled Size Big Small height width a
 -> PlainArray
      Unpacked Arbitrary Filled Filled Size Big Small height ShapeInt a)
-> Tall height width a -> Tall height ShapeInt a
forall packA propA lowerA upperA measA vertA horizA heightA widthA
       packB propB lowerB upperB measB vertB horizB heightB widthB a b.
(ToPlain
   packA propA lowerA upperA measA vertA horizA heightA widthA,
 FromPlain
   packB propB lowerB upperB measB vertB horizB heightB widthB) =>
(PlainArray
   packA propA lowerA upperA measA vertA horizA heightA widthA a
 -> PlainArray
      packB propB lowerB upperB measB vertB horizB heightB widthB b)
-> ArrayMatrix
     packA propA lowerA upperA measA vertA horizA heightA widthA a
-> ArrayMatrix
     packB propB lowerB upperB measB vertB horizB heightB widthB b
ArrMatrix.lift1 PlainArray
  Unpacked Arbitrary Filled Filled Size Big Small height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Small height ShapeInt a
forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
Plain.complement


{- |
> affineFrameFromFiber a b == (c,d)

Means:
An affine subspace is given implicitly by {x : a#*|x == b}.
Convert this into an explicit representation {c#*|y|+|d : y}.
Matrix @a@ must have full rank,
otherwise the explicit representation will miss dimensions
and we cannot easily determine the origin @d@ as a minimum norm solution.

The computation is like

> c = complement $ adjoint a
> d = minimumNorm a b

but the QR decomposition of 'a' is computed only once.
-}
affineFrameFromFiber ::
   (Shape.C width, Eq width, Shape.C height, Eq height, Class.Floating a) =>
   Wide height width a -> Vector height a ->
   (Tall width ShapeInt a, Vector width a)
affineFrameFromFiber :: Wide height width a
-> Vector height a -> (Tall width ShapeInt a, Vector width a)
affineFrameFromFiber Wide height width a
a Vector height a
b =
   let qr :: Householder Size Big Small width height a
qr = Full Size Big Small width height a
-> Householder Size Big Small width height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Householder meas vert horiz height width a
HH.fromMatrix (Full Size Big Small width height a
 -> Householder Size Big Small width height a)
-> Full Size Big Small width height a
-> Householder Size Big Small width height a
forall a b. (a -> b) -> a -> b
$ (PlainArray
   Unpacked Arbitrary Filled Filled Size Small Big height width a
 -> PlainArray
      Unpacked Arbitrary Filled Filled Size Big Small width height a)
-> Wide height width a -> Full Size Big Small width height a
forall packA propA lowerA upperA measA vertA horizA heightA widthA
       packB propB lowerB upperB measB vertB horizB heightB widthB a b.
(ToPlain
   packA propA lowerA upperA measA vertA horizA heightA widthA,
 FromPlain
   packB propB lowerB upperB measB vertB horizB heightB widthB) =>
(PlainArray
   packA propA lowerA upperA measA vertA horizA heightA widthA a
 -> PlainArray
      packB propB lowerB upperB measB vertB horizB heightB widthB b)
-> ArrayMatrix
     packA propA lowerA upperA measA vertA horizA heightA widthA a
-> ArrayMatrix
     packB propB lowerB upperB measB vertB horizB heightB widthB b
ArrMatrix.lift1 PlainArray
  Unpacked Arbitrary Filled Filled Size Small Big height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Small width height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.adjoint Wide height width a
a
   in (PlainArray
  Unpacked Arbitrary Filled Filled Size Big Small width ShapeInt a
-> Tall width ShapeInt a
forall pack prop lower upper meas vert horiz height width a.
FromPlain pack prop lower upper meas vert horiz height width =>
PlainArray pack prop lower upper meas vert horiz height width a
-> ArrayMatrix pack prop lower upper meas vert horiz height width a
ArrMatrix.lift0 (PlainArray
   Unpacked Arbitrary Filled Filled Size Big Small width ShapeInt a
 -> Tall width ShapeInt a)
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Small width ShapeInt a
-> Tall width ShapeInt a
forall a b. (a -> b) -> a -> b
$ Householder Size Big Small width height a -> Tall width ShapeInt a
forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
Plain.extractComplement Householder Size Big Small width height a
qr,
       Order
-> (General height () a -> General width () a)
-> Vector height a
-> Vector width a
forall height0 a height1 b.
Order
-> (General height0 () a -> General height1 () b)
-> Vector height0 a
-> Vector height1 b
ArrMatrix.unliftColumn Order
Layout.ColumnMajor (Householder Size Big Small width height a
-> General height () a -> General width () a
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 Eq width, C nrhs, Floating a) =>
Householder meas vert Small width height a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
HH.minimumNorm Householder Size Big Small width height a
qr) Vector height a
b)

{- |
This conversion is somehow inverse to 'affineFrameFromFiber'.
However, it is not precisely inverse in either direction.
This is because both 'affineFrameFromFiber' and 'affineFiberFromFrame'
accept non-orthogonal matrices but always return orthogonal ones.

In @affineFiberFromFrame c d@,
matrix @c@ should have full rank,
otherwise the implicit representation will miss dimensions.
-}
affineFiberFromFrame ::
   (Shape.C width, Eq width, Shape.C height, Eq height, Class.Floating a) =>
   Tall height width a -> Vector height a ->
   (Wide ShapeInt height a, Vector ShapeInt a)
affineFiberFromFrame :: Tall height width a
-> Vector height a -> (Wide ShapeInt height a, Vector ShapeInt a)
affineFiberFromFrame Tall height width a
c Vector height a
d =
   let a :: Full Size Small Big ShapeInt height a
a = Full Size Big Small height ShapeInt a
-> Full Size Small Big ShapeInt height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.adjoint (Full Size Big Small height ShapeInt a
 -> Full Size Small Big ShapeInt height a)
-> Full Size Big Small height ShapeInt a
-> Full Size Small Big ShapeInt height a
forall a b. (a -> b) -> a -> b
$ Tall height width a -> Full Size Big Small height ShapeInt a
forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
Plain.complement (Tall height width a -> Full Size Big Small height ShapeInt a)
-> Tall height width a -> Full Size Big Small height ShapeInt a
forall a b. (a -> b) -> a -> b
$ Tall height width a
-> PlainArray
     Unpacked Arbitrary Filled Filled Size Big Small height width a
forall pack property lower upper meas vert horiz height width a.
ToPlain pack property lower upper meas vert horiz height width =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> PlainArray
     pack property lower upper meas vert horiz height width a
ArrMatrix.toVector Tall height width a
c
   in (PlainArray
  Unpacked Arbitrary Filled Filled Size Small Big ShapeInt height a
-> Wide ShapeInt height a
forall pack prop lower upper meas vert horiz height width a.
FromPlain pack prop lower upper meas vert horiz height width =>
PlainArray pack prop lower upper meas vert horiz height width a
-> ArrayMatrix pack prop lower upper meas vert horiz height width a
ArrMatrix.lift0 Full Size Small Big ShapeInt height a
PlainArray
  Unpacked Arbitrary Filled Filled Size Small Big ShapeInt height a
a, Full Size Small Big ShapeInt height a
-> Vector height a -> Vector ShapeInt a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Eq width,
 Floating a) =>
Full meas vert horiz height width a
-> Vector width a -> Vector height a
Basic.multiplyVector Full Size Small Big ShapeInt height a
a Vector height a
d)