This is an alternative site for discovering Elm packages. You may be looking for the official Elm package site instead.

# Matrix

A matrix library written completely in Elm. This library aims to be a reasonably complete suite of linear algebra tools.

Some highlights are that this library has generic sized matrices, transposes, multiplication, and inversion.

``````import Matrix as Mt --and program away!
``````

# The Matrix Type

type Matrix = Mat Matnxn | Err String

The Matrix type. It can either be an actual matrix or an error string.

# Creating Matrices

fromList : ( Int, Int ) -> List Float -> Matrix

Create a (n rows x m columns) matrix with the list as the elements. Fails if dimension mismatch. Elements need to be specified in row-major order.

``````matrix =
Matrix.fromList ( 2, 3 ) [ 2, 2, 2, 3, 3, 3 ]
``````
from2DList : List (List Float) -> Matrix

Create a (n x m) matrix with inner lists being rows. The following is a 2 x 3 matrix:

``````matrix =
Matrix.from2DList
[ [ 2, 2, 2 ]
, [ 3, 3, 3 ]
]
``````
mat : List (List Float) -> Matrix

Create a (n x m) matrix with inner lists being rows. The following is a 2 x 3 matrix:

``````matrix =
Matrix.mat
[ [ 2, 2, 2 ]
, [ 3, 3, 3 ]
]
``````
fromString : String -> Matrix

Create a (n x m) matrix with inner lists being rows. In string format you can use Matlab/Julia-esque syntax. Spaces denote elements in a row. Semicolons denote elements in a column. The string must begin with [ and end with ].

The following is a 2 x 3 matrix:

``````matrix =
Matrix.fromString "[ 2 2 2; 3 3 3]"
``````

Any alpha/garbage characters will be set to zero.

mats : String -> Matrix

Shorthand for fromString

zeroes : ( Int, Int ) -> Matrix

Generate a matrix of zeroes.

``````lots_of_zeroes =
Matrix.zeroes ( 3, 4 )
``````
ones : ( Int, Int ) -> Matrix

Generate a matrix of ones.

``````lots_of_ones =
Matrix.ones ( 4, 3 )
``````
eye : Int -> Matrix

Create an nxn identity matrix.

``````identity =
Matrix.eye 3
``````
upper : Int -> Matrix

Create an nxn upper triangular matrix.

``````triangle =
Matrix.upper 4
``````
lower : Int -> Matrix

Create an nxn lower triangular matrix.

``````ltriangle =
Matrix.lower 4
``````
strictLower : Int -> Matrix

Create an nxn strict lower triangular matrix. This means that elements along the diagonal are zero.

``````sltriangle =
Matrix.strictLower 4
``````
strictUpper : Int -> Matrix

Create an nxn strict upper triangular matrix. This means that elements along the diagonal are zero.

``````striangle =
Matrix.strictUpper 4
``````

# Creating Vectors

cvecFromList : List Float -> Matrix

Create a column vector from a list.

``````column =
Matrix.cvecFromList [ 1, 2, 3, 4 ]
``````
rvecFromList : List Float -> Matrix

Create a row vector from a list.

cvec : List Float -> Matrix

Create a column vector from a list.

rvec : List Float -> Matrix

Create a row vector from a list.

vec : List Float -> Matrix

Create a column vector from a list.

# Vector Operations

cross : Matrix -> Matrix -> Matrix

Get the cross product of two 3d vectors. a >< b

``````x =
Matrix.vec [ 1, 0, 0 ]

y =
Matrix.vec [ 0, 1, 0 ]

z =
Matrix.cross x y
``````
dot : Matrix -> Matrix -> Maybe Float

Performs the dot product of two nxn vectors

``````x =
Matrix.vec [ 1, 0, 0 ]

y =
Matrix.vec [ 0, 1, 0 ]

zero =
Matrix.dot x y
``````

# Matrix Element-wise Operations

add : Matrix -> Matrix -> Matrix

Add two matrices of identical dimensions together

``````a =
Matrix.add (Matrix.zeroes ( 2, 2 )) (Matrix.ones ( 2, 2 ))
``````
equivalent : Float -> Matrix -> Matrix -> Bool

Checks if two matrices are equivalent within some epsilon.

``````epsilon =
10 ^ (-4)

is_equivalent =
equivalent epsilon a b
``````
sMul : Float -> Matrix -> Matrix

Perform scalar multiplication on a matrix.

sDiv : Float -> Matrix -> Matrix

Perform scalar division on a matrix.

map : (Float -> Float) -> Matrix -> Matrix

Map a function over all elements individually

map2 : (Float -> Float -> Float) -> Matrix -> Matrix -> Matrix

Map a function over elements of same index between matrices

eMul : Matrix -> Matrix -> Matrix

