# [AHK_L] Matrix class

16 replies to this topic

### #1 IsNull

IsNull
• Fellows
• 990 posts

Posted 13 May 2012 - 08:10 PM

Matrices are a very powerful mathematical thing, especially useful in programming for:

[*:fhte51gc] Transformations (2D/3D)
[*:fhte51gc] Color Transformations
[*:fhte51gc] Solving Linear Equations (Determinant or Gauss)
[*:fhte51gc] etc...
The Matrix Class offers several such Matrix related operations:

- Addition/Subtraction of two Matrices
- Multiplication of two Matrices
- Calculate the Determinant of any Matrix (laplace expansion) 8)
- Scalar-Multiplication of a Matrix
- Transposing a Matrix
- Generating Identity Matrices
- Generating filled Matrices of [n]
- Inverting a Matrix
- Transform a Matrix to the row echelon form (Gauss)
- Solve Linear Equation system (Gauss, using a pivot strategy)

All Methods can be used either as instance Methods on a Matrix-Object, or as static Helper Methods, eg:
or
A := new Matrix([1,2])
B := new Matrix([4,5])

Matrix Class

/**
* AHK static Matrix calculations
* by IsNull 2012
*
* Gauss/Pivot Strategies implemented by horst/Babba
*/
class MatrixStatic
{
/**
* Calculates the determinant (scalar) of the given Matrix
*
* Implementation Note:
* The given matrix will be reduced by laplace expansion
* until the matrix dimension is (2,2). The determinant of
* the 2,2 matrix is then directly calculated
*/
Det(m){

det := 0

colCnt := MatrixStatic.ColumnCount(m)
if(!MatrixStatic.IsSquare(m) ||  colCnt < 2)
throw new Exception("The matrix must be squared and its dimensions must be greater or equal than (2,2)!")

if(colCnt == 2){
; a 2,2 Matrix, we can calculate the determinant directly
det := m[1,1] * m[2,2] - m[2,1] * m[1,2]
}else{
; Laplace expansion
; @see: http://en.wikipedia.org/wiki/Laplace_expansion

i := 1 ;row
k := 1 ;col
det := 0
curVal := 0 ; current Value where [i,k] is pointing at
Loop, % colCnt
{
k := A_index

curVal := m[i,k]

if(curVal != 0) ; we can skip coz multiplication by zero
{
laplace := MatrixStatic.ExtractLaplace(m, [i,k]) ; extract a sub matrix by laplace
cofactor := curVal * (((-1)**(i+k)) * MatrixStatic.Det(laplace))
det += cofactor
}
}
}

return det
}

/**
* Extract the sub-matrix, by laplace expansion
*
* m    -   Matrix   (n,n)       source matrix
* pnt   -   Point   [row,col]   coord origin of Laplace
*
* returns Laplace-Matrix of [pnt]
*
* @see: http://en.wikipedia.org/wiki/Laplace_expansion
*/
ExtractLaplace(m, pnt){
laplace := {}

colCnt := MatrixStatic.ColumnCount(m)
rowCnt := MatrixStatic.RowCount(m)

pntRow := pnt[1]
pntCol := pnt[2]

laplaceRow := 1
laplaceCol := 1
bincRow := false

irow := 1
Loop, % colCnt
{
irow := a_index
bincRow := false
laplaceCol := 1
Loop, % colCnt
{
icol := a_index

if(pntCol != icol && pntRow != irow) ; // omit values in the range of the laplace origin
{
laplace[laplaceRow, laplaceCol] := m[irow, icol]
laplaceCol++
bincRow := true
}
}
if(bincRow)
laplaceRow++
}

return laplace
}

/**
* Multiply each element in the matrix whith the given scalar
*/
MultiplyScalar(A, n){

isColVector := (MatrixStatic.RowCount(A) == 1)

loop % MatrixStatic.RowCount(A)
{
rowi := a_index
Loop, % MatrixStatic.ColumnCount(A)
{
if(isColVector)
add[a_index] := A[A_Index] * n
else
add[rowi,a_index] := A[rowi,A_Index] * n
}
}
}

/**
* Addition of A and B
*/
if(MatrixStatic.ColumnCount(A) != MatrixStatic.ColumnCount(B)
|| MatrixStatic.RowCount(A) != MatrixStatic.RowCount(B))
{
throw Exception("Matrix Addition Error: All dimensions must agree!")
}

isColVector := (MatrixStatic.RowCount(A) == 1)

loop % MatrixStatic.RowCount(A)
{
rowi := a_index
Loop, % MatrixStatic.ColumnCount(A)
{
if(isColVector)
add[a_index] := A[A_Index] + B[A_Index]
else
add[rowi,a_index] := A[rowi,A_Index] + B[rowi,A_Index]
}
}
}

