{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

-----------------------------------------------------------------------------
{- |
Module      :  Internal.Util
Copyright   :  (c) Alberto Ruiz 2013
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

-}
-----------------------------------------------------------------------------

module Internal.Util(

    -- * Convenience functions
    vector, matrix,
    disp,
    formatSparse,
    approxInt,
    dispDots,
    dispBlanks,
    formatShort,
    dispShort,
    zeros, ones,
    diagl,
    row,
    col,
    (&), (¦), (|||), (——), (===),
    (?), (¿),
    Indexable(..), size,
    Numeric,
    rand, randn,
    cross,
    norm,
    ℕ,ℤ,ℝ,ℂ,iC,
    Normed(..), norm_Frob, norm_nuclear,
    magnit,
    normalize,
    mt,
    (~!~),
    pairwiseD2,
    rowOuters,
    null1,
    null1sym,
    -- * Convolution
    -- ** 1D
    corr, conv, corrMin,
    -- ** 2D
    corr2, conv2, separable,
    block2x2,block3x3,view1,unView1,foldMatrix,
    gaussElim_1, gaussElim_2, gaussElim,
    luST, luSolve', luSolve'', luPacked', luPacked'',
    invershur
) where

import Internal.Vector
import Internal.Matrix hiding (size)
import Internal.Numeric
import Internal.Element
import Internal.Container
import Internal.Vectorized
import Internal.IO
import Internal.Algorithms hiding (Normed,linearSolve',luSolve', luPacked')
import Numeric.Matrix()
import Numeric.Vector()
import Internal.Random
import Internal.Convolution
import Control.Monad(when,forM_)
import Text.Printf
import Data.List.Split(splitOn)
import Data.List(intercalate,sortBy,foldl')
import Control.Arrow((&&&),(***))
import Data.Complex
import Data.Function(on)
import Internal.ST
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

type ℝ = Double
type ℕ = Int
type ℤ = Int
type ℂ = Complex Double

-- | imaginary unit
iC :: C
iC :: C
iC = Double
0Double -> Double -> C
forall a. a -> a -> Complex a
:+Double
1

{- | Create a real vector.

>>> vector [1..5]
[1.0,2.0,3.0,4.0,5.0]
it :: Vector R

-}
vector :: [R] -> Vector R
vector :: [Double] -> Vector Double
vector = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList

{- | Create a real matrix.

>>> matrix 5 [1..15]
(3><5)
 [  1.0,  2.0,  3.0,  4.0,  5.0
 ,  6.0,  7.0,  8.0,  9.0, 10.0
 , 11.0, 12.0, 13.0, 14.0, 15.0 ]

-}
matrix
  :: Int -- ^ number of columns
  -> [R] -- ^ elements in row order
  -> Matrix R
matrix :: Int -> [Double] -> Matrix Double
matrix Int
c = Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
c (Vector Double -> Matrix Double)
-> ([Double] -> Vector Double) -> [Double] -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList


{- | print a real matrix with given number of digits after the decimal point

>>> disp 5 $ ident 2 / 3
2x2
0.33333  0.00000
0.00000  0.33333

-}
disp :: Int -> Matrix Double -> IO ()

disp :: Int -> Matrix Double -> IO ()
disp Int
n = String -> IO ()
putStr (String -> IO ())
-> (Matrix Double -> String) -> Matrix Double -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Matrix Double -> String
dispf Int
n


{- | create a real diagonal matrix from a list

>>> diagl [1,2,3]
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]

-}
diagl :: [Double] -> Matrix Double
diagl :: [Double] -> Matrix Double
diagl = Vector Double -> Matrix Double
forall a. (Num a, Element a) => Vector a -> Matrix a
diag (Vector Double -> Matrix Double)
-> ([Double] -> Vector Double) -> [Double] -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList

-- | a real matrix of zeros
zeros :: Int -- ^ rows
      -> Int -- ^ columns
      -> Matrix Double
zeros :: Int -> Int -> Matrix Double
zeros Int
r Int
c = Double -> (Int, Int) -> Matrix Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst Double
0 (Int
r,Int
c)

-- | a real matrix of ones
ones :: Int -- ^ rows
     -> Int -- ^ columns
     -> Matrix Double
ones :: Int -> Int -> Matrix Double
ones Int
r Int
c = Double -> (Int, Int) -> Matrix Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst Double
1 (Int
r,Int
c)

-- | concatenation of real vectors
infixl 3 &
(&) :: Vector Double -> Vector Double -> Vector Double
Vector Double
a & :: Vector Double -> Vector Double -> Vector Double
& Vector Double
b = [Vector Double] -> Vector Double
forall t. Storable t => [Vector t] -> Vector t
vjoin [Vector Double
a,Vector Double
b]

{- | horizontal concatenation

>>> ident 3 ||| konst 7 (3,4)
(3><7)
 [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
 , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
 , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]

-}
infixl 3 |||
(|||) :: Element t => Matrix t -> Matrix t -> Matrix t
Matrix t
a ||| :: Matrix t -> Matrix t -> Matrix t
||| Matrix t
b = [[Matrix t]] -> Matrix t
forall t. Element t => [[Matrix t]] -> Matrix t
fromBlocks [[Matrix t
a,Matrix t
b]]

-- | a synonym for ('|||') (unicode 0x00a6, broken bar)
infixl 3 ¦
(¦) :: Matrix Double -> Matrix Double -> Matrix Double
¦ :: Matrix Double -> Matrix Double -> Matrix Double
(¦) = Matrix Double -> Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t -> Matrix t
(|||)


-- | vertical concatenation
--
(===) :: Element t => Matrix t -> Matrix t -> Matrix t
infixl 2 ===
Matrix t
a === :: Matrix t -> Matrix t -> Matrix t
=== Matrix t
b = [[Matrix t]] -> Matrix t
forall t. Element t => [[Matrix t]] -> Matrix t
fromBlocks [[Matrix t
a],[Matrix t
b]]

-- | a synonym for ('===') (unicode 0x2014, em dash)
(——) :: Matrix Double -> Matrix Double -> Matrix Double
infixl 2 ——
—— :: Matrix Double -> Matrix Double -> Matrix Double
(——) = Matrix Double -> Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t -> Matrix t
(===)


-- | create a single row real matrix from a list
--
-- >>> row [2,3,1,8]
-- (1><4)
--  [ 2.0, 3.0, 1.0, 8.0 ]
--
row :: [Double] -> Matrix Double
row :: [Double] -> Matrix Double
row = Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
asRow (Vector Double -> Matrix Double)
-> ([Double] -> Vector Double) -> [Double] -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList

-- | create a single column real matrix from a list
--
-- >>> col [7,-2,4]
-- (3><1)
--  [  7.0
--  , -2.0
--  ,  4.0 ]
--
col :: [Double] -> Matrix Double
col :: [Double] -> Matrix Double
col = Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
asColumn (Vector Double -> Matrix Double)
-> ([Double] -> Vector Double) -> [Double] -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList

{- | extract rows

>>> (20><4) [1..] ? [2,1,1]
(3><4)
 [ 9.0, 10.0, 11.0, 12.0
 , 5.0,  6.0,  7.0,  8.0
 , 5.0,  6.0,  7.0,  8.0 ]

-}
infixl 9 ?
(?) :: Element t => Matrix t -> [Int] -> Matrix t
? :: Matrix t -> [Int] -> Matrix t
(?) = ([Int] -> Matrix t -> Matrix t) -> Matrix t -> [Int] -> Matrix t
forall a b c. (a -> b -> c) -> b -> a -> c
flip [Int] -> Matrix t -> Matrix t
forall t. Element t => [Int] -> Matrix t -> Matrix t
extractRows

{- | extract columns

(unicode 0x00bf, inverted question mark, Alt-Gr ?)

>>> (3><4) [1..] ¿ [3,0]
(3><2)
 [  4.0, 1.0
 ,  8.0, 5.0
 , 12.0, 9.0 ]

-}
infixl 9 ¿
(¿) :: Element t => Matrix t -> [Int] -> Matrix t
¿ :: Matrix t -> [Int] -> Matrix t
(¿)= ([Int] -> Matrix t -> Matrix t) -> Matrix t -> [Int] -> Matrix t
forall a b c. (a -> b -> c) -> b -> a -> c
flip [Int] -> Matrix t -> Matrix t
forall t. Element t => [Int] -> Matrix t -> Matrix t
extractColumns


cross :: Product t => Vector t -> Vector t -> Vector t
-- ^ cross product (for three-element vectors)
cross :: Vector t -> Vector t -> Vector t
cross Vector t
x Vector t
y | Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 Bool -> Bool -> Bool
&& Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 = [t] -> Vector t
forall a. Storable a => [a] -> Vector a
fromList [t
z1,t
z2,t
z3]
          | Bool
otherwise = String -> Vector t
forall a. HasCallStack => String -> a
error (String -> Vector t) -> String -> Vector t
forall a b. (a -> b) -> a -> b
$ String
"the cross product requires 3-element vectors (sizes given: "
                                String -> String -> String
forall a. [a] -> [a] -> [a]
++Int -> String
forall a. Show a => a -> String
show (Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
x)String -> String -> String
forall a. [a] -> [a] -> [a]
++String
" and "String -> String -> String
forall a. [a] -> [a] -> [a]
++Int -> String
forall a. Show a => a -> String
show (Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
y)String -> String -> String
forall a. [a] -> [a] -> [a]
++String
")"
  where
    [t
x1,t
x2,t
x3] = Vector t -> [t]
forall a. Storable a => Vector a -> [a]
toList Vector t
x
    [t
y1,t
y2,t
y3] = Vector t -> [t]
forall a. Storable a => Vector a -> [a]
toList Vector t
y
    z1 :: t
z1 = t
x2t -> t -> t
forall a. Num a => a -> a -> a
*t
y3t -> t -> t
forall a. Num a => a -> a -> a
-t
x3t -> t -> t
forall a. Num a => a -> a -> a
*t
y2
    z2 :: t
z2 = t
x3t -> t -> t
forall a. Num a => a -> a -> a
*t
y1t -> t -> t
forall a. Num a => a -> a -> a
-t
x1t -> t -> t
forall a. Num a => a -> a -> a
*t
y3
    z3 :: t
z3 = t
x1t -> t -> t
forall a. Num a => a -> a -> a
*t
y2t -> t -> t
forall a. Num a => a -> a -> a
-t
x2t -> t -> t
forall a. Num a => a -> a -> a
*t
y1

{-# SPECIALIZE cross :: Vector Double -> Vector Double -> Vector Double #-}
{-# SPECIALIZE cross :: Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) #-}

norm :: Vector Double -> Double
-- ^ 2-norm of real vector
norm :: Vector Double -> Double
norm = NormType -> Vector Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2

-- | p-norm for vectors, operator norm for matrices
class Normed a
  where
    norm_0   :: a -> R
    norm_1   :: a -> R
    norm_2   :: a -> R
    norm_Inf :: a -> R


instance Normed (Vector R)
  where
    norm_0 :: Vector Double -> Double
norm_0 Vector Double
v = Vector Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
sumElements (Vector Double -> Vector Double
forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
step (Vector Double -> Vector Double
forall a. Num a => a -> a
abs Vector Double
v Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
- Double -> Vector Double
forall (c :: * -> *) e. Container c e => e -> c e
scalar (Double
epsDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Vector Double -> RealOf Double
forall e. Product e => Vector e -> RealOf e
normInf Vector Double
v)))
    norm_1 :: Vector Double -> Double
norm_1 = NormType -> Vector Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1
    norm_2 :: Vector Double -> Double
norm_2 = NormType -> Vector Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2
    norm_Inf :: Vector Double -> Double
norm_Inf = NormType -> Vector Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity

instance Normed (Vector C)
  where
    norm_0 :: Vector C -> Double
norm_0 Vector C
v = Vector Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
sumElements (Vector Double -> Vector Double
forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
step ((Vector Double, Vector Double) -> Vector Double
forall a b. (a, b) -> a
fst (Vector C -> (Vector Double, Vector Double)
forall t (c :: * -> *).
(Convert t, Complexable c, RealElement t) =>
c (Complex t) -> (c t, c t)
fromComplex (Vector C -> Vector C
forall a. Num a => a -> a
abs Vector C
v)) Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
- Double -> Vector Double
forall (c :: * -> *) e. Container c e => e -> c e
scalar (Double
epsDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Vector C -> RealOf C
forall e. Product e => Vector e -> RealOf e
normInf Vector C
v)))
    norm_1 :: Vector C -> Double
norm_1 = NormType -> Vector C -> RealOf C
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1
    norm_2 :: Vector C -> Double
norm_2 = NormType -> Vector C -> RealOf C
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2
    norm_Inf :: Vector C -> Double
norm_Inf = NormType -> Vector C -> RealOf C
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity

instance Normed (Matrix R)
  where
    norm_0 :: Matrix Double -> Double
norm_0 = Vector Double -> Double
forall a. Normed a => a -> Double
norm_0 (Vector Double -> Double)
-> (Matrix Double -> Vector Double) -> Matrix Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
flatten
    norm_1 :: Matrix Double -> Double
norm_1 = NormType -> Matrix Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1
    norm_2 :: Matrix Double -> Double
norm_2 = NormType -> Matrix Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2
    norm_Inf :: Matrix Double -> Double
norm_Inf = NormType -> Matrix Double -> RealOf Double
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity

instance Normed (Matrix C)
  where
    norm_0 :: Matrix C -> Double
norm_0 = Vector C -> Double
forall a. Normed a => a -> Double
norm_0 (Vector C -> Double)
-> (Matrix C -> Vector C) -> Matrix C -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix C -> Vector C
forall t. Element t => Matrix t -> Vector t
flatten
    norm_1 :: Matrix C -> Double
norm_1 = NormType -> Matrix C -> RealOf C
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1
    norm_2 :: Matrix C -> Double
norm_2 = NormType -> Matrix C -> RealOf C
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2
    norm_Inf :: Matrix C -> Double
norm_Inf = NormType -> Matrix C -> RealOf C
forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity

instance Normed (Vector I)
  where
    norm_0 :: Vector I -> Double
norm_0 = I -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (I -> Double) -> (Vector I -> I) -> Vector I -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector I -> I
forall (c :: * -> *) e. Container c e => c e -> e
sumElements (Vector I -> I) -> (Vector I -> Vector I) -> Vector I -> I
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector I -> Vector I
forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
step (Vector I -> Vector I)
-> (Vector I -> Vector I) -> Vector I -> Vector I
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector I -> Vector I
forall a. Num a => a -> a
abs
    norm_1 :: Vector I -> Double
norm_1 = I -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (I -> Double) -> (Vector I -> I) -> Vector I -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector I -> I
forall e. Product e => Vector e -> RealOf e
norm1
    norm_2 :: Vector I -> Double
norm_2 Vector I
v = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> (I -> Double) -> I -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. I -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (I -> Double) -> I -> Double
forall a b. (a -> b) -> a -> b
$ Vector I -> Vector I -> I
forall t. Numeric t => Vector t -> Vector t -> t
dot Vector I
v Vector I
v
    norm_Inf :: Vector I -> Double
norm_Inf = I -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (I -> Double) -> (Vector I -> I) -> Vector I -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector I -> I
forall e. Product e => Vector e -> RealOf e
normInf

instance Normed (Vector Z)
  where
    norm_0 :: Vector Z -> Double
norm_0 = Z -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Z -> Double) -> (Vector Z -> Z) -> Vector Z -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Z -> Z
forall (c :: * -> *) e. Container c e => c e -> e
sumElements (Vector Z -> Z) -> (Vector Z -> Vector Z) -> Vector Z -> Z
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Z -> Vector Z
forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
step (Vector Z -> Vector Z)
-> (Vector Z -> Vector Z) -> Vector Z -> Vector Z
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Z -> Vector Z
forall a. Num a => a -> a
abs
    norm_1 :: Vector Z -> Double
norm_1 = Z -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Z -> Double) -> (Vector Z -> Z) -> Vector Z -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Z -> Z
forall e. Product e => Vector e -> RealOf e
norm1
    norm_2 :: Vector Z -> Double
norm_2 Vector Z
v = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> (Z -> Double) -> Z -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Z -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Z -> Double) -> Z -> Double
forall a b. (a -> b) -> a -> b
$ Vector Z -> Vector Z -> Z
forall t. Numeric t => Vector t -> Vector t -> t
dot Vector Z
v Vector Z
v
    norm_Inf :: Vector Z -> Double
norm_Inf = Z -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Z -> Double) -> (Vector Z -> Z) -> Vector Z -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Z -> Z
forall e. Product e => Vector e -> RealOf e
normInf

