matrix functions

Post your working scripts, libraries and tools for AHK v1.1 and older
User avatar
jeeswg
Posts: 6902
Joined: 19 Dec 2016, 01:58
Location: UK

matrix functions

12 May 2019, 12:48

I wanted some functions to do matrix multiplication and matrix transposition.

I also provide an explanation of matrix multiplication and some info regarding doing matrix calculations on the Wolfram Alpha website.

An explanation of matrix multiplication:

Code: Select all

;e.g. matrix multiplication:

;[1 2 3] . [7 10] = [50 68]
;[4 5 6]   [8 11]   [122 167]
;          [9 12]

;to visualise it, move the second matrix upwards:
;        [7 10]
;        [8 11]
;        [9 12]
;[1 2 3] [a b]
;[4 5 6] [c d]
;where a = (1*7 + 2*8 + 3*9) = 50
;where b = (1*10 + 2*11 + 3*12) = 68
;where c = (4*7 + 5*8 + 6*9) = 122
;where d = (4*10 + 5*11 + 6*12) = 167

;in Wolfram Alpha:
;{{1, 2, 3}, {4, 5, 6}} . {{7, 10}, {8, 11}, {9, 12}} - Wolfram|Alpha
;https://www.wolframalpha.com/input/?i=%7B%7B1,+2,+3%7D,+%7B4,+5,+6%7D%7D+.+%7B%7B7,+10%7D,+%7B8,+11%7D,+%7B9,+12%7D%7D

Some matrix functions:

Code: Select all

q:: ;matrix transposition and matrix multiplication
oMtx1 := [[1, 2, 3]
	, [4, 5, 6]]

oMtx2 := [[1, 4]
	, [2, 5]
	, [3, 6]]

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

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

;check if matrix is valid:
;oMtxInvalid := ObjClone(oMtx4)
;oMtxInvalid.Push([5])
;MsgBox, % JEE_MtxIsValid(oMtx4)
;MsgBox, % JEE_MtxIsValid(oMtxInvalid)