/**
* Subtraction of B from A
*/
Sub(A, B){
if(MatrixStatic.ColumnCount(A) != MatrixStatic.ColumnCount(B)
|| MatrixStatic.RowCount(A) != MatrixStatic.RowCount(B))
{
throw Exception("Matrix Addition Error: All dimensions must agree!")
}

isColVector := (MatrixStatic.RowCount(A) == 1)

loop % MatrixStatic.RowCount(A)
{
rowi := a_index
Loop, % MatrixStatic.ColumnCount(A)
{
if(isColVector)
add[a_index] := A[A_Index] - B[A_Index]
else
add[rowi,a_index] := A[rowi,A_Index] - B[rowi,A_Index]
}
}
}

/*
* Inverts the given Matrix
* so that: A * inv(A) = I
*/
Inverse(A){
if(!MatrixStatic.IsSquare(A))
throw Exception("Only square matrices have an inverse! Your given matrix is not square.")

size := MatrixStatic.ColumnCount(A)
identity := MatrixStatic.Eye(size)

inverse := MatrixStatic.Gauss(A, identity)
return inverse
}

/*
* Multiplies Matrix A with B.
*
* Note that A*B != B*A
*/
Multiply(A, B){
mul := {}
if(MatrixStatic.ColumnCount(A) != MatrixStatic.RowCount(B))
{
throw Exception("Matrix Multiplication Error: Inner dimensions must agree!")
}

Loop, % MatrixStatic.RowCount(A)
{
rowi := A_index
Loop, % MatrixStatic.ColumnCount(B)
{
aRow := MatrixStatic.RangeRow(A, rowi)
bCol := MatrixStatic.RangeCol(B, A_index)
bColT := MatrixStatic.Transpose(bCol)

mul[rowi,A_index] := MatrixStatic.dotP(aRow,bColT)
}
}
return mul
}

/*
* Dot-Product (scalar product)
*/
dotP(v1,v2){
dotp := 0
Loop, % v1.MaxIndex()
dotp += v1[A_index] * v2[A_index]
return dotp
}

/**
* Returns the Column-Vector at the given index
*/
RangeCol(m, colIndex){
col := {}

Loop, % MatrixStatic.RowCount(m)
{
rowi := A_index
col[rowi,1] := m[rowi,colIndex]
}

return col
}

/*
*  Returns the Row-Vector at the given index
*/
RangeRow(m, rowIndex){
return MatrixStatic.Clone(m[rowIndex])
}

/*
* Returns a transposed Matrix of the given
*/
Transpose(m){
mt := {} ; transposed matrix
i := 0
for each, row in m
{
i := A_index

if(MatrixStatic.ColumnCount(m) == 1){
mt[i] := row[1]
}else{
if(IsObject(row))
{
for each, item in row
mt[A_index,i] := item
}else{
mt[i,1] := row
}
}
}
return mt
}

IsSquare(m){
return IsObject(m) && MatrixStatic.ColumnCount(m) == MatrixStatic.RowCount(m)
}

/**
* Returns the count of columns in the given Matrix
*/
ColumnCount(m){
return IsObject(m[1]) ? m[1].MaxIndex() : m.MaxIndex()
}

/**
* Returns the count of rows in the given Matrix
*/
RowCount(m){
return IsObject(m[1]) ? m.MaxIndex() : 1
}

/**
* Clones the given Matrix
*/
Clone(m){
mc := {}
for each, row in m
{
if(IsObject(row))
{
rIndex := A_index
for each, item in row
mc[rIndex,A_index] := item
}else
mc[A_index] := row
}
return mc
}

/*
* Generates an quadratic identity matrix with a size of [n]
*/
Eye(n){
eye := {}
loop, % n
{
ri := A_Index
loop, % n
eye[ri,A_index] := (ri == a_index) ? 1 : 0
}

return eye
}

Zeros(n){
return this.Fill(n, 0)
}

/*
* Generates an quadratic matrix with a size of [n]
* and each element set to [fillNum]
*/
Fill(n, fillNum){
filled := {}
loop, % n
{
ri := A_Index
loop, % n
filled[ri,A_index] := fillNum
}

return filled
}

