16 replies to this topic

[IsNull]

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:

added := Matrix.Add([1,2],[4,5])

or

A := new Matrix([1,2])
B := new Matrix([4,5])
added := A.Add(B)

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){   
      add := {}

      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
         }
      }
      return add
   }
   
   /**
   * Addition of A and B
   */
   Add(A, B){
      if(MatrixStatic.ColumnCount(A) != MatrixStatic.ColumnCount(B) 
      || MatrixStatic.RowCount(A) != MatrixStatic.RowCount(B))
      {
         throw Exception("Matrix Addition Error: All dimensions must agree!")
      }
      
      add := {}

      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]
         }
      }
      return add
   }
   
   /**
   * 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!")
      }
      
      add := {}

      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]
         }
      }
      return add
   }
   
	/*
	* 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
   */
   Add(m, B=0){
      
      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)


AplusA := Matrix.Add(A,A)
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
R := Matrix.Rotate2D(radAngle)

;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

[IsNull]

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:

[infogulch]

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

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

[IsNull]

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.

[infogulch]

How about:

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.

[IsNull]

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.

[Uberi]

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

[IsNull]

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

[Uberi]

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 .

[IsNull]

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
   }

;...
}

[ironcates]

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

[IsNull]

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.

[Uberi]

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.

[IsNull]

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

[Uberi]

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.