;equivalent to oMtx3 definition above
oMtx3 := JEE_MtxFromString("
(
1 2 3
4 5 6
7 8 9
)")

MsgBox, % "matrix display:`r`n`r`n" JEE_MtxDisp(oMtx3) "`r`n`r`n" JEE_MtxDispWA(oMtx3)

MsgBox, % "matrix dimensions:`r`n"
. "`r`n" JEE_MtxWidth(oMtx1) "x" JEE_MtxHeight(oMtx1)
. "`r`n" JEE_MtxWidth(oMtx2) "x" JEE_MtxHeight(oMtx2)
. "`r`n" JEE_MtxWidth(oMtx3) "x" JEE_MtxHeight(oMtx3)
. "`r`n" JEE_MtxWidth(oMtx4) "x" JEE_MtxHeight(oMtx4)

vCount := 4
Loop, % vCount
{
	oMtxTemp := JEE_MtxTranspose(oMtx%A_Index%)
	vDisp1 := JEE_MtxDisp(oMtx%A_Index%)
	vDisp2 := JEE_MtxDisp(oMtxTemp)
	MsgBox, % "matrix transposition:`r`n`r`n" vDisp1 "`r`n`r`n" vDisp2
}

vCount := 4
Loop, % vCount
{
	vNum := A_Index
	Loop, % vCount
	{
		oMtxTemp := JEE_MtxMul(oMtx%vNum%, oMtx%A_Index%)
		vDisp1 := JEE_MtxDisp(oMtx%vNum%)
		vDisp2 := JEE_MtxDisp(oMtx%A_Index%)
		vDisp3 := IsObject(oMtxTemp) ? JEE_MtxDisp(oMtxTemp) : "`tINVALID"
		vDisp2 := StrReplace(vDisp2, "[", "`t[")
		vDisp3 := StrReplace(vDisp3, "[", "`t[")
		MsgBox, % "matrix multiplication:`r`n`r`n" vDisp2 "`r`n.`r`n" vDisp1 "`r`n=`r`n" vDisp3
	}
}
return

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

;note: checks that each row has the same number of values
JEE_MtxIsValid(oMtx)
{
	local
	vWidth := oMtx.1.Length()
	for _, oRow in oMtx
		if !(oRow.Length() = vWidth)
			return 0
	return 1
}

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

JEE_MtxWidth(oMtx)
{
	local
	return oMtx.1.Length()
}

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

JEE_MtxHeight(oMtx)
{
	local
	return oMtx.Length()
}

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

JEE_MtxMul(oMtx1, oMtx2)
{
	local
	vW1 := JEE_MtxWidth(oMtx1)
	vH2 := JEE_MtxHeight(oMtx2)
	if !(vW1 = vH2)
		return ""
	vW2 := JEE_MtxWidth(oMtx2)
	vH1 := JEE_MtxHeight(oMtx1)
	oMtx3 := []
	Loop, % vH1
	{
		oMtx3.Push([])
		vNum1 := A_Index
		Loop, % vW2
		{
			vSum := 0
			vNum2 := A_Index
			Loop, % vW1
				vSum += oMtx1[vNum1, A_Index] * oMtx2[A_Index, vNum2]
			oMtx3[vNum1, A_Index] := vSum
		}
	}
	return oMtx3
}

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

JEE_MtxDisp(oMtx, vSep:=", ", vPfx:="[", vSfx:="]", vSepLn:="`r`n")
{
	local
	vOutput := ""
	for vRow, oTemp in oMtx
	{
		vOutput .= vPfx
		for vCol, vValue in oTemp
			vOutput .= (A_Index=1?"":vSep) vValue
		vOutput .= vSfx vSepLn
	}
	return SubStr(vOutput, 1, -StrLen(vSepLn))
}

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

;Wolfram|Alpha Examples: Matrices
;https://www.wolframalpha.com/examples/mathematics/algebra/matrices/

;{{1, 2, 3}, {4, 5, 6}} . {{1, 4}, {2, 5}, {3, 6}} - Wolfram|Alpha
;https://www.wolframalpha.com/input/?i=%7B%7B1,+2,+3%7D,+%7B4,+5,+6%7D%7D+.+%7B%7B1,+4%7D,+%7B2,+5%7D,+%7B3,+6%7D%7D
;{{1, 4}, {2, 5}, {3, 6}} . {{1, 2, 3}, {4, 5, 6}} - Wolfram|Alpha
;https://www.wolframalpha.com/input/?i=%7B%7B1,+4%7D,+%7B2,+5%7D,+%7B3,+6%7D%7D+.+%7B%7B1,+2,+3%7D,+%7B4,+5,+6%7D%7D
;{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}^2 - Wolfram|Alpha
;https://www.wolframalpha.com/input/?i=%7B%7B1,+2,+3%7D,+%7B4,+5,+6%7D,+%7B7,+8,+9%7D%7D%5E2

;note: for Wolfram Alpha
JEE_MtxDispWA(oMtx)
{
	local
	return "{" JEE_MtxDisp(oMtx, ", ", "{", "}", ", ") "}"
}

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

JEE_MtxFromString(vText)
{
	local
	vText := StrReplace(vText, ",", " ")
	vText := RegExReplace(vText, " +", " ")
	oMtx := []
	for vRow, vRowText in StrSplit(vText, "`n", "`r")
		oMtx.Push(StrSplit(vRowText, " "))
	return oMtx
}

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

JEE_MtxTranspose(oMtx)
{
	local
	oMtx2 := []
	for vRow, oRow in oMtx
		for vCol, vValue in oRow
			oMtx2[vCol, vRow] := vValue
	return oMtx2
}

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

Links:
[AHK_L] Matrix class - Scripts and Functions - AutoHotkey Community
https://autohotkey.com/board/topic/80487-ahk-l-matrix-class/
Matrix multiplication - Rosetta Code
https://rosettacode.org/wiki/Matrix_multiplication#AutoHotkey
MS Paint: rotate by scaling/shearing (stretching/skewing) - AutoHotkey Community
https://autohotkey.com/boards/viewtopic.php?f=6&t=64481

See also:
Wolfram|Alpha: Computational Intelligence
https://www.wolframalpha.com/
complex numbers - AutoHotkey Community
https://autohotkey.com/boards/viewtopic.php?f=6&t=35809
homepage | tutorials | wish list | fun threads | donate
WARNING: copy your posts/messages before hitting Submit as you may lose them due to CAPTCHA
nnrxin
Posts: 8
Joined: 07 May 2018, 06:02

Re: matrix functions

05 Jan 2020, 09:36

[AHK_L] Matrix class - Scripts and Functions - AutoHotkey Community
https://autohotkey.com/board/topic/80487-ahk-l-matrix-class/

this Matrix class has a bug! :!:

Code: Select all


C := [[6]
       ,[5]
       ,[4]]

m := [3,2,1]

ddd := MatrixStatic.Multiply(c , m)
MsgBox % ddd[1,1] " " ddd[1,2] " " ddd[1,3] "`r`n" ddd[2,1] " " ddd[2,2] " " ddd[2,3] "`r`n" ddd[3,1] " " ddd[3,2] " " ddd[3,3]

ddd := MatrixStatic.Multiply(m , c)
MsgBox % ddd[1,1] " " ddd[1,2] " " ddd[1,3] "`r`n" ddd[2,1] " " ddd[2,2] " " ddd[2,3] "`r`n" ddd[3,1] " " ddd[3,2] " " ddd[3,3]


/**
* AHK static Matrix calculations
* by IsNull 2012
*
* Gauss/Pivot Strategies implemented by horst/Babba
*/
class MatrixStatic
{
	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  Broken Link for safety
         
         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  Broken Link for safety
   */
   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)
   }
}
after rewrite these functions, the bug disappeared :lol:

Code: Select all


	RangeCol(m, colIndex)
	{
		rowN := MatrixStatic.RowCount(m)
		col := []
		Loop % rowN
		{
			if (rowN = 1)
				col[1,1] := m[colIndex]
			else
				col[A_index,1] := m[A_index,colIndex]
		}
		return col
	}

	RangeRow(m, rowIndex)
	{
		rowN := MatrixStatic.RowCount(m)
		if (rowN = 1)
			return MatrixStatic.Clone(m)
		else
			return MatrixStatic.Clone(m[rowIndex])
	}



Return to “Scripts and Functions (v1)”

Who is online

Users browsing this forum: No registered users and 164 guests