/**
* Checks if the given two matrices are equal
*/
Equals(m,m2){
equal := false

mRowCnt := MatrixStatic.RowCount(m)
mColCnt := MatrixStatic.ColumnCount(m)

if(mRowCnt == MatrixStatic.RowCount(m2)
&& mColCnt == MatrixStatic.ColumnCount(m2)){

equal := true
isColVector := (mRowCnt == 1)

loop % mRowCnt
{
rowi := a_index
Loop, % mColCnt
{
if(isColVector)
{
MsgBox "it is" %rowi% %A_Index%
if(m[A_Index] != m2[A_Index])
{
equal := false
break, 2 ; break the outer loop
}
}else{
if(m[rowi,A_Index] != m2[rowi,A_Index])
{
equal := false
break, 2 ; break the outer loop
}
}
}
}
}
return equal
}

/*
* Returns a console/msgbox friendly string. Useful for debugging
*/
ToString(m){

if(!IsObject(m))
{
prnt := "No Matrix!"
}else{
prnt := "(" MatrixStatic.RowCount(m) "," MatrixStatic.ColumnCount(m) ") Matrix:`n---------`n"

if(MatrixStatic.RowCount(m) != 1)
{
for each, row in m
{
for each, item in row
prnt .= item " "
prnt .= "`n"
}
}else{
for each, val in m
prnt .= val " "
prnt .= "`n"
}
prnt .= "---------"
}

return prnt
}

/**
* Gets the mirror matrix for the straight line mirror y;
* y := m*x + b
*
*
* if b is NOT zero, you must move the vertices first with v = (0,-b)
* 1. move vertices along v = (0,-b)  | --> move the mirror axis to the origin
* 2. mirror with the matrix returned by this function |---> newVertex = M * vertex
* 3. move vertices back along -v = (0,b) | --> move the vertices back
*/
Mirror2D(m){
angle := this.Mirror2DByAngle(ATan(m))
}

/**
* Gets the mirror matrix for the straight line intersection the origin (0,0)
* angle     =    angle to x axis, in radians
*/
Mirror2DByAngle(angle){
angle *= 2
c := cos(angle)
s := sin(angle)

M := [[c, s]
,[s,-c]]

return M
}

Rotate2D(angle){
c := cos(angle)
s := sin(angle)

R := [[c,-s]
,[s, c]]

return R
}

Rotate3DZ(angle){
c := cos(angle)
s := sin(angle)

R := [[c,-s, 0]
,[s, c, 0]
,[0, 0, 1]]

return R
}

Rotate3DY(angle){
c := cos(angle)
s := sin(angle)

R := [[ c, 0, s]
,[ 0, 1, 0]
,[-s, 0, c]]

return R
}

Rotate3DX(angle){
c := cos(angle)
s := sin(angle)

R := [[ 1, 0, 0]
,[ 0, c,-s]
,[ 0, s, c]]

return R
}

Mirror3D(axis){

m1 := (axis == 1 || axis = "x") ? -1 : 1
m2 := (axis == 2 || axis = "y") ? -1 : 1
m3 := (axis == 3 || axis = "z") ? -1 : 1

M := [[ m1, 0,  0]
,[ 0, m2,  0]
,[ 0,  0, m3]]

return M
}

/*
* Given a matrix a and an (optional) index pair i,
* MatrixStatic.Pivot returns an index pair of a pivot located below and to the right of i if it finds some.
*
*/
Pivot(a, i="")
{
colnum:=a[1].maxindex()
rownum:=a.maxindex()

if(i == "")
i1 := i2 := 1
else
{
i1 := i[1] + 1
i2 := i[2] + 1
}

p2 := i2
while (p2 <= colnum)
{
p1 := i1
x := abs(a[i1,p2])

while (i1 + a_index <= rownum){
y := abs(a[i1+a_index, p2])
if(y > x)
{
x := y
p1 := i1+a_index
}
}
if(x)
return [p1,p2]
else
p2++
}
}

;=========================================================================