Perform element by element multiplication on a matrix.

# Matrix Operations

mul : Matrix -> Matrix -> Matrix

Perform matrix multiplication

``````A * B

a =
Matrix.eye (2, 2)

b =
Matrix.mats "[2, 3; 4 5]"

c =
mul a b
``````
vcat : Matrix -> Matrix -> Matrix

Concatenate two matrices vertically.

hcat : Matrix -> Matrix -> Matrix

Concatenate two matrices horizontally.

get : ( Int, Int ) -> Matrix -> Maybe Float

Get an item at index (row, column). Indices are 1-indexed.

``````a =
Matrix.mats "[3 4; 5 6]"

b =
-- equals 4
get ( 1, 2 ) a
``````
set : ( Int, Int ) -> Float -> Matrix -> Matrix

Get an item at index (row, column). Indices are 1-indexed.

transpose : Matrix -> Matrix

Transpose a matrix.

determinant : Matrix -> Maybe Float

Get the determinant of a square matrix.

``````a =
Matrix.mats "[1 2 3; 4 5 6; 7 8 9]"

is_singular =
if (determinant a) == 0 then
"Matrix is singular"
else
"Matrix is not singular"
``````
det : Matrix -> Maybe Float

Shorthand for determinant.

solveV : Matrix -> Matrix -> Matrix

Solve a system of equations of the form Ax = b. You provide A and b and get back x Where A is a matrix, and b and x are vectors

``````a =
Matrix.eye 3

b =
Matrix.vec [ 3, 2, 1 ]

x =
Matrix.solve a b
``````
solve : Matrix -> Matrix -> Matrix

Solve a system of equations of the form AX = B. You provide A and b and get back x Where A, B, X are matrices B is a matrix of solution vectors horizontally concatenated.

``````a =
Matrix.eye 3

b =
Matrix.hcat (Matrix.vec [ 3, 2, 1 ]) (Matrix.vec [ 1, 2, 3 ])

x =
Matrix.solve a b
``````
invert : Matrix -> Matrix

Invert a square matrix.

``````a =
"[ 2 5; 6 7]"

inva =
invert a

identity =
mul a inva
``````
inv : Matrix -> Matrix

Shorthand for invert.

luDecomp : Matrix -> ( Matrix, Matrix )

Get the lu decomposition of a matrix Since pivoting isn't implemented, watch out for numerical instability

getRows : Matrix -> List Matrix

Returns the rows of a matrix in a list.

getColumns : Matrix -> List Matrix

Returns the columns of a matrix in a list.

size : Matrix -> ( Int, Int )

Returns size of matrix, (rows, columns)

# Matrix Display

toString : Matrix -> String

Change matrix into string form, such as what would be displayed in the terminal.

``````Matrix.toString (Matrix.eye 3) == " 1 0 0\n 0 1 0\n 0 0 1"
``````
debugPrint : Matrix -> String

Helper to debug print. Most useful in repl.

``````> Matrix.debugPrint (Matrix.eye 3)
(3, 3) Matrix
1 0 0
0 1 0
0 0 1
: ""
""
``````

# Interop

to2DList : Matrix -> List (List Float)

Returns matrix as 2d list. Returns empty list if Matrix is in error.

``````Matrix.to2DList (Matrix.eye 2) == [ [1, 0], [0, 1] ]
``````
toFlatList : Matrix -> List Float

Returns matrix as flat list

``````Matrix.toFlatList (Matrix.eye 2) == [ 1, 0, 0, 1 ]
``````
``````module Matrix
exposing
( fromList
, from2DList
, fromString
, mat
, mats
, cvecFromList
, rvecFromList
, cvec
, rvec
, zeroes
, ones
, eye
, upper
, lower
, strictUpper
, strictLower
, vec
, dot
, cross
, sMul
, sDiv
, equivalent
, eMul
, get
, set
, map
, map2
, mul
, vcat
, hcat
, getRows
, getColumns
, transpose
, determinant
, det
, solveV
, solve
, toString
, debugPrint
, to2DList
, toFlatList
, luDecomp
, invert
, inv
, size
, Matrix
)

{-| A matrix library written completely in Elm.
This library aims to be a reasonably complete suite of linear algebra tools.

Some highlights are that this library has generic sized matrices,
transposes, multiplication, and inversion.

import Matrix as Mt --and program away!

# The Matrix Type

@docs Matrix

# Creating Matrices

@docs fromList, from2DList, mat, fromString, mats, zeroes, ones, eye, upper, lower, strictLower, strictUpper

# Creating Vectors

@docs cvecFromList, rvecFromList, cvec, rvec, vec

# Vector Operations

@docs cross, dot

# Matrix Element-wise Operations

@docs add, equivalent, sMul, sDiv, map, map2, eMul

# Matrix Operations

@docs mul, vcat, hcat, get, set, transpose, determinant, det, solveV, solve, invert, inv, luDecomp, getRows, getColumns, size

# Matrix Display

@docs toString, debugPrint

# Interop

@docs to2DList, toFlatList

-}

