{-# LANGUAGE PatternGuards #-}
-- |
-- Module    : Statistics.Matrix
-- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Basic matrix operations.
--
-- There isn't a widely used matrix package for Haskell yet, so
-- we implement the necessary minimum here.

module Statistics.Matrix
    ( -- * Data types
      Matrix(..)
    , Vector
      -- * Conversion from/to lists/vectors
    , fromVector
    , fromList
    , fromRowLists
    , fromRows
    , fromColumns
    , toVector
    , toList
    , toRows
    , toColumns
    , toRowLists
      -- * Other
    , generate
    , generateSym
    , ident
    , diag
    , dimension
    , center
    , multiply
    , multiplyV
    , transpose
    , power
    , norm
    , column
    , row
    , map
    , for
    , unsafeIndex
    , hasNaN
    , bounds
    , unsafeBounds
    ) where

import Prelude hiding (exponent, map)
import Control.Applicative ((<$>))
import Control.Monad.ST
import qualified Data.Vector.Unboxed as U
import           Data.Vector.Unboxed   ((!))
import qualified Data.Vector.Unboxed.Mutable as UM

import Numeric.Sum         (sumVector,kbn)
import Statistics.Matrix.Function
import Statistics.Matrix.Types
import Statistics.Matrix.Mutable  (unsafeNew,unsafeWrite,unsafeFreeze)



----------------------------------------------------------------
-- Conversion to/from vectors/lists
----------------------------------------------------------------

-- | Convert from a row-major list.
fromList :: Int                 -- ^ Number of rows.
         -> Int                 -- ^ Number of columns.
         -> [Double]            -- ^ Flat list of values, in row-major order.
         -> Matrix
fromList :: Int -> Int -> [Double] -> Matrix
fromList Int
r Int
c = Int -> Int -> Vector Double -> Matrix
fromVector Int
r Int
c (Vector Double -> Matrix)
-> ([Double] -> Vector Double) -> [Double] -> Matrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList

-- | create a matrix from a list of lists, as rows
fromRowLists :: [[Double]] -> Matrix
fromRowLists :: [[Double]] -> Matrix
fromRowLists = [Vector Double] -> Matrix
fromRows ([Vector Double] -> Matrix)
-> ([[Double]] -> [Vector Double]) -> [[Double]] -> Matrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> Vector Double) -> [[Double]] -> [Vector Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList

-- | Convert from a row-major vector.
fromVector :: Int               -- ^ Number of rows.
           -> Int               -- ^ Number of columns.
           -> U.Vector Double   -- ^ Flat list of values, in row-major order.
           -> Matrix
fromVector :: Int -> Int -> Vector Double -> Matrix
fromVector Int
r Int
c Vector Double
v
  | Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
c Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
len = [Char] -> Matrix
forall a. HasCallStack => [Char] -> a
error [Char]
"input size mismatch"
  | Bool
otherwise  = Int -> Int -> Vector Double -> Matrix
Matrix Int
r Int
c Vector Double
v
  where len :: Int
len    = Vector Double -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Double
v

-- | create a matrix from a list of vectors, as rows
fromRows :: [Vector] -> Matrix
fromRows :: [Vector Double] -> Matrix
fromRows [Vector Double]
xs
  | [] <- [Vector Double]
xs        = [Char] -> Matrix
forall a. HasCallStack => [Char] -> a
error [Char]
"Statistics.Matrix.fromRows: empty list of rows!"
  | (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/=Int
nCol) [Int]
ns = [Char] -> Matrix
forall a. HasCallStack => [Char] -> a
error [Char]
"Statistics.Matrix.fromRows: row sizes do not match"
  | Int
nCol Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0       = [Char] -> Matrix
forall a. HasCallStack => [Char] -> a
error [Char]
"Statistics.Matrix.fromRows: zero columns in matrix"
  | Bool
otherwise       = Int -> Int -> Vector Double -> Matrix
fromVector Int
nRow Int
nCol ([Vector Double] -> Vector Double
forall a. Unbox a => [Vector a] -> Vector a
U.concat [Vector Double]
xs)
  where
    Int
nCol:[Int]
ns = Vector Double -> Int
forall a. Unbox a => Vector a -> Int
U.length (Vector Double -> Int) -> [Vector Double] -> [Int]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Vector Double]
xs
    nRow :: Int
