commit b7322ea304ecf05f4dff9f230a3930168c204037 Author: entuland Date: Mon Jul 2 15:13:47 2018 +0200 Initial commit diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..dfe0770 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,2 @@ +# Auto detect text files and perform LF normalization +* text=auto diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..cce010c --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,35 @@ +LuaMatrix License +----------- + +LuaMatrix ( http://luamatrix.luaforge.net/ ) is licensed under the +same terms as Lua (MIT license) reproduced below. This means that +LuaMatrix is free software and can be used for both academic and +commercial purposes at absolutely no cost. + +For details and rationale, see http://www.lua.org/license.html . + +=============================================================================== + +Copyright (C) 2007-2010 Michael Lutz. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. + +=============================================================================== + +(end of COPYRIGHT) diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..2f42372 --- /dev/null +++ b/README.txt @@ -0,0 +1,2 @@ +https://github.com/entuland/lua-matrix is just a thin wrapper to make this library available as a mod inside the Minetest game. +Please see https://github.com/davidm/lua-matrix for the original version \ No newline at end of file diff --git a/init.lua b/init.lua new file mode 100644 index 0000000..7f4aa6a --- /dev/null +++ b/init.lua @@ -0,0 +1,8 @@ + +-- this whole mod is just a thin wrapper around https://github.com/davidm/lua-matrix + +local mod_path = minetest.get_modpath(minetest.get_current_modname()) + +matrix = dofile(mod_path .. "/matrix.lua") + +-- that's it diff --git a/matrix.lua b/matrix.lua new file mode 100644 index 0000000..83f7d57 --- /dev/null +++ b/matrix.lua @@ -0,0 +1,1251 @@ +--[[ + +LUA MODULE + + matrix v$(_VERSION) - matrix functions implemented with Lua tables + +SYNOPSIS + + local matrix = require 'matrix' + m1 = matrix{{8,4,1},{6,8,3}} + m2 = matrix{{-8,1,3},{5,2,1}} + assert(m1 + m2 == matrix{{0,5,4},{11,10,4}}) + +DESCRIPTION + + With simple matrices this script is quite useful, though for more + exact calculations, one would probably use a program like Matlab instead. + Matrices of size 100x100 can still be handled very well. + The error for the determinant and the inverted matrix is around 10^-9 + with a 100x100 matrix and an element range from -100 to 100. + + Characteristics: + + - functions called via matrix. should be able to handle + any table matrix of structure t[i][j] = value + - can handle a type of complex matrix + - can handle symbolic matrices. (Symbolic matrices cannot be + used with complex matrices.) + - arithmetic functions do not change the matrix itself + but build and return a new matrix + - functions are intended to be light on checks + since one gets a Lua error on incorrect use anyways + - uses mainly Gauss-Jordan elimination + - for Lua tables optimised determinant calculation (fast) + but not invoking any checks for special types of matrices + - vectors can be set up via vec1 = matrix{{ 1,2,3 }}^'T' or matrix{1,2,3} + - vectors can be multiplied to a scalar via num = vec1^'T' * vec2 + where num will be a matrix with the result in mtx[1][1], + or use num = vec1:scalar( vec2 ), where num is a number + +API + + matrix function list: + + matrix.add + matrix.columns + matrix.concath + matrix.concatv + matrix.copy + matrix.cross + matrix.det + matrix.div + matrix.divnum + matrix.dogauss + matrix.elementstostring + matrix.getelement + matrix.gsub + matrix.invert + matrix.ipairs + matrix.latex + matrix.len + matrix.mul + matrix.mulnum + matrix:new + matrix.normf + matrix.normmax + matrix.pow + matrix.print + matrix.random + matrix.replace + matrix.root + matrix.rotl + matrix.rotr + matrix.round + matrix.rows + matrix.scalar + matrix.setelement + matrix.size + matrix.solve + matrix.sqrt + matrix.sub + matrix.subm + matrix.tostring + matrix.transpose + matrix.type + + See code and test_matrix.lua. + +DEPENDENCIES + + None (other than Lua 5.1 or 5.2). May be used with complex.lua. + +HOME PAGE + + http://luamatrix.luaforge.net + http://lua-users.org/wiki/LuaMatrix + +DOWNLOAD/INSTALL + + ./util.mk + cd tmp/* + luarocks make + +LICENSE + + Licensed under the same terms as Lua itself. + + Developers: + Michael Lutz (chillcode) - original author + David Manura http://lua-users.org/wiki/DavidManura +--]] + +--//////////// +--// matrix // +--//////////// + +local matrix = {_TYPE='module', _NAME='matrix', _VERSION='0.2.11.20120416'} + +-- access to the metatable we set at the end of the file +local matrix_meta = {} + +--///////////////////////////// +--// Get 'new' matrix object // +--///////////////////////////// + +--// matrix:new ( rows [, columns [, value]] ) +-- if rows is a table then sets rows as matrix +-- if rows is a table of structure {1,2,3} then it sets it as a vector matrix +-- if rows and columns are given and are numbers, returns a matrix with size rowsxcolumns +-- if num is given then returns a matrix with given size and all values set to num +-- if rows is given as number and columns is "I", will return an identity matrix of size rowsxrows +function matrix:new( rows, columns, value ) + -- check for given matrix + if type( rows ) == "table" then + -- check for vector + if type(rows[1]) ~= "table" then -- expect a vector + return setmetatable( {{rows[1]},{rows[2]},{rows[3]}},matrix_meta ) + end + return setmetatable( rows,matrix_meta ) + end + -- get matrix table + local mtx = {} + local value = value or 0 + -- build identity matrix of given rows + if columns == "I" then + for i = 1,rows do + mtx[i] = {} + for j = 1,rows do + if i == j then + mtx[i][j] = 1 + else + mtx[i][j] = 0 + end + end + end + -- build new matrix + else + for i = 1,rows do + mtx[i] = {} + for j = 1,columns do + mtx[i][j] = value + end + end + end + -- return matrix with shared metatable + return setmetatable( mtx,matrix_meta ) +end + +--// matrix ( rows [, comlumns [, value]] ) +-- set __call behaviour of matrix +-- for matrix( ... ) as matrix.new( ... ) +setmetatable( matrix, { __call = function( ... ) return matrix.new( ... ) end } ) + + +-- functions are designed to be light on checks +-- so we get Lua errors instead on wrong input +-- matrix. should handle any table of structure t[i][j] = value +-- we always return a matrix with scripts metatable +-- cause its faster than setmetatable( mtx, getmetatable( input matrix ) ) + +--/////////////////////////////// +--// matrix 'matrix' functions // +--/////////////////////////////// + +--// for real, complex and symbolic matrices //-- + +-- note: real and complex matrices may be added, subtracted, etc. +-- real and symbolic matrices may also be added, subtracted, etc. +-- but one should avoid using symbolic matrices with complex ones +-- since it is not clear which metatable then is used + +--// matrix.add ( m1, m2 ) +-- Add two matrices; m2 may be of bigger size than m1 +function matrix.add( m1, m2 ) + local mtx = {} + for i = 1,#m1 do + local m3i = {} + mtx[i] = m3i + for j = 1,#m1[1] do + m3i[j] = m1[i][j] + m2[i][j] + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.sub ( m1 ,m2 ) +-- Subtract two matrices; m2 may be of bigger size than m1 +function matrix.sub( m1, m2 ) + local mtx = {} + for i = 1,#m1 do + local m3i = {} + mtx[i] = m3i + for j = 1,#m1[1] do + m3i[j] = m1[i][j] - m2[i][j] + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.mul ( m1, m2 ) +-- Multiply two matrices; m1 columns must be equal to m2 rows +-- e.g. #m1[1] == #m2 +function matrix.mul( m1, m2 ) + -- multiply rows with columns + local mtx = {} + for i = 1,#m1 do + mtx[i] = {} + for j = 1,#m2[1] do + local num = m1[i][1] * m2[1][j] + for n = 2,#m1[1] do + num = num + m1[i][n] * m2[n][j] + end + mtx[i][j] = num + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.div ( m1, m2 ) +-- Divide two matrices; m1 columns must be equal to m2 rows +-- m2 must be square, to be inverted, +-- if that fails returns the rank of m2 as second argument +-- e.g. #m1[1] == #m2; #m2 == #m2[1] +function matrix.div( m1, m2 ) + local rank; m2,rank = matrix.invert( m2 ) + if not m2 then return m2, rank end -- singular + return matrix.mul( m1, m2 ) +end + +--// matrix.mulnum ( m1, num ) +-- Multiply matrix with a number +-- num may be of type 'number' or 'complex number' +-- strings get converted to complex number, if that fails then to symbol +function matrix.mulnum( m1, num ) + local mtx = {} + -- multiply elements with number + for i = 1,#m1 do + mtx[i] = {} + for j = 1,#m1[1] do + mtx[i][j] = m1[i][j] * num + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.divnum ( m1, num ) +-- Divide matrix by a number +-- num may be of type 'number' or 'complex number' +-- strings get converted to complex number, if that fails then to symbol +function matrix.divnum( m1, num ) + local mtx = {} + -- divide elements by number + for i = 1,#m1 do + local mtxi = {} + mtx[i] = mtxi + for j = 1,#m1[1] do + mtxi[j] = m1[i][j] / num + end + end + return setmetatable( mtx, matrix_meta ) +end + + +--// for real and complex matrices only //-- + +--// matrix.pow ( m1, num ) +-- Power of matrix; mtx^(num) +-- num is an integer and may be negative +-- m1 has to be square +-- if num is negative and inverting m1 fails +-- returns the rank of matrix m1 as second argument +function matrix.pow( m1, num ) + assert(num == math.floor(num), "exponent not an integer") + if num == 0 then + return matrix:new( #m1,"I" ) + end + if num < 0 then + local rank; m1,rank = matrix.invert( m1 ) + if not m1 then return m1, rank end -- singular + num = -num + end + local mtx = matrix.copy( m1 ) + for i = 2,num do + mtx = matrix.mul( mtx,m1 ) + end + return mtx +end + +local function number_norm2(x) + return x * x +end + +--// matrix.det ( m1 ) +-- Calculate the determinant of a matrix +-- m1 needs to be square +-- Can calc the det for symbolic matrices up to 3x3 too +-- The function to calculate matrices bigger 3x3 +-- is quite fast and for matrices of medium size ~(100x100) +-- and average values quite accurate +-- here we try to get the nearest element to |1|, (smallest pivot element) +-- os that usually we have |mtx[i][j]/subdet| > 1 or mtx[i][j]; +-- with complex matrices we use the complex.abs function to check if it is bigger or smaller +function matrix.det( m1 ) + + -- check if matrix is quadratic + assert(#m1 == #m1[1], "matrix not square") + + local size = #m1 + + if size == 1 then + return m1[1][1] + end + + if size == 2 then + return m1[1][1]*m1[2][2] - m1[2][1]*m1[1][2] + end + + if size == 3 then + return ( m1[1][1]*m1[2][2]*m1[3][3] + m1[1][2]*m1[2][3]*m1[3][1] + m1[1][3]*m1[2][1]*m1[3][2] + - m1[1][3]*m1[2][2]*m1[3][1] - m1[1][1]*m1[2][3]*m1[3][2] - m1[1][2]*m1[2][1]*m1[3][3] ) + end + + --// no symbolic matrix supported below here + local e = m1[1][1] + local zero = type(e) == "table" and e.zero or 0 + local norm2 = type(e) == "table" and e.norm2 or number_norm2 + + --// matrix is bigger than 3x3 + -- get determinant + -- using Gauss elimination and Laplace + -- start eliminating from below better for removals + -- get copy of matrix, set initial determinant + local mtx = matrix.copy( m1 ) + local det = 1 + -- get det up to the last element + for j = 1,#mtx[1] do + -- get smallest element so that |factor| > 1 + -- and set it as last element + local rows = #mtx + local subdet,xrow + for i = 1,rows do + -- get element + local e = mtx[i][j] + -- if no subdet has been found + if not subdet then + -- check if element it is not zero + if e ~= zero then + -- use element as new subdet + subdet,xrow = e,i + end + -- check for elements nearest to 1 or -1 + elseif e ~= zero and math.abs(norm2(e)-1) < math.abs(norm2(subdet)-1) then + subdet,xrow = e,i + end + end + -- only cary on if subdet is found + if subdet then + -- check if xrow is the last row, + -- else switch lines and multiply det by -1 + if xrow ~= rows then + mtx[rows],mtx[xrow] = mtx[xrow],mtx[rows] + det = -det + end + -- traverse all fields setting element to zero + -- we don't set to zero cause we don't use that column anymore then anyways + for i = 1,rows-1 do + -- factor is the dividor of the first element + -- if element is not already zero + if mtx[i][j] ~= zero then + local factor = mtx[i][j]/subdet + -- update all remaining fields of the matrix, with value from xrow + for n = j+1,#mtx[1] do + mtx[i][n] = mtx[i][n] - factor * mtx[rows][n] + end + end + end + -- update determinant and remove row + if math.fmod( rows,2 ) == 0 then + det = -det + end + det = det * subdet + table.remove( mtx ) + else + -- break here table det is 0 + return det * 0 + end + end + -- det ready to return + return det +end + +--// matrix.dogauss ( mtx ) +-- Gauss elimination, Gauss-Jordan Method +-- this function changes the matrix itself +-- returns on success: true, +-- returns on failure: false,'rank of matrix' + +-- locals +-- checking here for the element nearest but not equal to zero (smallest pivot element). +-- This way the `factor` in `dogauss` will be >= 1, which +-- can give better results. +local pivotOk = function( mtx,i,j,norm2 ) + -- find min value + local iMin + local normMin = math.huge + for _i = i,#mtx do + local e = mtx[_i][j] + local norm = math.abs(norm2(e)) + if norm > 0 and norm < normMin then + iMin = _i + normMin = norm + end + end + if iMin then + -- switch lines if not in position. + if iMin ~= i then + mtx[i],mtx[iMin] = mtx[iMin],mtx[i] + end + return true + end + return false +end + +local function copy(x) + return type(x) == "table" and x.copy(x) or x +end + +-- note: in --// ... //-- we have a way that does no divison, +-- however with big number and matrices we get problems since we do no reducing +function matrix.dogauss( mtx ) + local e = mtx[1][1] + local zero = type(e) == "table" and e.zero or 0 + local one = type(e) == "table" and e.one or 1 + local norm2 = type(e) == "table" and e.norm2 or number_norm2 + + local rows,columns = #mtx,#mtx[1] + -- stairs left -> right + for j = 1,rows do + -- check if element can be setted to one + if pivotOk( mtx,j,j,norm2 ) then + -- start parsing rows + for i = j+1,rows do + -- check if element is not already zero + if mtx[i][j] ~= zero then + -- we may add x*otherline row, to set element to zero + -- tozero - x*mtx[j][j] = 0; x = tozero/mtx[j][j] + local factor = mtx[i][j]/mtx[j][j] + --// this should not be used although it does no division, + -- yet with big matrices (since we do no reducing and other things) + -- we get too big numbers + --local factor1,factor2 = mtx[i][j],mtx[j][j] //-- + mtx[i][j] = copy(zero) + for _j = j+1,columns do + --// mtx[i][_j] = mtx[i][_j] * factor2 - factor1 * mtx[j][_j] //-- + mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j] + end + end + end + else + -- return false and the rank of the matrix + return false,j-1 + end + end + -- stairs right <- left + for j = rows,1,-1 do + -- set element to one + -- do division here + local div = mtx[j][j] + for _j = j+1,columns do + mtx[j][_j] = mtx[j][_j] / div + end + -- start parsing rows + for i = j-1,1,-1 do + -- check if element is not already zero + if mtx[i][j] ~= zero then + local factor = mtx[i][j] + for _j = j+1,columns do + mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j] + end + mtx[i][j] = copy(zero) + end + end + mtx[j][j] = copy(one) + end + return true +end + +--// matrix.invert ( m1 ) +-- Get the inverted matrix or m1 +-- matrix must be square and not singular +-- on success: returns inverted matrix +-- on failure: returns nil,'rank of matrix' +function matrix.invert( m1 ) + assert(#m1 == #m1[1], "matrix not square") + local mtx = matrix.copy( m1 ) + local ident = setmetatable( {},matrix_meta ) + local e = m1[1][1] + local zero = type(e) == "table" and e.zero or 0 + local one = type(e) == "table" and e.one or 1 + for i = 1,#m1 do + local identi = {} + ident[i] = identi + for j = 1,#m1 do + identi[j] = copy((i == j) and one or zero) + end + end + mtx = matrix.concath( mtx,ident ) + local done,rank = matrix.dogauss( mtx ) + if done then + return matrix.subm( mtx, 1,(#mtx[1]/2)+1,#mtx,#mtx[1] ) + else + return nil,rank + end +end + +--// matrix.sqrt ( m1 [,iters] ) +-- calculate the square root of a matrix using "Denman Beavers square root iteration" +-- condition: matrix rows == matrix columns; must have a invers matrix and a square root +-- if called without additional arguments, the function finds the first nearest square root to +-- input matrix, there are others but the error between them is very small +-- if called with agument iters, the function will return the matrix by number of iterations +-- the script returns: +-- as first argument, matrix^.5 +-- as second argument, matrix^-.5 +-- as third argument, the average error between (matrix^.5)^2-inputmatrix +-- you have to determin for yourself if the result is sufficent enough for you +-- local average error +local function get_abs_avg( m1, m2 ) + local dist = 0 + local e = m1[1][1] + local abs = type(e) == "table" and e.abs or math.abs + for i=1,#m1 do + for j=1,#m1[1] do + dist = dist + abs(m1[i][j]-m2[i][j]) + end + end + -- norm by numbers of entries + return dist/(#m1*2) +end +-- square root function +function matrix.sqrt( m1, iters ) + assert(#m1 == #m1[1], "matrix not square") + local iters = iters or math.huge + local y = matrix.copy( m1 ) + local z = matrix(#y, 'I') + local dist = math.huge + -- iterate, and get the average error + for n=1,iters do + local lasty,lastz = y,z + -- calc square root + -- y, z = (1/2)*(y + z^-1), (1/2)*(z + y^-1) + y, z = matrix.divnum((matrix.add(y,matrix.invert(z))),2), + matrix.divnum((matrix.add(z,matrix.invert(y))),2) + local dist1 = get_abs_avg(y,lasty) + if iters == math.huge then + if dist1 >= dist then + return lasty,lastz,get_abs_avg(matrix.mul(lasty,lasty),m1) + end + end + dist = dist1 + end + return y,z,get_abs_avg(matrix.mul(y,y),m1) +end + +--// matrix.root ( m1, root [,iters] ) +-- calculate any root of a matrix +-- source: http://www.dm.unipi.it/~cortona04/slides/bruno.pdf +-- m1 and root have to be given;(m1 = matrix, root = number) +-- conditions same as matrix.sqrt +-- returns same values as matrix.sqrt +function matrix.root( m1, root, iters ) + assert(#m1 == #m1[1], "matrix not square") + local iters = iters or math.huge + local mx = matrix.copy( m1 ) + local my = matrix.mul(mx:invert(),mx:pow(root-1)) + local dist = math.huge + -- iterate, and get the average error + for n=1,iters do + local lastx,lasty = mx,my + -- calc root of matrix + --mx,my = ((p-1)*mx + my^-1)/p, + -- ((((p-1)*my + mx^-1)/p)*my^-1)^(p-2) * + -- ((p-1)*my + mx^-1)/p + mx,my = mx:mulnum(root-1):add(my:invert()):divnum(root), + my:mulnum(root-1):add(mx:invert()):divnum(root) + :mul(my:invert():pow(root-2)):mul(my:mulnum(root-1) + :add(mx:invert())):divnum(root) + local dist1 = get_abs_avg(mx,lastx) + if iters == math.huge then + if dist1 >= dist then + return lastx,lasty,get_abs_avg(matrix.pow(lastx,root),m1) + end + end + dist = dist1 + end + return mx,my,get_abs_avg(matrix.pow(mx,root),m1) +end + + +--// Norm functions //-- + +--// matrix.normf ( mtx ) +-- calculates the Frobenius norm of the matrix. +-- ||mtx||_F = sqrt(SUM_{i,j} |a_{i,j}|^2) +-- http://en.wikipedia.org/wiki/Frobenius_norm#Frobenius_norm +function matrix.normf(mtx) + local mtype = matrix.type(mtx) + local result = 0 + for i = 1,#mtx do + for j = 1,#mtx[1] do + local e = mtx[i][j] + if mtype ~= "number" then e = e:abs() end + result = result + e^2 + end + end + local sqrt = (type(result) == "number") and math.sqrt or result.sqrt + return sqrt(result) +end + +--// matrix.normmax ( mtx ) +-- calculates the max norm of the matrix. +-- ||mtx||_{max} = max{|a_{i,j}|} +-- Does not work with symbolic matrices +-- http://en.wikipedia.org/wiki/Frobenius_norm#Max_norm +function matrix.normmax(mtx) + local abs = (matrix.type(mtx) == "number") and math.abs or mtx[1][1].abs + local result = 0 + for i = 1,#mtx do + for j = 1,#mtx[1] do + local e = abs(mtx[i][j]) + if e > result then result = e end + end + end + return result +end + + +--// only for number and complex type //-- +-- Functions changing the matrix itself + +--// matrix.round ( mtx [, idp] ) +-- perform round on elements +local numround = function( num,mult ) + return math.floor( num * mult + 0.5 ) / mult +end +local tround = function( t,mult ) + for i,v in ipairs(t) do + t[i] = math.floor( v * mult + 0.5 ) / mult + end + return t +end +function matrix.round( mtx, idp ) + local mult = 10^( idp or 0 ) + local fround = matrix.type( mtx ) == "number" and numround or tround + for i = 1,#mtx do + for j = 1,#mtx[1] do + mtx[i][j] = fround(mtx[i][j],mult) + end + end + return mtx +end + +--// matrix.random( mtx [,start] [, stop] [, idip] ) +-- fillmatrix with random values +local numfill = function( _,start,stop,idp ) + return math.random( start,stop ) / idp +end +local tfill = function( t,start,stop,idp ) + for i in ipairs(t) do + t[i] = math.random( start,stop ) / idp + end + return t +end +function matrix.random( mtx,start,stop,idp ) + local start,stop,idp = start or -10,stop or 10,idp or 1 + local ffill = matrix.type( mtx ) == "number" and numfill or tfill + for i = 1,#mtx do + for j = 1,#mtx[1] do + mtx[i][j] = ffill( mtx[i][j], start, stop, idp ) + end + end + return mtx +end + + +--////////////////////////////// +--// Object Utility Functions // +--////////////////////////////// + +--// for all types and matrices //-- + +--// matrix.type ( mtx ) +-- get type of matrix, normal/complex/symbol or tensor +function matrix.type( mtx ) + local e = mtx[1][1] + if type(e) == "table" then + if e.type then + return e:type() + end + return "tensor" + end + return "number" +end + +-- local functions to copy matrix values +local num_copy = function( num ) + return num +end +local t_copy = function( t ) + local newt = setmetatable( {}, getmetatable( t ) ) + for i,v in ipairs( t ) do + newt[i] = v + end + return newt +end + +--// matrix.copy ( m1 ) +-- Copy a matrix +-- simple copy, one can write other functions oneself +function matrix.copy( m1 ) + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + local mtx = {} + for i = 1,#m1[1] do + mtx[i] = {} + for j = 1,#m1 do + mtx[i][j] = docopy( m1[i][j] ) + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.transpose ( m1 ) +-- Transpose a matrix +-- switch rows and columns +function matrix.transpose( m1 ) + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + local mtx = {} + for i = 1,#m1[1] do + mtx[i] = {} + for j = 1,#m1 do + mtx[i][j] = docopy( m1[j][i] ) + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.subm ( m1, i1, j1, i2, j2 ) +-- Submatrix out of a matrix +-- input: i1,j1,i2,j2 +-- i1,j1 are the start element +-- i2,j2 are the end element +-- condition: i1,j1,i2,j2 are elements of the matrix +function matrix.subm( m1,i1,j1,i2,j2 ) + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + local mtx = {} + for i = i1,i2 do + local _i = i-i1+1 + mtx[_i] = {} + for j = j1,j2 do + local _j = j-j1+1 + mtx[_i][_j] = docopy( m1[i][j] ) + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.concath( m1, m2 ) +-- Concatenate two matrices, horizontal +-- will return m1m2; rows have to be the same +-- e.g.: #m1 == #m2 +function matrix.concath( m1,m2 ) + assert(#m1 == #m2, "matrix size mismatch") + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + local mtx = {} + local offset = #m1[1] + for i = 1,#m1 do + mtx[i] = {} + for j = 1,offset do + mtx[i][j] = docopy( m1[i][j] ) + end + for j = 1,#m2[1] do + mtx[i][j+offset] = docopy( m2[i][j] ) + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.concatv ( m1, m2 ) +-- Concatenate two matrices, vertical +-- will return m1 +-- m2 +-- columns have to be the same; e.g.: #m1[1] == #m2[1] +function matrix.concatv( m1,m2 ) + assert(#m1[1] == #m2[1], "matrix size mismatch") + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + local mtx = {} + for i = 1,#m1 do + mtx[i] = {} + for j = 1,#m1[1] do + mtx[i][j] = docopy( m1[i][j] ) + end + end + local offset = #mtx + for i = 1,#m2 do + local _i = i + offset + mtx[_i] = {} + for j = 1,#m2[1] do + mtx[_i][j] = docopy( m2[i][j] ) + end + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.rotl ( m1 ) +-- Rotate Left, 90 degrees +function matrix.rotl( m1 ) + local mtx = matrix:new( #m1[1],#m1 ) + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + for i = 1,#m1 do + for j = 1,#m1[1] do + mtx[#m1[1]-j+1][i] = docopy( m1[i][j] ) + end + end + return mtx +end + +--// matrix.rotr ( m1 ) +-- Rotate Right, 90 degrees +function matrix.rotr( m1 ) + local mtx = matrix:new( #m1[1],#m1 ) + local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy + for i = 1,#m1 do + for j = 1,#m1[1] do + mtx[j][#m1-i+1] = docopy( m1[i][j] ) + end + end + return mtx +end + +local function tensor_tostring( t,fstr ) + if not fstr then return "["..table.concat(t,",").."]" end + local tval = {} + for i,v in ipairs( t ) do + tval[i] = string.format( fstr,v ) + end + return "["..table.concat(tval,",").."]" +end +local function number_tostring( e,fstr ) + return fstr and string.format( fstr,e ) or e +end + +--// matrix.tostring ( mtx, formatstr ) +-- tostring function +function matrix.tostring( mtx, formatstr ) + local ts = {} + local mtype = matrix.type( mtx ) + local e = mtx[1][1] + local tostring = mtype == "tensor" and tensor_tostring or + type(e) == "table" and e.tostring or number_tostring + for i = 1,#mtx do + local tstr = {} + for j = 1,#mtx[1] do + tstr[j] = tostring(mtx[i][j],formatstr) + end + ts[i] = table.concat(tstr, "\t") + end + return table.concat(ts, "\n") +end + +--// matrix.print ( mtx [, formatstr] ) +-- print out the matrix, just calls tostring +function matrix.print( ... ) + print( matrix.tostring( ... ) ) +end + +--// matrix.latex ( mtx [, align] ) +-- LaTeX output +function matrix.latex( mtx, align ) + -- align : option to align the elements + -- c = center; l = left; r = right + -- \usepackage{dcolumn}; D{.}{,}{-1}; aligns number by . replaces it with , + local align = align or "c" + local str = "$\\left( \\begin{array}{"..string.rep( align, #mtx[1] ).."}\n" + local getstr = matrix.type( mtx ) == "tensor" and tensor_tostring or number_tostring + for i = 1,#mtx do + str = str.."\t"..getstr(mtx[i][1]) + for j = 2,#mtx[1] do + str = str.." & "..getstr(mtx[i][j]) + end + -- close line + if i == #mtx then + str = str.."\n" + else + str = str.." \\\\\n" + end + end + return str.."\\end{array} \\right)$" +end + + +--// Functions not changing the matrix + +--// matrix.rows ( mtx ) +-- return number of rows +function matrix.rows( mtx ) + return #mtx +end + +--// matrix.columns ( mtx ) +-- return number of columns +function matrix.columns( mtx ) + return #mtx[1] +end + +--// matrix.size ( mtx ) +-- get matrix size as string rows,columns +function matrix.size( mtx ) + if matrix.type( mtx ) == "tensor" then + return #mtx,#mtx[1],#mtx[1][1] + end + return #mtx,#mtx[1] +end + +--// matrix.getelement ( mtx, i, j ) +-- return specific element ( row,column ) +-- returns element on success and nil on failure +function matrix.getelement( mtx,i,j ) + if mtx[i] and mtx[i][j] then + return mtx[i][j] + end +end + +--// matrix.setelement( mtx, i, j, value ) +-- set an element ( i, j, value ) +-- returns 1 on success and nil on failure +function matrix.setelement( mtx,i,j,value ) + if matrix.getelement( mtx,i,j ) then + -- check if value type is number + mtx[i][j] = value + return 1 + end +end + +--// matrix.ipairs ( mtx ) +-- iteration, same for complex +function matrix.ipairs( mtx ) + local i,j,rows,columns = 1,0,#mtx,#mtx[1] + local function iter() + j = j + 1 + if j > columns then -- return first element from next row + i,j = i + 1,1 + end + if i <= rows then + return i,j + end + end + return iter +end + +--/////////////////////////////// +--// matrix 'vector' functions // +--/////////////////////////////// + +-- a vector is defined as a 3x1 matrix +-- get a vector; vec = matrix{{ 1,2,3 }}^'T' + +--// matrix.scalar ( m1, m2 ) +-- returns the Scalar Product of two 3x1 matrices (vectors) +function matrix.scalar( m1, m2 ) + return m1[1][1]*m2[1][1] + m1[2][1]*m2[2][1] + m1[3][1]*m2[3][1] +end + +--// matrix.cross ( m1, m2 ) +-- returns the Cross Product of two 3x1 matrices (vectors) +function matrix.cross( m1, m2 ) + local mtx = {} + mtx[1] = { m1[2][1]*m2[3][1] - m1[3][1]*m2[2][1] } + mtx[2] = { m1[3][1]*m2[1][1] - m1[1][1]*m2[3][1] } + mtx[3] = { m1[1][1]*m2[2][1] - m1[2][1]*m2[1][1] } + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.len ( m1 ) +-- returns the Length of a 3x1 matrix (vector) +function matrix.len( m1 ) + return math.sqrt( m1[1][1]^2 + m1[2][1]^2 + m1[3][1]^2 ) +end + + +--// matrix.replace (mtx, func, ...) +-- for each element e in the matrix mtx, replace it with func(mtx, ...). +function matrix.replace( m1, func, ... ) + local mtx = {} + for i = 1,#m1 do + local m1i = m1[i] + local mtxi = {} + for j = 1,#m1i do + mtxi[j] = func( m1i[j], ... ) + end + mtx[i] = mtxi + end + return setmetatable( mtx, matrix_meta ) +end + +--// matrix.remcomplex ( mtx ) +-- set the matrix elements to strings +-- IMPROVE: tostring v.s. tostringelements confusing +function matrix.elementstostrings( mtx ) + local e = mtx[1][1] + local tostring = type(e) == "table" and e.tostring or tostring + return matrix.replace(mtx, tostring) +end + +--// matrix.solve ( m1 ) +-- solve; tries to solve a symbolic matrix to a number +function matrix.solve( m1 ) + assert( matrix.type( m1 ) == "symbol", "matrix not of type 'symbol'" ) + local mtx = {} + for i = 1,#m1 do + mtx[i] = {} + for j = 1,#m1[1] do + mtx[i][j] = tonumber( loadstring( "return "..m1[i][j][1] )() ) + end + end + return setmetatable( mtx, matrix_meta ) +end + +--////////////////////////-- +--// METATABLE HANDLING //-- +--////////////////////////-- + +--// MetaTable +-- as we declaired on top of the page +-- local/shared metatable +-- matrix_meta + +-- note '...' is always faster than 'arg1,arg2,...' if it can be used + +-- Set add "+" behaviour +matrix_meta.__add = function( ... ) + return matrix.add( ... ) +end + +-- Set subtract "-" behaviour +matrix_meta.__sub = function( ... ) + return matrix.sub( ... ) +end + +-- Set multiply "*" behaviour +matrix_meta.__mul = function( m1,m2 ) + if getmetatable( m1 ) ~= matrix_meta then + return matrix.mulnum( m2,m1 ) + elseif getmetatable( m2 ) ~= matrix_meta then + return matrix.mulnum( m1,m2 ) + end + return matrix.mul( m1,m2 ) +end + +-- Set division "/" behaviour +matrix_meta.__div = function( m1,m2 ) + if getmetatable( m1 ) ~= matrix_meta then + return matrix.mulnum( matrix.invert(m2),m1 ) + elseif getmetatable( m2 ) ~= matrix_meta then + return matrix.divnum( m1,m2 ) + end + return matrix.div( m1,m2 ) +end + +-- Set unary minus "-" behavior +matrix_meta.__unm = function( mtx ) + return matrix.mulnum( mtx,-1 ) +end + +-- Set power "^" behaviour +-- if opt is any integer number will do mtx^opt +-- (returning nil if answer doesn't exist) +-- if opt is 'T' then it will return the transpose matrix +-- only for complex: +-- if opt is '*' then it returns the complex conjugate matrix + local option = { + -- only for complex + ["*"] = function( m1 ) return matrix.conjugate( m1 ) end, + -- for both + ["T"] = function( m1 ) return matrix.transpose( m1 ) end, + } +matrix_meta.__pow = function( m1, opt ) + return option[opt] and option[opt]( m1 ) or matrix.pow( m1,opt ) +end + +-- Set equal "==" behaviour +matrix_meta.__eq = function( m1, m2 ) + -- check same type + if matrix.type( m1 ) ~= matrix.type( m2 ) then + return false + end + -- check same size + if #m1 ~= #m2 or #m1[1] ~= #m2[1] then + return false + end + -- check elements equal + for i = 1,#m1 do + for j = 1,#m1[1] do + if m1[i][j] ~= m2[i][j] then + return false + end + end + end + return true +end + +-- Set tostring "tostring( mtx )" behaviour +matrix_meta.__tostring = function( ... ) + return matrix.tostring( ... ) +end + +-- set __call "mtx( [formatstr] )" behaviour, mtx [, formatstr] +matrix_meta.__call = function( ... ) + matrix.print( ... ) +end + +--// __index handling +matrix_meta.__index = {} +for k,v in pairs( matrix ) do + matrix_meta.__index[k] = v +end + + +--///////////////////////////////// +--// symbol class implementation +--///////////////////////////////// + +-- access to the symbolic metatable +local symbol_meta = {}; symbol_meta.__index = symbol_meta +local symbol = symbol_meta + +function symbol_meta.new(o) + return setmetatable({tostring(o)}, symbol_meta) +end +symbol_meta.to = symbol_meta.new + +-- symbol( arg ) +-- same as symbol.to( arg ) +-- set __call behaviour of symbol +setmetatable( symbol_meta, { __call = function( _,s ) return symbol_meta.to( s ) end } ) + + +-- Converts object to string, optionally with formatting. +function symbol_meta.tostring( e,fstr ) + return string.format( fstr,e[1] ) +end + +-- Returns "symbol" if object is a symbol type, else nothing. +function symbol_meta:type() + if getmetatable(self) == symbol_meta then + return "symbol" + end +end + +-- Performs string.gsub on symbol. +-- for use in matrix.replace +function symbol_meta:gsub(from, to) + return symbol.to( string.gsub( self[1],from,to ) ) +end + +-- creates function that replaces one letter by something else +-- makereplacer( "a",4,"b",7, ... )(x) +-- will replace a with 4 and b with 7 in symbol x. +-- for use in matrix.replace +function symbol_meta.makereplacer( ... ) + local tosub = {} + local args = {...} + for i = 1,#args,2 do + tosub[args[i]] = args[i+1] + end + local function func( a ) return tosub[a] or a end + return function(sym) + return symbol.to( string.gsub( sym[1], "%a", func ) ) + end +end + +-- applies abs function to symbol +function symbol_meta.abs(a) + return symbol.to("(" .. a[1] .. "):abs()") +end + +-- applies sqrt function to symbol +function symbol_meta.sqrt(a) + return symbol.to("(" .. a[1] .. "):sqrt()") +end + +function symbol_meta.__add(a,b) + return symbol.to(a .. "+" .. b) +end + +function symbol_meta.__sub(a,b) + return symbol.to(a .. "-" .. b) +end + +function symbol_meta.__mul(a,b) + return symbol.to("(" .. a .. ")*(" .. b .. ")") +end + +function symbol_meta.__div(a,b) + return symbol.to("(" .. a .. ")/(" .. b .. ")") +end + +function symbol_meta.__pow(a,b) + return symbol.to("(" .. a .. ")^(" .. b .. ")") +end + +function symbol_meta.__eq(a,b) + return a[1] == b[1] +end + +function symbol_meta.__tostring(a) + return a[1] +end + +function symbol_meta.__concat(a,b) + return tostring(a) .. tostring(b) +end + +matrix.symbol = symbol + + +-- return matrix +return matrix + +--///////////////-- +--// chillcode //-- +--///////////////-- diff --git a/mod.conf b/mod.conf new file mode 100644 index 0000000..5a69441 --- /dev/null +++ b/mod.conf @@ -0,0 +1 @@ +name = matrix