import Array

type alias Matnxn =
{ dimensions : ( Int, Int )
, elements : Array.Array Float
}

{-| The Matrix type. It can either be an actual matrix or an error string.
-}
type Matrix
= Mat Matnxn
| Err String

-- Matrix Creation

{-| Create a (n rows x m columns) matrix with the list as the elements.
Fails if dimension mismatch. Elements need to be specified in row-major order.

matrix =
Matrix.fromList ( 2, 3 ) [ 2, 2, 2, 3, 3, 3 ]

-}
fromList : ( Int, Int ) -> List Float -> Matrix
fromList dimensions elements =
fromArray dimensions <| Array.fromList elements

{-| Create a (n rows x m columns) matrix with the list as the elements.
Fails if dimension mismatch. Elements need to be specified in row-major order.
-}
fromArray : ( Int, Int ) -> Array.Array Float -> Matrix
fromArray ( rows, columns ) elements =
if rows * columns == Array.length elements then
Mat { dimensions = ( rows, columns ), elements = elements }
else
let
dimensions =
Basics.toString (rows * columns)

numelements =
Basics.toString <| Array.length elements
in
Err <| "The dimensions, row * columns: " ++ dimensions ++ ", do not match the number of elements: " ++ numelements

{-| Create a (n x m) matrix with inner lists being rows.
The following is a 2 x 3 matrix:

matrix =
Matrix.from2DList
[ [ 2, 2, 2 ]
, [ 3, 3, 3 ]
]

-}
from2DList : List (List Float) -> Matrix
from2DList a =
Just data ->
let
columns =
List.length <| data

rows =
List.length a

-- check if all sublists are same length
is_correct =
List.all ((==) columns) <| List.map List.length a
in
if is_correct then
fromList ( rows, columns ) <| List.concat a
else
Err <| "One or more of the sublist rows are malformed"

Nothing ->
fromList ( 0, 0 ) []

{-| Create a (n x m) matrix with inner lists being rows.
The following is a 2 x 3 matrix:

matrix =
Matrix.mat
[ [ 2, 2, 2 ]
, [ 3, 3, 3 ]
]

-}
mat : List (List Float) -> Matrix
mat =
from2DList

{-| Create a (n x m) matrix with inner lists being rows.
In string format you can use Matlab/Julia-esque syntax.
Spaces denote elements in a row.
Semicolons denote elements in a column.
The string must begin with [ and end with ].

The following is a 2 x 3 matrix:

matrix =
Matrix.fromString "[ 2 2 2; 3 3 3]"

Any alpha/garbage characters will be set to zero.

-}
fromString : String -> Matrix
fromString string =
from2DList <| matParser string

{-| Shorthand for fromString
-}
mats : String -> Matrix
mats =
fromString

{-| Create a column vector from a list.

column =
Matrix.cvecFromList [ 1, 2, 3, 4 ]

-}
cvecFromList : List Float -> Matrix
cvecFromList a =
fromList ( List.length a, 1 ) a

{-| Create a column vector from a list.
-}
cvec : List Float -> Matrix
cvec =
cvecFromList

{-| Create a column vector from a list.
-}
vec : List Float -> Matrix
vec =
cvec

{-| Create a row vector from a list.
-}
rvecFromList : List Float -> Matrix
rvecFromList a =
fromList ( 1, List.length a ) a

{-| Create a row vector from a list.
-}
rvec : List Float -> Matrix
rvec =
rvecFromList

{-| Generate a random matrix
-}
rand : ( Int, Int ) -> Matrix
rand a =
Err "Not Implemented"

{-| Generate a matrix of ones.

lots_of_ones =
Matrix.ones ( 4, 3 )

-}
ones : ( Int, Int ) -> Matrix
ones ( rows, columns ) =
fromArray ( rows, columns ) <| Array.repeat (rows * columns) 1.0

{-| Generate a matrix of zeroes.

lots_of_zeroes =
Matrix.zeroes ( 3, 4 )

-}
zeroes : ( Int, Int ) -> Matrix
zeroes ( rows, columns ) =
fromArray ( rows, columns ) <| Array.repeat (rows * columns) 0.0

{-| Create an nxn identity matrix.

identity =
Matrix.eye 3

-}
eye : Int -> Matrix
eye diagonal =
let
gen x =
if x % (diagonal + 1) == 0 then
1
else
0
in
Array.initialize (diagonal * diagonal) gen
|> fromArray ( diagonal, diagonal )