nRow    = [Vector Double] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Vector Double]
xs


-- | create a matrix from a list of vectors, as columns
fromColumns :: [Vector] -> Matrix
fromColumns :: [Vector Double] -> Matrix
fromColumns = Matrix -> Matrix
transpose (Matrix -> Matrix)
-> ([Vector Double] -> Matrix) -> [Vector Double] -> Matrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Vector Double] -> Matrix
fromRows

-- | Convert to a row-major flat vector.
toVector :: Matrix -> U.Vector Double
toVector :: Matrix -> Vector Double
toVector (Matrix Int
_ Int
_ Vector Double
v) = Vector Double
v

-- | Convert to a row-major flat list.
toList :: Matrix -> [Double]
toList :: Matrix -> [Double]
toList = Vector Double -> [Double]
forall a. Unbox a => Vector a -> [a]
U.toList (Vector Double -> [Double])
-> (Matrix -> Vector Double) -> Matrix -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix -> Vector Double
toVector

-- | Convert to a list of lists, as rows
toRowLists :: Matrix -> [[Double]]
toRowLists :: Matrix -> [[Double]]
toRowLists (Matrix Int
_ Int
nCol Vector Double
v)
  = [Double] -> [[Double]]
forall a. [a] -> [[a]]
chunks ([Double] -> [[Double]]) -> [Double] -> [[Double]]
forall a b. (a -> b) -> a -> b
$ Vector Double -> [Double]
forall a. Unbox a => Vector a -> [a]
U.toList Vector Double
v
  where
    chunks :: [a] -> [[a]]
chunks [] = []
    chunks [a]
