%D \module
%D   [       file=m-matrix,
%D        version=2014.11.04, % already a year older
%D          title=\CONTEXT\ Extra Modules,
%D       subtitle=Matrices,
%D         author={Jeong Dalyoung \& Hans Hagen},
%D           date=\currentdate,
%D      copyright={PRAGMA ADE \& \CONTEXT\ Development Team}]
%C
%C This module is part of the \CONTEXT\ macro||package and is
%C therefore copyrighted by \PRAGMA. See mreadme.pdf for
%C details.

%D This code is based on a post by Dalyoung on the context list. After that
%D we turned it into a module and improved the code a bit. Feel free to ask
%D us for more. Once we're satisfied, a more general helper l-matrix could
%D be made. Dalyoung does the clever bits, and Hans only cleanes up and
%D optimizes a bit.

% \registerctxluafile{l-matrix}{} % not yet

\startmodule[matrix]

\startluacode

local tonumber, type = tonumber, type

local settings_to_hash = utilities.parsers.settings_to_hash
local formatters       = string.formatters
local copy             = table.copy
local insert           = table.insert
local remove           = table.remove
local random           = math.random

local context          = context

local matrix      = { }
moduledata.matrix = matrix

local f_matrix_slot = formatters["%s_{%s%s}"]

function matrix.symbolic(sym, x, y, nx ,ny) -- symMatrix("a", "m", "n")
    local nx = nx or 2
    local ny = ny or nx
    local function filled(i,y)
        local mrow = { }
        for j=1,nx do
            mrow[#mrow+1] = f_matrix_slot(sym,i,j)
        end
        mrow[#mrow+1] = "\\cdots"
        mrow[#mrow+1] = f_matrix_slot(sym,i,y)
        return mrow
    end
    local function dummy()
        local mrow = { }
        for j=1,nx do
            mrow[#mrow+1] = "\\vdots"
        end
        mrow[#mrow+1] = "\\ddots"
        mrow[#mrow+1] = "\\vdots"
        return mrow
    end
    --
    local mm = { }
    for i=1,ny do
        mm[i] = filled(i,y)
    end
    mm[#mm+1] = dummy()
    mm[#mm+1] = filled(x,y)
    return mm
end

-- todo: define a matrix at the tex end so that we have more control

-- local fences = {
--     parentheses = { left = "\\left(\\,", right = "\\,\\right)" },
--     brackets    = { left = "\\left[\\,", right = "\\,\\right]" },
--     bars        = { left = "\\left|\\,", right = "\\,\\right|" },
-- }

local fences = {
    parentheses = { "matrix:parentheses" },
    brackets    = { "matrix:brackets" },
    bars        = { "matrix:bars" },
}

-- one can add more fences

fences.bar         = fences.bars
fences.parenthesis = fences.parentheses
fences.bracket     = fences.brackets

-- one can set the template

matrix.template    = "%0.3F"

function matrix.typeset(m,options)
    if type(m) == "table" then
        local options  = settings_to_hash(options or "")
        local whatever = options.determinant == "yes" and fences.bars or fences.parentheses
        if options.fences then
            whatever = fences[options.fences] or whatever
        elseif options.determinant then
         -- whatever = fences.brackets
            whatever = fences.bars
        end
        local template = options.template or matrix.template
        if template == "yes" then
            template = matrix.template
        elseif template == "no" then
            template = false
        elseif tonumber(template) then
            template = "%0." .. template .. "F"
        end
        context.startnamedmatrix(whatever)
            if type(m[1]) ~= "table" then
                m = { copy(m) }
            end
            for i=1,#m do
                local mi = m[i]
                for j=1,#mi do
                    context.NC()
                    local n = mi[j]
                    if template and tonumber(n) then
                        context(template,n)
                    else
                        context(mi[j])
                    end
                end
                context.NR()
            end
        context.stopnamedmatrix()
    elseif m then
        context(m)
    end
end

function matrix.swaprows(t,i,j)
    local ti = t[i]
    if not ti then
        return "error: no row"
    end
    local tj = t[j]
    if not tj then
        return "error: no row"
    end
    t[i], t[j] = tj, ti
    return t
end

-- interchange two columns (i-th, j-th)

function matrix.swapcolumns(t, i, j)
    local t1 = t[1]
    if not t1 then
        return "error: no rows"
    end
    local n = #t1
    if i > n or j > n then
        return "error: no column"
    end
    for k = 1, #t do
        local tk = t[k]
        tk[i], tk[j] = tk[j], tk[i]
    end
    return t
end

matrix.swapC       = matrix.swapcolumns
matrix.swapR       = matrix.swaprows
matrix.swapcolumns = matrix.swapcolumns
matrix.swaprows    = matrix.swaprows
matrix.swap        = matrix.swaprows

-- replace i-th row with factor * (i-th row)

function matrix.multiply(m,i,factor)
    local mi = m[i]
    for k=1,#mi do
        mi[k] = factor * mi[k]
    end
    return m
end

-- scalar product "factor * m"

function matrix.scalar(m, factor)
    for i=1,#m do
        local mi = m[i]
        for j=1,#mi do
            mi[j] = factor * mi[j]
        end
    end
    return m
end

-- replace i-th row with i-th row + factor * (j-th row)

function matrix.sumrow(m,i,j,factor)
    local mi = m[i]
    local mj = m[j]
    for k=1,#mi do
        mi[k] = mi[k] + factor * mj[k]
    end
end

-- transpose of a matrix

function matrix.transpose(m)
    local t = { }
    for j=1,#m[1] do
        local r = { }
        for i=1,#m do
            r[i] = m[i][j]
        end
        t[j] = r
    end
    return t
end

-- inner product of two vectors

function matrix.inner(u,v)
    local nu = #u
    if nu == 0 then
        return 0
    end
    local nv = #v
    if nv ~= nu then
       return "error: size mismatch"
    end
    local result = 0
    for i=1,nu do
        result = result + u[i] * v[i]
    end
    return result
end

-- product of two matrices

function matrix.product(m1,m2)
    if #m1[1] == #m2 then
        local product = { }
        for i=1,#m1 do
            local m1i  = m1[i]
            local mrow = { }
            for j=1,#m2[1] do
                local temp = 0
                for k=1,#m1[1] do
                    temp = temp + m1i[k] * m2[k][j]
                end
                mrow[j] = temp
            end
            product[i] = mrow
        end
        return product
    else
        return "error: size mismatch"
    end
end

local function uppertri(m,sign)
    local temp = copy(m)
    for i=1,#temp-1 do
        local pivot = temp[i][i]
        if pivot == 0 then
            local pRow = i +1
            while temp[pRow][i] == 0 do
                pRow = pRow + 1
                if pRow > #temp then -- if there is no nonzero number
                    return temp
                end
            end
            temp[i], temp[pRow] = temp[pRow], temp[i]
            if sign then
                sign = -sign
            end
        end
        local mi = temp[i]
        for k=i+1, #temp do
            local factor = -temp[k][i]/mi[i]
            local mk = temp[k]
            for l=i,#mk do
                mk[l] = mk[l] + factor * mi[l]
            end
        end
    end
    if sign then
        return temp, sign
    else
        return temp
    end
end

matrix.uppertri = uppertri

local function determinant(m)
    if #m == #m[1] then
        local d = 1
        local t, s = uppertri(m,1)
        for i=1,#t do
            d = d * t[i][i]
        end
        return s*d
    else
        return "error: not a square matrix"
    end
end

matrix.determinant = determinant

local function rowechelon(m,r)
    local temp = copy(m)
    local pRow = 1
    local pCol = 1
    while pRow <=  #temp  do
        local pivot = temp[pRow][pCol]
        if  pivot == 0 then
            local i = pRow
            local n = #temp
            while temp[i][pCol] == 0 do
                i = i + 1
                if i > n then
                    -- no nonzero number in a column
                    pCol = pCol + 1
                    if pCol > #temp[pRow] then
                        -- there is no nonzero number in a row
                        return temp
                    end
                    i = pRow
                end
            end
            temp[pRow], temp[i] = temp[i], temp[pRow]
        end
        local row = temp[pRow]
        pivot = row[pCol]
        for l=pCol,#row do
            row[l] = row[l]/pivot
        end

        if r == 1 then
            -- make the "reduced row echelon form"
            local row = temp[pRow]
            for k=1,pRow-1 do
                local current = temp[k]
                local factor  = -current[pCol]
                local mk      = current
                for l=pCol,#mk do
                    mk[l] = mk[l] + factor * row[l]
                end
            end
        end
        -- just make the row echelon form
        local row = temp[pRow]
        for k=pRow+1, #temp do
            local current = temp[k]
            local factor  = -current[pCol]
            local mk      = current
             for l=pCol,#mk do
                mk[l] = mk[l] + factor * row[l]
            end
        end
        pRow = pRow + 1
        pCol = pCol + 1

        if pRow > #temp or pCol > #temp[1] then
            pRow = #temp + 1
        end
    end
    return temp
end

matrix.rowechelon = rowechelon
matrix.rowEchelon = rowechelon

-- make matrices until its determinant is not 0

local function make(m,n,low,high) -- m and n swapped
    if not n then
        n = 10
    end
    if not m then
        m = 10
    end
    if not low then
        low = 0
    end
    if not high then
        high = 10
    end
    local t = { }
    for i=1,m do
        t[i] = { }
    end
    while true do
        for i=1,m do
            local ti = t[i]
            for j=1,n do
                ti[j] = random(low,high)
            end
        end
        if n ~= m or determinant(t,1) ~= 0 then
            return t
        end
    end
end

matrix.make  = make
matrix.makeR = matrix.make

-- extract submatrix by using

local function submatrix(t,i,j)
    local rows    = #t
    local columns = #t[1]
    local sign    = 1 -- not used
    if i <= rows and j <= columns then
        local c = copy(t)
        remove(c,i)
        for k=1,rows-1 do
            remove(c[k],j)
        end
        return c
    else
        return "error: out of bound"
    end
end

matrix.submatrix = submatrix
matrix.subMatrix = submatrix

-- calculating determinant using Laplace Expansion

local function laplace(t) -- not sure if this is the most effient but
    local factors = { 1 }  -- it's not used for number crunching anyway
    local data    = copy(t)
    local det     = 0
    while #data > 0 do
        local mat = { }
        local siz = #data[1]
        if siz == 0 then
            return "error: no determinant"
        elseif siz == 1 then
            det = data[1][1]
            return det
        end
        for i=1,siz do
            mat[i] = data[1]
            remove(data,1)
        end
        local factor = remove(factors,1)
        local m1 = mat[1]
        if siz == 2 then
            local m2 = mat[2]
            det = det + factor * (m1[1]*m2[2] - m1[2]*m2[1])
        else
            for j=1,#m1 do
                local m1j = m1[j]
                if m1j ~= 0 then
                    insert(factors, (-1)^(j+1) * factor * m1j)
                    local m = submatrix(mat,1,j)
                    for k, v in next, m do
                         insert(data,v)
                    end
                end
            end
        end
    end
    return det
end

matrix.laplace = laplace

--  solve the linear equation m X = c

local function solve(m,c)
    local n = #m
    if n ~= #c then
     -- return "error: size mismatch"
        return nil
    end
    local newm = copy(m)
    local temp = copy(c)
    local solution = copy(c)
    for i=1,n do
        insert(newm[i],temp[i])
    end
    newm = uppertri(newm, 0)
    for k = n,1,-1 do
        local val = 0
        local new = newm[k]
        for j = k+1, n do
            val = val + new[j] * solution[j]
        end
        if new[k] == 0 then
         -- return "error: no unique solution"
            return nil
        else
            solution[k] =  (new[n+1] - val)/new[k]
        end
    end
    return solution
end

matrix.solve = solve

-- find the inverse matrix of m

local function inverse(m)
    local n = #m
    local temp = copy(m)
    if n ~= #m[1] then
        return temp
    end
    for i=1,n do
        for j=1,n do
            insert(temp[i],j == i and 1 or 0)
        end
    end
    temp = rowechelon(temp,1)
    for i=1,n do
        for j=1,n do
            remove(temp[i], 1)
        end
    end
    return temp
end

matrix.inverse = inverse

-- create zero and identity matrix

local function makeM(k,v)
    local tt = { }
    for i=1,k do
        local temp = { }
        for j=1,k do
            temp[j] = 0
        end
        tt[i] = temp
    end
    if v and v > 0 then
        for i=1,k do
            tt[i][i] = 1
        end
    end
    return tt
end

matrix.makeM        = makeM
matrix.makeidentity = makeM
matrix.makezero     = makeM

-- append the rows of the second matrix  to the bottom of the first matrix

local function joinrows(t1, t2)
    local nt1 = #t1
    local nt2 = #t2
    if (nt1*nt2 > 0) and (#t1[1] ~= #t2[1]) then
        return "error: different number of columns"
    else
        for i=1,nt2 do
            t1[nt1+i] = t2[i]
        end
        return t1
    end
end

matrix.joinrows = joinrows
matrix.joinRows = joinrows

-- append the columns of the second matrix  to the right of the first matrix

local function joincolumns(t1, t2)
    local nt1 = #t1
    local nt2 = #t2
    if nt1 == 0 then
        return t2
    end
    if nt2 == 0 then
        return t1
    end
    if nt1 ~= nt2 then
        return "error: different number of rows"
    end
    nt3 = #t2[1]
    for i=1,nt1 do
        local temp = t2[i]
        for j= 1, nt3 do
            insert(t1[i],temp[j])
        end
    end
    return t1
end

matrix.joincolumns = joincolumns
matrix.joinColumns = joincolumns

\stopluacode

\stopmodule

\unexpanded\def\ctxmodulematrix#1{\ctxlua{moduledata.matrix.#1}}

\continueifinputfile{m-matrix.mkiv}

\usemodule[m-matrix]
\usemodule[art-01]

\starttext

\startluacode
document.DemoMatrixA = {
    { 0, 2,  4, -4, 1 },
    { 0, 0,  2,  3, 4 },
    { 2, 2, -6,  2, 4 },
    { 2, 0, -6,  9, 7 },
    { 2, 3,  4,  5, 6 },
    { 6, 6, -6,  6, 6 },
}

document.DemoMatrixB = {
    { 0, 2,  4, -4, 1 },
    { 0, 0,  2,  3, 4 },
    { 2, 2, -6,  3, 4 },
    { 2, 0, -6,  9, 7 },
    { 2, 2, -6,  2, 4 },
}

document.DemoMatrixC = {
    {  3, 3, -1,  3 },
    { -1, 4,  1,  3 },
    {  5, 4,  0,  2 },
    {  2, 4,  0, -1 },
}
\stopluacode

\startbuffer[demo]
\typebuffer
\startalignment[middle]
    \dontleavehmode\inlinebuffer
\stopalignment
\stopbuffer

\setuphead[section][before={\testpage[5]\blank[2*big]}]

\startsubject[title={A symbolic matrix}]

\startbuffer
\ctxmodulematrix{typeset(moduledata.matrix.symbolic("a", "m", "n"))}
$\qquad\qquad$
\ctxmodulematrix{typeset(moduledata.matrix.symbolic("a", "m", "n", 4, 8))}
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Generate a new $m \times n$ matrix}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(moduledata.matrix.makeR(4,3, 0,5))
    context.qquad()
    context("\\qquad")
    moduledata.matrix.typeset(moduledata.matrix.makeR(5,5,-1,5))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Swap two rows (ex: 2 and 4)}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixA)
    context("$\\qquad \\Rightarrow \\qquad$")
    moduledata.matrix.typeset(moduledata.matrix.swaprows(document.DemoMatrixA,2,4))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Swap two columns (ex: 1 and 3)}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixA)
    context("$\\qquad \\Rightarrow \\qquad$")
    moduledata.matrix.typeset(moduledata.matrix.swapcolumns(document.DemoMatrixA,1, 3))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Multiply 3 to row 2($3 \times r_2$)}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixA)
    context("$\\qquad \\Rightarrow \\qquad$")
    moduledata.matrix.typeset(moduledata.matrix.multiply(document.DemoMatrixA,2,3))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Add 4 times of row 3 to row 2($r_2 + 4 \times r_3$)}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixA)
    context("$\\qquad \\Rightarrow \\qquad$")
    moduledata.matrix.sumrow(document.DemoMatrixA,2,3,4)
    moduledata.matrix.typeset(document.DemoMatrixA)
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Transpose a matrix}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixA)
    context("$\\qquad \\Rightarrow \\qquad$")
    moduledata.matrix.typeset(moduledata.matrix.transpose(document.DemoMatrixA))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={The inner product of two vectors}]