/***************************************
*
* Given matrices a and (optional) b   (b is allowed to be a vector), mat_rowechtrafo reduces a to row echelon form.
* On b the same (only a-dependent) operations take place.
****************************************
*/
ToRowEchelonForm(aorig, b=""){

a := MatrixStatic.Clone(aorig)

rownum := a.maxindex()
colnum := a[1].maxindex()

b_is_2d := IsObject(b[1])

if (b_is_2d)
colnum_b := b[1].maxindex()

k := 1
while(k < rownum){
pivot := MatrixStatic.Pivot(a, pivot)
p1 := pivot[1]
p2 := pivot[2]

if (!pivot)
break

if (p1 != k){
swap  := a[k]
a[k]  := a[p1]
a[p1] := swap
swap  := b[k]
b[k]  := b[p1]
b[p1] := swap

pivot[1] := p1 := k
}

l1 := p1 + 1
while(l1 <= rownum){
l2 := p2 + 1
fact := a[l1,p2] / a[p1,p2]

while(l2 <= colnum){
a[l1,l2] -= fact * a[p1,l2]
l2 += 1
}
a[l1,p2] := 0

if (b_is_2d)
while (a_index <= colnum_b)
b[l1,a_index] -= fact * b[p1,a_index]
else
b[l1] -= fact * b[p1]
l1++
}
k++
}
return a
}

;=========================================================================

/**
* Given a matrix a in row echelon form and a vector   b   (b is allowed to be a matrix),   mat_RowEchSol   returns a solution if it finds one.
* In case b is matrix,   mat_RowEchSol   interprets the columns of   b   to be vectors and searches solutions to those vectors. If a solution to the ith column is found, it is inserted in the   output array   as ith column. Otherwise the   output array   lacks the ith column.
*/
RowEchelonSolve(a, b, pivot_row2col=""){

rownum := a.maxindex()
colnum := a[1].maxindex()

if (!pivot_row2col){
pivot_row2col := {}
while (pivot := MatrixStatic.Pivot(a,pivot))
pivot_row2col[pivot[1]] := pivot[2]
}

rnk := pivot_row2col.maxindex()
if (rnk="")
rnk := 0

if (!IsObject(b[1])){
j1 := rownum
while(j1 >= 1){
if (b[j1] != 0)
{
if (j1 <= rnk)
break
else
return
}
j1--
}

sol_vec:={}
loop, % colnum
sol_vec[a_index] := 0

while (rnk>=1)
{
sol_vec[pivot_row2col[rnk]] := b[rnk] / a[rnk,pivot_row2col[rnk]]
while (a_index<rnk)
b[a_index] -= a[a_index,pivot_row2col[rnk]] * sol_vec[pivot_row2col[rnk]]
rnk--
}
}else{

valid_ix := {}
colnum_b := b[1].maxindex()
some_index := 1

loop, % colnum_b
{
j2 := a_index
j1 := rownum

while (j1 >= 1){
if (b[j1,j2] != 0){
if (j1 <= rnk){
valid_ix[some_index] := j2
some_index += 1
}
break
}
j1--
}
}

rnk_BkUp := rnk
sol_vec := {}

for count,i2 in valid_ix
{
loop, % colnum
sol_vec[a_index,i2] := 0

rnk := rnk_BkUP
while(rnk >= 1){
sol_vec[pivot_row2col[rnk],i2] := b[rnk,i2] / a[rnk,pivot_row2col[rnk]]
while(a_index < rnk)
b[a_index,i2] -= a[a_index,pivot_row2col[rnk]] * sol_vec[pivot_row2col[rnk],i2]
rnk--
}
}
}
return sol_vec
}

/**
* Gauss solve the System Ax = B.
* The System A is first transformed into row echolon form (gaussian elimination)
* Afterwards, the System is solved by substitiution.
*
* returns x, the solution of Ax = B
*/
Gauss(A, B)
{
workerB := this.Clone(B)
rowechelon := MatrixStatic.ToRowEchelonForm(A, workerB)
return MatrixStatic.RowEchelonSolve(rowechelon, workerB)
}

}

/**
* AHK Matrix class
*
*/
class Matrix extends MatrixStatic
{
__new(m)
{
this.Prototype(m)
}

/**
* Calculates the determinant (scalar) of the given Matrix
*
* Implementation Note:
* The given matrix will be reduced by laplace expansion
* until the matrix dimension is (2,2). The determinant of
* the 2,2 matrix is then directly calculated
*/
Det(m=0){
if (this != Matrix)
m := this
return base.Det(m)
}

/**
* Multiply each element in this matrix whith the given scalar [n] and returns the new matrix
*/
MultiplyScalar(m, n=0){

if (this != Matrix)
{
n := m
m := this
}

return new Matrix(base.MultiplyScalar(m, n))
}

/**
* Add the given matrix to this one and return the new matrix
*/

if (this != Matrix)
{
B := m
m := this
}

return new Matrix(base.Add(m, B))
}

