{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.Function (
   SqRt, sqrt, sqrtSchur, sqrtDenmanBeavers,
   Exp, exp, expRealHermitian,
   Log, log, logUnipotentUpper,
   LiftReal, liftReal,
   ) where

import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Matrix.Lazy.UpperTriangular as LazyUpper
import qualified Numeric.LAPACK.Matrix.Array.Mosaic as ArrMosaic
import qualified Numeric.LAPACK.Matrix.Mosaic.Basic as Mosaic
import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian
import qualified Numeric.LAPACK.Matrix.Symmetric as Symmetric
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Quadratic as Quad
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Diagonal as Diagonal
import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Scalar as Scalar
import qualified Numeric.LAPACK.Shape.Private as ShapePriv
import qualified Numeric.LAPACK.Shape as ExtShape
import Numeric.LAPACK.Matrix.Array.Mosaic
         (UnitUpperP, LowerP, UpperP, FlexUpperP,
          HermitianP, HermitianPosSemidefP)
import Numeric.LAPACK.Matrix.Array.Private (Quadratic, ArrayMatrix, packTag)
import Numeric.LAPACK.Matrix.Square (Square)
import Numeric.LAPACK.Matrix
         ((#!), (#+#), (#-#), (#\##), (#*\), (#*##), (##*#))
import Numeric.LAPACK.Vector ((|+|), (|-|))

import qualified Numeric.Netlib.Class as Class

import qualified Type.Data.Bool as Bool
import Type.Data.Bool (False, True)

import Foreign.Storable (Storable)

import qualified Data.Array.Comfort.Storable.Unchecked as Array
import qualified Data.Array.Comfort.Shape as Shape

import qualified Data.NonEmpty as NonEmpty
import qualified Data.Complex as Complex
import qualified Data.List.HT as ListHT
import qualified Data.Stream as Stream
import Data.Complex (Complex)
import Data.Semigroup ((<>))
import Data.Function.HT (nest)
import Data.Tuple.HT (mapFst, swap)
import Data.Stream (Stream)

import qualified Prelude as P
import Prelude hiding (sqrt, exp, log)


class (MatrixShape.Property property) => SqRt property where
   sqrt ::
      (Layout.Packing pack,
       MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
       Shape.C sh, Class.Real a) =>
      Quadratic pack property lower upper sh a ->
      Quadratic pack property lower upper sh a

instance SqRt Omni.Unit where
   sqrt a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerIdentity -> a
         Omni.PowerUpperTriangular -> sqrtUpper (const Scalar.one) a
         Omni.PowerLowerTriangular ->
            Triangular.transpose .
            sqrtUpper (const Scalar.one) .
            Triangular.transpose $ a
         Omni.PowerDiagonal -> error "non-unit diagonal impossible"
         Omni.PowerFull -> error "unit full matrix impossible"

{- |
For Full matrices:
Explicit solution for matrices up to size 2.
Solution via 'sqrtDenmanBeavers' for larger sizes.
-}
instance SqRt Omni.Arbitrary where
   sqrt m =
      case Omni.powerSingleton $ ArrMatrix.shape m of
         Omni.PowerDiagonal -> liftDiagonal P.sqrt m
         Omni.PowerUpperTriangular -> sqrtUpper P.sqrt m
         Omni.PowerLowerTriangular ->
            Triangular.transpose . sqrtUpper P.sqrt . Triangular.transpose $ m
         Omni.PowerFull ->
            if Shape.size (Quad.size m) <= 2
               then
                  ArrMatrix.Array $ Array.fromList (ArrMatrix.shape m) $
                  case Array.toList $ ArrMatrix.unwrap m of
                     [qA,qB,qC,qD] ->
                        let ((a,b),(c,d)) = sqrt2 ((qA,qB),(qC,qD))
                        in [a,b,c,d]
                     [qA] -> [P.sqrt qA]
                     _ -> []
               else
                  case Scalar.precisionOfFunctor m of
                     Scalar.Float -> sqrtDenmanBeavers m
                     Scalar.Double -> sqrtDenmanBeavers m

instance SqRt Omni.Symmetric where
   sqrt a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerSymmetric ->
            Symmetric.fromHermitian . sqrtHermitian . Hermitian.fromSymmetric
               $ a
         _ -> error "Symmetric.sqrt: impossible shape"

instance
   (neg ~ False, Bool.C zero, Bool.C pos) =>
      SqRt (Omni.Hermitian neg zero pos) where
   sqrt a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerDiagonal -> liftHermitianDiagonal P.sqrt a
         Omni.PowerHermitian -> sqrtHermitian a
         _ -> error "Hermitian.sqrt: impossible shape"

sqrtHermitian ::
   (Layout.Packing pack,
    Bool.C neg, Bool.C zero, Bool.C pos, Omni.Hermitian neg zero pos ~ herm,
    Shape.C sh, Class.Real a) =>
   Quadratic pack herm Layout.Filled Layout.Filled sh a ->
   Quadratic pack herm Layout.Filled Layout.Filled sh a
sqrtHermitian a =
   case Scalar.precisionOfFunctor a of
      Scalar.Float -> liftHermitian P.sqrt a
      Scalar.Double -> liftHermitian P.sqrt a

sqrtUpper ::
   (Layout.Packing pack, Omni.TriDiag diag, Shape.C sh, Class.Floating a) =>
    (a -> a) -> FlexUpperP pack diag sh a -> FlexUpperP pack diag sh a
sqrtUpper sqrtF a =
   Triangular.adaptOrder a $ ArrMatrix.lift1 Mosaic.reunpack $
   LazyUpper.toStorable $ LazyUpper.sqrt sqrtF $ LazyUpper.fromStorable a

{-
/A B\ = /a b\. /a b\ = /a*a+b*c a*b+b*d\ = /a²+b*c  b*(a+d)\
\C D/   \c d/  \c d/   \c*a+d*c c*b+d*d/   \c*(a+d) d²+b*c /

b/c = B/C   B*c = C*b
b*C/B = c
A=a²+b²*C/B
D=d²+b²*C/B

a²-d² = A-D
B = b*(a+d)
B*(a-d) = b*(a²-d²) = b*(A-D)
C*(a-d) = c*(a²-d²) = c*(A-D)

(4*B*C-(A+D)^2)*b^4 + 2*B^2*(A-D)*b^2 - B^4 = 0

with x = B/b = C/c = a+d:
maxima:
   poly_buchberger([a²+b*c-A,a*b+b*d-B,c*a+d*c-C,b*c+d²-D,x-(a+d)],[a,d,b,c,x]);

x^4-2*(A+D)*x^2+(A-D)^2+4*B*C = 0
c = C/x, b = B/x
a² = A - b*c = A - B*C/x²
-}
sqrt2 :: (Floating a) => ((a,a),(a,a)) -> ((a,a),(a,a))
sqrt2 ((a,b),(c,d)) =
   let x = P.sqrt $ a+d + 2 * P.sqrt (a*d - b*c)
       y = (d-a)/x
   in (((x-y)/2, b/x), (c/x, (x+y)/2))


{- |
Square root solver that works on the Schur decomposition.

Schur decomposition enables computing the square root
of (some) singular matrices like @((1,0),(0,0))@.
However, the Schur decomposition might emit small negative values
on the diagonal, where exact computation would yield zeros.
This would let the square root solver fail.
And there are singular matrices that have no square root, at all,
e.g. @((0,1),(0,0))@.

The solver is restricted to a real triangular Schur matrix.
The check for non-real eigenvalues may exclude matrices
that actually have a real-valued square root.
E.g. @sqrt ((0,-2),(2,0)) = ((1,-1),(1,-1))@

In the future we might fix this
by solving 2x2 blocks at the diagonal using 'sqrt2'.
-}
sqrtSchur :: (Shape.C sh, Class.Real a) => Square sh a -> Square sh a
sqrtSchur = applyUnchecked $ applyPermutable $ \a ->
   let (q,u) = Square.schur $ checkZeroOffDiagonal a
   in q <> sqrtUpper P.sqrt (Triangular.takeUpper u) #*## Square.adjoint q

checkZeroOffDiagonal ::
   (Shape.Indexed sh, Class.Floating a) =>
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
checkZeroOffDiagonal a =
   if all (\ij -> Scalar.isZero (a#!ij)) $
      ListHT.mapAdjacent (flip (,)) $ Shape.indices $ Quad.size a
      then a
      else error "sqrtSchur: non-real eigenvalues"

{- |
Iterative square root solver, similar to Newton iteration.

Eigenvalues must all be positive,
otherwise, the iteration might loop forever,
or if an eigenvalue is zero, the computation of matrix inverse will fail.
-}
sqrtDenmanBeavers ::
   (ArrMatrix.Homogeneous prop, ArrMatrix.Additive prop) =>
   (Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
   (Shape.C sh, Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar) =>
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
sqrtDenmanBeavers = applyUnchecked $ \a ->
   limit (Scalar.selectReal 1e-6 1e-14) $ fmap fst $
   Stream.iterate
      (\(b,c) ->
         (Matrix.scaleReal 0.5 (b #+# Matrix.inverse c),
          Matrix.scaleReal 0.5 (Matrix.inverse b #+# c)))
      (a, Matrix.identityFrom a)

limit ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width,
    Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar,
    ArrMatrix.ArrayMatrix
      pack prop lower upper meas vert horiz height width ~ m) =>
   ar -> Stream (m a) -> m a
limit eps =
   ArrMatrix.Array . snd . Stream.head .
   Stream.dropWhile
      (\(b0,b1) ->
         Vector.normInf (b0|-|b1) > 0.5*eps * Vector.normInf (b0|+|b1)) .
   streamMapAdjacent (,) . fmap ArrMatrix.unwrap

streamMapAdjacent :: (a -> a -> b) -> Stream a -> Stream b
streamMapAdjacent f xs = Stream.zipWith f xs (Stream.tail xs)


class (MatrixShape.Property property) => Exp property where
   exp ::
      (Layout.Packing pack,
       MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
       Shape.C sh, Class.Floating a) =>
      Quadratic pack property lower upper sh a ->
      Quadratic pack property lower upper sh a

instance Exp Omni.Arbitrary where
   exp a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerDiagonal ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real -> liftDiagonal P.exp a
               Scalar.Complex -> liftDiagonal P.exp a
         Omni.PowerFull ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real ->
                  ArrMatrix.liftUnpacked1 id $
                  applyUnchecked (expScaledPade (#\##)) $ Matrix.toFull a
               Scalar.Complex ->
                  ArrMatrix.liftUnpacked1 id $
                  applyUnchecked (expScaledPade (#\##)) $ Matrix.toFull a
         Omni.PowerUpperTriangular ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real    -> applyUnchecked (expScaledPade solveUpper) a
               Scalar.Complex -> applyUnchecked (expScaledPade solveUpper) a
         Omni.PowerLowerTriangular ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real    -> applyUnchecked (expScaledPade solveLower) a
               Scalar.Complex -> applyUnchecked (expScaledPade solveLower) a

solveUpper ::
   (Layout.Packing pack, Shape.C sh, Eq sh, Class.Floating a) =>
   UpperP pack sh a -> UpperP pack sh a -> UpperP pack sh a
solveUpper d n =
   ArrMosaic.assureMirrored $ d #\## Triangular.toSquare n

solveLower ::
   (Layout.Packing pack, Shape.C sh, Eq sh, Class.Floating a) =>
   LowerP pack sh a -> LowerP pack sh a -> LowerP pack sh a
solveLower d n =
   ArrMosaic.assureMirrored $ d #\## Triangular.toSquare n


instance
   (neg ~ True, zero ~ True, pos ~ True) =>
      Exp (Omni.Hermitian neg zero pos) where
   exp a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerDiagonal ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real -> liftHermitianDiagonal P.exp a
               Scalar.Complex -> liftHermitianDiagonal P.exp a
         Omni.PowerHermitian ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real -> liftHermitian P.exp a
               Scalar.Complex -> liftHermitian P.exp a
         _ -> error "Hermitian.exp: impossible shape"

instance Exp Omni.Symmetric where
   exp a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerSymmetric ->
            case Scalar.complexSingletonOfFunctor a of
               Scalar.Real ->
                  Symmetric.fromHermitian . liftHermitian P.exp .
                  Hermitian.fromSymmetric $ a
               Scalar.Complex ->
                  applyUnchecked
                     (expScaledPade
                        (\d n ->
                           Symmetric.assureSymmetry $
                              d #\## Symmetric.toSquare n))
                     a
         _ -> error "Symmetric.exp: impossible shape"

{-
"Nineteen Dubious Ways to Compute the Exponential of a Matrix,
Twenty-Five Years Later."
by Cleve Moler and Charles Van Loan
11. Matrix Exponential in MATLAB.
-}
expScaledPade ::
   (ArrMatrix.Homogeneous prop, ArrMatrix.Subtractive prop) =>
   (Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
   (Shape.C sh, Eq sh,
    Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar) =>
   (Quadratic pack prop lower upper sh a ->
    Quadratic pack prop lower upper sh a ->
    Quadratic pack prop lower upper sh a) ->
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
expScaledPade solve a0 =
   let s = max 0 $ (1+) $ ceiling $ logBase 2 $ norm1 a0
       a = Matrix.scaleReal ((1/2)^s) a0
       (odds,evens) = deinterleave $ NonEmpty.tail expPadeCoefficients
       Stream.Cons eye as = Matrix.powers a
       (oddPowers,evenPowers) = deinterleave $ Stream.toList as
       v = foldl (#+#) eye $ zipWith Matrix.scaleReal evens evenPowers
       u =
         fmap (NonEmpty.foldl1 (#+#)) $ NonEmpty.fetch $
         zipWith Matrix.scaleReal odds oddPowers
       op vm f um = maybe vm (f vm) um
   in nest s Matrix.square $ solve (op v (#-#) u) (op v (#+#) u)

{-
Move to utility-ht.
Cf. storable-vector, synthesizer-core.
However, they differ in details.
-}
deinterleave :: [a] -> ([a],[a])
deinterleave =
   let go [] = ([],[])
       go (x:xs) = mapFst (x:) $ swap $ go xs
   in go

expPadeCoefficients :: (Class.Real a) => NonEmpty.T [] a
expPadeCoefficients =
   let eps = Scalar.selectReal 1e-8 1e-16
       q = expPadeOrder eps
       coeff k = fromIntegral (q-k+1) / fromIntegral (k*(2*q-k+1))
   in NonEmpty.scanl (*) (1 `asTypeOf` eps) $ map coeff [1..q]

expPadeOrder :: (Class.Real a) => a -> Int
expPadeOrder eps =
   let factorials = scanl (*) 1 $ iterate (1+) 1
   in subtract 1 $ length $ takeWhile id $
      zipWith3 (\num den twoPower -> num > den*twoPower*eps)
         factorials (ListHT.sieve 2 factorials) (iterate (2*) 1)


{- |
Mathematically the name 'expRealSymmetric' would be more common,
but we support definiteness tags only for the 'Hermitian' type.
Formally the result is always positive definite,
but negative eigenvalues easily yield numerically singular matrices as result.
-}
expRealHermitian ::
   (Layout.Packing pack, Shape.C sh, Class.Real a) =>
   HermitianP pack sh a -> HermitianPosSemidefP pack sh a
expRealHermitian =
   Hermitian.relaxSemidefinite .
   Hermitian.assurePositiveDefiniteness .
   exp


class (MatrixShape.Property property) => Log property where
   log ::
      (Layout.Packing pack,
       MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
       Shape.C sh, Class.Real a) =>
      Quadratic pack property lower upper sh a ->
      Quadratic pack property lower upper sh a

instance Log Omni.Arbitrary where
   log = liftReal P.log

instance
   (neg ~ True, zero ~ True, pos ~ True) =>
      Log (Omni.Hermitian neg zero pos) where
   log = liftReal P.log

instance Log Omni.Symmetric where
   log = liftReal P.log


logUnipotentUpper ::
   (Layout.Packing pack, Shape.C sh, Class.Floating a) =>
   UnitUpperP pack sh a -> UpperP pack sh a
logUnipotentUpper = logUnipotent . Triangular.relaxUnitDiagonal

{- |
Logarithm for unipotent matrices.
It is an unchecked error, if the matrix is not unipotent.
-}
logUnipotent ::
   (Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
   (ArrMatrix.Scale prop, ArrMatrix.Subtractive prop) =>
   (Shape.C sh, Class.Floating a) =>
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
logUnipotent a =
   case Scalar.complexSingletonOfFunctor a of
      Scalar.Real    -> logUnipotentAux a
      Scalar.Complex -> logUnipotentAux a

-- cf. numeric-prelude:PowerSeries
logUnipotentAux ::
   (Layout.Packing pack, Omni.PowerStrip lower, Omni.PowerStrip upper) =>
   (ArrMatrix.Scale prop, ArrMatrix.Subtractive prop) =>
   (Shape.C sh, Class.Floating a, Scalar.RealOf a ~ ar, Class.Real ar) =>
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
logUnipotentAux = applyUnchecked $ \a ->
   let b = a #-# Matrix.identityFrom a
   in foldl (#+#) (ArrMatrix.zero $ ArrMatrix.shape b) $
      zipWith ArrMatrix.scale
         (zipWith (/) (cycle [1,-1]) (iterate (1+) 1))
         (Stream.takeWhile ((0<) . normInf) $ Matrix.powers1 b)


class (MatrixShape.Property property) => LiftReal property where
   {- |
   Lift any function with a Taylor expansion
   to a diagonalizable matrix.
   -}
   liftReal ::
      (Layout.Packing pack,
       MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
       Shape.C sh, Class.Real a) =>
      (a -> a) ->
      Quadratic pack property lower upper sh a ->
      Quadratic pack property lower upper sh a

{- |
Generic algorithm that applies a scalar function
to the elements of the diagonal factor
of a full, triangular or diagonal matrix with distinct eigenvalues.
It is not checked whether the matrix has distinct eigenvalues.
-}
instance LiftReal Omni.Arbitrary where
   liftReal f a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerDiagonal -> liftDiagonal f a
         Omni.PowerUpperTriangular -> flip applyUnchecked a $ \b ->
            let (vr,d,vlAdj) = Triangular.eigensystem b
                scal = Triangular.takeDiagonal $ vlAdj <> vr
            in ArrMosaic.assureMirrored $
               Triangular.toSquare vr
                  #*\ Vector.divide (Array.map f d) scal
                  ##*# vlAdj
         Omni.PowerLowerTriangular -> flip applyUnchecked a $ \b ->
            let (vr,d,vlAdj) = Triangular.eigensystem b
                scal = Triangular.takeDiagonal $ vlAdj <> vr
            in ArrMosaic.assureMirrored $
               Triangular.toSquare vr
                  #*\ Vector.divide (Array.map f d) scal
                  ##*# vlAdj
         Omni.PowerFull ->
            case Scalar.precisionOfFunctor a of
               Scalar.Float  -> liftRealFull f a
               Scalar.Double -> liftRealFull f a

liftRealFull ::
   (MatrixShape.Property prop,
    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
    Shape.C sh, Class.Real a, Scalar.RealOf a ~ a) =>
   (a -> a) ->
   Quadratic Layout.Unpacked prop lower upper sh a ->
   Quadratic Layout.Unpacked prop lower upper sh a
liftRealFull f = applyPermutable $ applyUnchecked $ \a ->
   let (vr,d,vlAdj) = Square.eigensystem $ Matrix.toFull a
       vrR = matrixRealPart vr
       vlAdjR = matrixRealPart vlAdj
       dR = Array.map Complex.realPart d
       scal = Square.takeDiagonal $ vlAdjR <> vrR
   in if Scalar.isZero $ Vector.normInf $ Array.map Complex.imagPart d
         then ArrMatrix.liftUnpacked1 id $
              vrR #*\ Vector.divide (Array.map f dR) scal ##*# vlAdjR
         else error "liftReal: non-real eigenvalues"

matrixRealPart ::
   (ArrayMatrix pack property lower upper meas vert horiz height width ~ matrix,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Class.Real a) =>
   matrix (Complex a) -> matrix a
matrixRealPart =
   ArrMatrix.Array . Array.map Complex.realPart . ArrMatrix.unwrap


instance
   (neg ~ True, zero ~ True, pos ~ True) =>
      LiftReal (Omni.Hermitian neg zero pos) where
   liftReal f a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerDiagonal -> liftHermitianDiagonal f a
         Omni.PowerHermitian -> liftHermitianReal f a
         _ -> error "Hermitian.liftReal: impossible shape"

instance LiftReal Omni.Symmetric where
   liftReal f a =
      case Omni.powerSingleton $ ArrMatrix.shape a of
         Omni.PowerSymmetric ->
            Symmetric.fromHermitian .
            liftHermitianReal f . Hermitian.fromSymmetric $ a
         _ -> error "Symmetric.liftReal: impossible shape"

liftHermitianReal ::
   (Layout.Packing pack, Shape.C sh, Class.Real a) =>
   (a -> a) -> HermitianP pack sh a -> HermitianP pack sh a
liftHermitianReal f a =
   case Scalar.precisionOfFunctor a of
      Scalar.Float -> liftHermitian f a
      Scalar.Double -> liftHermitian f a


norm1 ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Class.Floating a) =>
   ArrMatrix.ArrayMatrix pack prop lower upper meas vert horiz height width a ->
   Scalar.RealOf a
norm1 = Vector.norm1 . ArrMatrix.unwrap

normInf ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Class.Floating a) =>
   ArrMatrix.ArrayMatrix pack prop lower upper meas vert horiz height width a ->
   Scalar.RealOf a
normInf = Vector.normInf . ArrMatrix.unwrap


liftDiagonal ::
   (Layout.Packing pack, Shape.C sh, Class.Floating a, Class.Floating b) =>
   (a -> b) ->
   Quadratic pack Omni.Arbitrary Layout.Empty Layout.Empty sh a ->
   Quadratic pack Omni.Arbitrary Layout.Empty Layout.Empty sh b
liftDiagonal f = Diagonal.lift $ Array.map f

liftHermitianDiagonal ::
   (Layout.Packing pack,
    Bool.C neg, Bool.C zero, Bool.C pos, Omni.Hermitian neg zero pos ~ herm,
    Shape.C sh, Class.Floating a, Class.Floating b) =>
   (a -> b) ->
   Quadratic pack herm Layout.Empty Layout.Empty sh a ->
   Quadratic pack herm Layout.Empty Layout.Empty sh b
liftHermitianDiagonal f a =
   case packTag a of
      Layout.Packed -> ArrMatrix.lift1 (Array.map f) a
      Layout.Unpacked ->
         let b = Square.diagonal $ Array.map f $ Quad.takeDiagonal a
         in ArrMatrix.liftUnpacked1 id $
            if ArrMatrix.order a == ArrMatrix.order b
               then b
               else Square.transpose b

liftHermitian ::
   (Layout.Packing pack,
    Bool.C neg, Bool.C zero, Bool.C pos, Omni.Hermitian neg zero pos ~ herm,
    Shape.C sh, Class.Floating a, Scalar.RealOf a ~ ar, Storable ar) =>
   (ar -> ar) ->
   Quadratic pack herm Layout.Filled Layout.Filled sh a ->
   Quadratic pack herm Layout.Filled Layout.Filled sh a
liftHermitian f = applyPermutable $ applyUnchecked $ \a ->
   let (q,d) = Hermitian.eigensystem a
   in ArrMatrix.lift1 id $
      Hermitian.congruenceDiagonalAdjoint (Square.toFull q) (Array.map f d)


applyPermutable ::
   (Shape.C sh, Class.Floating a) =>
   (Quadratic pack prop lower upper (ExtShape.IntIndexed sh) a ->
    Quadratic pack prop lower upper (ExtShape.IntIndexed sh) a) ->
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
applyPermutable f =
   Quad.mapSize ExtShape.deconsIntIndexed .
   f .
   Quad.mapSize ExtShape.IntIndexed

applyUnchecked ::
   (Shape.C sh, Class.Floating a) =>
   (Quadratic pack prop lower upper (ShapePriv.Unchecked sh) a ->
    Quadratic pack prop lower upper (ShapePriv.Unchecked sh) a) ->
   Quadratic pack prop lower upper sh a ->
   Quadratic pack prop lower upper sh a
applyUnchecked f =
   Quad.mapSize ShapePriv.deconsUnchecked .
   f .
   Quad.mapSize ShapePriv.Unchecked