{-# 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 :: Complex Double
iC = Double
0forall 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 = 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 = forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
c forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 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 = forall a. (Num a, Element a) => Vector a -> Matrix a
diag forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = 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 = 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 = 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 ||| :: forall t. Element t => Matrix t -> Matrix t -> Matrix t
||| Matrix t
b = 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
(¦) = 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 === :: forall t. Element t => Matrix t -> Matrix t -> Matrix t
=== Matrix t
b = 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
(——) = 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 = forall a. Storable a => Vector a -> Matrix a
asRow forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall a. Storable a => Vector a -> Matrix a
asColumn forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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
? :: forall t. Element t => Matrix t -> [Int] -> Matrix t
(?) = forall a b c. (a -> b -> c) -> b -> a -> c
flip 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
¿ :: forall t. Element t => Matrix t -> [Int] -> Matrix t
(¿)= forall a b c. (a -> b -> c) -> b -> a -> c
flip 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 :: forall t. Product t => Vector t -> Vector t -> Vector t
cross Vector t
x Vector t
y | forall t. Storable t => Vector t -> Int
dim Vector t
x forall a. Eq a => a -> a -> Bool
== Int
3 Bool -> Bool -> Bool
&& forall t. Storable t => Vector t -> Int
dim Vector t
y forall a. Eq a => a -> a -> Bool
== Int
3 = forall a. Storable a => [a] -> Vector a
fromList [t
z1,t
z2,t
z3]
          | Bool
otherwise = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"the cross product requires 3-element vectors (sizes given: "
                                forall a. [a] -> [a] -> [a]
++forall a. Show a => a -> String
show (forall t. Storable t => Vector t -> Int
dim Vector t
x)forall a. [a] -> [a] -> [a]
++String
" and "forall a. [a] -> [a] -> [a]
++forall a. Show a => a -> String
show (forall t. Storable t => Vector t -> Int
dim Vector t
y)forall a. [a] -> [a] -> [a]
++String
")"
  where
    [t
x1,t
x2,t
x3] = forall a. Storable a => Vector a -> [a]
toList Vector t
x
    [t
y1,t
y2,t
y3] = forall a. Storable a => Vector a -> [a]
toList Vector t
y
    z1 :: t
z1 = t
x2forall a. Num a => a -> a -> a
*t
y3forall a. Num a => a -> a -> a
-t
x3forall a. Num a => a -> a -> a
*t
y2
    z2 :: t
z2 = t
x3forall a. Num a => a -> a -> a
*t
y1forall a. Num a => a -> a -> a
-t
x1forall a. Num a => a -> a -> a
*t
y3
    z3 :: t
z3 = t
x1forall a. Num a => a -> a -> a
*t
y2forall a. Num a => a -> a -> a
-t
x2forall 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 = 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 = forall (c :: * -> *) e. Container c e => c e -> e
sumElements (forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
step (forall a. Num a => a -> a
abs Vector Double
v forall a. Num a => a -> a -> a
- forall (c :: * -> *) e. Container c e => e -> c e
scalar (Double
epsforall a. Num a => a -> a -> a
*forall e. Product e => Vector e -> RealOf e
normInf Vector Double
v)))
    norm_1 :: Vector Double -> Double
norm_1 = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1
    norm_2 :: Vector Double -> Double
norm_2 = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2
    norm_Inf :: Vector Double -> Double
norm_Inf = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity

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

instance Normed (Matrix C)
  where
    norm_0 :: Matrix (Complex Double) -> Double
norm_0 = forall a. Normed a => a -> Double
norm_0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Vector t
flatten
    norm_1 :: Matrix (Complex Double) -> Double
norm_1 = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1
    norm_2 :: Matrix (Complex Double) -> Double
norm_2 = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2
    norm_Inf :: Matrix (Complex Double) -> Double
norm_Inf = 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 = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (c :: * -> *) e. Container c e => c e -> e
sumElements forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
step forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs
    norm_1 :: Vector I -> Double
norm_1 = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall e. Product e => Vector e -> RealOf e
norm1
    norm_2 :: Vector I -> Double
norm_2 Vector I
v = forall a. Floating a => a -> a
sqrt forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall t. Numeric t => Vector t -> Vector t -> t
dot Vector I
v Vector I
v
    norm_Inf :: Vector I -> Double
norm_Inf = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall e. Product e => Vector e -> RealOf e
normInf

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

instance Normed (Vector Float)
  where
    norm_0 :: Vector Float -> Double
norm_0 = forall a. Normed a => a -> Double
norm_0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_1 :: Vector Float -> Double
norm_1 = forall a. Normed a => a -> Double
norm_1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_2 :: Vector Float -> Double
norm_2 = forall a. Normed a => a -> Double
norm_2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_Inf :: Vector Float -> Double
norm_Inf = forall a. Normed a => a -> Double
norm_Inf forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall a. Normed a => a -> Double
norm_0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_1 :: Vector (Complex Float) -> Double
norm_1 = forall a. Normed a => a -> Double
norm_1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_2 :: Vector (Complex Float) -> Double
norm_2 = forall a. Normed a => a -> Double
norm_2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    norm_Inf :: Vector (Complex Float) -> Double
norm_Inf = forall a. Normed a => a -> Double
norm_Inf forall b c a. (b -> c) -> (a -> b) -> a -> 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 :: forall t. (Normed (Vector t), Element t) => Matrix t -> Double
norm_Frob = forall a. Normed a => a -> Double
norm_2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall t. Field t => Matrix t -> Double
norm_nuclear = forall (c :: * -> *) e. Container c e => c e -> e
sumElements forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall t. (Element t, Normed (Vector t)) => Double -> t -> Bool
magnit Double
e t
x = forall a. Normed a => a -> Double
norm_1 (forall a. Storable a => [a] -> Vector a
fromList [t
x]) 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 :: forall t.
(Normed (Vector t), Num (Vector t), Field t) =>
Vector t -> Vector t
normalize Vector t
v = Vector t
v forall a. Fractional a => a -> a -> a
/ forall t (c :: * -> *).
(Convert t, Complexable c) =>
c (RealOf t) -> c t
real (forall (c :: * -> *) e. Container c e => e -> c e
scalar (forall a. Normed a => a -> Double
norm_2 Vector t
v))


-- | trans . inv
mt :: Matrix Double -> Matrix Double
mt :: Matrix Double -> Matrix Double
mt = forall t. Matrix t -> Matrix t
trans forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 :: forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size = forall (c :: * -> *) t. Container c t => c t -> 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 (Complex Double) -> Int -> Complex Double
(!) = 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 = forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector (Int
jforall a. Num a => a -> a -> a
*Int
c) Int
c (forall t. Element t => Matrix t -> Vector t
flatten Matrix t
m)
      where
        c :: Int
c = 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 forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
oy forall a. Num a => a -> a -> a
+ Vector Double
ox forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
y2 forall a. Num a => a -> a -> a
- Matrix Double
2forall a. Num a => a -> a -> a
* Matrix Double
x forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> forall t. Matrix t -> Matrix t
trans Matrix Double
y
               | Bool
otherwise = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"pairwiseD2 with different number of columns: "
                                   forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Matrix Double
x) forall a. [a] -> [a] -> [a]
++ String
", " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Matrix Double
y)
  where
    ox :: Vector Double