\startbuffer
\startluacode
    context("$<1,2,3> \\cdot <3,1,2> \\ =\\ $ ")
    context( moduledata.matrix.inner({ 1, 2, 3 }, { 3, 1, 2 }))
\stopluacode
\stopbuffer

\getbuffer[demo]

\startbuffer
\startluacode
context("$<1,2,3> \\cdot <3,1,2, 4> \\ =\\ $ ")
context(moduledata.matrix.inner({ 1, 2, 3 }, { 3, 1, 2, 4 }))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={The product of two matrices}]

\startbuffer
\startluacode
    context("$\\ $")
    moduledata.matrix.typeset(document.DemoMatrixB)
    context("$\\cdot$")
    moduledata.matrix.typeset(document.DemoMatrixA)
    context("$ = $")
    moduledata.matrix.typeset(moduledata.matrix.product
        (document.DemoMatrixB,document.DemoMatrixB))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={An Upper Triangular Matrix}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixB)
    context("$\\qquad \\Rightarrow \\qquad$")
    moduledata.matrix.typeset(moduledata.matrix.uppertri(document.DemoMatrixB))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Determinant: using triangulation}]

\startbuffer
\startluacode
    local m = {
        { 1, 2,  4 },
        { 0, 0,  2 },
        { 2, 2, -6 },
        { 2, 2, -6 },
    }
    moduledata.matrix.typeset(m, {fences="bars"})
    context("$\\qquad = \\qquad$")
    context(moduledata.matrix.determinant(m))
