{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE FlexibleInstances #-} module Numeric.LAPACK.Matrix.Plain.Format where 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.Output as Output import Numeric.LAPACK.Output (Output, formatRow, (<+>)) import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor, ColumnMajor), UnaryProxy) import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Scalar (conjugate) import Numeric.LAPACK.Private (caseRealComplexFunc) import qualified Numeric.Netlib.Class as Class import qualified Type.Data.Num.Unary.Literal as TypeNum import qualified Type.Data.Num.Unary as Unary import Type.Data.Num (integralFromProxy) import qualified Data.Array.Comfort.Storable.Unchecked as Array import qualified Data.Array.Comfort.Boxed as BoxedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable (Array) import Data.Array.Comfort.Shape ((::+)) import Text.Printf (PrintfArg, printf) import qualified Data.List.Match as Match import qualified Data.List.HT as ListHT import Data.Functor.Compose (Compose(Compose, getCompose)) import Data.Foldable (Foldable, fold) import Data.List (mapAccumL, transpose) import Data.Complex (Complex((:+))) import Data.Ix (Ix) deflt :: String deflt = "%.4g" class (Shape.C sh) => FormatArray sh where {- We use constraint @(Class.Floating a)@ and not @(Format a)@ because it allows us to align the components of complex numbers. -} formatArray :: (Class.Floating a, Output out) => String -> Array sh a -> out instance (Integral i) => FormatArray (Shape.ZeroBased i) where formatArray = formatVector instance (Integral i) => FormatArray (Shape.OneBased i) where formatArray = formatVector instance (Ix i) => FormatArray (Shape.Range i) where formatArray = formatVector instance (Integral i) => FormatArray (Shape.Shifted i) where formatArray = formatVector instance (Enum enum, Bounded enum) => FormatArray (Shape.Enumeration enum) where formatArray = formatVector instance (FormatArray sh) => FormatArray (Shape.Deferred sh) where formatArray fmt = formatArray fmt . Array.mapShape (\(Shape.Deferred sh) -> sh) instance (FormatArray sh0, FormatArray sh1) => FormatArray (sh0::+sh1) where formatArray fmt v = formatArray fmt (Vector.takeLeft v) <+> formatArray fmt (Vector.takeRight v) formatVector :: (Shape.C sh, Class.Floating a, Output out) => String -> Array sh a -> out formatVector fmt = formatRow . map (Output.text . fold . printfFloating fmt) . Array.toList instance (Extent.Measure meas, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) => FormatArray (Layout.Full meas vert horiz height width) where formatArray = formatFull formatFull :: (Extent.Measure meas, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width, Output out, Class.Floating a) => String -> Full meas vert horiz height width a -> out formatFull fmt m = let Layout.Full order extent = Array.shape m in formatAligned (printfFloating fmt) $ splitRows order (Extent.dimensions extent) $ Array.toList m instance (Eq lower, Extent.Measure meas, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) => FormatArray (Layout.Split lower meas vert horiz height width) where formatArray = formatHouseholder formatHouseholder :: (Eq lower, Extent.Measure meas, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width, Class.Floating a, Output out) => String -> Array (Layout.Split lower meas vert horiz height width) a -> out formatHouseholder fmt m = let Layout.Split _ order extent = Array.shape m in formatSeparateTriangle (printfFloating fmt) $ splitRows order (Extent.dimensions extent) $ Array.toList m instance (Shape.C size) => FormatArray (Layout.Hermitian size) where formatArray = formatMirrored conjugate instance (Shape.C size) => FormatArray (Layout.Symmetric size) where formatArray = formatMirrored id formatMirrored :: (Shape.C size, Class.Floating a, Output out) => (a -> a) -> String -> Array (Layout.Mosaic Layout.Packed mirror Shape.Upper size) a -> out formatMirrored adapt fmt m = let Layout.Mosaic Layout.Packed _mirror Layout.Upper order size = Array.shape m in formatSeparateTriangle (printfFloating fmt) $ complementTriangle adapt order (Shape.size size) $ Array.toList m complementTriangle :: (Class.Floating a) => (a -> a) -> Order -> Int -> [a] -> [[a]] complementTriangle adapt order n xs = let mergeTriangles lower upper = zipWith (++) (map (map adapt . init) lower) upper in case order of RowMajor -> let tri = slice (take n $ iterate pred n) xs trans = reverse $ transpose $ map reverse tri in mergeTriangles trans tri ColumnMajor -> let tri = slice (take n [1..]) xs in mergeTriangles tri (transpose tri) formatDiagonal :: (Shape.C size, Class.Floating a, Output out) => String -> Order -> size -> [a] -> out formatDiagonal fmt order size xs = let n0 = Unary.unary TypeNum.u0 in formatAligned (printfFloatingMaybe fmt) $ padBanded (n0,n0) order (size,size) xs instance (Layout.UpLo uplo, Shape.C size) => FormatArray (Layout.Triangular uplo size) where formatArray = formatTriangular formatTriangular :: (Layout.UpLo uplo, Shape.C size, Class.Floating a, Output out) => String -> Array (Layout.Triangular uplo size) a -> out formatTriangular fmt m = let Layout.Mosaic Layout.Packed Layout.NoMirror uplo order size = Array.shape m in formatAligned (printfFloatingMaybe fmt) $ padTriangle uplo order (Shape.size size) $ Array.toList m padTriangle :: Layout.UpLoSingleton uplo -> Order -> Int -> [a] -> [[Maybe a]] padTriangle uplo = case uplo of Layout.Lower -> padLowerTriangle Layout.Upper -> padUpperTriangle padUpperTriangle :: Order -> Int -> [a] -> [[Maybe a]] padUpperTriangle order n xs = let mxs = map Just xs nothings = iterate (Nothing:) [] in case order of RowMajor -> zipWith (++) nothings (slice (take n $ iterate pred n) mxs) ColumnMajor -> transpose $ zipWith (++) (slice (take n [1..]) mxs) (reverse $ take n nothings) padLowerTriangle :: Order -> Int -> [a] -> [[Maybe a]] padLowerTriangle order n xs = map (map Just) $ case order of RowMajor -> slice (take n [1..]) xs ColumnMajor -> foldr (\(y:ys) zs -> [y] : zipWith (:) ys zs) [] (slice (take n $ iterate pred n) xs) slice :: [Int] -> [a] -> [[a]] slice ns xs = snd $ mapAccumL (\ys n -> let (vs,ws) = splitAt n ys in (ws,vs)) xs ns instance (Unary.Natural sub, Unary.Natural super, Extent.Measure meas, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) => FormatArray (Layout.Banded sub super meas vert horiz height width) where formatArray fmt m = let Layout.Banded offDiag order extent = Array.shape m in formatAligned (printfFloatingMaybe fmt) $ padBanded offDiag order (Extent.dimensions extent) $ Array.toList m padBanded :: (Shape.C height, Shape.C width, Unary.Natural sub, Unary.Natural super) => (UnaryProxy sub, UnaryProxy super) -> Order -> (height, width) -> [a] -> [[Maybe a]] padBanded (sub,super) order (height,width) xs = let slices = ListHT.sliceVertical (Layout.bandedBreadth (sub,super)) xs m = Shape.size height n = Shape.size width in case order of RowMajor -> map (take n) $ zipWith (shiftRow Nothing) (iterate (1+) (- integralFromProxy sub)) (map (map Just) slices) ColumnMajor -> let ku = integralFromProxy super in take m $ drop ku $ foldr (\col mat -> zipWith (:) (map Just col ++ repeat Nothing) ([]:mat)) (replicate (ku + m - n) []) slices instance (Unary.Natural offDiag, Shape.C size) => FormatArray (Layout.BandedHermitian offDiag size) where formatArray fmt m = let Layout.BandedHermitian offDiag order size = Array.shape m in formatSeparateTriangle (printfFloatingMaybe fmt) $ padBandedHermitian offDiag order size $ Array.toList m padBandedHermitian :: (Unary.Natural offDiag, Shape.C size, Class.Floating a) => UnaryProxy offDiag -> Order -> size -> [a] -> [[Maybe a]] padBandedHermitian offDiag order _size xs = let k = integralFromProxy offDiag slices = ListHT.sliceVertical (k + 1) xs in case order of RowMajor -> foldr (\row square -> Match.take ([]:square) (map Just row) : zipWith (:) (tail $ map (Just . conjugate) row ++ repeat Nothing) square) [] slices ColumnMajor -> zipWith (shiftRow Nothing) (iterate (1+) (-k)) $ map (map Just) $ zipWith (++) (map (map conjugate . init) slices) (drop k $ foldr (\column band -> zipWith (++) (map (:[]) column ++ repeat []) ([]:band)) (replicate k []) slices) shiftRow :: a -> Int -> [a] -> [a] shiftRow pad k = if k<=0 then drop (-k) else (replicate k pad ++) splitRows :: (Shape.C height, Shape.C width) => Order -> (height, width) -> [a] -> [[a]] splitRows order (height,width) = case order of RowMajor -> ListHT.sliceVertical (Shape.size width) ColumnMajor -> ListHT.sliceHorizontal (Shape.size height) formatAligned :: (Functor f, Foldable f) => Output out => (a -> f String) -> [[a]] -> out formatAligned printFmt = Output.formatAligned . map (map (fmap Output.text . printFmt)) formatSeparateTriangle :: (Functor f, Foldable f) => Output out => (a -> f String) -> [[a]] -> out formatSeparateTriangle printFmt = Output.formatSeparateTriangle . map (map (fmap Output.text . printFmt)) data TupleShape a = TupleShape instance (Class.Floating a) => Shape.C (TupleShape a) where size sh = caseRealComplexFunc sh 1 2 type Tuple a = BoxedArray.Array (TupleShape a) fillTuple :: (Class.Floating a) => b -> Tuple a b fillTuple = BoxedArray.replicate TupleShape newtype ToTuple a = ToTuple {getToTuple :: a -> Tuple a String} printfFloating :: (Class.Floating a) => String -> a -> Tuple a String printfFloating fmt = getToTuple $ Class.switchFloating (ToTuple $ fillTuple . printf fmt) (ToTuple $ fillTuple . printf fmt) (ToTuple $ printfComplex fmt) (ToTuple $ printfComplex fmt) printfFloatingMaybe :: (Class.Floating a) => String -> Maybe a -> Tuple a String printfFloatingMaybe = maybe (fillTuple "") . printfFloating printfComplex :: (Class.Real a) => String -> Complex a -> Tuple (Complex a) String printfComplex fmt = getToTuple $ getCompose $ Class.switchReal (Compose $ ToTuple $ printfComplexAux fmt) (Compose $ ToTuple $ printfComplexAux fmt) printfComplexAux :: (PrintfArg a, Class.Real a) => String -> Complex a -> Tuple (Complex a) String printfComplexAux fmt (r:+i) = if i<0 || isNegativeZero i then complexTuple (printf (fmt ++ "-") r) (printf (fmt ++ "i") (-i)) else complexTuple (printf (fmt ++ "+") r) (printf (fmt ++ "i") i) complexTuple :: (Class.Real a) => b -> b -> Tuple (Complex a) b complexTuple b0 b1 = BoxedArray.fromList TupleShape [b0,b1]