ox = forall {e} {d} {c :: * -> *}. (Konst e d c, Num e) => d -> c e
one (forall t. Matrix t -> Int
rows Matrix Double
x)
    oy :: Vector Double
oy = forall {e} {d} {c :: * -> *}. (Konst e d c, Num e) => d -> c e
one (forall t. Matrix t -> Int
rows Matrix Double
y)
    oc :: Vector Double
oc = forall {e} {d} {c :: * -> *}. (Konst e d c, Num e) => d -> c e
one (forall t. Matrix t -> Int
cols Matrix Double
x)
    one :: d -> c e
one d
k = forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst e
1 d
k
    x2 :: Vector Double
x2 = Matrix Double
x forall a. Num a => a -> a -> a
* Matrix Double
x 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 forall a. Num a => a -> a -> a
* Matrix Double
y forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<> Vector Double
oc
    ok :: Bool
ok = forall t. Matrix t -> Int
cols Matrix Double
x forall a. Eq a => a -> a -> Bool
== 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' forall a. Num a => a -> a -> a
* Matrix Double
b'
  where
    a' :: Matrix Double
a' = forall t. Product t => Matrix t -> Matrix t -> Matrix t
kronecker Matrix Double
a (Int -> Int -> Matrix Double
ones Int
1 (forall t. Matrix t -> Int
cols Matrix Double
b))
    b' :: Matrix Double