instance Normed (Vector Float)
  where
    norm_0 :: Vector Float -> Double
norm_0 = Vector Double -> Double
forall a. Normed a => a -> Double
norm_0 (Vector Double -> Double)
-> (Vector Float -> Vector Double) -> Vector Float -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Float -> Vector Double
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_1 :: Vector Float -> Double
norm_1 = Vector Double -> Double
forall a. Normed a => a -> Double
norm_1 (Vector Double -> Double)
-> (Vector Float -> Vector Double) -> Vector Float -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Float -> Vector Double
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_2 :: Vector Float -> Double
norm_2 = Vector Double -> Double
forall a. Normed a => a -> Double
norm_2 (Vector Double -> Double)
-> (Vector Float -> Vector Double) -> Vector Float -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Float -> Vector Double
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_Inf :: Vector Float -> Double
norm_Inf = Vector Double -> Double
forall a. Normed a => a -> Double
norm_Inf (Vector Double -> Double)
-> (Vector Float -> Vector Double) -> Vector Float -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Float -> Vector Double
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double

instance Normed (Vector (Complex Float))
  where
    norm_0 :: Vector (Complex Float) -> Double
norm_0 = Vector C -> Double
forall a. Normed a => a -> Double
norm_0 (Vector C -> Double)
-> (Vector (Complex Float) -> Vector C)
-> Vector (Complex Float)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Complex Float) -> Vector C
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_1 :: Vector (Complex Float) -> Double
norm_1 = Vector C -> Double
forall a. Normed a => a -> Double
norm_1 (Vector C -> Double)
-> (Vector (Complex Float) -> Vector C)
-> Vector (Complex Float)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Complex Float) -> Vector C
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_2 :: Vector (Complex Float) -> Double
norm_2 = Vector C -> Double
forall a. Normed a => a -> Double
norm_2 (Vector C -> Double)
-> (Vector (Complex Float) -> Vector C)
-> Vector (Complex Float)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Complex Float) -> Vector C
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_Inf :: Vector (Complex Float) -> Double
norm_Inf = Vector C -> Double
forall a. Normed a => a -> Double
norm_Inf (Vector C -> Double)
-> (Vector (Complex Float) -> Vector C)
-> Vector (Complex Float)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Complex Float) -> Vector C
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double

-- | Frobenius norm (Schatten p-norm with p=2)
norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
norm_Frob :: Matrix t -> Double
norm_Frob = Vector t -> Double
forall a. Normed a => a -> Double
norm_2 (Vector t -> Double)
-> (Matrix t -> Vector t) -> Matrix t -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten

-- | Sum of singular values (Schatten p-norm with p=1)
norm_nuclear :: Field t => Matrix t -> R
norm_nuclear :: Matrix t -> Double
norm_nuclear = Vector Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
sumElements (Vector Double -> Double)
-> (Matrix t -> Vector Double) -> Matrix t -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Vector Double
forall t. Field t => Matrix t -> Vector Double
singularValues

{- | Check if the absolute value or complex magnitude is greater than a given threshold

>>> magnit 1E-6 (1E-12 :: R)
False
>>> magnit 1E-6 (3+iC :: C)
True
>>> magnit 0 (3 :: I ./. 5)
True

-}
magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool
magnit :: Double -> t -> Bool
magnit Double
e t
x = Vector t -> Double
forall a. Normed a => a -> Double
norm_1 ([t] -> Vector t
forall a. Storable a => [a] -> Vector a
fromList [t
x]) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
e


-- | Obtains a vector in the same direction with 2-norm=1
normalize :: (Normed (Vector t), Num (Vector t), Field t) => Vector t -> Vector t
normalize :: Vector t -> Vector t
normalize Vector t
v = Vector t
v Vector t -> Vector t -> Vector t
forall a. Fractional a => a -> a -> a
/ Vector (RealOf t) -> Vector t
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c (RealOf t) -> c t
real (Double -> Vector Double
forall (c :: * -> *) e. Container c e => e -> c e
scalar (Vector t -> Double
forall a. Normed a => a -> Double
norm_2 Vector t
v))


-- | trans . inv
mt :: Matrix Double -> Matrix Double
mt :: Matrix Double -> Matrix Double
mt = Matrix Double -> Matrix Double
forall t. Matrix t -> Matrix t
trans (Matrix Double -> Matrix Double)
-> (Matrix Double -> Matrix Double)
-> Matrix Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> Matrix Double
forall t. Field t => Matrix t -> Matrix t
inv

--------------------------------------------------------------------------------
{- |

>>> size $ vector [1..10]
10
>>> size $ (2><5)[1..10::Double]
(2,5)

-}
size :: Container c t => c t -> IndexOf c
size :: c t -> IndexOf c
size = c t -> IndexOf c
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
size'