/**
* Subtract B from this matrix and returns the new matrix
*/
Sub(B){

if (this != Matrix)
{
B := m
m := this
}

return new Matrix(base.Sub(m, B))
}

/*
* Returns the Inverse of this Matrix
* so that: A * inv(A) = I
*/
Inverse(A=0){
if (this != Matrix)
A := this
return new Matrix(base.Inverse(A))
}

/*
* Multiplies this Matrix with B from right and returns the new matrix
*
*/
MultiplyRight(B){
return new Matrix(base.Multiply(this, B))
}

/*
* Multiplies this Matrix with B from left and returns the new matrix
*
*/
MultiplyLeft(B){
return new Matrix(base.Multiply(B, this))
}

/**
* Returns the Column-Vector at the given index
*/
RangeCol(m, colIndex=0){

if (this != Matrix)
{
colIndex := m
m := this
}

return new Matrix(base.RangeCol(m, colIndex))
}

/*
*  Returns the Row-Vector at the given index
*/
RangeRow(m,rowIndex=0){

if (this != Matrix)
{
rowIndex := m
m := this
}

return new Matrix(base.RangeRow(m, rowIndex))
}

/*
* Returns a transposed Matrix of this matrix
*/
Transpose(m=0){
if (this != Matrix)
m := this
return new Matrix(base.Transpose(m))
}

IsSquare(m=0){
if (this != Matrix)
m := this
return base.IsSquare(m)
}

/**
* Returns the count of columns in the given Matrix
*/
ColumnCount(m=0){
if (this != Matrix)
m := this
return base.ColumnCount(m)
}

/**
* Returns the count of rows in the given Matrix
*/
RowCount(m=0){
if (this != Matrix)
m := this
return base.RowCount(m)
}

/**
* Clones the given Matrix
*/
Clone(m=0){
if (this != Matrix)
m := this
return new Matrix(m)
}

Prototype(m){
for each, row in m
{
if(IsObject(row))
{
rIndex := A_index
for each, item in row
this[rIndex,A_index] := item
}else
this[A_index] := row
}
return mc
}

Equals(m,m2=0){
if (this != Matrix)
{
m2 := m
m := this
}
return base.Equals(m,m2)
}

/*
* Returns a console/msgbox friendly string. Useful for debugging
*/
ToString(m=0){
if (this != Matrix)
m := this
return base.ToString(m)
}

/*
* Generates an quadratic identity matrix with a size of [n]
*/
Eye(n){
return new Matrix(base.Eye(n))
}

Zeros(n){
return new Matrix(base.Zeros(n))
}

/*
* Generates an quadratic matrix with a size of [n]
* and each element set to [fillNum]
*/
Fill(n, fillNum){
return  new Matrix(base.Fill(n, fillNum))
}

Mirror2D(m){
return new Matrix(base.Mirror2D(m))
}

/**
* Gets the mirror matrix for the straight line intersection the origin (0,0)
* angle     =    angle to x axis, in radians
*/
Mirror2DByAngle(angle){
return new Matrix(base.Mirror2DByAngle(angle))
}

Rotate2D(angle){
return new Matrix(base.Rotate2D(angle))
}

Rotate3DZ(angle){
return new Matrix(base.Rotate3DZ(angle))
}

Rotate3DY(angle){
return new Matrix(base.Rotate3DY(angle))
}

Rotate3DX(angle){
return new Matrix(base.Rotate3DX(angle))
}

Mirror3D(axis){
return new Matrix(base.Mirror3D(axis))
}

Gauss(A, B=""){
if(B="")
{
B := A
A := this
}
return new Matrix(base.Gauss(A,B))
}

ToRowEchelonForm(a, b=""){
if (this != Matrix)
{
b := a
a := this
}
return base.ToRowEchelonForm(a, b)
}
}
Notes:
Performance optimization may follow too, but actually you should use a native implementation of matrix arithmetic if you're after super performance.

Show-case Script:

A := [[1,0,2]
,[2,3,1]]

B := [[1,2]
,[3,4]
,[5,6]
,[7,8]]

C := [[1]
,[4]]

D := [1,2,3,4]

E := [[2,3]
,[1,0]]

F := [[1,0,2]
,[2,3,1]
,[8,7,6]]

H := [[1,3,5,-1]
,[2,2,3, 0]
,[0,1,2, 3]
,[4,5,3, 2]]

zeros := Matrix.Zeros(6)
MsgBox % zeros.ToString()

