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

   determinant,
   determinantAbsolute,
   complement,

   householder,
   householderTall,
   ) where

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

import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import Numeric.LAPACK.Matrix.Array (Full, Tall, Square)
import Numeric.LAPACK.Matrix.Private (ShapeInt)
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.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full horiz Extent.Small height width a ->
   Full vert horiz height nrhs a ->
   Full vert horiz width nrhs a
leastSquares = ArrMatrix.lift2 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.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full Extent.Small vert height width a ->
   Full vert horiz height nrhs a ->
   Full vert horiz width nrhs a
minimumNorm = ArrMatrix.lift2 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.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   RealOf a ->
   Full horiz vert height width a ->
   Full vert horiz height nrhs a ->
   (Int, Full vert horiz width nrhs a)
leastSquaresMinimumNormRCond rcond a b =
   mapSnd ArrMatrix.lift0 $
   Plain.leastSquaresMinimumNormRCond
      rcond (ArrMatrix.toVector a) (ArrMatrix.toVector b)


pseudoInverseRCond ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   RealOf a ->
   Full vert horiz height width a ->
   (Int, Full horiz vert width height a)
pseudoInverseRCond rcond =
   mapSnd (ArrMatrix.lift0 . Basic.recheck) .
   Plain.pseudoInverseRCond rcond .
   Basic.uncheck . ArrMatrix.toVector


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

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


determinant :: (Shape.C sh, Class.Floating a) => Square sh a -> a
determinant = HH.determinant . 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.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   Full vert horiz height width a -> RealOf a
determinantAbsolute = Plain.determinantAbsolute . 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 = ArrMatrix.lift1 Plain.complement