%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 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_p = {
    left  = "\\left(\\,",
    right = "\\,\\right)",
}

local fences_b = {
    left  = "\\left[\\,",
    right = "\\,\\right]",
}

function matrix.typeset(m,options)
    local options = settings_to_hash(options or "")
    context.startmatrix(options.determinant and fences_b or fences_p)
        for i=1, #m do
            local mi = m[i]
            for j=1,#mi do
                context.NC(mi[j])
            end
            context.NR()
        end
    context.stopmatrix()
end

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

function matrix.swap(t,i,j)
    t[i], t[j] = t[j], t[i]
end

-- 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 0
    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)
    local product = { }
    if #m1[1] == #m2 then
        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
    end
    return product
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

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

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

--  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 + $3 \times r_4$}]

\startluacode
moduledata.matrix.typeset(document.DemoMatrixA)
context.blank()
moduledata.matrix.sumrow(document.DemoMatrixA, 2, 3, 4)
context.blank()
moduledata.matrix.typeset(document.DemoMatrixA)
\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))
\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)
moduledata.matrix.typeset(moduledata.matrix.rowechelon(m,1))
\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))
\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