identity := Matrix.Eye(6)
MsgBox % identity.ToString()

MsgBox % identity.Equals(zeros) ? "the matrices are equal" : "the matrices are NOT equal"

At := Matrix.Transpose(A)
MsgBox % Matrix.ToString(A) "`n`nTransposed`n" At.ToString()

P := Matrix.Multiply(E, A)
MsgBox % Matrix.ToString(E) "`nmultiplicated by `n" Matrix.ToString(A) "`nyields:" Matrix.ToString(P)

MsgBox % "A added to A" Matrix.ToString(AplusA)

Amul3 := Matrix.MultiplyScalar(A,3)
MsgBox % "A multiplicated by 3" Matrix.ToString(Amul3)

det := Matrix.Det(H) ; calculate the determinant of matrix H
MsgBox % Matrix.ToString(H) "`nthe determinant is:`n" det

Hm := new Matrix(H) ; you can also use instances of the matrix class
MsgBox % "using the matrix object:`n" Hm.ToString() "`nthe determinant is:`n" Hm.Det()

;
; 1. define a mirror matrix
; 2. define a rotation matrix
; 3. combine them to a single transformation matrix by multiplicating them
; 4. use it to transform a single 2D vertex
;

;1.
;y = 2*x
m := 2
M := Matrix.Mirror2D(m)

;2.
angleDegree := 90
radAngle :=  angleDegree * (45 / atan(1)) ; 90° degree to radians

;3. T = R * M
T := R.MultiplyRight(M)

;4. v --> vertex

v := [-2,5]

; v' = T * v
vTransformed := T.MultiplyRight(v)

; gaussian calculations (code provided by Babba/horst)
a:=[[2,0,1]
,[0,3,0]
,[10,0,-4]]

; calc the inverse of an matrix
msgbox % Matrix.Inverse(a).ToString()

/*
Solve a linear equation system like:

3x + 2y - z   =  1
2x - 2y + 4z  = -2
-x + 1/2y - z =  0

expected result:
x = 1
y = -2
z = -2

*/
s := [[ 3, 2, -1]
,[ 2,-2,  4]
,[-1,1/2,-1]]

b := [1
,-2
,0]

x := Matrix.Gauss(s,b)

MsgBox %  x.ToString()

ExitApp

Hope some of you can use it

### #2 IsNull

IsNull
• Fellows
• 990 posts

Posted 14 May 2012 - 04:54 PM

I've added support for calculating the determinant of any matrix. This was achieved whit the Laplace expansion.

The determinant is the first step (and probably the last) for a linear equation solver.

Btw: I've to learn this in my math class, so the best way to understand something is to automate it. Thus I write a Lib for it. :mrgreen:

### #3 infogulch

infogulch
• Moderators
• 717 posts

Posted 14 May 2012 - 05:22 PM

Nice! Is it possible to support N-dimensional matrices?

Also, why only static references to the Matrix class? (I.e. no constructors.)

### #4 IsNull

IsNull
• Fellows
• 990 posts

Posted 14 May 2012 - 05:53 PM

Thanks for your feedback! Glad to see that those strange math excursions take up some attention. :mrgreen:

Nice! Is it possible to support N-dimensional matrices?

What do you mean exactly by N-Dimensional? The lib supports any dimension Matrices, but some Operations are simply not defined for non-square Matrices, such as the determinant.
All others support any (n,m) Matrices when they are in a valid combination for the given operation.

Also, why only static references to the Matrix class? (I.e. no constructors.)

This is something I thought about myself. I'll add support for real matrix objects, but the definition of a matrix would be more complicated. Eg.:

B := [[1,2]
,[3,4]
,[5,6]
,[7,8]]
..would become...
B := new Matrix([[1,2]
,[3,4]
,[5,6]
,[7,8]])

However, for the sake of "good" OOD, I'll support both usages by a little trick soon.

### #5 infogulch

infogulch
• Moderators
• 717 posts

Posted 14 May 2012 - 05:59 PM

B := [[1,2]
,[3,4]
,[5,6]
,[7,8]]

B.base := Matrix

