{-# 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.Private as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Type.Private 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.Private (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)
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 = ArrMatrix.lift2 Plain.leastSquares
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 = ArrMatrix.lift2 Plain.minimumNorm
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 rcond a b =
mapSnd ArrMatrix.lift0 $
Plain.leastSquaresMinimumNormRCond
rcond (ArrMatrix.toVector a) (ArrMatrix.toVector 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 rcond =
mapSnd (ArrMatrix.lift0 . Basic.recheck) .
Plain.pseudoInverseRCond rcond .
Basic.uncheck . ArrMatrix.toVector
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 b d x =
x
|-|
ArrMatrix.unliftColumn Layout.ColumnMajor
(minimumNorm b) (Multiply.matrixVector b x |-| d)
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 a c b d =
Plain.leastSquaresConstraint
(ArrMatrix.toVector a) c
(ArrMatrix.toVector b) d
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 a b d =
if Shape.size (Matrix.height a) == 0
then (Vector.zero (Matrix.width a), Vector.zero (Matrix.width b))
else Plain.gaussMarkovLinearModel
(ArrMatrix.toVector a) (ArrMatrix.toVector b) d
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 a =
let hh = HH.fromMatrix a
in (HH.extractQ hh, HH.extractR 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 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
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 = Plain.determinantAbsolute . ArrMatrix.toVector
complement ::
(Shape.C height, Shape.C width, Class.Floating a) =>
Tall height width a -> Tall height ShapeInt a
complement = ArrMatrix.lift1 Plain.complement
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 a b =
let qr = HH.fromMatrix $ ArrMatrix.lift1 Basic.adjoint a
in (ArrMatrix.lift0 $ Plain.extractComplement qr,
ArrMatrix.unliftColumn Layout.ColumnMajor (HH.minimumNorm qr) b)
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 c d =
let a = Basic.adjoint $ Plain.complement $ ArrMatrix.toVector c
in (ArrMatrix.lift0 a, Basic.multiplyVector a d)