\stopluacode
\stopbuffer

\getbuffer[demo]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixC, { fences = "bars" })
    context("$\\qquad = \\qquad$")
    context(moduledata.matrix.determinant(document.DemoMatrixC))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Determinant: using Laplace Expansion}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(document.DemoMatrixC, { fences = "bars" })
    context("$\\qquad = \\qquad$")
    context(moduledata.matrix.laplace(document.DemoMatrixC))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Example of Laplace Expansion using submatrix function}]

\startbuffer
\startluacode
    local m = {
        { 1, 5, 4, 2 },
        { 5, 2, 0, 4 },
        { 2, 2, 1, 1 },
        { 1, 0, 0, 5 },
    }
    local options = {fences = "bars"}

    moduledata.matrix.typeset(m,options)
    context("\\par $=$")
    for j = 1, #m[1] do
        local mm =  moduledata.matrix.submatrix(m, 1, j)
        local factor = (-1)^(1+j) *(m[1][j])
        context("\\ ($%d$) \\cdot ", factor)
        moduledata.matrix.typeset(mm, options)
        if j < #m[1] then
            context("\\ $+$ ")
        end
    end
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Row echelon form}]

\startbuffer
\startluacode
    local m = {
        { 1, 3, -2,  0, 2,  0,  0 },
        { 2, 6, -5, -2, 4, -3, -1 },
        { 0, 0,  5, 10, 0, 15,  5 },
        { 2, 6,  0,  8, 4, 18,  6 },
    }
    moduledata.matrix.typeset(m)
    context("$\\Rightarrow$")
    moduledata.matrix.typeset(moduledata.matrix.rowechelon(m,1))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Solving linear equation}]