{-| Create an nxn upper triangular matrix.

triangle =
Matrix.upper 4

-}
upper : Int -> Matrix
upper diagonal =
let
range =
List.range 0 (diagonal - 1)

f i =
List.append (List.repeat i 0.0) (List.repeat (diagonal - i) 1.0)

list =
List.map f range
in
from2DList list

{-| Create an nxn strict upper triangular matrix.
This means that elements along the diagonal are zero.

striangle =
Matrix.strictUpper 4

-}
strictUpper : Int -> Matrix
strictUpper diagonal =
map2 (-) (upper diagonal) (eye diagonal)

{-| Create an nxn lower triangular matrix.

ltriangle =
Matrix.lower 4

-}
lower : Int -> Matrix
lower diagonal =
transpose <| upper diagonal

{-| Create an nxn strict lower triangular matrix.
This means that elements along the diagonal are zero.

sltriangle =
Matrix.strictLower 4

-}
strictLower : Int -> Matrix
strictLower diagonal =
map2 (-) (transpose <| upper diagonal) (eye diagonal)

-- Operations

{-| Perform matrix multiplication

A * B

a =
Matrix.eye (2, 2)

b =
Matrix.mats "[2, 3; 4 5]"

c =
mul a b

-}
mul : Matrix -> Matrix -> Matrix
mul a b =
forwardError2 "[in mul]" mulList a b

mulArray : Matnxn -> Matnxn -> Matrix
mulArray a_ b_ =
Err "Not Implemented"

{-| Implementation of multiplication using lists.
Probably slow.
-}
mulList : Matnxn -> Matnxn -> Matrix
mulList a_ b_ =
if numColumns a_ /= numRows b_ then
let
acolumns =
Basics.toString <| numColumns a_

brows =
Basics.toString <| numRows b_
in
Err <| "Dimension mismatch in a*b: a.columns = " ++ acolumns ++ " b.rows = " ++ brows ++ "."
else
let
a_list =
to2DListBase a_

b_list =
to2DList <| transposeBase b_

--collapse a row from the left and column from the right into scalar
collapse x y =
List.sum <| List.map2 (*) x y

--element level
constructList x y =
case x of
[] ->
[]

m :: ms ->
List.map (collapse m) y :: constructList ms y
in
from2DList <| constructList a_list b_list

{-| Get the lu decomposition of a matrix
Since pivoting isn't implemented, watch out for numerical instability
-}
luDecomp : Matrix -> ( Matrix, Matrix )
luDecomp a =
case a of
Mat a_ ->
if numRows a_ == numColumns a_ then
luSplit a_
else
( Err "Must be a square matrix", Err "" )

Err string ->
( Err string, Err "" )

{-| Splits the lu factorized single matrix into two
-}
luSplit : Matnxn -> ( Matrix, Matrix )
luSplit a =
let
single =
luNoPivotSingle a

dim =
numColumns a

l =
eMul single (strictLower dim)

u =
eMul single (upper dim)
in
( l, u )

{-| Performs lu factorization
-}
luNoPivotSingle : Matnxn -> Matrix
luNoPivotSingle a =
let
lu =
zeroes a.dimensions

--i is inner loop j outer loop
-- bind computelem to a
luCompute index dest =
luComputeElem index a dest

indices =
genIndices a.dimensions
in
case lu of
Mat lu_ ->
Mat <| List.foldl luCompute lu_ indices

_ ->
lu

genIndices : ( Int, Int ) -> List ( Int, Int )
genIndices ( i_range, j_range ) =
let
listj =
List.concat <| List.map (List.repeat i_range) <| List.range 1 j_range

listi =
List.concat <| List.repeat j_range (List.range 1 i_range)
in
List.map2 (,) listi listj

luComputeElem : ( Int, Int ) -> Matnxn -> Matnxn -> Matnxn
luComputeElem ( i, j ) original lu =
let
tiny =
10 ^ (-40)

getWithDefault index mat =
Maybe.withDefault 0.0 (getBase index mat)

aij =
getWithDefault ( i, j ) original

bjj =
Maybe.withDefault tiny <| getBase ( j, j ) lu

compute index =
getWithDefault ( i, index ) lu * getWithDefault ( index, j ) lu
in
if i > j then
let
k =
List.range 1 (j - 1)
in
-- NRIC 2.3.13
setBase ( i, j ) ((aij - (List.sum <| List.map compute k)) / bjj) lu
else
let
k =
List.range 1 (i - 1)
in
-- NRIC 2.3.12
setBase ( i, j ) (aij - (List.sum <| List.map compute k)) lu

getScaling : List (List Float) -> List Float
getScaling mat =
List.map (\x -> List.map abs x) mat
|> List.map (List.maximum)
|> List.map (Maybe.withDefault 0)
|> List.map (\x -> 1 / x)

