%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}{1.001} % not yet

\startmodule[matrix]

\startluacode

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|" },
}

-- 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.startmatrix(whatever)
            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.stopmatrix()
    elseif m then
        context(m)
    end
end

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

-- function matrix.swaprows(t,i,j)
--     if i <= #t and j <= #t then
--         t[i], t[j] = t[j], t[i]
--         return t
--     else
--         return "error: out of bound"
--     end
-- end

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

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

-- function matrix.swapcolumns(t, i, j)
--     if i <= #t[1] and j <= #t[1] then
--         for k = 1, #t do
--             t[k][i], t[k][j] = t[k][j], t[k][i]
--         end
--         return t
--     else
--         return "error: out of bound"
--     end
-- end

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 then
        return "error: no row i"
    end
    if j <= n then
        return "error: no row j"
    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.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

function matrix.make(n,m,low,high)
    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 = 100
    end
    local t = { }  -- make an empty n1 x n2 array
    local again = true
    for i=1,n do
        t[i] = { }
    end
    while true do
        for i=1,n do
            local ti = t[i]
            for j=1,m do
                ti[j] = random(low,high)
            end
        end
        if n ~= m or determinant(t,1) ~= 0 then
            return t
        end
    end
end

-- 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

-- calculating determinant using Laplace Expansion

function matrix.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

--  solve the linear equation m X = c

local function solve(m,c)
    local n = #m
    if n ~= #c then
        return copy(m)
    end
    local newm = copy(m)
    local temp = copy(c)
    for i=1,n do
        insert(newm[i],temp[i])
    end
    return rowechelon(newm,1)
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

\stopluacode

\stopmodule

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

\continueifinputfile{m-matrix.mkiv}

\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,  2, 4 },
    { 2, 0, -6,  9, 7 },
    { 2, 2, -6,  2, 4 },
    { 2, 2, -6,  2, 4 },
}
\stopluacode

\startsubject[title={A symbolic matrix}]

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

\stopsubject

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

\startluacode
moduledata.matrix.typeset(document.DemoMatrixA)
context.blank()
moduledata.matrix.swap(document.DemoMatrixA, 2, 4)
context.blank()
moduledata.matrix.typeset(document.DemoMatrixA)
\stopluacode

\stopsubject

\startsubject[title={Multiply $3 \times r_2$}]

\startluacode
moduledata.matrix.typeset(document.DemoMatrixA)
context.blank()
moduledata.matrix.typeset(moduledata.matrix.multiply(document.DemoMatrixA, 2, 3))
\stopluacode

\stopsubject

\startsubject[title={Row 2 + $4 \times r_3$}]

\startluacode
moduledata.matrix.typeset(document.DemoMatrixA)
context.blank()
moduledata.matrix.sumrow(document.DemoMatrixA, 2, 3, 4)
context.blank()
moduledata.matrix.typeset(document.DemoMatrixA,{ fences = "bars" } )
\stopluacode

\stopsubject

\startsubject[title={Transpose a matrix}]

\startluacode
moduledata.matrix.typeset(document.DemoMatrixA)
context.blank()
moduledata.matrix.typeset(moduledata.matrix.transpose(document.DemoMatrixA))
\stopluacode

\stopsubject

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

\startluacode
context(moduledata.matrix.inner({ 1, 2, 3 }, { 3, 1, 2 }))
context.blank()
context(moduledata.matrix.inner({ 1, 2, 3 }, { 3, 1, 2, 4 }))
\stopluacode

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

\startluacode
moduledata.matrix.typeset(document.DemoMatrixA)
context.blank()
moduledata.matrix.typeset(moduledata.matrix.product(document.DemoMatrixA,document.DemoMatrixA))
\stopluacode

\stopsubject

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

\ctxmodulematrix{typeset(moduledata.matrix.uppertri(document.DemoMatrixB))}

\startsubject[title={A determinant}]

\startluacode
local m = {
    { 1, 2,  4 },
    { 0, 0,  2 },
    { 2, 2, -6 },
}
context(moduledata.matrix.determinant(m, "determinant=yes" ))
\stopluacode

\stopsubject

\startsubject[title={Row echelon form}]

\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.blank()
moduledata.matrix.typeset(moduledata.matrix.rowechelon(m,1), { determinant = "yes" })
\stopluacode

\stopsubject

\startsubject[title={Solving linear equation}]

\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

\stopsubject

\startsubject[title={Inverse matrix}]

\startcombination[2*1]
    {\ctxlua{moduledata.matrix.typeset { { 1, 1, 1 }, { 0, 2, 3 }, { 3, 2, 1 } }}} {}
    {\ctxlua{moduledata.matrix.typeset(moduledata.matrix.inverse { { 1, 1, 1 }, { 0, 2, 3 }, { 3, 2, 1 } })}} {}
\stopcombination

\stopsubject

\stoptext