\startbuffer
\startluacode
    local m = {
        {  1,  3, -2,  0 },
        {  2,  0,  1,  2 },
        {  6, -5, -2,  4 },
        { -3, -1,  5, 10 },
    }

    local c = { 5, 2, 6, 8 }

    moduledata.matrix.typeset(moduledata.matrix.solve(m,c))
    context.blank()
    moduledata.matrix.typeset(moduledata.matrix.solve(m,c), { template = 6 })
    context.blank()
    moduledata.matrix.typeset(moduledata.matrix.solve(m,c), { template = "no" })
    context.blank()
    moduledata.matrix.typeset(moduledata.matrix.solve(m,c), { template = "%0.3f" })
    context.blank()
    moduledata.matrix.typeset(moduledata.matrix.solve(m,c), { template = "%0.4F" })
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={Inverse matrix}]

\startbuffer
\startluacode
    local m = {
        { 1, 1, 1 },
        { 0, 2, 3 },
        { 3, 2, 1 },
    }
    context("$A =\\quad$")
    moduledata.matrix.typeset(m)
    context("$\\qquad A^{-1} = \\quad$")
    moduledata.matrix.typeset(moduledata.matrix.inverse(m))
    context("\\blank\\ ")
    moduledata.matrix.typeset(m)
    context("$\\cdot$")
    moduledata.matrix.typeset(moduledata.matrix.inverse(m))
    context("$ = $")
    moduledata.matrix.typeset(moduledata.matrix.product(m, moduledata.matrix.inverse(m)))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={make matrices(zero, identiry, random}]