About that trick, were you thinking of something like this?
MultiplyScalar(A, n = "") {
if (this != Matrix) ; i.e. it was not invoked with Matrix.MultiplyScalar(x,y)
n := A, A := this
...
I think that should do it.

### #6 IsNull

IsNull
• Fellows
• 990 posts

Posted 14 May 2012 - 07:04 PM

Done. Supports now Instance and Static usage.

About that trick, were you thinking of something like this?

Yup, exactly. Additionally, I switch the method parameters for consistency, when the last was omitted.

### #7 Uberi

Uberi
• Fellows
• 1050 posts

Posted 19 May 2012 - 01:40 AM

Very useful. Would it be possible to directly support importing OpenGL or DirectX matrixes? That would make certain things such as quaternions much easier when working in 3D.

I ought to brush up on my matrix math now

### #8 IsNull

IsNull
• Fellows
• 990 posts

Posted 19 May 2012 - 09:52 AM

Hey uberi,

Would it be possible to directly support importing OpenGL or DirectX matrixes?

In which format/data structure are they?

Parsing Strings / Data structures into the Matrix shouldn't be a problem. If it gets very specific to a certain problem domain, I may move the parser/wrapper into a separate MatrixFactory/Builder class.

Btw, I can also offer to implement the common translation base matrices such as:

[*:123afw8f] 2D: rotation
[*:123afw8f] 2D: mirroring along a given straight line (y=mx+
[*:123afw8f] 3D: mirroring matrices for x, y and z
[*:123afw8f] 3D: rotation matrices for x, y and z
[*:123afw8f] ...

### #9 Uberi

Uberi
• Fellows
• 1050 posts

Posted 19 May 2012 - 02:21 PM

Hey uberi,

Would it be possible to directly support importing OpenGL or DirectX matrixes?

In which format/data structure are they?

I believe OpenGL matrices are 4x4 column major, stored as a single array of floats. That makes it an array of 16 floats (not sure if they support doubles). I seem to recall that Direct3D matrices are row major, but it should still be the same layout in memory.

Btw, I can also offer to implement the common translation base matrices such as:

[*:1pmyhmlw] 2D: rotation
[*:1pmyhmlw] 2D: mirroring along a given straight line (y=mx+
[*:1pmyhmlw] 3D: mirroring matrices for x, y and z
[*:1pmyhmlw] 3D: rotation matrices for x, y and z
[*:1pmyhmlw] ...

That would be really useful .

### #10 IsNull

IsNull
• Fellows
• 990 posts

Posted 23 May 2012 - 01:19 PM

I've added several static Transformation helper Methods to the the Matrix class. They basically generate transformation Matrixes which can be applied (Multiplied) to vertices to make the actual transformation.

A snippet from the new code:
class Matrix
{
;...

/**
* Gets the mirror matrix for the straight line mirror y;
* y := m*x + b
*
*
* if b is NOT zero, you must move the vertices first with v = (0,-b)
* 1. move vertices along v = (0,-b)  | --> move the mirror axis to the origin
* 2. mirror with the matrix returned by this function |---> newVertex = M * vertex
* 3. move vertices back along -v = (0,b) | --> move the vertices back
*/
Mirror2D(m){
angle := this.Mirror2DByAngle(ATan(m))
}

/**
* Gets the mirror matrix for the straight line intersection the origin (0,0)
* angle     =    angle to x axis, in radians
*/
Mirror2DByAngle(angle){
angle *= 2
c := cos(angle)
s := sin(angle)

M := [[c, s]
,[s,-c]]

return M
}

Rotate2D(angle){
c := cos(angle)
s := sin(angle)

R := [[c,-s]
,[s, c]]

return R
}

Rotate3DZ(angle){
c := cos(angle)
s := sin(angle)

R := [[c,-s, 0]
,[s, c, 0]
,[0, 0, 1]]

return R
}

Rotate3DY(angle){
c := cos(angle)
s := sin(angle)

R := [[ c, 0, s]
,[ 0, 1, 0]
,[-s, 0, c]]

return R
}

Rotate3DX(angle){
c := cos(angle)
s := sin(angle)

R := [[ 1, 0, 0]
,[ 0, c,-s]
,[ 0, s, c]]

return R
}

Mirror3D(axis){

m1 := (axis == 1 || axis = "x") ? -1 : 1
m2 := (axis == 2 || axis = "y") ? -1 : 1
m3 := (axis == 3 || axis = "z") ? -1 : 1

M := [[ m1, 0,  0]
,[ 0, m2,  0]
,[ 0,  0, m3]]

return M
}

;...
}

### #11 ironcates

ironcates
• Members
• 32 posts

Posted 11 June 2012 - 05:08 AM

You probably already know this, but the adj(A) would be an excellent step toward the inv(A) since inv(A)=adj(A)/det(A).

<!-- m -->http://en.wikipedia....Adjugate_matrix<!-- m -->

### #12 IsNull

IsNull
• Fellows
• 990 posts

Posted 11 June 2012 - 09:11 AM

AFAIK, Gauss is the fastest possibility to solve equations. Generating an inverse of a matrix is nothing other than solving a big equation system.
In the German Forum, there are examples of it: <!-- m -->http://de.autohotkey...opic.php?t=9753<!-- m -->

I may take this up to create the inverse - its more about time restrictions than about missing opportunities.

### #13 Uberi

Uberi
• Fellows
• 1050 posts

Posted 18 June 2012 - 12:57 AM

Just want to mention, I've implemented Guassian elimination for row-major matrices:

LinearSystemSolve(Matrix)
{
Loop, % Matrix[1].MaxIndex()
{
Column := A_Index
Loop, % Matrix.MaxIndex() - Column
{
Row := Column + A_Index
Values := []
For Index, Value In Matrix[Column]
Values[Index] := Value * -(Matrix[Row][Column] / Matrix[Column][Column])
For Index, Value In Values
Matrix[Row][Index] += Value
}
}

Output := []
Index := 0
RowIndex := Matrix.MaxIndex()
Loop, %RowIndex%
{
;substitute in all known coefficients
Row := Matrix[RowIndex - Index]
ColumnIndex := Row.MaxIndex()
Inner := 0
Loop, %Index%
Inner += Output[A_Index] * Row[ColumnIndex - A_Index]
Index ++
Output.Insert((Row[ColumnIndex] - Inner) / Row[ColumnIndex - Index])
}

;reverse the matrix rows
Temp1 := []
For Index, Value In Output
Temp1.Insert(1,Value)
Return, Temp1
}

Feel free to adapt it to your needs.

It's a lightly optimized translation of this code.

### #14 IsNull

IsNull
• Fellows
• 990 posts

Posted 19 June 2012 - 04:30 PM

Hey Uberi,

Thanks for your solution. It seems very similiar to my "compact" solution here:
<!-- m -->http://de.autohotkey... ... 6308#76308<!-- m -->

But like my version, you are using divison of the elements whithout any pre-analysis. When the equation has Zeros it wont work; Eg:
s := [[0, 1, 1, 15]
,  [2, 3, 1, 28]
,  [3,-4, 2,  3]]

Pivot strategies[1] solve the issue. A very simple approach is to change the order of the equations so that the ones whit very small numbers (especially Zeros) are at the bottom in our matrix.
I'll post a slightly adapted version soon

[1] <!-- m -->http://www.math.okst... ... 13-l12.pdf<!-- m -->

### #15 Uberi

Uberi
• Fellows
• 1050 posts

Posted 19 June 2012 - 04:48 PM

Here's what I've done since posting the above:

LinearSystemSolve(Matrix)
{
;wip: check input

Unknowns := Matrix.MaxIndex()
Column := 0
Loop, % Unknowns - 1
{
Column ++

;find the row with the largest value in the current column
RowMaximum := Column
RowIndex := Column
While, RowIndex < Unknowns ;iterate through each row after the current one
{
RowIndex ++ ;obtain the row index
If Abs(Matrix[RowIndex][Column]) > Abs(Matrix[RowMaximum][Column]) ;value greater than previous maximum
RowMaximum := RowIndex
}

;swap the maximum row and the current row
Row := Matrix[Column]
Matrix[Column] := Matrix[RowMaximum]
Matrix[RowMaximum] := Row

;check if the matrix has a singular unique solution
If Abs(Matrix[Column][Column]) < 0.000001
throw Exception("Could not find singular unique solution.")

;eliminate the unknown in the current column
RowIndex := Column
While, RowIndex < Unknowns ;iterate through each row after the current one
{
RowIndex ++ ;obtain the row index
Index := Unknowns + 1
While, Index >= Column
{
Matrix[RowIndex][Index] -= (Matrix[Column][Index] * Matrix[RowIndex][Column]) / Matrix[Column][Column]
Index --
}
}
}

;substitute known coefficients for unknowns
Result := []
RowIndex := Unknowns
While, RowIndex > 0 ;loop through each unknown backwards
{
Row := Matrix[RowIndex]
Value := 0
Index := Unknowns
While, Index > RowIndex
Value += Result[Index] * Row[Index], Index --
Result[RowIndex] := (Row[Unknowns + 1] - Value) / Row[RowIndex]
RowIndex --
}
Return, Result
}

Significantly faster than before and has partial pivoting included.