xs = case Int -> [a] -> ([a], [a])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
nCol [a]
xs of
      ([a]
rowE,[a]
rest) -> [a]
rowE [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [a] -> [[a]]
chunks [a]
rest


-- | Convert to a list of vectors, as rows
toRows :: Matrix -> [Vector]
toRows :: Matrix -> [Vector Double]
toRows (Matrix Int
_ Int
nCol Vector Double
v) = Vector Double -> [Vector Double]
forall a. Unbox a => Vector a -> [Vector a]
chunks Vector Double
v
  where
    chunks :: Vector a -> [Vector a]
chunks Vector a
xs
      | Vector a -> Bool
forall a. Unbox a => Vector a -> Bool
U.null Vector a
xs = []
      | Bool
otherwise = case Int -> Vector a -> (Vector a, Vector a)
forall a. Unbox a => Int -> Vector a -> (Vector a, Vector a)
U.splitAt Int
nCol Vector a
xs of
          (Vector a
rowE,Vector a
rest) -> Vector a
rowE Vector a -> [Vector a] -> [Vector a]
forall a. a -> [a] -> [a]
: Vector a -> [Vector a]
chunks Vector a
rest

-- | Convert to a list of vectors, as columns
toColumns :: Matrix -> [Vector]
toColumns :: Matrix -> [Vector Double]
toColumns = Matrix -> [Vector Double]
toRows (Matrix -> [Vector Double])
-> (Matrix -> Matrix) -> Matrix -> [Vector Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix -> Matrix
transpose



----------------------------------------------------------------
-- Other
----------------------------------------------------------------

-- | Generate matrix using function
generate :: Int                 -- ^ Number of rows
         -> Int                 -- ^ Number of columns
         -> (Int -> Int -> Double)
            -- ^ Function which takes /row/ and /column/ as argument.
         -> Matrix
generate :: Int -> Int -> (Int -> Int -> Double) -> Matrix
generate Int
nRow Int
nCol Int -> Int -> Double
f
  = Int -> Int -> Vector Double -> Matrix
Matrix Int
nRow Int
nCol (Vector Double -> Matrix) -> Vector Double -> Matrix
forall a b. (a -> b) -> a -> b
$ Int -> (Int -> Double) -> Vector Double
forall a. Unbox a => Int -> (Int -> a) -> Vector a
U.generate (Int
nRowInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nCol) ((Int -> Double) -> Vector Double)
-> (Int -> Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ \Int
i ->
      let (Int
r,Int
c) = Int
i Int -> Int -> (Int, Int)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Int
nCol in Int -> Int -> Double
f Int
r Int
c

-- | Generate symmetric square matrix using function
generateSym
  :: Int                 -- ^ Number of rows and columns
  -> (Int -> Int -> Double)
     -- ^ Function which takes /row/ and /column/ as argument. It must
     --   be symmetric in arguments: @f i j == f j i@
  -> Matrix
generateSym :: Int -> (Int -> Int -> Double) -> Matrix
generateSym Int
n Int -> Int -> Double
f = (forall s. ST s Matrix) -> Matrix
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s Matrix) -> Matrix)
-> (forall s. ST s Matrix) -> Matrix
forall a b. (a -> b) -> a -> b
$ do
  MMatrix s
m <- Int -> Int -> ST s (MMatrix s)
forall s. Int -> Int -> ST s (MMatrix s)
unsafeNew Int
n Int
n
  Int -> Int -> (Int -> ST s ()) -> ST s ()
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for Int
0 Int
n ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
r -> do
    MMatrix s -> Int -> Int -> Double -> ST s ()
forall s. MMatrix s -> Int -> Int -> Double -> ST s ()
unsafeWrite MMatrix s
m Int
r Int
r (Int -> Int -> Double
f Int
r Int
r)
    Int -> Int -> (Int -> ST s ()) -> ST s ()
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for (Int
rInt -> 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
c -> do
      let x :: Double
x = Int -> Int -> Double
f Int
r Int
c
      MMatrix s -> Int -> Int -> Double -> ST s ()
forall s. MMatrix s -> Int -> Int -> Double -> ST s ()
unsafeWrite MMatrix s
m Int
r Int
c Double
x
      MMatrix s -> Int -> Int -> Double -> ST s ()
forall s. MMatrix s -> Int -> Int -> Double -> ST s ()
unsafeWrite MMatrix s
m Int
c Int
r Double
x
  MMatrix s -> ST s Matrix
forall s. MMatrix s -> ST s Matrix
unsafeFreeze MMatrix s
m


-- | Create the square identity matrix with given dimensions.
ident :: Int -> Matrix
ident :: Int -> Matrix
ident Int
n = Vector Double -> Matrix
diag (Vector Double -> Matrix) -> Vector Double -> Matrix
forall a b. (a -> b) -> a -> b
$ Int -> Double -> Vector Double
forall a. Unbox a => Int -> a -> Vector a
U.replicate Int
n Double
1.0

-- | Create a square matrix with given diagonal, other entries default to 0
diag :: Vector -> Matrix
diag :: Vector Double -> Matrix
diag Vector Double
v
  = Int -> Int -> Vector Double -> Matrix
Matrix Int
n Int
n (Vector Double -> Matrix) -> Vector Double -> Matrix
forall a b. (a -> b) -> a -> b
$ (forall s. ST s (MVector s Double)) -> Vector Double
forall a. Unbox a => (forall s. ST s (MVector s a)) -> Vector a
U.create ((forall s. ST s (MVector s Double)) -> Vector Double)
-> (forall s. ST s (MVector s Double)) -> Vector Double
forall a b. (a -> b) -> a -> b
$ do
      MVector s Double
arr <- Int -> Double -> ST s (MVector (PrimState (ST s)) Double)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
UM.replicate (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) Double
0
      Int -> Int -> (Int -> ST s ()) -> ST s ()
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
for Int
0 Int
n ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i ->
        MVector (PrimState (ST s)) Double -> Int -> Double -> ST s ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
UM.unsafeWrite MVector s Double
MVector (PrimState (ST s)) Double
arr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i) (Vector Double
v Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
! Int
i)
      MVector s Double -> ST s (MVector s Double)
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Double
arr
  where
    n :: Int
n = Vector Double -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Double
v

-- | Return the dimensions of this matrix, as a (row,column) pair.
dimension :: Matrix -> (Int, Int)
dimension :: Matrix -> (Int, Int)
dimension (Matrix Int
r Int
c Vector Double
_) = (Int
r, Int
c)

-- | Matrix-matrix multiplication. Matrices must be of compatible
-- sizes (/note: not checked/).
multiply :: Matrix -> Matrix -> Matrix
multiply :: Matrix -> Matrix -> Matrix
multiply m1 :: Matrix
m1@(Matrix Int
r1 Int
_ Vector Double
_) m2 :: Matrix
m2@(Matrix Int
_ Int
c2 Vector Double
_) =
  Int -> Int -> Vector Double -> Matrix
Matrix Int
r1 Int
c2 (Vector Double -> Matrix) -> Vector Double -> Matrix
forall a b. (a -> b) -> a -> b
$ Int -> (Int -> Double) -> Vector Double
forall a. Unbox a => Int -> (Int -> a) -> Vector a
U.generate (Int
r1Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
c2) Int -> Double
go
  where
    go :: Int -> Double
go Int
t = (KBNSum -> Double) -> Vector Double -> Double
forall (v :: * -> *) s.
(Vector v Double, Summation s) =>
(s -> Double) -> v Double -> Double
sumVector KBNSum -> Double
kbn (Vector Double -> Double) -> Vector Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
U.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) (Matrix -> Int -> Vector Double
row Matrix
m1 Int
i) (Matrix -> Int -> Vector Double
column Matrix
m2 Int
j)
      where (Int
i,Int
j) = Int
t Int -> Int -> (Int, Int)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Int
c2