\startbuffer
\startluacode
    moduledata.matrix.typeset(moduledata.matrix.makeM(3, 0))
    context.qquad()
    moduledata.matrix.typeset(moduledata.matrix.makeM(3, 1))
    context.qquad()
    moduledata.matrix.typeset(moduledata.matrix.makeR(4,3))
\stopluacode
\stopbuffer

\getbuffer[demo]

\stopsubject

\startsubject[title={join rows, join columns}]

\startbuffer
\startluacode
    local mat1 = moduledata.matrix.makeR(3, 4)
    local mat2 = moduledata.matrix.makeR(4, 3)

    context("Appending as columns: ")
    context.blank()
    moduledata.matrix.typeset(mat1)
    context("$\\&$")
    moduledata.matrix.typeset(mat1)
    context("\\quad $\\Rightarrow$ \\quad")
    moduledata.matrix.joinColumns(mat1, mat1)
    moduledata.matrix.typeset(mat1)
    context.blank()
    context("Appending as rows: ")
    context.blank()
    moduledata.matrix.typeset(mat2)
    context("$\\&$")
    moduledata.matrix.typeset(mat2)
    context("\\quad $\\Rightarrow$ \\quad")
    moduledata.matrix.joinRows(mat2, mat2)
    moduledata.matrix.typeset(mat2)
\stopluacode
\stopbuffer

\getbuffer[demo]

\stoptext