{-| Get an item at index (row, column). Indices are 1-indexed.

a =
Matrix.mats "[3 4; 5 6]"

b =
-- equals 4
get ( 1, 2 ) a

-}
get : ( Int, Int ) -> Matrix -> Maybe Float
get indices a =
forwardErrorF "[in get]" (getBase indices) a

getBase : ( Int, Int ) -> Matnxn -> Maybe Float
getBase ( r_index, c_index ) a_ =
let
check_r_bounds =
0 < r_index && r_index <= numRows a_

check_c_bounds =
0 < c_index && c_index <= numColumns a_
in
if check_r_bounds && check_c_bounds then
Array.get ((r_index - 1) * (numColumns a_) + c_index - 1)
a_.elements
else
Nothing

{-| Get an item at index (row, column). Indices are 1-indexed.
-}
set : ( Int, Int ) -> Float -> Matrix -> Matrix
set indices data a =
forwardError "[in set]" (\x -> (Mat <| setBase indices data x)) a

setBase : ( Int, Int ) -> Float -> Matnxn -> Matnxn
setBase ( r_index, c_index ) data a_ =
let
check_r_bounds =
0 < r_index && r_index <= numRows a_

check_c_bounds =
0 < c_index && c_index <= numColumns a_
in
if check_r_bounds && check_c_bounds then
{ a_
| elements =
Array.set ((r_index - 1) * (numColumns a_) + c_index - 1)
data
a_.elements
}
else
a_

{-| Add two matrices of identical dimensions together

a =
Matrix.add (Matrix.zeroes ( 2, 2 )) (Matrix.ones ( 2, 2 ))

-}
add : Matrix -> Matrix -> Matrix

addBase : Matnxn -> Matnxn -> Matrix
if equalSize a_ b_ then
fromArray ( numRows a_, numColumns a_ ) <| arraymap2 (+) a_.elements b_.elements
else
let
dimToString a_

bdims =
dimToString b_
in
Err <| "Matrices not equal size. a: " ++ adims ++ ", b: " ++ bdims

{-| Map a function over all elements individually
-}
map : (Float -> Float) -> Matrix -> Matrix
map f a =
forwardError "[in map]"
(\a_ -> fromArray ( numRows a_, numColumns a_ ) <| Array.map f a_.elements)
a

{-| Map a function over elements of same index between matrices
-}
map2 : (Float -> Float -> Float) -> Matrix -> Matrix -> Matrix
map2 f a b =
forwardError2 "[in map2]" (map2Base f) a b

map2Base : (Float -> Float -> Float) -> Matnxn -> Matnxn -> Matrix
map2Base f a b =
if equalSize a b then
fromArray a.dimensions <| arraymap2 f a.elements b.elements
else
Err "Unequal Matrix sizes"

{-| Perform scalar multiplication on a matrix.
-}
sMul : Float -> Matrix -> Matrix
sMul a b =
map ((*) a) b

{-| Perform element by element multiplication on a matrix.
-}
eMul : Matrix -> Matrix -> Matrix
eMul a b =
map2 (*) a b

{-| Perform scalar division on a matrix.
-}
sDiv : Float -> Matrix -> Matrix
sDiv a b =
sMul (1 / a) b

{-| Invert a square matrix.

a =
"[ 2 5; 6 7]"

inva =
invert a

identity =
mul a inva

-}
invert : Matrix -> Matrix
invert a =
case a of
Mat a_ ->
solve a (eye (numRows a_))

Err string ->
Err string

{-| Shorthand for invert.
-}
inv : Matrix -> Matrix
inv =
invert

{-| Solve a system of equations of the form
AX = B. You provide A and b and get back x
Where A, B, X are matrices
B is a matrix of solution vectors horizontally concatenated.

a =
Matrix.eye 3

b =
Matrix.hcat (Matrix.vec [ 3, 2, 1 ]) (Matrix.vec [ 1, 2, 3 ])

x =
Matrix.solve a b

-}
solve : Matrix -> Matrix -> Matrix
solve a b =
case ( a, b ) of
( Mat a_, Mat b_ ) ->
if numRows a_ == numRows b_ then
let
bs =
getColumns b

single =
luNoPivotSingle a_

boundApply vec =
case vec of
Mat vec_ ->
applySubstitution single vec_

_ ->
vec
in
List.foldr hcat (fromList ( numRows a_, 0 ) []) <| List.map boundApply bs
else
Err "Dimension mismatch"

( _, _ ) ->
Err "One of the inputs was malformed"

{-| Solve a system of equations of the form
Ax = b. You provide A and b and get back x
Where A is a matrix, and b and x are vectors

a =
Matrix.eye 3

b =
Matrix.vec [ 3, 2, 1 ]

x =
Matrix.solve a b

-}
solveV : Matrix -> Matrix -> Matrix
solveV a b =
forwardError2 "[in solveV]" solveVBase a b