b' = forall t. Product t => Matrix t -> Matrix t -> Matrix t
kronecker (Int -> Int -> Matrix Double
ones Int
1 (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 = forall a. [a] -> a
last forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> [Vector t]
toColumns forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall a. [a] -> a
last forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> [Vector t]
toColumns forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Field t => Herm t -> (Vector Double, Matrix t)
eigSH

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

infixl 0 ~!~
Bool
c ~!~ :: Bool -> String -> f ()
~!~ String
msg = forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
c (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) = forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
sep 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 = forall r. PrintfType r => String -> r
printf String
"%.0f" t
x

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

approxInt :: Matrix Double -> Maybe (Matrix Double)
approxInt Matrix Double
m
    | forall a. Normed a => a -> Double
norm_Inf (Vector Double
v forall a. Num a => a -> a -> a
- Vector Double
vi) forall a. Ord a => a -> a -> Bool
< Double
2forall a. Num a => a -> a -> a
*forall x. RealFloat x => x
peps forall a. Num a => a -> a -> a
* forall a. Normed a => a -> Double
norm_Inf Vector Double
v = forall a. a -> Maybe a
Just (forall t. Storable t => Int -> Vector t -> Matrix t
reshape (forall t. Matrix t -> Int
cols Matrix Double
m) Vector Double
vi)
    | Bool
otherwise = forall a. Maybe a
Nothing
  where
    v :: Vector Double
v = 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 forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> String -> String -> Int -> Matrix Double -> String
formatSparse String
"." (forall a. Int -> a -> [a]
replicate Int
n Char
' ') String
"  " Int
n

dispBlanks :: Int -> Matrix Double -> IO ()
dispBlanks Int
n = String -> IO ()
putStr 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) = forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Matrix t
m
    (Int
r1,Int
r2,Int
r3)
        | Int
rm forall a. Ord a => a -> a -> Bool
<= Int
maxr = (Int
rm,Int
0,Int
0)
        | Bool
otherwise  = (Int
maxrforall a. Num a => a -> a -> a
-Int
3,Int
rmforall a. Num a => a -> a -> a
-Int
maxrforall a. Num a => a -> a -> a
+Int
1,Int
2)
    (Int
c1,Int
c2,Int
c3)
        | Int
cm forall a. Ord a => a -> a -> Bool
<= Int
maxc = (Int
cm,Int
0,Int
0)
        | Bool
otherwise  = (Int
maxcforall a. Num a => a -> a -> a
-Int
3,Int
cmforall a. Num a => a -> a -> a
-Int
maxcforall 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]] = 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 = 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 forall a. Ord a => a -> a -> Bool
> Int
maxc = forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
"|" t -> String
fmt Matrix t
auxm
        | Bool
otherwise = forall t.
Element t =>
String -> (t -> String) -> Matrix t -> String
format String
sep t -> String
fmt Matrix t
auxm
    auxm3 :: [String]
auxm3
        | Int