-- | Matrix-vector multiplication.
multiplyV :: Matrix -> Vector -> Vector
multiplyV :: Matrix -> Vector Double -> Vector Double
multiplyV Matrix
m Vector Double
v
  | Matrix -> Int
cols Matrix
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
c = Int -> (Int -> Double) -> Vector Double
forall a. Unbox a => Int -> (Int -> a) -> Vector a
U.generate (Matrix -> Int
rows Matrix
m) ((KBNSum -> Double) -> Vector Double -> Double
forall (v :: * -> *) s.
(Vector v Double, Summation s) =>
(s -> Double) -> v Double -> Double
sumVector KBNSum -> Double
kbn (Vector Double -> Double)
-> (Int -> Vector Double) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
U.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) Vector Double
v (Vector Double -> Vector Double)
-> (Int -> Vector Double) -> Int -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix -> Int -> Vector Double
row Matrix
m)
  | Bool
otherwise   = [Char] -> Vector Double
forall a. HasCallStack => [Char] -> a
error ([Char] -> Vector Double) -> [Char] -> Vector Double
forall a b. (a -> b) -> a -> b
$ [Char]
"matrix/vector unconformable " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> [Char]
forall a. Show a => a -> [Char]
show (Matrix -> Int
cols Matrix
m,Int
c)
  where c :: Int
c = Vector Double -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Double
v

-- | Raise matrix to /n/th power. Power must be positive
-- (/note: not checked).
power :: Matrix -> Int -> Matrix
power :: Matrix -> Int -> Matrix
power Matrix
mat Int
1 = Matrix
mat
power Matrix
mat Int
n = Matrix
res
  where
    mat2 :: Matrix
mat2 = Matrix -> Int -> Matrix
power Matrix
mat (Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2)
    pow :: Matrix
pow  = Matrix -> Matrix -> Matrix
multiply Matrix
mat2 Matrix
mat2
    res :: Matrix
res | Int -> Bool
forall a. Integral a => a -> Bool
odd Int
n     = Matrix -> Matrix -> Matrix
multiply Matrix
pow Matrix
mat
        | Bool
otherwise = Matrix
pow

-- | Element in the center of matrix (not corrected for exponent).
center :: Matrix -> Double
center :: Matrix -> Double
center mat :: Matrix
mat@(Matrix Int
r Int
c Vector Double
_) =
    (Vector Double -> Int -> Double) -> Matrix -> Int -> Int -> Double
forall r. (Vector Double -> Int -> r) -> Matrix -> Int -> Int -> r
unsafeBounds Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Matrix
mat (Int
r Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2) (Int
c Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2)