{- | Alternative indexing function.

>>> vector [1..10] ! 3
4.0

On a matrix it gets the k-th row as a vector:

>>> matrix 5 [1..15] ! 1
[6.0,7.0,8.0,9.0,10.0]
it :: Vector Double

>>> matrix 5 [1..15] ! 1 ! 3
9.0

-}
class Indexable c t | c -> t , t -> c
  where
    infixl 9 !
    (!) :: c -> Int -> t

instance Indexable (Vector Double) Double
  where
    (!) = Vector Double -> Int -> Double
forall t. Storable t => Vector t -> Int -> t
(@>)

instance Indexable (Vector Float) Float
  where
    (!) = Vector Float -> Int -> Float
forall t. Storable t => Vector t -> Int -> t
(@>)

instance Indexable (Vector I) I
  where
    (!) = Vector I -> Int -> I
forall t. Storable t => Vector t -> Int -> t
(@>)

instance Indexable (Vector Z) Z
  where
    (!) = Vector Z -> Int -> Z
forall t. Storable t => Vector t -> Int -> t
(@>)

instance Indexable (Vector (Complex Double)) (Complex Double)
  where
    (!) = Vector C -> Int -> C
forall t. Storable t => Vector t -> Int -> t
(@>)

instance Indexable (Vector (Complex Float)) (Complex Float)
  where
    (!) = Vector (Complex Float) -> Int -> Complex Float
forall t. Storable t => Vector t -> Int -> t
(@>)

instance Element t => Indexable (Matrix t) (Vector t)
  where
    Matrix t
m ! :: Matrix t -> Int -> Vector t
! Int
j = Int -> Int -> Vector t -> Vector t
forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
c) Int
c (Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten Matrix t
m)
      where
        c :: Int
c = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
m

--------------------------------------------------------------------------------

-- | Matrix of pairwise squared distances of row vectors
-- (using the matrix product trick in blog.smola.org)
pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
pairwiseD2 Matrix Double
x Matrix Double
y | Bool
ok = Vector Double
x2 Vector Double -> Vector Double -> Matrix Double
forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
oy Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Vector Double
ox Vector Double -> Vector Double -> Matrix Double
forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
y2 Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
- Matrix Double
2Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
* Matrix Double
x Matrix Double -> Matrix Double -> Matrix Double
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix Double -> Matrix Double
forall t. Matrix t -> Matrix t
trans Matrix Double
y
               | Bool
otherwise = String -> Matrix Double
forall a. HasCallStack => String -> a
error (String -> Matrix Double) -> String -> Matrix Double
forall a b. (a -> b) -> a -> b
$ String
"pairwiseD2 with different number of columns: "
                                   String -> String -> String
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> String
forall a. Show a => a -> String
show (Matrix Double -> IndexOf Matrix
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
size Matrix Double
x) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
", " String -> String -> String
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> String
forall a. Show a => a -> String
show (Matrix Double -> IndexOf Matrix
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
size Matrix Double
y)
  where
    ox :: Vector Double
ox = Int -> Vector Double
forall e d (c :: * -> *). (Konst e d c, Num e) => d -> c e
one (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
x)
    oy :: Vector Double
oy = Int -> Vector Double
forall e d (c :: * -> *). (Konst e d c, Num e) => d -> c e
one (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
y)
    oc :: Vector Double
oc = Int -> Vector Double
forall e d (c :: * -> *). (Konst e d c, Num e) => d -> c e
one (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
x)
    one :: d -> c e
one d
k = e -> d -> c e
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst e
1 d
k
    x2 :: Vector Double
x2 = Matrix Double
x Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
* Matrix Double
x Matrix Double -> Vector Double -> Vector Double
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Vector Double
oc
    y2 :: Vector Double
y2 = Matrix Double
y Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
* Matrix Double
y Matrix Double -> Vector Double -> Vector Double
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Vector Double
oc
    ok :: Bool
ok = Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
y

--------------------------------------------------------------------------------

{- | outer products of rows

>>> a
(3><2)
 [   1.0,   2.0
 ,  10.0,  20.0
 , 100.0, 200.0 ]
>>> b
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]

>>> rowOuters a (b ||| 1)
(3><8)
 [   1.0,   2.0,   3.0,   1.0,    2.0,    4.0,    6.0,   2.0
 ,  40.0,  50.0,  60.0,  10.0,   80.0,  100.0,  120.0,  20.0
 , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ]

-}
rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
rowOuters Matrix Double
a Matrix Double
b = Matrix Double
a' Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
* Matrix Double
b'
  where
    a' :: Matrix Double
a' = Matrix Double -> Matrix Double -> Matrix Double
forall t. Product t => Matrix t -> Matrix t -> Matrix t
kronecker Matrix Double
a (Int -> Int -> Matrix Double
ones Int
1 (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
b))
    b' :: Matrix Double
b' = Matrix Double -> Matrix Double -> Matrix Double
forall t. Product t => Matrix t -> Matrix t -> Matrix t
kronecker (Int -> Int -> Matrix Double
ones Int
1 (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
a)) Matrix Double
b

--------------------------------------------------------------------------------

-- | solution of overconstrained homogeneous linear system
null1 :: Matrix R -> Vector R
null1 :: Matrix Double -> Vector Double
null1 = [Vector Double] -> Vector Double
forall a. [a] -> a
last ([Vector Double] -> Vector Double)
-> (Matrix Double -> [Vector Double])
-> Matrix Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
toColumns (Matrix Double -> [Vector Double])
-> (Matrix Double -> Matrix Double)
-> Matrix Double
-> [Vector Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double, Matrix Double) -> Matrix Double
forall a b. (a, b) -> b
snd ((Vector Double, Matrix Double) -> Matrix Double)
-> (Matrix Double -> (Vector Double, Matrix Double))
-> Matrix Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> (Vector Double, Matrix Double)
forall t. Field t => Matrix t -> (Vector Double, Matrix t)
rightSV

-- | solution of overconstrained homogeneous symmetric linear system
null1sym :: Herm R -> Vector R
null1sym :: Herm Double -> Vector Double
null1sym = [Vector Double] -> Vector Double
forall a. [a] -> a
last ([Vector Double] -> Vector Double)
-> (Herm Double -> [Vector Double]) -> Herm Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
toColumns (Matrix Double -> [Vector Double])
-> (Herm Double -> Matrix Double) -> Herm Double -> [Vector Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double, Matrix Double) -> Matrix Double
forall a b. (a, b) -> b
snd ((Vector Double, Matrix Double) -> Matrix Double)
-> (Herm Double -> (Vector Double, Matrix Double))
-> Herm Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Herm Double -> (Vector Double, Matrix Double)
forall t. Field t => Herm t -> (Vector Double, Matrix t)
eigSH

--------------------------------------------------------------------------------

infixl 0 ~!~
Bool
c ~!~ :: Bool -> String -> f ()
~!~ String
msg = Bool -> f () -> f ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
c (String -> f ()
forall a. HasCallStack => String -> a
error String
msg)

--------------------------------------------------------------------------------

formatSparse :: String -> String -> String -> Int -> Matrix Double -> String

formatSparse :: String -> String -> String -> Int -> Matrix Double -> String
formatSparse String
zeroI String
_zeroF String
sep Int
_ (Matrix Double -> Maybe (Matrix Double)
approxInt -> Just Matrix Double
m) = String -> (Double -> String) -> Matrix Double -> String
forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
sep Double -> String
forall t. (Eq t, Num t, PrintfArg t) => t -> String
f Matrix Double
m
  where
    f :: t -> String
f t
0 = String
zeroI
    f t
x = String -> t -> String
forall r. PrintfType r => String -> r
printf String
"%.0f" t
x

formatSparse String
zeroI String
zeroF String
sep Int
n Matrix Double
m = String -> (Double -> String) -> Matrix Double -> String
forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
sep Double -> String
f Matrix Double
m
  where
    f :: Double -> String
f Double
x | Double -> Double
forall a. Num a => a -> a
abs (Double
x::Double) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall x. RealFloat x => x
peps = String
zeroIString -> String -> String
forall a. [a] -> [a] -> [a]
++String
zeroF
        | Double -> Double