cm forall a. Ord a => a -> a -> Bool
> Int
maxc = forall a b. (a -> b) -> [a] -> [b]
map ([String] -> String
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall a. [a] -> [[a]] -> [a]
intercalate String
sep (forall a. Int -> [a] -> [a]
take (Int
maxcforall a. Num a => a -> a -> a
-Int
3) [String]
items) forall a. [a] -> [a] -> [a]
++ String
"  .. " forall a. [a] -> [a] -> [a]
++
              forall a. [a] -> [[a]] -> [a]
intercalate String
sep (forall a. Int -> [a] -> [a]
drop (Int
maxcforall a. Num a => a -> a -> a
-Int
3) [String]
items)
    auxm4 :: String
auxm4
        | Int
rm forall a. Ord a => a -> a -> Bool
> Int
maxr = [String] -> String
unlines (forall a. Int -> [a] -> [a]
take (Int
maxrforall a. Num a => a -> a -> a
-Int
3) [String]
auxm3 forall a. [a] -> [a] -> [a]
++ String
vsep forall a. a -> [a] -> [a]
: forall a. Int -> [a] -> [a]
drop (Int
maxrforall a. Num a => a -> a -> a
-Int
3) [String]
auxm3)
        | Bool
otherwise = [String] -> String
unlines [String]
auxm3
    vsep :: String
vsep = forall a b. (a -> b) -> [a] -> [b]
map Char -> Char
g (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 =
    forall r. PrintfType r => String -> r
printf String
"%dx%d\n%s" (forall t. Matrix t -> Int
rows Matrix Double
m) (forall t. Matrix t -> Int
cols Matrix Double
m) (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 = forall r. PrintfType r => String -> r
printf (String
"%."forall a. [a] -> [a] -> [a]
++forall a. Show a => a -> String
show Int
dec 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 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 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 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 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 forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? ([Extractor]
er forall a. [a] -> Int -> a
!! Int
i, [Extractor]
ec 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
rforall a. Num a => a -> a -> a
-Int
1), Int -> Int -> Int -> Extractor
Range Int
r Int
1 (Int
rforall a. Num a => a -> a -> a
+Int
nrforall a. Num a => a -> a -> a
-Int
1), Int -> Extractor
Drop (Int
nrforall a. Num a => a -> a -> a
+Int
r) ]
    ec :: [Extractor]
ec = [ Int -> Int -> Int -> Extractor
Range Int
0 Int
1 (Int
cforall a. Num a => a -> a -> a
-Int
1), Int -> Int -> Int -> Extractor
Range Int
c Int
1 (Int
cforall a. Num a => a -> a -> a
+Int
ncforall a. Num a => a -> a -> a
-Int
1), Int -> Extractor
Drop (Int
ncforall a. Num a => a -> a -> a
+Int
c) ]

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

unView1 :: Numeric t => View1 t -> Matrix t
unView1 :: forall t. Numeric t => View1 t -> Matrix t
unView1 (t
e,Vector t
r,Vector t
c,Matrix t
m) = forall t. Element t => [[Matrix t]] -> Matrix t
fromBlocks [[forall (c :: * -> *) e. Container c e => e -> c e
scalar t
e, forall a. Storable a => Vector a -> Matrix a
asRow Vector t
r],[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 :: 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 ( (View1 t -> View1 t
f forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Numeric t => Matrix t -> Maybe (View1 t)
view1 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)) = forall t. Numeric t => View1 t -> Matrix t
unView1 (t
e, Vector t
r, Vector t
c, 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 e -> (IndexOf Vector, Matrix e)
swapMax Int
k Matrix e
m
    | forall t. Matrix t -> Int
rows Matrix e
m forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& IndexOf Vector
jforall a. Ord a => a -> a -> Bool
>IndexOf Vector
0 = (IndexOf Vector
j, Matrix e
m 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 e
m)
  where
    j :: IndexOf Vector
j = forall (c :: * -> *) t. Container c t => c t -> IndexOf c
maxIndex forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
abs (forall m mt. Transposable m mt => m -> mt
tr Matrix e
m forall c t. Indexable c t => c -> Int -> t
! Int
k)
    swapped :: [Int]
swapped = IndexOf Vector
jforall a. a -> [a] -> [a]
:[Int
1..IndexOf Vector
jforall a. Num a => a -> a -> a
-Int
1] forall a. [a] -> [a] -> [a]
++ Int
0forall a. a -> [a] -> [a]
:[IndexOf Vector
jforall a. Num a => a -> a -> a
+Int
1..forall t. Matrix t -> Int
rows Matrix e
mforall 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 = forall t.
Numeric t =>
(Matrix t -> Matrix t)
-> (View1 t -> View1 t) -> Matrix t -> Matrix t
foldMatrix Matrix t -> Matrix t
g forall {t} {a} {c}.
(Eq t, Num a, Num c, Product t, Container Vector t, Fractional t,
 Num (Vector t)) =>
(t, Vector t, Vector t, Matrix t) -> (a, Vector t, c, Matrix t)
f Matrix t
a
  where
    f :: (t, Vector t, Vector t, Matrix t) -> (a, Vector t, c, Matrix t)
f (t
e,Vector t
r,Vector t
c,Matrix t
m)
        | t
e forall a. Eq a => a -> a -> Bool
/= t
0    = (a
1, Vector t
r', c
0, Matrix t
m forall a. Num a => a -> a -> a
- forall t. Product t => Vector t -> Vector t -> Matrix t
outer Vector t
c Vector t
r')
        | Bool
otherwise = forall a. HasCallStack => String -> a
error String
"singular!"
      where
        r' :: Vector t
r' = Vector t
r forall a. Fractional a => a -> a -> a
/ forall (c :: * -> *) e. Container c e => e -> c e
scalar t
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 :: forall t.
(Eq t, Fractional t, Num (Vector t), Numeric t) =>
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 = forall t. Element t => Matrix t -> Matrix t
flipud forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Matrix t
fliprl
    splitColsAt :: Int -> Matrix t -> (Matrix t, Matrix t)
splitColsAt Int
n = (forall t. Element t => Int -> Matrix t -> Matrix t
takeColumns Int
n forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& 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 = forall {t}. Element t => Int -> Matrix t -> (Matrix t, Matrix t)
splitColsAt (forall t. Matrix t -> Int
cols Matrix t
a) (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 forall a b. (a -> b) -> a -> b
$ Matrix t
x forall t. Element t => Matrix t -> Matrix t -> Matrix t
||| Matrix t
y)
    (Matrix t
a1,Matrix t
b1) = forall {t}.
(Numeric t, Eq t, Fractional t, Num (Vector t)) =>
(Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
go (forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {e}.
(Container Vector e, Num (Vector e), CTrans e) =>
Int -> Matrix e -> (Int, Matrix e)
swapMax Int
0) Matrix t
a Matrix t
b
    ( Matrix t
_, Matrix t
r) = forall {t}.
(Numeric t, Eq t, Fractional t, Num (Vector t)) =>
(Matrix t -> Matrix t)
-> Matrix t -> Matrix t -> (Matrix t, Matrix t)
go forall a. a -> a
id (Matrix t -> Matrix t
flipudrl forall a b. (a -> b) -> a -> b
$ Matrix t
a1) (Matrix t -> Matrix t
flipudrl 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 :: forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Matrix t -> Matrix t -> Matrix t
gaussElim_1 Matrix t
x Matrix t
y = forall t. Element t => Int -> Matrix t -> Matrix t
dropColumns (forall t. Matrix t -> Int
rows Matrix t
x) (forall t. Element t => Matrix t -> Matrix t
flipud forall a b. (a -> b) -> a -> b
$ forall t. Element t => [Vector t] -> Matrix t
fromRows [Vector t]
s2)
  where
    rs :: [Vector t]
rs = forall t. Element t => Matrix t -> [Vector t]
toRows forall a b. (a -> b) -> a -> b
$ Matrix t
x forall t. Element t => Matrix t -> Matrix t -> Matrix t
||| Matrix t
y
    s1 :: Matrix t
s1 = forall t. Element t => [Vector t] -> Matrix t
fromRows forall a b. (a -> b) -> a -> b
$ forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> Int -> [Vector t] -> [Vector t]
pivotDown (forall t. Matrix t -> Int
rows Matrix t
x) Int
0 [Vector t]
rs      -- interesting
    s2 :: [Vector t]
s2 = forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> [Vector t] -> [Vector t]
pivotUp (forall t. Matrix t -> Int
rows Matrix t
xforall a. Num a => a -> a -> a
-Int
1) (forall t. Element t => Matrix t -> [Vector t]
toRows forall a b. (a -> b) -> a -> b
$ 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 :: 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
n [Vector t]
xs
    | Int
t forall a. Eq a => a -> a -> Bool
== Int
n    = []
    | Bool
otherwise = Vector t
y forall a. a -> [a] -> [a]
: 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
nforall a. Num a => a -> a -> a
+Int
1) [Vector t]
ys
  where
    Vector t
y:[Vector t]
ys = (Int, [Vector t]) -> [Vector t]
redu (forall {b} {a}.
(Ord b, Num b, Indexable a b) =>
Int -> [a] -> (Int, [a])
pivot Int
n [Vector t]
xs)

    pivot :: Int -> [a] -> (Int, [a])
pivot Int
k = (forall a b. a -> b -> a
const Int
k forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& forall a. a -> a
id)
            forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Ord a => a -> a -> Ordering
compare forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (forall a. Num a => a -> a
absforall b c a. (b -> c) -> (a -> b) -> a -> c
. (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 forall a. Eq a => a -> a -> Bool
== t
0 = forall a. HasCallStack => String -> a
error String
"gauss: singular!"  -- FIXME
        | Bool
otherwise = Vector t
u forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map Vector t -> Vector t
f [Vector t]
zs
      where
        p :: t
p = Vector t
xforall c t. Indexable c t => c -> Int -> t
!Int
k
        u :: Vector t
u = forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (forall a. Fractional a => a -> a
recip (Vector t
xforall 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 forall a. Num a => a -> a -> a
- forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (Vector t
zforall 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 :: forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> [Vector t] -> [Vector t]
pivotUp Int
n [Vector t]
xs
    | Int
n forall a. Eq a => a -> a -> Bool
== -Int
1 = []
    | Bool
otherwise = Vector t
y forall a. a -> [a] -> [a]
: forall t.
(Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
 Numeric t) =>
Int -> [Vector t] -> [Vector t]
pivotUp (Int
nforall 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 forall a. a -> [a] -> [a]
: 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 forall a. Num a => a -> a -> a
- forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (Vector t
zforall 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 = forall t. Element t => Int -> Matrix t -> Matrix t
dropColumns (forall t. Matrix t -> Int
rows Matrix t
a) forall a b. (a -> b) -> a -> b
$ forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable 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 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
rforall 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 = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (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   = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (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   = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (forall t. t -> RowRange -> ColRange -> RowOper t
SCAL t
a (Int -> RowRange
Row Int
i) ColRange
AllCols) STMatrix s t
m
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
n] forall a b. (a -> b) -> a -> b
$ \Int
i -> do
        Int
c <- forall (c :: * -> *) t. Container c t => c t -> IndexOf c
maxIndex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Vector t
flatten forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> 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)
        forall {t} {s}.
(Num t, Element t) =>
STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
x Int
i (Int
iforall a. Num a => a -> a -> a
+Int
c)
        t
a <- forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
i Int
i
        forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (t
a forall a. Eq a => a -> a -> Bool
== t
0) forall a b. (a -> b) -> a -> b
$ forall a. HasCallStack => String -> a
error String
"singular!"
        forall {t} {s}.
(Num t, Element t) =>
STMatrix s t -> t -> Int -> ST s ()
scal STMatrix s t
x (forall a. Fractional a => a -> a
recip t
a) Int
i
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
iforall a. Num a => a -> a -> a
+Int
1..Int
n] forall a b. (a -> b) -> a -> b
$ \Int
j -> do
            t
b <- forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
j Int
i
            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
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
n,Int
nforall a. Num a => a -> a -> a
-Int
1..Int
1] forall a b. (a -> b) -> a -> b
$ \Int
i -> do
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
iforall a. Num a => a -> a -> a
-Int
1,Int
iforall a. Num a => a -> a -> a
-Int
2..Int
0] forall a b. (a -> b) -> a -> b
$ \Int
j -> do
            t
b <- forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
j Int
i
            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 = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (forall t. t -> Int -> Int -> ColRange -> RowOper t
AXPY t
a Int
i Int
j (Int -> ColRange
FromCol (Int
iforall 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   = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
j ColRange
AllCols)           STMatrix s t
m
    STVector s Int
p <- forall t s. Storable t => Int -> ST s (STVector s t)
newUndefinedVector Int
r
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
rforall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \Int
i -> do
        Int
k <- forall (c :: * -> *) t. Container c t => c t -> IndexOf c
maxIndex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Vector t
flatten forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> 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)
        forall t s. Storable t => STVector s t -> Int -> t -> ST s ()
writeVector STVector s Int
p Int
i (Int
kforall a. Num a => a -> a -> a
+Int
i)
        forall {t} {s}.
(Num t, Element t) =>
STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
x Int
i (Int
iforall a. Num a => a -> a -> a
+Int
k)
        t
a <- forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
i Int
i
        forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (t -> Bool
ok t
a) forall a b. (a -> b) -> a -> b
$ do
            forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
iforall a. Num a => a -> a -> a
+Int
1..Int
rforall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \Int
j -> do
                t
b <- (forall a. Fractional a => a -> a -> a
/t
a) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s t
x Int
j Int
i
                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
                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 <- forall t s. Storable t => STVector s t -> ST s (Vector t)
unsafeFreezeVector STVector s Int
p
    forall (m :: * -> *) a. Monad m => a -> m a
return (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 = forall t. Matrix t -> [Int] -> LU t
LU Matrix t
m [Int]
p
  where
    (Matrix t
m,[Int]
p) = forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable (forall {t} {b} {s}.
(Container Vector t, Num (Vector t), Fractional t) =>
(t -> Bool) -> (Int, b) -> STMatrix s t -> ST s [Int]
luST (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) = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (forall t. t -> RowRange -> ColRange -> RowOper t
SCAL t
a (Int -> Int -> RowRange
RowRange Int
r0 (Int
r0forall a. Num a => a -> a -> a
+Int
nrforall a. Num a => a -> a -> a
-Int
1)) (Int -> Int -> ColRange
ColRange Int
c0 (Int
c0forall a. Num a => a -> a -> a
+Int
ncforall 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 <- 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
rforall a. Num a => a -> a -> a
-Int
1forall a. Num a => a -> a -> a
-Int
k
        o :: Int
o  = if Int
k forall a. Ord a => a -> a -> Bool
< Int
rforall a. Num a => a -> a -> a
-Int
1 then Int
1 else Int
0
        s :: Slice s a
s = forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
x (Int
kforall a. Num a => a -> a -> a
+Int
1) (Int
kforall a. Num a => a -> a -> a
+Int
1) Int
rr Int
rr
        u :: Slice s a
u = forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
x Int
k     (Int
kforall a. Num a => a -> a -> a
+Int
1) Int
o  Int
rr
        l :: Slice s a
l = forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
Slice STMatrix s a
x (Int
kforall a. Num a => a -> a -> a
+Int
1) Int
k     Int
rr Int
o
    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 <- 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 <- forall t s. Storable t => STVector s t -> ST s (Vector t)
unsafeFreezeVector STVector s t
p
    forall (m :: * -> *) a. Monad m => a -> m a
return Vector t
v


luPacked'' :: Matrix t -> (Matrix t, [Int])
luPacked'' Matrix t
m = (forall a. a -> a
id forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** forall a. Storable a => Vector a -> [a]
toList) (forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable (forall {t} {t} {t} {s} {a}.
Storable t =>
Int
-> (t -> t -> STVector s t -> ST s a) -> t -> t -> ST s (Vector t)
withVec (forall t. Matrix t -> Int
rows Matrix t
m) 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
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0..Int
rforall a. Num a => a -> a -> a
-Int
1] forall a b. (a -> b) -> a -> b
$ \Int
k -> do
            forall {e} {s}.
(Container Vector e, Num (Vector e), Num e) =>
STMatrix s e -> 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) <- 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
            forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (forall t. (Element t, Normed (Vector t)) => Double -> t -> Bool
magnit Double
0 t
d) forall a b. (a -> b) -> a -> b
$ do
                forall {t} {s}. (Num t, Element t) => t -> Slice s t -> ST s ()
scalS (forall a. Fractional a => a -> a
recip t
d) Slice s t
l
                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 e -> STVector s Int -> Int -> ST s ()
pivot STMatrix s e
x STVector s Int
p Int
k = do
        Int
j <- forall (c :: * -> *) t. Container c t => c t -> IndexOf c
maxIndex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Vector t
flatten forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a t s.
Element a =>
STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
extractMatrix STMatrix s e
x (Int -> RowRange
FromRow Int
k) (Int -> ColRange
Col Int
k)
        forall t s. Storable t => STVector s t -> Int -> t -> ST s ()
writeVector STVector s Int
p Int
k (Int
jforall a. Num a => a -> a -> a
+Int
k)
        Int -> Int -> ST s ()
swap Int
k (Int
kforall a. Num a => a -> a -> a
+Int
j)
      where
        swap :: Int -> Int -> ST s ()
swap Int
i Int
j = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
j ColRange
AllCols) STMatrix s e
x

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

rowRange :: Matrix t -> [Int]
rowRange Matrix t
m = [Int
0..forall t. Matrix t -> Int
rows Matrix t
m 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 = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' 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
rhsforall t. Element t => Matrix t -> [Int] -> 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 <- forall {t}. Matrix t -> [Int]
rowRange Matrix t
lup ]
      where
        d :: Int -> Matrix t
d Int
k = Matrix t
lup 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 forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
at Int
k, Int -> Extractor
Drop (Int
kforall a. Num a => a -> a -> a
+Int
1))
        b :: Int -> Matrix t
b Int
k = Matrix t
rhs 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 forall a. Num a => a -> a -> a
- a t
uforall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
(Mul a b c, Product t) =>
a t -> b t -> c t
<>Matrix t
x) forall a. Fractional a => a -> a -> a
/ Matrix t
d
                       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 = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' 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
rhsforall 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 <- forall {t}. Matrix t -> [Int]
rowRange Matrix t
lup ]
      where
        l :: Int -> Matrix t
l Int
k = Matrix t
lup 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 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
                   forall t. Element t => Matrix t -> Matrix t -> Matrix t
===
                (Matrix t
b forall a. Num a => a -> a -> a
- a t
lforall (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 = forall {t}.
(Container Vector t, Fractional t, Product t, Num (Vector t)) =>
Matrix t -> Matrix t -> Matrix t
backSust' Matrix t
lup (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 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 = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ 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 <- forall t s. Storable t => Matrix t -> ST s (STMatrix s t)
unsafeThawMatrix Matrix t
lup
        let go :: Int -> ST s ()
go Int
k = forall t s.
Element t =>
t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
gemmm t
1 (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) (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) (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)
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Int -> ST s ()
go [Int
0..Int
rforall a. Num a => a -> a -> a
-Int
1]


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

            go :: Int -> ST s ()
go Int
k = forall t s.
Element t =>
t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
gemmm t
1 (Int -> Slice s t
b Int
k) (-t
1) (Int -> Slice s t
u Int
k) (Int -> Slice s t
x Int
k) forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Int -> ST s ()
scal Int
k
        forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Int -> ST s ()
go [Int
rforall a. Num a => a -> a -> a
-Int
1,Int
rforall 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 = forall {t}.
(Fractional t, Container Vector t) =>
Matrix t -> Matrix t -> Matrix t
backSust Matrix t
lup (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 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
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) forall a. Eq a => a -> a -> Bool
== (Int
1,Int
1) = forall t b. t -> MatrixView t b
Elem (forall t. Storable t => Matrix t -> Int -> Int -> t
atM' Matrix t
m Int
0 Int
0)
    | Bool
otherwise        = 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) = forall (c :: * -> *) t. Container c t => c t -> 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
ctforall 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
rtforall 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
rtforall a. Num a => a -> a -> a
-Int
r,Int
ctforall a. Num a => a -> a -> a
-Int
c) Matrix t
m
    subm :: (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
subm = forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix

viewBlock :: Matrix t -> MatrixView t (Matrix t)
viewBlock Matrix t
m = 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 = forall t. Matrix t -> Int
rows Matrix t
m forall a. Integral a => a -> a -> a
`div` Int
2

invershur :: Matrix t -> Matrix t
invershur (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) = 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 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 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 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
r4forall 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 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 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 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
r1forall a. Num a => a -> a -> a
-Matrix t
r7
    d' :: Matrix t
d' = -Matrix t
r6

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