solveVBase : Matnxn -> Matnxn -> Matrix
solveVBase a b =
let
b_is_vector =
numColumns b == 1

dimensions_match =
numRows b == numColumns a

a_is_square =
numColumns a == numRows a
in
case ( b_is_vector, dimensions_match, a_is_square ) of
( True, True, True ) ->
applySubstitution (luNoPivotSingle a) b

( False, _, _ ) ->
Err "b is not a vector."

( _, False, _ ) ->
Err "Dimensions of A and b do not match"

( _, _, False ) ->
Err "A is not square"

{- Applies forward and backward substitution, decoupling substitution from
computing lu decomp
-}

applySubstitution : Matrix -> Matnxn -> Matrix
applySubstitution single b =
let
brows =
numRows b

z =
zeroes ( brows, 1 )
in
case ( z, single ) of
( Mat z_, Mat single_ ) ->
let
is =
List.range 1 (numRows b)

fsBound i y =
forwardSubstitution i single_ b y

y =
List.foldl fsBound z_ is

bsBound i x =
backSubstitution i single_ y x

x =
List.foldr bsBound z_ is
in
Mat x

_ ->
z

{-| Perform forward substitution remembering that the diagonal of the lower lu
matrix is all ones.
l is technically the lower lu matrix
b is the original passed in vector of b's
y is the solution, updated piecewise, of L * y = b

Numerical Recipes in C 2.3.6

-}
forwardSubstitution : Int -> Matnxn -> Matnxn -> Matnxn -> Matnxn
forwardSubstitution i l b y =
let
getWithDefault index mat =
Maybe.withDefault 0.0 (getBase index mat)

getVecDefault index vec =
getWithDefault ( index, 1 ) vec

yi =
case i of
1 ->
getVecDefault i b

_ ->
let
bi =
getVecDefault i b

j =
List.range 1 (i - 1)

lijs =
List.map (\curr_j -> getWithDefault ( i, curr_j ) l) j

yjs =
List.map (\curr_j -> getVecDefault curr_j y) j

sum =
List.sum <| List.map2 (*) lijs yjs
in
bi - sum
in
setBase ( i, 1 ) yi y

{-| Perform forward substitution remembering that the diagonal of the lower lu
matrix is all ones.
u is technically the upper lu matrix
y is the original passed in vector of b's
x is the solution, updated piecewise of BOTH Ux = y and the original Ax = b

Numerical Recipes in C 2.3.7

-}
backSubstitution : Int -> Matnxn -> Matnxn -> Matnxn -> Matnxn
backSubstitution i u y x =
let
getWithDefault index mat =
Maybe.withDefault 0.0 (getBase index mat)

getVecDefault index vec =
getWithDefault ( index, 1 ) vec

uii =
getWithDefault ( i, i ) u

xrows =
numRows x

xi =
if i == xrows then
(getVecDefault i y) / (uii)
else
let
--aij = getWithDefault ( i, j ) original
yi =
getVecDefault i y

j =
List.range (i + 1) (xrows)

uijs =
List.map (\curr_j -> getWithDefault ( i, curr_j ) u) j

xjs =
List.map (\curr_j -> getVecDefault curr_j x) j

sum =
List.sum <| List.map2 (*) uijs xjs
in
(yi - sum) / uii
in
setBase ( i, 1 ) xi x

{-| Transpose a matrix.
-}
transpose : Matrix -> Matrix
transpose a =
forwardError "[in transpose]" transposeBase a