forall a. Num a => a -> a
abs (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round Double
x::Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall x. RealFloat x => x
peps
            = String -> Double -> String
forall r. PrintfType r => String -> r
printf (String
"%.0f."String -> String -> String
forall a. [a] -> [a] -> [a]
++Int -> Char -> String
forall a. Int -> a -> [a]
replicate Int
n Char
' ') Double
x
        | Bool
otherwise = String -> Double -> String
forall r. PrintfType r => String -> r
printf (String
"%."String -> String -> String
forall a. [a] -> [a] -> [a]
++Int -> String
forall a. Show a => a -> String
show Int
nString -> String -> String
forall a. [a] -> [a] -> [a]
++String
"f") Double
x

approxInt :: Matrix Double -> Maybe (Matrix Double)
approxInt Matrix Double
m
    | Vector Double -> Double
forall a. Normed a => a -> Double
norm_Inf (Vector Double
v Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
- Vector Double
vi) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall x. RealFloat x => x
peps Double -> Double -> Double
forall a. Num a => a -> a -> a
* Vector Double -> Double
forall a. Normed a => a -> Double
norm_Inf Vector Double
v = Matrix Double -> Maybe (Matrix Double)
forall a. a -> Maybe a
Just (Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
reshape (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m) Vector Double
vi)
    | Bool
otherwise = Maybe (Matrix Double)
forall a. Maybe a
Nothing
  where
    v :: Vector Double
v = Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
flatten Matrix Double
m
    vi :: Vector Double
vi = Vector Double -> Vector Double
roundVector Vector Double
v

dispDots :: Int -> Matrix Double -> IO ()
dispDots Int
n = String -> IO ()
putStr (String -> IO ())
-> (Matrix Double -> String) -> Matrix Double -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> String -> String -> Int -> Matrix Double -> String
formatSparse String
"." (Int -> Char -> String
forall a. Int -> a -> [a]
replicate Int
n Char
' ') String
"  " Int
n

dispBlanks :: Int -> Matrix Double -> IO ()
dispBlanks Int
n = String -> IO ()
putStr (String -> IO ())
-> (Matrix Double -> String) -> Matrix Double -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> String -> String -> Int -> Matrix Double -> String
formatSparse String
"" String
"" String
"  " Int
n

formatShort :: String -> (t -> String) -> Int -> Int -> Matrix t -> String
formatShort String
sep t -> String
fmt Int
maxr Int
maxc Matrix t
m = String
auxm4
  where
    (Int
rm,Int
cm) = Matrix t -> IndexOf Matrix
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
size Matrix t
m
    (Int
r1,Int
r2,Int
r3)
        | Int
rm Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
maxr = (Int
rm,Int
0,Int
0)
        | Bool
otherwise  = (Int
maxrInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3,Int
rmInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
maxrInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
2)
    (Int
c1,Int
c2,Int
c3)
        | Int
cm Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
maxc = (Int
cm,Int
0,Int
0)
        | Bool
otherwise  = (Int
maxcInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3,Int
cmInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
maxcInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
2)
    [ [Matrix t
a,Matrix t
_,Matrix t
b]
     ,[Matrix t
_,Matrix t
_,Matrix t
_]
     ,[Matrix t
c,Matrix t
_,Matrix t
d]] = [Int] -> [Int] -> Matrix t -> [[Matrix t]]
forall t. Element t => [Int] -> [Int] -> Matrix t -> [[Matrix t]]
toBlocks [Int
r1,Int
r2,Int
r3]
                          [Int
c1,Int
c2,Int
c3] Matrix t
m
    auxm :: Matrix t
auxm = [[Matrix t]] -> Matrix t
forall t. Element t => [[Matrix t]] -> Matrix t
fromBlocks [[Matrix t
a,Matrix t
b],[Matrix t
c,Matrix t
d]]
    auxm2 :: String
auxm2
        | Int
cm Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
maxc = String -> (t -> String) -> Matrix t -> String
forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
"|" t -> String
fmt Matrix t
auxm
        | Bool
otherwise = String -> (t -> String) -> Matrix t -> String
forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
sep t -> String
fmt Matrix t
auxm
    auxm3 :: [String]
auxm3
        | Int
cm Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
maxc = (String -> String) -> [String] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map ([String] -> String
f ([String] -> String) -> (String -> [String]) -> String -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> String -> [String]
forall a. Eq a => [a] -> [a] -> [[a]]
splitOn String
"|") (String -> [String]
lines String
auxm2)
        | Bool
otherwise = (String -> [String]
lines String
auxm2)
    f :: [String] -> String
f [String]
items = String -> [String] -> String
forall a. [a] -> [[a]] -> [a]
intercalate String
sep (Int -> [String] -> [String]
forall a. Int -> [a] -> [a]
take (Int
maxcInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) [String]
items) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"  .. " String -> String -> String
forall a. [a] -> [a] -> [a]
++
              String -> [String] -> String
forall a. [a] -> [[a]] -> [a]
intercalate String
sep (Int -> [String] -> [String]
forall a. Int -> [a] -> [a]
drop (Int
maxcInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) [String]
items)
    auxm4 :: String
auxm4
        | Int
rm Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
maxr = [String] -> String
unlines (Int -> [String] -> [String]
forall a. Int -> [a] -> [a]
take (Int
maxrInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) [String]
auxm3 [String] -> [String] -> [String]
forall a. [a] -> [a] -> [a]
++ String
vsep String -> [String] -> [String]
forall a. a -> [a] -> [a]
: Int -> [String] -> [String]
forall a. Int -> [a] -> [a]
drop (Int
maxrInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) [String]
auxm3)
        | Bool
otherwise = [String] -> String
unlines [String]
auxm3
    vsep :: String
vsep = (Char -> Char) -> String -> String
forall a b. (a -> b) -> [a] -> [b]
map Char -> Char
g ([String] -> String
forall a. [a] -> a
head [String]
auxm3)
    g :: Char -> Char
g Char
'.' = Char
':'
    g Char
_ = Char
' '


dispShort :: Int -> Int -> Int -> Matrix Double -> IO ()
dispShort :: Int -> Int -> Int -> Matrix Double -> IO ()
dispShort Int
maxr Int
maxc Int
dec Matrix Double
m =
    String -> Int -> Int -> String -> IO ()
forall r. PrintfType r => String -> r
printf String
"%dx%d\n%s" (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m) (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m) (String
-> (Double -> String) -> Int -> Int -> Matrix Double -> String
forall t.
(Num t, Container Vector t) =>
String -> (t -> String) -> Int -> Int -> Matrix t -> String
formatShort String
"  " Double -> String
fmt Int
maxr Int
maxc Matrix Double
m)
  where
    fmt :: Double -> String
fmt = String -> Double -> String
forall r. PrintfType r => String -> r
printf (String
"%."String -> String -> String
forall a. [a] -> [a] -> [a]
++Int -> String
forall a. Show a => a -> String
show Int
dec String -> String -> String
forall a. [a] -> [a] -> [a]
++String
"f")

--------------------------------------------------------------------------------

-- matrix views

block2x2 :: Int -> Int -> Matrix t -> [[Matrix t]]
block2x2 Int
r Int
c Matrix t
m = [[Matrix t
m11,Matrix t
m12],[Matrix t
m21,Matrix t
m22]]
  where
    m11 :: Matrix t
m11 = Matrix t
m Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Take Int
r, Int -> Extractor
Take Int
c)
    m12 :: Matrix t
m12 = Matrix t
m Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Take Int
r, Int -> Extractor
Drop Int
c)
    m21 :: Matrix t
m21 = Matrix t
m Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Drop Int
r, Int -> Extractor
Take Int
c)
    m22 :: Matrix t
m22 = Matrix t
m Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Drop Int
r, Int -> Extractor
Drop Int
c)

block3x3 :: Int -> Int -> Int -> Int -> Matrix t -> [[Matrix t]]
block3x3 Int
r Int
nr Int
c Int
nc Matrix t
m = [[Matrix t
m Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? ([Extractor]
er [Extractor] -> Int -> Extractor
forall a. [a] -> Int -> a
!! Int
i, [Extractor]
ec [Extractor] -> Int -> Extractor
forall a. [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0..Int
2] ] | Int
i <- [Int
0..Int
2] ]
  where
    er :: [Extractor]
er = [ Int -> Int -> Int -> Extractor
Range Int
0 Int
1 (Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), Int -> Int -> Int -> Extractor
Range Int
r Int
1 (Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nrInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), Int -> Extractor
Drop (Int
nrInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
r) ]
    ec :: [Extractor]
ec = [ Int -> Int -> Int -> Extractor
Range Int
0 Int
1 (Int
cInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), Int -> Int -> Int -> Extractor
Range Int
c Int
1 (Int
cInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ncInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), Int -> Extractor
Drop (Int
ncInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
c) ]

view1 :: Numeric t => Matrix t -> Maybe (View1 t)
view1 :: Matrix t -> Maybe (View1 t)
view1 Matrix t
m
    | Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 = View1 t -> Maybe (View1 t)
forall a. a -> Maybe a
Just (t
e, Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten Matrix t
m12, Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten Matrix t
m21 , Matrix t
m22)
    | Bool
otherwise = Maybe (View1 t)
forall a. Maybe a
Nothing
  where
    [[Matrix t
m11,Matrix t
m12],[Matrix t
m21,Matrix t
m22]] = Int -> Int -> Matrix t -> [[Matrix t]]
forall t. Element t => Int -> Int -> Matrix t -> [[Matrix t]]
block2x2 Int
1 Int
1 Matrix t
m
    e :: t
e = Matrix t
m11 Matrix t -> IndexOf Matrix -> t
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
0, Int
0)

unView1 :: Numeric t => View1 t -> Matrix t
unView1 :: View1 t -> Matrix t
unView1 (t
e,Vector t
r,Vector t
c,Matrix t
m) = [[Matrix t]] -> Matrix t
forall t. Element t => [[Matrix t]] -> Matrix t
fromBlocks [[t -> Matrix t
forall (c :: * -> *) e. Container c e => e -> c e
scalar t
e, Vector t -> Matrix t
forall a. Storable a => Vector a -> Matrix a
asRow Vector t
r],[Vector t -> Matrix t
forall a. Storable a => Vector a -> Matrix a
asColumn Vector t
c, Matrix t
m]]

type View1 t = (t, Vector t, Vector t, Matrix t)

foldMatrix :: Numeric t => (Matrix t -> Matrix t) -> (View1 t -> View1 t) -> (Matrix t -> Matrix t)
foldMatrix :: (Matrix t -> Matrix t)
-> (View1 t -> View1 t) -> Matrix t -> Matrix t
foldMatrix Matrix t -> Matrix t
g View1 t -> View1 t
f ( (View1 t -> View1 t
f (View1 t -> View1 t) -> Maybe (View1 t) -> Maybe (View1 t)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>) (Maybe (View1 t) -> Maybe (View1 t))
-> (Matrix t -> Maybe (View1 t)) -> Matrix t -> Maybe (View1 t)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Maybe (View1 t)
forall t. Numeric t => Matrix t -> Maybe (View1 t)
view1 (Matrix t -> Maybe (View1 t))
-> (Matrix t -> Matrix t) -> Matrix t -> Maybe (View1 t)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Matrix t
g -> Just (t
e,Vector t
r,Vector t
c,Matrix t
m)) = View1 t -> Matrix t
forall t. Numeric t => View1 t -> Matrix t
unView1 (t
e, Vector t
r, Vector t
c, (Matrix t -> Matrix t)
-> (View1 t -> View1 t) -> Matrix t -> Matrix t
forall t.
Numeric t =>
(Matrix t -> Matrix t)
-> (View1 t -> View1 t) -> Matrix t -> Matrix t
foldMatrix Matrix t -> Matrix t
g View1 t -> View1 t
f Matrix t
m)
foldMatrix Matrix t -> Matrix t
_ View1 t -> View1 t
_ Matrix t
m = Matrix t
m


swapMax :: Int -> Matrix t -> (IndexOf Vector, Matrix t)
swapMax Int
k Matrix t
m
    | Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& IndexOf Vector
jIndexOf Vector -> IndexOf Vector -> Bool
forall a. Ord a => a -> a -> Bool
>IndexOf Vector
0 = (IndexOf Vector
j, Matrix t
m Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Vector I -> Extractor
Pos ([Int] -> Vector I
idxs [Int]
swapped), Extractor
All))
    | Bool
otherwise  = (IndexOf Vector
0,Matrix t
m)
  where
    j :: IndexOf Vector
j = Vector t -> IndexOf Vector
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
maxIndex (Vector t -> IndexOf Vector) -> Vector t -> IndexOf Vector
forall a b. (a -> b) -> a -> b
$ Vector t -> Vector t
forall a. Num a => a -> a
abs (Matrix t -> Matrix t
forall m mt. Transposable m mt => m -> mt
tr Matrix t
m Matrix t -> Int -> Vector t
forall c t. Indexable c t => c -> Int -> t
! Int
k)
    swapped :: [Int]
swapped = Int
IndexOf Vector
jInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int
1..Int
IndexOf Vector
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ Int
0Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int
IndexOf Vector
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

down :: (Matrix t -> Matrix t) -> Matrix t -> Matrix t
down Matrix t -> Matrix t
g Matrix t
a = (Matrix t -> Matrix t)
-> (View1 t -> View1 t) -> Matrix t -> Matrix t
forall t.
Numeric t =>
(Matrix t -> Matrix t)
-> (View1 t -> View1 t) -> Matrix t -> Matrix t
foldMatrix Matrix t -> Matrix t
g View1 t -> View1 t
forall a a c.
(Eq a, Num a, Num c, Product a, Fractional a, Container Vector a,
 Num (Vector a)) =>