-- | Calculate the Euclidean norm of a vector.
norm :: Vector -> Double
norm :: Vector Double -> Double
norm = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double)
-> (Vector Double -> Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (KBNSum -> Double) -> Vector Double -> Double
forall (v :: * -> *) s.
(Vector v Double, Summation s) =>
(s -> Double) -> v Double -> Double
sumVector KBNSum -> Double
kbn (Vector Double -> Double)
-> (Vector Double -> Vector Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
U.map Double -> Double
square

-- | Return the given column.
column :: Matrix -> Int -> Vector
column :: Matrix -> Int -> Vector Double
column (Matrix Int
r Int
c Vector Double
v) Int
i = Vector Double -> Vector Int -> Vector Double
forall a. Unbox a => Vector a -> Vector Int -> Vector a
U.backpermute Vector Double
v (Vector Int -> Vector Double) -> Vector Int -> Vector Double
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Int -> Vector Int
forall a. (Unbox a, Num a) => a -> a -> Int -> Vector a
U.enumFromStepN Int
i Int
c Int
r
{-# INLINE column #-}

-- | Return the given row.
row :: Matrix -> Int -> Vector
row :: Matrix -> Int -> Vector Double
row (Matrix Int
_ Int
c Vector Double
v) Int
i = Int -> Int -> Vector Double -> Vector Double
forall a. Unbox a => Int -> Int -> Vector a -> Vector a
U.slice (Int
cInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i) Int
c Vector Double
v

unsafeIndex :: Matrix
            -> Int              -- ^ Row.
            -> Int              -- ^ Column.
            -> Double
unsafeIndex :: Matrix -> Int -> Int -> Double
unsafeIndex = (Vector Double -> Int -> Double) -> Matrix -> Int -> Int -> Double
forall r. (Vector Double -> Int -> r) -> Matrix -> Int -> Int -> r
unsafeBounds Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex

-- | Apply function to every element of matrix
map :: (Double -> Double) -> Matrix -> Matrix
map :: (Double -> Double) -> Matrix -> Matrix
map Double -> Double
f (Matrix Int
r Int
c Vector Double
v) = Int -> Int -> Vector Double -> Matrix
Matrix Int
r Int
c ((Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
U.map Double -> Double
f Vector Double
v)

-- | Indicate whether any element of the matrix is @NaN@.
hasNaN :: Matrix -> Bool
hasNaN :: Matrix -> Bool
hasNaN = (Double -> Bool) -> Vector Double -> Bool
forall a. Unbox a => (a -> Bool) -> Vector a -> Bool
U.any Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN (Vector Double -> Bool)
-> (Matrix -> Vector Double) -> Matrix -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix -> Vector Double
toVector

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector.
bounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
bounds :: (Vector Double -> Int -> r) -> Matrix -> Int -> Int -> r
bounds Vector Double -> Int -> r
k (Matrix Int
rs Int
cs Vector Double
v) Int
r Int
c
  | Int
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
rs = [Char] -> r
forall a. HasCallStack => [Char] -> a
error [Char]
"row out of bounds"
  | Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
cs = [Char] -> r
forall a. HasCallStack => [Char] -> a
error [Char]
"column out of bounds"
  | Bool
otherwise        = Vector Double -> Int -> r
k Vector Double
v (Int -> r) -> Int -> r
forall a b. (a -> b) -> a -> b
$! Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
cs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c
{-# INLINE bounds #-}

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector, without checking.
unsafeBounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
unsafeBounds :: (Vector Double -> Int -> r) -> Matrix -> Int -> Int -> r
unsafeBounds Vector Double -> Int -> r
k (Matrix Int
_ Int
cs Vector Double
v) Int
r Int
c = Vector Double -> Int -> r
k Vector Double
v (Int -> r) -> Int -> r
forall a b. (a -> b) -> a -> b
$! Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
cs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c
{-# INLINE unsafeBounds #-}

transpose :: Matrix -> Matrix
transpose :: Matrix -> Matrix
transpose m :: Matrix
m@(Matrix Int
r0 Int
c0 Vector Double
_) = Int -> Int -> Vector Double -> Matrix
Matrix Int
c0 Int
r0 (Vector Double -> Matrix)
-> ((Int -> Double) -> Vector Double) -> (Int -> Double) -> Matrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> (Int -> Double) -> Vector Double
forall a. Unbox a => Int -> (Int -> a) -> Vector a
U.generate (Int
r0Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
c0) ((Int -> Double) -> Matrix) -> (Int -> Double) -> Matrix
forall a b. (a -> b) -> a -> b
$ \Int
i ->
  let (Int
r,Int
c) = Int
i Int -> Int -> (Int, Int)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Int
r0
  in Matrix -> Int -> Int -> Double
unsafeIndex Matrix
m Int
c Int
r