transposeBase : Matnxn -> Matrix
transposeBase a_ =
let
f index =
let
mappedindex =
(index % numRows a_)
* (numColumns a_)
+ (index // numRows a_)
in
Maybe.withDefault 0.0 (Array.get mappedindex a_.elements)
in
Array.initialize (numColumns a_ * numRows a_) f
|> fromArray ( numColumns a_, numRows a_ )

{-| Get the determinant of a square matrix.

a =
Matrix.mats "[1 2 3; 4 5 6; 7 8 9]"

is_singular =
if (determinant a) == 0 then
"Matrix is singular"
else
"Matrix is not singular"

-}
determinant : Matrix -> Maybe Float
determinant a =
let
detBase a_ =
if numRows a_ == numColumns a_ then
let
single =
luNoPivotSingle a_
in
List.range 1 (numRows a_)
|> List.map (\x -> ( x, x ))
|> List.map (\x -> Maybe.withDefault 0.0 (get x single))
|> List.product
|> Just
else
Nothing
in
forwardErrorF "[in determinant]" detBase a

{-| Shorthand for determinant.
-}
det : Matrix -> Maybe Float
det =
determinant

{-| Performs the dot product of two nxn vectors

x =
Matrix.vec [ 1, 0, 0 ]

y =
Matrix.vec [ 0, 1, 0 ]

zero =
Matrix.dot x y

-}
dot : Matrix -> Matrix -> Maybe Float
dot a b =
case ( a, b ) of
( Mat a_, Mat b_ ) ->
let
a_is_vector =
numColumns a_ == 1

b_is_vector =
numColumns b_ == 1

same_length =
(numRows a_) == (numRows b_)
in
if a_is_vector && b_is_vector && same_length then
arraymap2 (*) a_.elements b_.elements
|> Array.foldr (+) 0
|> Just
else
Nothing

( _, _ ) ->
Nothing

-- is there a way to do good scalar error handling?

{-| Get the cross product of two 3d vectors. a >< b

x =
Matrix.vec [ 1, 0, 0 ]

y =
Matrix.vec [ 0, 1, 0 ]

z =
Matrix.cross x y

-}
cross : Matrix -> Matrix -> Matrix
cross a b =
forwardError2 "[in cross]" crossBase a b

crossBase : Matnxn -> Matnxn -> Matrix
crossBase a_ b_ =
if check3dVec a_ && check3dVec b_ then
let
getDefault x vec =
Maybe.withDefault 0.0 <| getBase x vec

ax =
getDefault ( 1, 1 ) a_

ay =
getDefault ( 2, 1 ) a_

az =
getDefault ( 3, 1 ) a_

bx =
getDefault ( 1, 1 ) b_

by =
getDefault ( 2, 1 ) b_

bz =
getDefault ( 3, 1 ) b_
in
let
i =
(ay * bz) - (by * az)

j =
-1 * ((ax * bz) - (bx * az))

k =
(ax * by) - (bx * ay)
in
vec <| [ i, j, k ]
else
Err "One or both vectors are malformed."

{-| Checks if two matrices are of equal size
-}
equalSize : Matnxn -> Matnxn -> Bool
equalSize a b =
(numRows a == numRows b) && (numColumns a == numColumns b)

{-| Checks if two matrices are equivalent within some epsilon.

epsilon =
10 ^ (-4)

is_equivalent =
equivalent epsilon a b

-}
equivalent : Float -> Matrix -> Matrix -> Bool
equivalent epsilon a b =
case ( a, b ) of
( Mat a_, Mat b_ ) ->
let
equal_size =
equalSize a_ b_

equal_members =
arraymap2 (-) a_.elements b_.elements
|> Array.map (\x -> abs x < epsilon)
|> Array.toList
|> List.all ((==) True)
in
if equal_size && equal_members then
True
else
False

_ ->
False

{-| Concatenate two matrices vertically.
-}
vcat : Matrix -> Matrix -> Matrix
vcat a b =
forwardError2 "[in vcat]" vcatBase a b

vcatBase : Matnxn -> Matnxn -> Matrix
vcatBase a_ b_ =
let
acols =
numColumns a_

bcols =
numColumns b_

arows =
numRows a_

brows =
numRows b_
in
if acols == bcols then
fromArray ( arows + brows, acols ) (Array.append a_.elements b_.elements)
else
Err <|
"Number of columns are not equal: a: "
++ Basics.toString acols
++ " b: "
++ Basics.toString bcols

{-| Concatenate two matrices horizontally.
-}
hcat : Matrix -> Matrix -> Matrix
hcat a b =
forwardError2 "[in hcat]" hcatBase a b

hcatBase : Matnxn -> Matnxn -> Matrix
hcatBase a b =
let
arows =
numRows a

brows =
numRows b
in
if arows == brows then
vcat (transposeBase a) (transposeBase b)
|> transpose
else
Err <|
"Number of rows are not equal: a: "
++ Basics.toString arows
++ " b: "
++ Basics.toString brows

{-| Returns matrix as flat list

Matrix.toFlatList (Matrix.eye 2) == [ 1, 0, 0, 1 ]

-}
toFlatList : Matrix -> List Float
toFlatList n =
case n of
Mat n_ ->
Array.toList n_.elements

_ ->
[]

{-| Returns matrix as 2d list.
Returns empty list if Matrix is in error.

Matrix.to2DList (Matrix.eye 2) == [ [1, 0], [0, 1] ]

-}
to2DList : Matrix -> List (List Float)
to2DList n =
case n of
Mat n_ ->
to2DListBase n_

_ ->
[ [] ]

to2DListBase : Matnxn -> List (List Float)
to2DListBase z =
make2D (numColumns z) (Array.toList z.elements)

{-| Returns size of matrix, (rows, columns)
-}
size : Matrix -> ( Int, Int )
size n =
case n of
Mat n_ ->
n_.dimensions

_ ->
( 0, 0 )

-- Auxiliary Functions

{-| Helper to catch errors in functions of two variables (Matrix -> Matrix) ->
Matrix
-}
forwardError2 : String -> (Matnxn -> Matnxn -> Matrix) -> Matrix -> Matrix -> Matrix
forwardError2 error f a b =
case ( a, b ) of
( Err string, Mat _ ) ->
Err <| "\n" ++ error ++ " Matrix a: " ++ string

( Mat _, Err string ) ->
Err <| "\n" ++ error ++ " Matrix b: " ++ string

( Err string, Err string2 ) ->
Err <| "\n" ++ error ++ " Matrix a: " ++ string ++ "\n Matrix b: " ++ string2

( Mat a_, Mat b_ ) ->
f a_ b_

forwardError : String -> (Matnxn -> Matrix) -> Matrix -> Matrix
forwardError error f a =
case a of
Mat a_ ->
f a_

Err string ->
Err <| "\n" ++ error ++ string

forwardErrorF : String -> (Matnxn -> Maybe Float) -> Matrix -> Maybe Float
forwardErrorF error f a =
case a of
Mat a_ ->
f a_

Err string ->
Nothing

{-| Helper to take a Matrix and stringify its dimensions
-}
dimToString : Matnxn -> String
dimToString a =
let
arows =
Basics.toString <| numRows a

acols =
Basics.toString <| numColumns a
in
"(" ++ arows ++ "," ++ acols ++ ")"

{-| Change matrix into string form, such as what would be displayed in the terminal.

Matrix.toString (Matrix.eye 3) == " 1 0 0\n 0 1 0\n 0 0 1"

-}
toString : Matrix -> String
toString a =
case a of
Mat mat ->
toStringBasic mat

Err string ->
string

{-| Change correctly formed matrix into string form
-}
toStringBasic : Matnxn -> String
toStringBasic a =
let
strings =
a.elements
|> Array.toList
|> List.map Basics.toString
|> List.map ((++) " ")

structured_strings =
make2D (numColumns a) strings

matrix_string =
structured_strings
|> List.intersperse [ "\n" ]
|> List.concat
|> List.foldr (++) ""
in
matrix_string

{-| Helper to re-2dify a flat matrix
-}
make2D : Int -> List a -> List (List a)
make2D num_row_elem list =
case list of
[] ->
[]

items ->
if List.length list < num_row_elem then
[ list ]
else
List.take num_row_elem list :: make2D num_row_elem (List.drop num_row_elem list)

{-| Returns the rows of a matrix in a list.
-}
getRows : Matrix -> List Matrix
getRows a =
List.map rvec <| to2DList a

{-| Returns the columns of a matrix in a list.
-}
getColumns : Matrix -> List Matrix
getColumns a =
List.map vec <| to2DList (transpose a)

{-| Helper to debug print. Most useful in repl.

> Matrix.debugPrint (Matrix.eye 3)
(3, 3) Matrix
1 0 0
0 1 0
0 0 1
: ""
""

-}
debugPrint : Matrix -> String
debugPrint a =
let
description_string =
Basics.toString (size a) ++ " Matrix\n"

final_string =
description_string ++ toString a ++ "\n"
in
Debug.log (final_string) ""

{-| Get the number of rows in a matrix
-}
numRows : Matnxn -> Int
numRows a =
Tuple.first a.dimensions

{-| Get the number of columns in a matrix
-}
numColumns : Matnxn -> Int
numColumns a =
Tuple.second a.dimensions

{-| Check if the matrix is a 3d vector
-}
check3dVec : Matnxn -> Bool
check3dVec a =
if numRows a == 3 && numColumns a == 1 then
True
else
False

{-| Map2 for arrays
-}
arraymap2 : (a -> b -> c) -> Array.Array a -> Array.Array b -> Array.Array c
arraymap2 f a b =
if Array.isEmpty a || Array.isEmpty b then
Array.fromList []
else
let
dropArr a =
Array.slice 1 (Array.length a) a

result =
case ( Array.get 0 a, Array.get 0 b ) of
( Just aval, Just bval ) ->
Array.fromList [ f aval bval ]

_ ->
Array.fromList []
in
Array.append result (arraymap2 f (dropArr a) (dropArr b))

--string must be enclosed by brackets

matParser : String -> List (List Float)
matParser string =
if String.startsWith "[" string && String.endsWith "]" string then
String.dropLeft 1 string
|> String.dropRight 1
|> String.split ";"
|> List.map (String.split " ")
|> List.map (List.filter (\x -> not <| String.isEmpty x))
|> List.map (List.map (\x -> Result.withDefault 0 (String.toFloat x)))
else
[]

-- Operators for convenience

{-| Matrix multiply
-}
(**) : Matrix -> Matrix -> Matrix
(**) =
mul