(a, Vector a, Vector a, Matrix a) -> (a, Vector a, c, Matrix a)
f Matrix t
a
  where
    f :: (a, Vector a, Vector a, Matrix a) -> (a, Vector a, c, Matrix a)
f (a
e,Vector a
r,Vector a
c,Matrix a
m)
        | a
e a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0    = (a
1, Vector a
r', c
0, Matrix a
m Matrix a -> Matrix a -> Matrix a
forall a. Num a => a -> a -> a
- Vector a -> Vector a -> Matrix a
forall t. Product t => Vector t -> Vector t -> Matrix t
outer Vector a
c Vector a
r')
        | Bool
otherwise = String -> (a, Vector a, c, Matrix a)
forall a. HasCallStack => String -> a
error String
"singular!"
      where
        r' :: Vector a
r' = Vector a
r Vector a -> Vector a -> Vector a
forall a. Fractional a => a -> a -> a
/ a -> Vector a
forall (c :: * -> *) e. Container c e => e -> c e
scalar a
e


-- | generic reference implementation of gaussian elimination
--
-- @a <> gaussElim a b = b@
--
gaussElim_2
  :: (Eq t, Fractional t, Num (Vector t), Numeric t)
  => Matrix t -> Matrix t -> Matrix t

gaussElim_2 :: Matrix t -> Matrix t -> Matrix t
gaussElim_2 Matrix t
a Matrix t
b = Matrix t -> Matrix t
flipudrl Matrix t
r
  where
    flipudrl :: Matrix t -> Matrix t
flipudrl = Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t
flipud (Matrix t -> Matrix t)
-> (Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t
fliprl
    splitColsAt :: Int -> Matrix t -> (Matrix t, Matrix t)
splitColsAt Int
n = (Int -> Matrix t -> Matrix t
forall t. Element t => Int -> Matrix t -> Matrix t
takeColumns Int
n (Matrix t -> Matrix t)
-> (Matrix t -> Matrix t) -> Matrix t -> (Matrix t, Matrix t)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& Int -> Matrix t -> Matrix t
forall t. Element t => Int -> Matrix t -> Matrix t
dropColumns Int
n)
    go :: (Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
go Matrix t -> Matrix t
f Matrix t
x Matrix t
y = Int -> Matrix t -> (Matrix t, Matrix t)
forall t. Element t => Int -> Matrix t -> (Matrix t, Matrix t)
splitColsAt (Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a) ((Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall t.
(Numeric t, Eq t, Num (Vector t), Fractional t) =>
(Matrix t -> Matrix t) -> Matrix t -> Matrix t
down Matrix t -> Matrix t
f (Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall a b. (a -> b) -> a -> b
$ Matrix t
x Matrix t -> Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t -> Matrix t
||| Matrix t
y)
    (Matrix t
a1,Matrix t
b1) = (Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
forall t.
(Numeric t, Eq t, Fractional t, Num (Vector t)) =>
(Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
go ((Int, Matrix t) -> Matrix t
forall a b. (a, b) -> b
snd ((Int, Matrix t) -> Matrix t)
-> (Matrix t -> (Int, Matrix t)) -> Matrix t -> Matrix t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Matrix t -> (Int, Matrix t)
forall t.
(CTrans t, Container Vector t, Num (Vector t)) =>
Int -> Matrix t -> (Int, Matrix t)
swapMax Int
0) Matrix t
a Matrix t
b
    ( Matrix t
_, Matrix t
r) = (Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
forall t.
(Numeric t, Eq t, Fractional t, Num (Vector t)) =>
(Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
go Matrix t -> Matrix t
forall a. a -> a
id (Matrix t -> Matrix t
flipudrl (Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall a b. (a -> b) -> a -> b
$ Matrix t
a1) (Matrix t -> Matrix t
flipudrl (Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall a b. (a -> b) -> a -> b
$ Matrix t
b1)

--------------------------------------------------------------------------------

gaussElim_1
  :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
  => Matrix t -> Matrix t -> Matrix t

gaussElim_1 :: Matrix t -> Matrix t -> Matrix t
gaussElim_1 Matrix t
x Matrix t
y = Int -> Matrix t -> Matrix t
forall t. Element t => Int -> Matrix t -> Matrix t
dropColumns (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x) (Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t
flipud (Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall a b. (a -> b) -> a -> b
$ [Vector t] -> Matrix t
forall t. Element t => [Vector t] -> Matrix t
fromRows [Vector t]
s2)
  where
    rs :: [Vector t]
rs = Matrix t -> [Vector t]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix t -> [Vector t]) -> Matrix t -> [Vector t]
forall a b. (a -> b) -> a -> b
$ Matrix t
x Matrix t -> Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t -> Matrix t
||| Matrix t
y
    s1 :: Matrix t
s1 = [Vector t] -> Matrix t
forall t. Element t => [Vector t] -> Matrix t
fromRows ([Vector t] -> Matrix t) -> [Vector t] -> Matrix t
forall a b. (a -> b) -> a -> b
$ Int -> Int -> [Vector t] -> [Vector t]
forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> Int -> [Vector t] -> [Vector t]
pivotDown (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x) Int
0 [Vector t]
rs      -- interesting
    s2 :: [Vector t]
s2 = Int -> [Vector t] -> [Vector t]
forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> [Vector t] -> [Vector t]
pivotUp (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Matrix t -> [Vector t]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix t -> [Vector t]) -> Matrix t -> [Vector t]
forall a b. (a -> b) -> a -> b
$ Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t
flipud Matrix t
s1)

pivotDown
  :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
  => Int -> Int -> [Vector t] -> [Vector t]
pivotDown :: Int -> Int -> [Vector t] -> [Vector t]
pivotDown Int
t Int
n [Vector t]
xs
    | Int
t Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n    = []
    | Bool
otherwise = Vector t
y Vector t -> [Vector t] -> [Vector t]
forall a. a -> [a] -> [a]
: Int -> Int -> [Vector t] -> [Vector t]
forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> Int -> [Vector t] -> [Vector t]
pivotDown Int
t (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) [Vector t]
ys
  where
    Vector t
y:[Vector t]
ys = (Int, [Vector t]) -> [Vector t]
redu (Int -> [Vector t] -> (Int, [Vector t])
forall a c.
(Ord a, Num a, Indexable c a) =>
Int -> [c] -> (Int, [c])
pivot Int
n [Vector t]
xs)

    pivot :: Int -> [c] -> (Int, [c])
pivot Int
k = (Int -> [c] -> Int
forall a b. a -> b -> a
const Int
k ([c] -> Int) -> ([c] -> [c]) -> [c] -> (Int, [c])
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& [c] -> [c]
forall a. a -> a
id)
            ([c] -> (Int, [c])) -> ([c] -> [c]) -> [c] -> (Int, [c])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (c -> c -> Ordering) -> [c] -> [c]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((a -> a -> Ordering) -> a -> a -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (a -> a -> Ordering) -> (c -> a) -> c -> c -> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (a -> a
forall a. Num a => a -> a
abs(a -> a) -> (c -> a) -> c -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (c -> Int -> a
forall c t. Indexable c t => c -> Int -> t
! Int
k)))

    redu :: (Int, [Vector t]) -> [Vector t]
    redu :: (Int, [Vector t]) -> [Vector t]
redu (Int
k,Vector t
x:[Vector t]
zs)
        | t
p t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
0 = String -> [Vector t]
forall a. HasCallStack => String -> a
error String
"gauss: singular!"  -- FIXME
        | Bool
otherwise = Vector t
u Vector t -> [Vector t] -> [Vector t]
forall a. a -> [a] -> [a]
: (Vector t -> Vector t) -> [Vector t] -> [Vector t]
forall a b. (a -> b) -> [a] -> [b]
map Vector t -> Vector t
f [Vector t]
zs
      where
        p :: t
p = Vector t
xVector t -> Int -> t
forall c t. Indexable c t => c -> Int -> t
!Int
k
        u :: Vector t
u = t -> Vector t -> Vector t
forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (t -> t
forall a. Fractional a => a -> a
recip (Vector t
xVector t -> Int -> t
forall c t. Indexable c t => c -> Int -> t
!Int
k)) Vector t
x
        f :: Vector t -> Vector t
f Vector t
z = Vector t
z Vector t -> Vector t -> Vector t
forall a. Num a => a -> a -> a
- t -> Vector t -> Vector t
forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (Vector t
zVector t -> Int -> t
forall c t. Indexable c t => c -> Int -> t
!Int
k) Vector t
u
    redu (Int
_,[]) = []


pivotUp
  :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
  => Int -> [Vector t] -> [Vector t]
pivotUp :: Int -> [Vector t] -> [Vector t]
pivotUp Int
n [Vector t]
xs
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
1 = []
    | Bool
otherwise = Vector t
y Vector t -> [Vector t] -> [Vector t]
forall a. a -> [a] -> [a]
: Int -> [Vector t] -> [Vector t]
forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> [Vector t] -> [Vector t]
pivotUp (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [Vector t]
ys
  where
    Vector t
y:[Vector t]
ys = (Int, [Vector t]) -> [Vector t]
redu' (Int
n,[Vector t]
xs)

    redu' :: (Int, [Vector t]) -> [Vector t]
    redu' :: (Int, [Vector t]) -> [Vector t]
redu' (Int
k,Vector t
x:[Vector t]
zs) = Vector t
u Vector t -> [Vector t] -> [Vector t]
forall a. a -> [a] -> [a]
: (Vector t -> Vector t) -> [Vector t] -> [Vector t]
forall a b. (a -> b) -> [a] -> [b]
map Vector t -> Vector t
f [Vector t]
zs
      where
        u :: Vector t
u = Vector t
x
        f :: Vector t -> Vector t
f Vector t
z = Vector t
z Vector t -> Vector t -> Vector t
forall a. Num a => a -> a -> a
- t -> Vector t -> Vector t
forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (Vector t
zVector t -> Int -> t
forall c t. Indexable c t => c -> Int -> t
!Int
k) Vector t
u
    redu' (Int
_,[]) = []

--------------------------------------------------------------------------------

gaussElim :: Matrix t -> Matrix t -> Matrix t
gaussElim Matrix t
a Matrix t
b = Int -> Matrix t -> Matrix t
forall t. Element t => Int -> Matrix t -> Matrix t
dropColumns (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a) (Matrix t -> Matrix t) -> Matrix t -> Matrix t
forall a b. (a -> b) -> a -> b
$ (Matrix t, ()) -> Matrix t
forall a b. (a, b) -> a
fst ((Matrix t, ()) -> Matrix t) -> (Matrix t, ()) -> Matrix t
forall a b. (a -> b) -> a -> b
$ (forall s. (Int, Int) -> STMatrix s t -> ST s ())
-> Matrix t -> (Matrix t, ())
forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable forall s. (Int, Int) -> STMatrix s t -> ST s ()
forall t b s.
(Container Vector t, Num (Vector t), Eq t, Fractional t) =>
(Int, b) -> STMatrix s t -> ST s ()
gaussST (Matrix t
a Matrix t -> Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t -> Matrix t
||| Matrix t
b)

gaussST :: (Int, b) -> STMatrix s t -> ST s ()
gaussST (Int
r,b
_) STMatrix s t
x = do
    let n :: Int
n = Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1
        axpy :: STMatrix s t -> t -> Int -> Int -> ST s ()
axpy STMatrix s t
m t
a Int
i Int
j = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (t -> Int -> Int -> ColRange -> RowOper t
forall t. t -> Int -> Int -> ColRange -> RowOper t
AXPY t
a Int
i Int
j ColRange
AllCols)     STMatrix s t
m
        swap :: STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
m Int
i Int
j   = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (Int -> Int -> ColRange -> RowOper t
forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
j ColRange
AllCols)       STMatrix s t
m
        scal :: STMatrix s t -> t -> Int -> ST s ()
scal STMatrix s t
m t
a Int
i   = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (t -> RowRange -> ColRange -> RowOper t
forall t. t -> RowRange -> ColRange -> RowOper t
SCAL t
a (Int -> RowRange
Row Int
i) ColRange
AllCols) STMatrix s t
m
    [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
n] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
        Int
c <- Vector t -> Int
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
maxIndex (Vector t -> Int) -> (Matrix t -> Vector t) -> Matrix t -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector t -> Vector t
forall a. Num a => a -> a
abs (Vector t -> Vector t)
-> (Matrix t -> Vector t) -> Matrix t -> Vector t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten (Matrix t -> Int) -> ST s (Matrix t) -> ST s Int
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> STMatrix s t -> RowRange -> ColRange -> ST s (Matrix t)
forall a t s.
Element a =>
STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
extractMatrix STMatrix s t
x (Int -> RowRange
FromRow Int
i) (Int -> ColRange
Col Int
i)
        STMatrix s t -> Int -> Int -> ST s ()
forall t s.
(Num t, Element t) =>
STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
x Int
i (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
c)
        t
a <- STMatrix s t -> Int -> Int -> ST s t
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
i Int
i
        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (t
a t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
0) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ String -> ST s ()
forall a. HasCallStack => String -> a
error String
"singular!"
        STMatrix s t -> t -> Int -> ST s ()
forall t s.
(Num t, Element t) =>
STMatrix s t -> t -> Int -> ST s ()
scal STMatrix s t
x (t -> t
forall a. Fractional a => a -> a
recip t
a) Int
i
        [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
n] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> do
            t
b <- STMatrix s t -> Int -> Int -> ST s t
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
j Int
i
            STMatrix s t -> t -> Int -> Int -> ST s ()
forall t s.
(Num t, Element t) =>
STMatrix s t -> t -> Int -> Int -> ST s ()
axpy STMatrix s t
x (-t
b) Int
i Int
j
    [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
n,Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1..Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
        [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2..Int
0] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> do
            t
b <- STMatrix s t -> Int -> Int -> ST s t
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
j Int
i
            STMatrix s t -> t -> Int -> Int -> ST s ()
forall t s.
(Num t, Element t) =>
STMatrix s t -> t -> Int -> Int -> ST s ()
axpy STMatrix s t
x (-t
b) Int
i Int
j



luST :: (t -> Bool) -> (Int, b) -> STMatrix s t -> ST s [Int]
luST t -> Bool
ok (Int
r,b
_) STMatrix s t
x = do
    let axpy :: STMatrix s t -> t -> Int -> Int -> ST s ()
axpy STMatrix s t
m t
a Int
i Int
j = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (t -> Int -> Int -> ColRange -> RowOper t
forall t. t -> Int -> Int -> ColRange -> RowOper t
AXPY t
a Int
i Int
j (Int -> ColRange
FromCol (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))) STMatrix s t
m
        swap :: STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
m Int
i Int
j   = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (Int -> Int -> ColRange -> RowOper t
forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
j ColRange
AllCols)           STMatrix s t
m
    STVector s Int
p <- Int -> ST s (STVector s Int)
forall t s. Storable t => Int -> ST s (STVector s t)
newUndefinedVector Int
r
    [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
        Int
k <- Vector t -> Int
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
maxIndex (Vector t -> Int) -> (Matrix t -> Vector t) -> Matrix t -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector t -> Vector t
forall a. Num a => a -> a
abs (Vector t -> Vector t)
-> (Matrix t -> Vector t) -> Matrix t -> Vector t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten (Matrix t -> Int) -> ST s (Matrix t) -> ST s Int
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> STMatrix s t -> RowRange -> ColRange -> ST s (Matrix t)
forall a t s.
Element a =>
STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
extractMatrix STMatrix s t
x (Int -> RowRange
FromRow Int
i) (Int -> ColRange
Col Int
i)
        STVector s Int -> Int -> Int -> ST s ()
forall t s. Storable t => STVector s t -> Int -> t -> ST s ()
writeVector STVector s Int
p Int
i (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i)
        STMatrix s t -> Int -> Int -> ST s ()
forall t s.
(Num t, Element t) =>
STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
x Int
i (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)
        t
a <- STMatrix s t -> Int -> Int -> ST s t
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
i Int
i
        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (t -> Bool
ok t
a) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
            [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> do
                t
b <- (t -> t -> t
forall a. Fractional a => a -> a -> a
/t
a) (t -> t) -> ST s t -> ST s t
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> STMatrix s t -> Int -> Int -> ST s t
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
j Int
i
                STMatrix s t -> t -> Int -> Int -> ST s ()
forall t s.
(Num t, Element t) =>
STMatrix s t -> t -> Int -> Int -> ST s ()
axpy STMatrix s t
x (-t
b) Int
i Int
j
                STMatrix s t -> Int -> Int -> t -> ST s ()
forall t s.
Storable t =>
STMatrix s t -> Int -> Int -> t -> ST s ()
writeMatrix STMatrix s t
x Int
j Int
i t
b
    Vector Int
v <- STVector s Int -> ST s (Vector Int)
forall t s. Storable t => STVector s t -> ST s (Vector t)
unsafeFreezeVector STVector s Int
p
    [Int] -> ST s [Int]
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector Int -> [Int]
forall a. Storable a => Vector a -> [a]
toList Vector Int
v)

{- | Experimental implementation of 'luPacked'
     for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.

>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
(5><5)
 [  1,  1,  2,  3,  4
 ,  5,  7,  7,  8,  9
 , 10, 11, 13, 13, 14
 , 15, 16,  0,  2,  2
 ,  3,  4,  5,  6,  8 ]

>>> let (l,u,p,s) = luFact $ luPacked' m
>>> l
(5><5)
 [  1,  0, 0,  0, 0
 ,  6,  1, 0,  0, 0
 , 12,  7, 1,  0, 0
 ,  7, 10, 7,  1, 0
 ,  8,  2, 6, 11, 1 ]
>>> u
(5><5)
 [ 15, 16,  0,  2,  2
 ,  0, 13,  7, 13, 14
 ,  0,  0, 15,  0, 11
 ,  0,  0,  0, 15, 15
 ,  0,  0,  0,  0,  1 ]

-}
luPacked' :: Matrix t -> LU t
luPacked' Matrix t
x = Matrix t -> [Int] -> LU t
forall t. Matrix t -> [Int] -> LU t
LU Matrix t
m [Int]
p
  where
    (Matrix t
m,[Int]
p) = (forall s. (Int, Int) -> STMatrix s t -> ST s [Int])
-> Matrix t -> (Matrix t, [Int])
forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable ((t -> Bool) -> (Int, Int) -> STMatrix s t -> ST s [Int]
forall t b s.
(Container Vector t, Num (Vector t), Fractional t) =>
(t -> Bool) -> (Int, b) -> STMatrix s t -> ST s [Int]
luST (Double -> t -> Bool
forall t. (Element t, Normed (Vector t)) => Double -> t -> Bool
magnit Double
0)) Matrix t
x

--------------------------------------------------------------------------------

scalS :: t -> Slice s t -> ST s ()
scalS t
a (Slice STMatrix s t
x Int
r0 Int
c0 Int
nr Int
nc) = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (t -> RowRange -> ColRange -> RowOper t
forall t. t -> RowRange -> ColRange -> RowOper t
SCAL t
a (Int -> Int -> RowRange
RowRange Int
r0 (Int
r0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
nrInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) (Int -> Int -> ColRange
ColRange Int
c0 (Int
c0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ncInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))) STMatrix s t
x

view :: STMatrix s a
-> Int -> Int -> ST s (a, Slice s a, Slice s a, Slice s a)
view STMatrix s a
x Int
k Int
r = do
    a
d <- STMatrix s a -> Int -> Int -> ST s a
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s a
x Int
k Int
k
    let rr :: Int
rr = Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k
        o :: Int
o  = if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 then Int
1 else Int
0
        s :: Slice s a
s = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
x (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
rr Int
rr
        u :: Slice s a
u = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
x Int
k     (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
o  Int
rr
        l :: Slice s a
l = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
x (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
k     Int
rr Int
o
    (a, Slice s a, Slice s a, Slice s a)
-> ST s (a, Slice s a, Slice s a, Slice s a)
forall (m :: * -> *) a. Monad m => a -> m a
return (a
d,Slice s a
u,Slice s a
l,Slice s a
s)

withVec :: Int
-> (t -> t -> STVector s t -> ST s a) -> t -> t -> ST s (Vector t)
withVec Int
r t -> t -> STVector s t -> ST s a
f = \t
s t
x -> do
    STVector s t
p <- Int -> ST s (STVector s t)
forall t s. Storable t => Int -> ST s (STVector s t)
newUndefinedVector Int
r
    a
_ <- t -> t -> STVector s t -> ST s a
f t
s t
x STVector s t
p
    Vector t
v <- STVector s t -> ST s (Vector t)
forall t s. Storable t => STVector s t -> ST s (Vector t)
unsafeFreezeVector STVector s t
p
    Vector t -> ST s (Vector t)
forall (m :: * -> *) a. Monad m => a -> m a
return Vector t
v


luPacked'' :: Matrix t -> (Matrix t, [Int])
luPacked'' Matrix t
m = (Matrix t -> Matrix t
forall a. a -> a
id (Matrix t -> Matrix t)
-> (Vector Int -> [Int])
-> (Matrix t, Vector Int)
-> (Matrix t, [Int])
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** Vector Int -> [Int]
forall a. Storable a => Vector a -> [a]
toList) ((forall s. (Int, Int) -> STMatrix s t -> ST s (Vector Int))
-> Matrix t -> (Matrix t, Vector Int)
forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable (Int
-> ((Int, Int) -> STMatrix s t -> STVector s Int -> ST s ())
-> (Int, Int)
-> STMatrix s t
-> ST s (Vector Int)
forall t t t s a.
Storable t =>
Int
-> (t -> t -> STVector s t -> ST s a) -> t -> t -> ST s (Vector t)
withVec (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m) (Int, Int) -> STMatrix s t -> STVector s Int -> ST s ()
forall t b s.
(Container Vector t, Num (Vector t), Normed (Vector t),
 Fractional t) =>
(Int, b) -> STMatrix s t -> STVector s Int -> ST s ()
lu2) Matrix t
m)
  where
    lu2 :: (Int, b) -> STMatrix s t -> STVector s Int -> ST s ()
lu2 (Int
r,b
_) STMatrix s t
x STVector s Int
p = do
        [Int] -> (Int -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
k -> do
            STMatrix s t -> STVector s Int -> Int -> ST s ()
forall t s.
(Container Vector t, Num (Vector t), Num t) =>
STMatrix s t -> STVector s Int -> Int -> ST s ()
pivot STMatrix s t
x STVector s Int
p Int
k
            (t
d,Slice s t
u,Slice s t
l,Slice s t
s) <- STMatrix s t
-> Int -> Int -> ST s (t, Slice s t, Slice s t, Slice s t)
forall a s.
Storable a =>
STMatrix s a
-> Int -> Int -> ST s (a, Slice s a, Slice s a, Slice s a)
view STMatrix s t
x Int
k Int
r
            Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double -> t -> Bool
forall t. (Element t, Normed (Vector t)) => Double -> t -> Bool
magnit Double
0 t
d) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                t -> Slice s t -> ST s ()
forall t s. (Num t, Element t) => t -> Slice s t -> ST s ()
scalS (t -> t
forall a. Fractional a => a -> a
recip t
d) Slice s t
l
                t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
forall t s.
Element t =>
t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
gemmm t
1 Slice s t
s (-t
1) Slice s t
l Slice s t
u

    pivot :: STMatrix s t -> STVector s Int -> Int -> ST s ()
pivot STMatrix s t
x STVector s Int
p Int
k = do
        Int
j <- Vector t -> Int
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
maxIndex (Vector t -> Int) -> (Matrix t -> Vector t) -> Matrix t -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector t -> Vector t
forall a. Num a => a -> a
abs (Vector t -> Vector t)
-> (Matrix t -> Vector t) -> Matrix t -> Vector t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> Vector t
forall t. Element t => Matrix t -> Vector t
flatten (Matrix t -> Int) -> ST s (Matrix t) -> ST s Int
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> STMatrix s t -> RowRange -> ColRange -> ST s (Matrix t)
forall a t s.
Element a =>
STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
extractMatrix STMatrix s t
x (Int -> RowRange
FromRow Int
k) (Int -> ColRange
Col Int
k)
        STVector s Int -> Int -> Int -> ST s ()
forall t s. Storable t => STVector s t -> Int -> t -> ST s ()
writeVector STVector s Int
p Int
k (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)
        Int -> Int -> ST s ()
swap Int
k (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j)
      where
        swap :: Int -> Int -> ST s ()
swap Int
i Int
j = RowOper t -> STMatrix s t -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (Int -> Int -> ColRange -> RowOper t
forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
j ColRange
AllCols) STMatrix s t
x

--------------------------------------------------------------------------------

rowRange :: Matrix t -> [Int]
rowRange Matrix t
m = [Int
0..Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

at :: Int -> Extractor
at Int
k = Vector I -> Extractor
Pos ([Int] -> Vector I
idxs[Int
k])

backSust' :: Matrix t -> Matrix t -> Matrix t
backSust' Matrix t
lup Matrix t
rhs = (Matrix t -> (Matrix t, Matrix t, Matrix t) -> Matrix t)
-> Matrix t -> [(Matrix t, Matrix t, Matrix t)] -> Matrix t
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Matrix t -> (Matrix t, Matrix t, Matrix t) -> Matrix t
forall t (a :: * -> *).
(Container Vector t, Fractional t, Num (Vector t),
 Mul a Matrix Matrix, Product t) =>
Matrix t -> (Matrix t, a t, Matrix t) -> Matrix t
f (Matrix t
rhsMatrix t -> [Int] -> Matrix t
forall t. Element t => Matrix t -> [Int] -> Matrix t
?[]) ([(Matrix t, Matrix t, Matrix t)]
-> [(Matrix t, Matrix t, Matrix t)]
forall a. [a] -> [a]
reverse [(Matrix t, Matrix t, Matrix t)]
ls)
  where
    ls :: [(Matrix t, Matrix t, Matrix t)]
ls  = [ (Int -> Matrix t
d Int
k , Int -> Matrix t
u Int
k , Int -> Matrix t
b Int
k) | Int
k <- Matrix t -> [Int]
forall t. Matrix t -> [Int]
rowRange Matrix t
lup ]
      where
        d :: Int -> Matrix t
d Int
k = Matrix t
lup Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
at Int
k, Int -> Extractor
at Int
k)
        u :: Int -> Matrix t
u Int
k = Matrix t
lup Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
at Int
k, Int -> Extractor
Drop (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
        b :: Int -> Matrix t
b Int
k = Matrix t
rhs Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
at Int
k, Extractor
All)

    f :: Matrix t -> (Matrix t, a t, Matrix t) -> Matrix t
f Matrix t
x (Matrix t
d,a t
u,Matrix t
b) = (Matrix t
b Matrix t -> Matrix t -> Matrix t
forall a. Num a => a -> a -> a
- a t
ua t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<>Matrix t
x) Matrix t -> Matrix t -> Matrix t
forall a. Fractional a => a -> a -> a
/ Matrix t
d
                       Matrix t -> Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t -> Matrix t
===
                        Matrix t
x


forwSust' :: Matrix t -> Matrix t -> Matrix t
forwSust' Matrix t
lup Matrix t
rhs = (Matrix t -> (Matrix t, Matrix t) -> Matrix t)
-> Matrix t -> [(Matrix t, Matrix t)] -> Matrix t
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Matrix t -> (Matrix t, Matrix t) -> Matrix t
forall t (a :: * -> *).
(Container Vector t, Num (Vector t), Mul a Matrix Matrix,
 Product t) =>
Matrix t -> (a t, Matrix t) -> Matrix t
f (Matrix t
rhsMatrix t -> [Int] -> Matrix t
forall t. Element t => Matrix t -> [Int] -> Matrix t
?[]) [(Matrix t, Matrix t)]
ls
  where
    ls :: [(Matrix t, Matrix t)]
ls  = [ (Int -> Matrix t
l Int
k , Int -> Matrix t
b Int
k) | Int
k <- Matrix t -> [Int]
forall t. Matrix t -> [Int]
rowRange Matrix t
lup ]
      where
        l :: Int -> Matrix t
l Int
k = Matrix t
lup Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
at Int
k, Int -> Extractor
Take Int
k)
        b :: Int -> Matrix t
b Int
k = Matrix t
rhs Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
at Int
k, Extractor
All)

    f :: Matrix t -> (a t, Matrix t) -> Matrix t
f Matrix t
x (a t
l,Matrix t
b) =     Matrix t
x
                   Matrix t -> Matrix t -> Matrix t
forall t. Element t => Matrix t -> Matrix t -> Matrix t
===
                (Matrix t
b Matrix t -> Matrix t -> Matrix t
forall a. Num a => a -> a -> a
- a t
la t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<>Matrix t
x)


luSolve'' :: LU t -> Matrix t -> Matrix t
luSolve'' (LU Matrix t
lup [Int]
p) Matrix t
b = Matrix t -> Matrix t -> Matrix t
forall t.
(Container Vector t, Fractional t, Product t, Num (Vector t)) =>
Matrix t -> Matrix t -> Matrix t
backSust' Matrix t
lup (Matrix t -> Matrix t -> Matrix t
forall t.
(Container Vector t, Product t, Num (Vector t)) =>
Matrix t -> Matrix t -> Matrix t
forwSust' Matrix t
lup Matrix t
pb)
  where
    pb :: Matrix t
pb = Matrix t
b Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Vector I -> Extractor
Pos ([Int] -> Vector I
fixPerm' [Int]
p), Extractor
All)

--------------------------------------------------------------------------------

forwSust :: Matrix t -> Matrix t -> Matrix t
forwSust Matrix t
lup Matrix t
rhs = (Matrix t, ()) -> Matrix t
forall a b. (a, b) -> a
fst ((Matrix t, ()) -> Matrix t) -> (Matrix t, ()) -> Matrix t
forall a b. (a -> b) -> a -> b
$ (forall s. (Int, Int) -> STMatrix s t -> ST s ())
-> Matrix t -> (Matrix t, ())
forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable forall s. (Int, Int) -> STMatrix s t -> ST s ()
f Matrix t
rhs
  where
    f :: (Int, Int) -> STMatrix s t -> ST s ()
f (Int
r,Int
c) STMatrix s t
x = do
        STMatrix s t
l <- Matrix t -> ST s (STMatrix s t)
forall t s. Storable t => Matrix t -> ST s (STMatrix s t)
unsafeThawMatrix Matrix t
lup
        let go :: Int -> ST s ()
go Int
k = t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
forall t s.
Element t =>
t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
gemmm t
1 (STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s t
x Int
k Int
0 Int
1 Int
c) (-t
1) (STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s t
l Int
k Int
0 Int
1 Int
k) (STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s t
x Int
0 Int
0 Int
k Int
c)
        (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Int -> ST s ()
go [Int
0..Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]


backSust :: Matrix a -> Matrix a -> Matrix a
backSust Matrix a
lup Matrix a
rhs = (Matrix a, ()) -> Matrix a
forall a b. (a, b) -> a
fst ((Matrix a, ()) -> Matrix a) -> (Matrix a, ()) -> Matrix a
forall a b. (a -> b) -> a -> b
$ (forall s. (Int, Int) -> STMatrix s a -> ST s ())
-> Matrix a -> (Matrix a, ())
forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable forall s. (Int, Int) -> STMatrix s a -> ST s ()
f Matrix a
rhs
  where
    f :: (Int, Int) -> STMatrix s a -> ST s ()
f (Int
r,Int
c) STMatrix s a
m = do
        STMatrix s a
l <- Matrix a -> ST s (STMatrix s a)
forall t s. Storable t => Matrix t -> ST s (STMatrix s t)
unsafeThawMatrix Matrix a
lup
        let d :: Int -> a
d Int
k = a -> a
forall a. Fractional a => a -> a
recip (Matrix a
lup Matrix a -> IndexOf Matrix -> a
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
k,Int
k))
            u :: Int -> Slice s a
u Int
k = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
l Int
k (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
1 (Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)
            b :: Int -> Slice s a
b Int
k = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
m Int
k Int
0 Int
1 Int
c
            x :: Int -> Slice s a
x Int
k = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
m (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
0 (Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k) Int
c
            scal :: Int -> ST s ()
scal Int
k = RowOper a -> STMatrix s a -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (a -> RowRange -> ColRange -> RowOper a
forall t. t -> RowRange -> ColRange -> RowOper t
SCAL (Int -> a
d Int
k) (Int -> RowRange
Row Int
k) ColRange
AllCols) STMatrix s a
m

            go :: Int -> ST s ()
go Int
k = a -> Slice s a -> a -> Slice s a -> Slice s a -> ST s ()
forall t s.
Element t =>
t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
gemmm a
1 (Int -> Slice s a
b Int
k) (-a
1) (Int -> Slice s a
u Int
k) (Int -> Slice s a
x Int
k) ST s () -> ST s () -> ST s ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Int -> ST s ()
scal Int
k
        (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Int -> ST s ()
go [Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2..Int
0]


{- | Experimental implementation of 'luSolve' for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.

>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
(2><2)
 [ 1, 2
 , 3, 5 ]
>>> b
(2><3)
 [ 5, 1, 3
 , 8, 6, 3 ]

>>> luSolve' (luPacked' a) b
(2><3)
 [ 4,  7, 4
 , 7, 10, 6 ]

-}
luSolve' :: LU t -> Matrix t -> Matrix t
luSolve' (LU Matrix t
lup [Int]
p) Matrix t
b = Matrix t -> Matrix t -> Matrix t
forall a.
(Fractional a, Container Vector a) =>
Matrix a -> Matrix a -> Matrix a
backSust Matrix t
lup (Matrix t -> Matrix t -> Matrix t
forall t. (Element t, Num t) => Matrix t -> Matrix t -> Matrix t
forwSust Matrix t
lup Matrix t
pb)
  where
    pb :: Matrix t
pb = Matrix t
b Matrix t -> (Extractor, Extractor) -> Matrix t
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Vector I -> Extractor
Pos ([Int] -> Vector I
fixPerm' [Int]
p), Extractor
All)


--------------------------------------------------------------------------------

data MatrixView t b
    = Elem t
    | Block b b b b
  deriving Int -> MatrixView t b -> String -> String
[MatrixView t b] -> String -> String
MatrixView t b -> String
(Int -> MatrixView t b -> String -> String)
-> (MatrixView t b -> String)
-> ([MatrixView t b] -> String -> String)
-> Show (MatrixView t b)
forall a.
(Int -> a -> String -> String)
-> (a -> String) -> ([a] -> String -> String) -> Show a
forall t b.
(Show t, Show b) =>
Int -> MatrixView t b -> String -> String
forall t b.
(Show t, Show b) =>
[MatrixView t b] -> String -> String
forall t b. (Show t, Show b) => MatrixView t b -> String
showList :: [MatrixView t b] -> String -> String
$cshowList :: forall t b.
(Show t, Show b) =>
[MatrixView t b] -> String -> String
show :: MatrixView t b -> String
$cshow :: forall t b. (Show t, Show b) => MatrixView t b -> String
showsPrec :: Int -> MatrixView t b -> String -> String
$cshowsPrec :: forall t b.
(Show t, Show b) =>
Int -> MatrixView t b -> String -> String
Show


viewBlock' :: Int -> Int -> Matrix t -> MatrixView t (Matrix t)
viewBlock' Int
r Int
c Matrix t
m
    | (Int
rt,Int
ct) (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
1,Int
1) = t -> MatrixView t (Matrix t)
forall t b. t -> MatrixView t b
Elem (Matrix t -> Int -> Int -> t
forall t. Storable t => Matrix t -> Int -> Int -> t
atM' Matrix t
m Int
0 Int
0)
    | Bool
otherwise        = Matrix t
-> Matrix t -> Matrix t -> Matrix t -> MatrixView t (Matrix t)
forall t b. b -> b -> b -> b -> MatrixView t b
Block Matrix t
m11 Matrix t
m12 Matrix t
m21 Matrix t
m22
  where
    (Int
rt,Int
ct) = Matrix t -> IndexOf Matrix
forall (c :: * -> *) e. Container c e => c e -> IndexOf c
size Matrix t
m
    m11 :: Matrix t
m11 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
subm (Int
0,Int
0) (Int
r,Int
c)       Matrix t
m
    m12 :: Matrix t
m12 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
subm (Int
0,Int
c) (Int
r,Int
ctInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
c)    Matrix t
m
    m21 :: Matrix t
m21 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
subm (Int
r,Int
0) (Int
rtInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r,Int
c)    Matrix t
m
    m22 :: Matrix t
m22 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
subm (Int
r,Int
c) (Int
rtInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r,Int
ctInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
c) Matrix t
m
    subm :: (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
subm = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix

viewBlock :: Matrix t -> MatrixView t (Matrix t)
viewBlock Matrix t
m = Int -> Int -> Matrix t -> MatrixView t (Matrix t)
forall t.
(Num t, Container Vector t) =>
Int -> Int -> Matrix t -> MatrixView t (Matrix t)
viewBlock' Int
n Int
n Matrix t
m
  where
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2

invershur :: Matrix t -> Matrix t
invershur (Matrix t -> MatrixView t (Matrix t)
forall t.
(Num t, Container Vector t) =>
Matrix t -> MatrixView t (Matrix t)
viewBlock -> Block Matrix t
a Matrix t
b Matrix t
c Matrix t
d) = [[Matrix t]] -> Matrix t
forall t. Element t => [[Matrix t]] -> Matrix t
fromBlocks [[Matrix t
a',Matrix t
b'],[Matrix t
c',Matrix t
d']]
  where
    r1 :: Matrix t
r1 = Matrix t -> Matrix t
invershur Matrix t
a
    r2 :: Matrix t
r2 = Matrix t
c Matrix t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix t
r1
    r3 :: Matrix t
r3 = Matrix t
r1 Matrix t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix t
b
    r4 :: Matrix t
r4 = Matrix t
c Matrix t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix t
r3
    r5 :: Matrix t
r5 = Matrix t
r4Matrix t -> Matrix t -> Matrix t
forall a. Num a => a -> a -> a
-Matrix t
d
    r6 :: Matrix t
r6 = Matrix t -> Matrix t
invershur Matrix t
r5
    b' :: Matrix t
b' = Matrix t
r3 Matrix t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix t
r6
    c' :: Matrix t
c' = Matrix t
r6 Matrix t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix t
r2
    r7 :: Matrix t
r7 = Matrix t
r3 Matrix t -> Matrix t -> Matrix t
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix t
c'
    a' :: Matrix t
a' = Matrix t
r1Matrix t -> Matrix t -> Matrix t
forall a. Num a => a -> a -> a
-Matrix t
r7
    d' :: Matrix t
d' = -Matrix t
r6

invershur Matrix t
x = Matrix t -> Matrix t
forall a. Fractional a => a -> a
recip Matrix t
x

--------------------------------------------------------------------------------

instance Testable (Matrix I) where
   checkT :: Matrix I -> (Bool, IO ())
checkT Matrix I
_ = (Bool, IO ())
test

test :: (Bool, IO())
test :: (Bool, IO ())
test = ([Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Bool]
ok, () -> IO ()
forall (m :: * -> *) a. Monad m => a -> m a
return ())
  where
    m :: Matrix I
m  = (Int
3Int -> Int -> [I] -> Matrix I
forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
4) [I
1..I
12] :: Matrix I
    r :: Matrix I
r  = (Int
2Int -> Int -> [I] -> Matrix I
forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
3) [I
1,I
2,I
3,I
4,I
3,I
2]
    c :: Matrix I
c  = (Int
3Int -> Int -> [I] -> Matrix I
forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
2) [I
0,I
4,I
4,I
1,I
2,I
3]
    p :: Matrix I
p  = (Int
9Int -> Int -> [I] -> Matrix I
forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
10) [I
0..I
89] :: Matrix I
    ep :: Matrix I
ep = (Int
2Int -> Int -> [I] -> Matrix I
forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
3) [I
10,I
24,I
32,I
44,I
31,I
23]
    md :: Matrix Double
md = Matrix I -> Matrix Double
forall (c :: * -> *) e. Container c e => c I -> c e
fromInt Matrix I
m      :: Matrix Double
    ok :: [Bool]
ok = [ Matrix I -> Matrix I
forall m mt. Transposable m mt => m -> mt
tr Matrix I
m Matrix I -> Matrix I -> Matrix I
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix I
m Matrix I -> Matrix I -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix Double -> Matrix I
forall (c :: * -> *) e. Container c e => c e -> c I
toInt (Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
tr Matrix Double
md Matrix Double -> Matrix Double -> Matrix Double
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix Double
md)
         , Matrix I
m Matrix I -> Matrix I -> Matrix I
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix I -> Matrix I
forall m mt. Transposable m mt => m -> mt
tr Matrix I
m Matrix I -> Matrix I -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix Double -> Matrix I
forall (c :: * -> *) e. Container c e => c e -> c I
toInt (Matrix Double
md Matrix Double -> Matrix Double -> Matrix Double
forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
tr Matrix Double
md)
         , Matrix I
m Matrix I -> (Extractor, Extractor) -> Matrix I
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Take Int
2, Int -> Extractor
Take Int
3) Matrix I -> Matrix I -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix I -> Matrix I -> Matrix I -> Matrix I
forall t. Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
remap (Vector I -> Matrix I
forall a. Storable a => Vector a -> Matrix a
asColumn (Int -> Vector I
range Int
2)) (Vector I -> Matrix I
forall a. Storable a => Vector a -> Matrix a
asRow (Int -> Vector I
range Int
3)) Matrix I
m
         , Matrix I -> Matrix I -> Matrix I -> Matrix I
forall t. Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
remap Matrix I
r (Matrix I -> Matrix I
forall m mt. Transposable m mt => m -> mt
tr Matrix I
c) Matrix I
p Matrix I -> Matrix I -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix I
ep
         , Matrix I -> Matrix I
forall m mt. Transposable m mt => m -> mt
tr Matrix I
p Matrix I -> (Extractor, Extractor) -> Matrix I
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Vector I -> Extractor
PosCyc ([Int] -> Vector I
idxs[-Int
5,Int
13]), Vector I -> Extractor
Pos ([Int] -> Vector I
idxs[Int
3,Int
7,Int
1])) Matrix I -> Matrix I -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
2Int -> Int -> [I] -> Matrix I
forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
3) [I
35,I
75,I
15,I
33,I
73,I
13]
         ]