/** This is an incomplete class that implements matrix operations. Please feel free to implement missing functions! */ class Matrix { /** The array that will actually store the matrix data. Note that this is a raw 2-dimensional Frink array, and the indices are zero-based in contrast with most matrix notation which is one-based. The class should try to enforce that this is a rectangular 2-dimensional array, but the code doesn't currently enforce that with some constructors. */ var array /** The number of rows in the matrix. */ var rows /** The number of columns in the matrix. */ var cols /** Construct a new Matrix with the specified number of rows and columns and the specified default value for each element. */ new[rowCount, colCount, default=0] := { rows = rowCount cols = colCount array = makeArray[[rows, cols], default] } /** Construct a new Matrix from a rectangular 2-dimensional array. This currently does not verify that the array is rectangular. */ new[a is array] := { rows = length[a] cols = length[a@0] array = a } /** Sets the value of the specified element, using 1-based coordinates, as is common in matrix notation. Stupid mathematicians. */ set[row, col, val] := array@(row-1)@(col-1) = val /** Gets the value of the specified element, using 1-based coordinates, as is common in matrix notation. */ get[row, col] := array@(row-1)@(col-1) /** Gets the number of columns. */ getCols[] := cols /** Gets the number of rows. */ getRows[] := rows /** Multiply two matrices. The number of columns in matrix a should equal the number of rows in matrix b. The resulting matrix will have the same number of rows as matrix a and the same number of columns as matrix b. TODO: Use a faster algorithm for large arrays? This is O(n^3). TODO: Optimize this for small arrays by precalculating. */ class multiply[a, b] := { if a.cols != b.rows return undef else items = a.cols resRows = a.rows resCols = b.cols resultArray = makeArray[[resRows, resCols], 0] aa = a.array bb = b.array for row = 0 to resRows-1 for col = 0 to resCols-1 { sum = aa@row@0 * bb@0@col // Make units come out right for k = 1 to items-1 sum = sum + aa@row@k * bb@k@col resultArray@row@col = sum } return new Matrix[resultArray] } /** Multiply this matrix by the specified matrix and return the result. */ multiply[b] := multiply[this, b] /** Multiplies all elements of a matrix by a scalar. */ multiplyByScalar[s] := { a = makeArray[[rows, cols]] for rowNum = 0 to rows-1 { row = array@rowNum for colNum = 0 to cols-1 a@rowNum@colNum = row@colNum * s } return new Matrix[a] } /** Multiplies a row in place by the specified value. (indices are 1-based.) */ multiplyRowByScalar[row, val] := { r = row-1 array@r = mul[array@r, val] } /** Transposes the Matrix and returns the new transposed Matrix. This uses the built-in array.transpose[] method. */ transpose[] := new Matrix[array.transpose[]] /** Returns true if this Matrix is square, false otherwise. */ isSquare[] := rows == cols /** Exchange the specified-numbered rows in the matrix. The rows are 1-indexed to match standard matrix notation. */ exchangeRows[row1, row2] := { r1 = row1-1 r2 = row2-1 [array@r1, array@r2] = [array@r2, array@r1] } /** Exchange the specified-numbered columns in the matrix. The columns are 1-indexed to match standard matrix notation. */ exchangeCols[c1, c2] := { for r = 0 to rows-1 { temp = array@r@(c2-1) array@r@(c2-1) = array@r@(c1-1) array@r@(c1-1) = temp } } /** Returns a matrix as a string with rows separated with newlines. */ toString[] := { result = "" for r = 0 to rows-1 result = result + " " + array@r + "\n" return result } /** Returns the matrix formatted by formatTable with rows separated with newlines. */ toTable[] := { return formatTable[array] } /** Formats a matrix with Unicode box drawing characters. The result is a single string. */ formatMatrix[] := { return formatMatrix[array] } /** Formats a matrix with Unicode box drawing characters. The result is fairly compact with blank lines separating rows and 2 spaces seperating columns. The result is a single string. */ formatMatrixCompact[] := { return formatMatrixCompact[array] } /** Formats a matrix with Unicode box drawing characters. The result is as compact as possible with no blank lines separating rows and 1 space separating columns. The result is a single string. */ formatMatrixVeryCompact[] := { return formatMatrixVeryCompact[array] } /** Returns the matrix formatted as a matrix with Unicode box drawing characters. The result is an array of lines that can be further formatted. */ formatMatrixLines[] := { m = formatTableLines[array] rows = length[m] width = length[m@0] result = new array result.push["\u250c" + repeat[" ", width] + "\u2510"] for n=0 to rows-1 result.push["\u2502" + m@n + "\u2502"] result.push["\u2514" + repeat[" ", width] + "\u2518"] return result } /** Gets the specified (1-based) column as a row array. */ getColumnAsArray[col] := { return deepCopy[array.getColumn[col-1]] } /** Gets the specified (1-based) column as a 1-column matrix */ getColumnAsMatrix[col] := { return new Matrix[deepCopy[array.getColumn[col-1].transpose[]]] } /** Gets the specified (1-based) row as a row array. */ getRowAsArray[row] := { return deepCopy[array@(row-1)] } /** Creates the LU (lower-upper) decomposition of the matrix. This creates two matrices, L and U, where L has the value 1 on the diagonal from top left to bottom right, like: 1 0 0 L = L21 1 0 L31 L32 1 U11 U12 U13 U = 0 U22 U23 0 0 U33 This is basically like solving equations by Gaussian elimination. */ LUDecompose[] := { return LUDecomposeCrout[] } /** This uses Crout's method to decompose a square matrix into an lower and upper triangular matrix. See: https://en.wikipedia.org/wiki/Crout_matrix_decomposition This creates two matrices, L and U, where U has the value 1 on the diagonal from top left to bottom right, like: L11 0 0 L = L21 L22 0 L31 L32 L33 1 U12 U13 U = 0 1 U23 0 0 1 returns [L, U] The original matrix can be obtained as L.multiply[U] Note that the Crout algorithm fails on certain matrices, for example m = new Matrix[[[1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[-1,1,1,1]]] See: https://semath.info/src/inverse-cofactor-ex4.html for more on that matrix. https://en.wikipedia.org/wiki/Crout_matrix_decomposition */ LUDecomposeCrout[] := { n = rows L = new array[[rows, cols], 0] U = new array[[rows, cols], 0] sum = 0 for i=0 to n-1 U@i@i = 1 for j=0 to n-1 { for i=j to n-1 { sum = 0 for k=0 to j-1 sum = sum + L@i@k * U@k@j L@i@j = array@i@j - sum } for i=j to n-1 { sum = 0 for k=0 to j-1 sum = sum + L@j@k * U@k@i if L@j@j == 0 { println["Matrix.LUDecomposeCrout: det(L) is zero. Can't divide by zero."] return undef } U@j@i = (array@j@i - sum) / L@j@j } } return [new Matrix[L], new Matrix[U]] } /** Another algorithm for Crout LU decomposition. See: http://www.mosismath.com/Matrix_LU/Martix_LU.html This returns [L, U] as matrices which may be transposes of the matrices returned by LUDecomposeCrout[] . See its notes for more about this decomposition and its properties. */ LUDecomposeCrout2[] := { if ! isSquare[] { println["Matrix.LUDecomposeCrout2 only works on square arrays."] return undef } n = rows L = new array[[n,n], 0] U = new array[[n,n], 0] a = 0 for i = 0 to n-1 L@i@i = 1 for j = 0 to n-1 { for i = 0 to n-1 { if i >= j { U@j@i = array@j@i for k = 0 to j-1 { U@j@i = U@j@i - U@k@i * L@j@k } } if i > j { L@i@j = array@i@j for k = 0 to j-1 L@i@j = L@i@j - U@k@j * L@i@k diag = U@j@j if diag == 0 { println["Matrix.LUDecomposeCrout2: det(L) is zero. Can't divide by zero."] println[formatMatrix[L]] println[formatMatrix[U]] return undef } L@i@j = L@i@j / diag } } } return [new Matrix[L], new Matrix[U]] } /** This uses the Cholesky-Banachiewicz algorithm to decompose a square Hermitian matrix into a lower triangular matrix. If you transpose the lower triangular matrix L by L.transpose[], you get an upper triangular matrix which is symmetrical to the lower triangular matrix. Multiplying the lower triangular matrix by the upper triangular matrix, that is, L.multiply[L.transpose[]] gives you back the original matrix! See: https://en.m.wikipedia.org/wiki/Cholesky_decomposition */ CholeskyB[] := { if ! isHermitian[] { println["Matrix.CholeskyB only works on square Hermitian matrices."] return undef } n = rows L = new array[[n,n], 0] for i = 0 to n-1 for j = 0 to i { sum = 0 for k = 0 to j-1 sum = sum + L@i@k * L@j@k if i == j L@i@j = sqrt[array@i@i - sum] else L@i@j = (1 / L@j@j * (array@i@j -sum)) } return new Matrix[L] } /** This uses the Cholesky-Crout algorithm to decompose a square Hermitian matrix into a lower triangular matrix. If you transpose the lower triangular matrix L by L.transpose[], you get an upper triangular matrix which is symmetrical to the lower triangular matrix. Multiplying the lower triangular matrix by the upper triangular matrix, that is, L.multiply[L.transpose[]] gives you back the original matrix! See: https://en.m.wikipedia.org/wiki/Cholesky_decomposition */ CholeskyCrout[] := { if ! isHermitian[] { println["Matrix.CholeskyCrout only works on square Hermitian matrices."] return undef } n = rows L = new array[[n,n], 0] for j = 0 to n-1 { sum = 0 for k = 0 to j-1 sum = sum + (L@j@k)^2 L@j@j = sqrt[array@j@j - sum] for i = j+1 to n-1 { sum = 0 for k = 0 to j-1 sum = sum + L@i@k * L@j@k L@i@j = (1 / L@j@j * (array@i@j -sum)) } } return new Matrix[L] } /** This uses Doolittle's method to decompose a square matrix into an lower and upper triangular matrix. This creates two matrices, L and U, where L has the value 1 on the diagonal from top left to bottom right, like: 1 0 0 L = L21 1 0 L31 L32 1 U11 U12 U13 U = 0 U22 U23 0 0 U33 The original matrix can be obtained as L.multiply[U] Note that the Doolittle algorithm fails on certain matrices, for example m = new Matrix[[[1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[-1,1,1,1]]] See: https://semath.info/src/inverse-cofactor-ex4.html for more on that matrix. See: https://www.geeksforgeeks.org/doolittle-algorithm-lu-decomposition/ */ LUDecomposeDoolittle[] := { if ! isSquare[] { println["Matrix.LUDecomposeDoolittle only works on square arrays."] return undef } n = rows L = new array[[n,n], 0] U = new array[[n,n], 0] // Decomposing matrix into Upper and Lower // triangular matrix for i = 0 to n-1 { // Upper Triangular for k = i to n-1 { // Summation of L(i, j) * U(j, k) sum = 0 for j = 0 to i-1 sum = sum + (L@i@j * U@j@k) // Evaluating U(i, k) U@i@k = array@i@k - sum } // Lower Triangular for k = i to n-1 { if i == k L@i@i = 1 // Diagonal as 1 else { // Summation of L(k, j) * U(j, i) sum = 0 for j = 0 to i-1 sum = sum + L@k@j * U@j@i diag = U@i@i if (diag == 0) { println["Matrix.LUDecomposeDoolittle: U@$i@$i is zero when solving for L@$k@$i. Can't divide by zero."] println[formatMatrix[L]] println[formatMatrix[U]] return undef } // Evaluating L(k, i) L@k@i = (array@k@i - sum) / diag } } } return [new Matrix[L], new Matrix[U]] } /** Returns the determinant of a square matrix. */ det[] := { if ! isSquare[] { println["Determinants are only defined for square matrices! Size was $rows x $cols"] return undef } // It's much faster to use precalculated determinant formulas for small // matrices than, to, say LUDecompose the things. if rows == 1 return array@0@0 if rows == 2 return array@0@0 array@1@1 - array@0@1 array@1@0 // ad - bc if rows == 3 return -array@0@2 array@1@1 array@2@0 + array@0@1 array@1@2 array@2@0 + array@0@2 array@1@0 array@2@1 - array@0@0 array@1@2 array@2@1 - array@0@1 array@1@0 array@2@2 + array@0@0 array@1@1 array@2@2 if rows == 4 return array@0@1 array@1@3 array@2@2 array@3@0 - array@0@1 array@1@2 array@2@3 array@3@0 - array@0@0 array@1@3 array@2@2 array@3@1 + array@0@0 array@1@2 array@2@3 array@3@1 - array@0@1 array@1@3 array@2@0 array@3@2 + array@0@0 array@1@3 array@2@1 array@3@2 + array@0@1 array@1@0 array@2@3 array@3@2 - array@0@0 array@1@1 array@2@3 array@3@2 + array@0@3 (array@1@2 array@2@1 array@3@0 - array@1@1 array@2@2 array@3@0 - array@1@2 array@2@0 array@3@1 + array@1@0 array@2@2 array@3@1 + array@1@1 array@2@0 array@3@2 - array@1@0 array@2@1 array@3@2) + (array@0@1 array@1@2 array@2@0 - array@0@0 array@1@2 array@2@1 - array@0@1 array@1@0 array@2@2 + array@0@0 array@1@1 array@2@2) array@3@3 + array@0@2 (-array@1@3 array@2@1 array@3@0 + array@1@1 array@2@3 array@3@0 + array@1@3 array@2@0 array@3@1 - array@1@0 array@2@3 array@3@1 - array@1@1 array@2@0 array@3@3 + array@1@0 array@2@1 array@3@3) if rows == 5 return array@0@2 array@1@4 array@2@3 array@3@1 array@4@0 - array@0@2 array@1@3 array@2@4 array@3@1 array@4@0 - array@0@1 array@1@4 array@2@3 array@3@2 array@4@0 + array@0@1 array@1@3 array@2@4 array@3@2 array@4@0 - array@0@2 array@1@4 array@2@1 array@3@3 array@4@0 + array@0@1 array@1@4 array@2@2 array@3@3 array@4@0 + array@0@2 array@1@1 array@2@4 array@3@3 array@4@0 - array@0@1 array@1@2 array@2@4 array@3@3 array@4@0 + array@0@2 array@1@3 array@2@1 array@3@4 array@4@0 - array@0@1 array@1@3 array@2@2 array@3@4 array@4@0 - array@0@2 array@1@1 array@2@3 array@3@4 array@4@0 + array@0@1 array@1@2 array@2@3 array@3@4 array@4@0 - array@0@2 array@1@4 array@2@3 array@3@0 array@4@1 + array@0@2 array@1@3 array@2@4 array@3@0 array@4@1 + array@0@0 array@1@4 array@2@3 array@3@2 array@4@1 - array@0@0 array@1@3 array@2@4 array@3@2 array@4@1 + array@0@2 array@1@4 array@2@0 array@3@3 array@4@1 - array@0@0 array@1@4 array@2@2 array@3@3 array@4@1 - array@0@2 array@1@0 array@2@4 array@3@3 array@4@1 + array@0@0 array@1@2 array@2@4 array@3@3 array@4@1 - array@0@2 array@1@3 array@2@0 array@3@4 array@4@1 + array@0@0 array@1@3 array@2@2 array@3@4 array@4@1 + array@0@2 array@1@0 array@2@3 array@3@4 array@4@1 - array@0@0 array@1@2 array@2@3 array@3@4 array@4@1 + array@0@1 array@1@4 array@2@3 array@3@0 array@4@2 - array@0@1 array@1@3 array@2@4 array@3@0 array@4@2 - array@0@0 array@1@4 array@2@3 array@3@1 array@4@2 + array@0@0 array@1@3 array@2@4 array@3@1 array@4@2 - array@0@1 array@1@4 array@2@0 array@3@3 array@4@2 + array@0@0 array@1@4 array@2@1 array@3@3 array@4@2 + array@0@1 array@1@0 array@2@4 array@3@3 array@4@2 - array@0@0 array@1@1 array@2@4 array@3@3 array@4@2 + array@0@1 array@1@3 array@2@0 array@3@4 array@4@2 - array@0@0 array@1@3 array@2@1 array@3@4 array@4@2 - array@0@1 array@1@0 array@2@3 array@3@4 array@4@2 + array@0@0 array@1@1 array@2@3 array@3@4 array@4@2 + array@0@2 array@1@4 array@2@1 array@3@0 array@4@3 - array@0@1 array@1@4 array@2@2 array@3@0 array@4@3 - array@0@2 array@1@1 array@2@4 array@3@0 array@4@3 + array@0@1 array@1@2 array@2@4 array@3@0 array@4@3 - array@0@2 array@1@4 array@2@0 array@3@1 array@4@3 + array@0@0 array@1@4 array@2@2 array@3@1 array@4@3 + array@0@2 array@1@0 array@2@4 array@3@1 array@4@3 - array@0@0 array@1@2 array@2@4 array@3@1 array@4@3 + array@0@1 array@1@4 array@2@0 array@3@2 array@4@3 - array@0@0 array@1@4 array@2@1 array@3@2 array@4@3 - array@0@1 array@1@0 array@2@4 array@3@2 array@4@3 + array@0@0 array@1@1 array@2@4 array@3@2 array@4@3 + array@0@2 array@1@1 array@2@0 array@3@4 array@4@3 - array@0@1 array@1@2 array@2@0 array@3@4 array@4@3 - array@0@2 array@1@0 array@2@1 array@3@4 array@4@3 + array@0@0 array@1@2 array@2@1 array@3@4 array@4@3 + array@0@1 array@1@0 array@2@2 array@3@4 array@4@3 - array@0@0 array@1@1 array@2@2 array@3@4 array@4@3 + array@0@4 (array@1@1 array@2@3 array@3@2 array@4@0 - array@1@1 array@2@2 array@3@3 array@4@0 - array@1@0 array@2@3 array@3@2 array@4@1 + array@1@0 array@2@2 array@3@3 array@4@1 - array@1@1 array@2@3 array@3@0 array@4@2 + array@1@0 array@2@3 array@3@1 array@4@2 + array@1@1 array@2@0 array@3@3 array@4@2 - array@1@0 array@2@1 array@3@3 array@4@2 + array@1@3 (array@2@2 array@3@1 array@4@0 - array@2@1 array@3@2 array@4@0 - array@2@2 array@3@0 array@4@1 + array@2@0 array@3@2 array@4@1 + array@2@1 array@3@0 array@4@2 - array@2@0 array@3@1 array@4@2) + (array@1@1 array@2@2 array@3@0 - array@1@0 array@2@2 array@3@1 - array@1@1 array@2@0 array@3@2 + array@1@0 array@2@1 array@3@2) array@4@3 + array@1@2 (-array@2@3 array@3@1 array@4@0 + array@2@1 array@3@3 array@4@0 + array@2@3 array@3@0 array@4@1 - array@2@0 array@3@3 array@4@1 - array@2@1 array@3@0 array@4@3 + array@2@0 array@3@1 array@4@3)) + (array@0@2 (-array@1@3 array@2@1 array@3@0 + array@1@1 array@2@3 array@3@0 + array@1@3 array@2@0 array@3@1 - array@1@0 array@2@3 array@3@1 - array@1@1 array@2@0 array@3@3 + array@1@0 array@2@1 array@3@3) + array@0@1 (array@1@3 array@2@2 array@3@0 - array@1@2 array@2@3 array@3@0 - array@1@3 array@2@0 array@3@2 + array@1@0 array@2@3 array@3@2 + array@1@2 array@2@0 array@3@3 - array@1@0 array@2@2 array@3@3) + array@0@0 (-array@1@3 array@2@2 array@3@1 + array@1@2 array@2@3 array@3@1 + array@1@3 array@2@1 array@3@2 - array@1@1 array@2@3 array@3@2 - array@1@2 array@2@1 array@3@3 + array@1@1 array@2@2 array@3@3)) array@4@4 + array@0@3 (-array@1@1 array@2@4 array@3@2 array@4@0 + array@1@1 array@2@2 array@3@4 array@4@0 + array@1@0 array@2@4 array@3@2 array@4@1 - array@1@0 array@2@2 array@3@4 array@4@1 + array@1@1 array@2@4 array@3@0 array@4@2 - array@1@0 array@2@4 array@3@1 array@4@2 - array@1@1 array@2@0 array@3@4 array@4@2 + array@1@0 array@2@1 array@3@4 array@4@2 + array@1@4 (-array@2@2 array@3@1 array@4@0 + array@2@1 array@3@2 array@4@0 + array@2@2 array@3@0 array@4@1 - array@2@0 array@3@2 array@4@1 - array@2@1 array@3@0 array@4@2 + array@2@0 array@3@1 array@4@2) - array@1@1 array@2@2 array@3@0 array@4@4 + array@1@0 array@2@2 array@3@1 array@4@4 + array@1@1 array@2@0 array@3@2 array@4@4 - array@1@0 array@2@1 array@3@2 array@4@4 + array@1@2 (array@2@4 array@3@1 array@4@0 - array@2@1 array@3@4 array@4@0 - array@2@4 array@3@0 array@4@1 + array@2@0 array@3@4 array@4@1 + array@2@1 array@3@0 array@4@4 - array@2@0 array@3@1 array@4@4)) // If we fell through to here, we have a larger matrix and have to try // to find its determinant through more inefficent methods. // TODO: Calculate determinant in terms of other determinants instead of // using LUDecompose when LUDecompose doesn't work. product = 1 [L, U] = LUDecomposeDoolittle[] // This will happen if the matrix is singular. if U == undef return undef // Multiply diagonals of the lower triangular matrix. The // upper triangular matrix has all ones on the diagonal. // Should there be a permutation matrix here to get the signs // right? for i=0 to rows-1 product = product * (U.array)@i@i return product } /** Create a square identity matrix with the specified number of rows and columns. The elements on the diagonal will be set to 1, the rest to zero. This requires Frink 2020-04-19 or later. */ class makeIdentity[dimension] := { return new Matrix[makeArray[[dimension, dimension], {|a,b| a==b ? 1 : 0}]] } /** Makes a diagonal matrix. This is passed in an array of elements (e.g. [1,2,3] that will make up the diagonal elements. This requires Frink 2020-04-22 or later. */ class makeDiagonal[arr] := { d = length[arr] return new Matrix[makeArray[[d,d], {|a,b,data| a==b ? data@a : 0}, arr]] } /** Makes a new one-column Matrix from an array. */ class makeColumn[arr] := { return new Matrix[arr.transpose[]] } /** Returns a matrix with the specified row and column removed. Row and column numbers are 1-indexed. This is used in many matrix operations including calculation of determinant or of the adjugate/adjoint matrix. */ removeRowColumn[rowToRemove, colToRemove] := { a = makeArray[[rows-1, cols-1]] newRow = 0 ROW: for origRow = 0 to rows-1 { if origRow == rowToRemove-1 next ROW newCol = 0 COL: for origCol = 0 to rows-1 { if origCol == colToRemove-1 next COL a@newRow@newCol = array@origRow@origCol newCol = newCol + 1 } newRow = newRow + 1 } return new Matrix[a] } /** Returns the adjugate or adjoint matrix. See: https://semath.info/src/inverse-cofactor-ex4.html TODO: Calculate hardcoded adjugate matrices for small matrices because this is expensive */ adjugate[] := { a = makeArray[[rows,cols]] for i = 0 to rows-1 for j = 0 to cols-1 a@i@j = (-1)^((i+1)+(j+1)) * removeRowColumn[j+1, i+1].det[] return new Matrix[a] } /** Returns the inverse of a matrix. See: https://semath.info/src/inverse-cofactor-ex4.html TODO: Calculate hardcoded inverse matrices for small matrices because this is expensive */ inverse[] := { if ! isSquare[] { println["Matrix.inverse only works on square arrays."] return undef } // It's much faster to calculate inverse matrices with hard-coded // equations for small matrices. if rows == 1 return 1 / array@0@0 if rows == 2 { invdet = 1/det[] // [ d -b ] // [ -c a ] / det return (new Matrix[[[array@1@1, -array@0@1], [-array@1@0, array@0@0]]]).multiplyByScalar[invdet] } if rows == 3 { invdet = 1/det[] return (new Matrix[[[-array@1@2 array@2@1 + array@1@1 array@2@2, array@0@2 array@2@1 - array@0@1 array@2@2, -array@0@2 array@1@1 + array@0@1 array@1@2], [ array@1@2 array@2@0 - array@1@0 array@2@2, -array@0@2 array@2@0 + array@0@0 array@2@2, array@0@2 array@1@0 - array@0@0 array@1@2], [-array@1@1 array@2@0 + array@1@0 array@2@1, array@0@1 array@2@0 - array@0@0 array@2@1, -array@0@1 array@1@0 + array@0@0 array@1@1]]]).multiplyByScalar[invdet] } if rows == 4 { invdet = 1/det[] return (new Matrix[[[-array@1@3 array@2@2 array@3@1 + array@1@2 array@2@3 array@3@1 + array@1@3 array@2@1 array@3@2 - array@1@1 array@2@3 array@3@2 - array@1@2 array@2@1 array@3@3 + array@1@1 array@2@2 array@3@3, array@0@3 array@2@2 array@3@1 - array@0@2 array@2@3 array@3@1 - array@0@3 array@2@1 array@3@2 + array@0@1 array@2@3 array@3@2 + array@0@2 array@2@1 array@3@3 - array@0@1 array@2@2 array@3@3, -array@0@3 array@1@2 array@3@1 + array@0@2 array@1@3 array@3@1 + array@0@3 array@1@1 array@3@2 - array@0@1 array@1@3 array@3@2 - array@0@2 array@1@1 array@3@3 + array@0@1 array@1@2 array@3@3, array@0@3 array@1@2 array@2@1 - array@0@2 array@1@3 array@2@1 - array@0@3 array@1@1 array@2@2 + array@0@1 array@1@3 array@2@2 + array@0@2 array@1@1 array@2@3 - array@0@1 array@1@2 array@2@3], [array@1@3 array@2@2 array@3@0 - array@1@2 array@2@3 array@3@0 - array@1@3 array@2@0 array@3@2 + array@1@0 array@2@3 array@3@2 + array@1@2 array@2@0 array@3@3 - array@1@0 array@2@2 array@3@3, -array@0@3 array@2@2 array@3@0 + array@0@2 array@2@3 array@3@0 + array@0@3 array@2@0 array@3@2 - array@0@0 array@2@3 array@3@2 - array@0@2 array@2@0 array@3@3 + array@0@0 array@2@2 array@3@3, array@0@3 array@1@2 array@3@0 - array@0@2 array@1@3 array@3@0 - array@0@3 array@1@0 array@3@2 + array@0@0 array@1@3 array@3@2 + array@0@2 array@1@0 array@3@3 - array@0@0 array@1@2 array@3@3, -array@0@3 array@1@2 array@2@0 + array@0@2 array@1@3 array@2@0 + array@0@3 array@1@0 array@2@2 - array@0@0 array@1@3 array@2@2 - array@0@2 array@1@0 array@2@3 + array@0@0 array@1@2 array@2@3], [-array@1@3 array@2@1 array@3@0 + array@1@1 array@2@3 array@3@0 + array@1@3 array@2@0 array@3@1 - array@1@0 array@2@3 array@3@1 - array@1@1 array@2@0 array@3@3 + array@1@0 array@2@1 array@3@3, array@0@3 array@2@1 array@3@0 - array@0@1 array@2@3 array@3@0 - array@0@3 array@2@0 array@3@1 + array@0@0 array@2@3 array@3@1 + array@0@1 array@2@0 array@3@3 - array@0@0 array@2@1 array@3@3, -array@0@3 array@1@1 array@3@0 + array@0@1 array@1@3 array@3@0 + array@0@3 array@1@0 array@3@1 - array@0@0 array@1@3 array@3@1 - array@0@1 array@1@0 array@3@3 + array@0@0 array@1@1 array@3@3, array@0@3 array@1@1 array@2@0 - array@0@1 array@1@3 array@2@0 - array@0@3 array@1@0 array@2@1 + array@0@0 array@1@3 array@2@1 + array@0@1 array@1@0 array@2@3 - array@0@0 array@1@1 array@2@3], [array@1@2 array@2@1 array@3@0 - array@1@1 array@2@2 array@3@0 - array@1@2 array@2@0 array@3@1 + array@1@0 array@2@2 array@3@1 + array@1@1 array@2@0 array@3@2 - array@1@0 array@2@1 array@3@2, -array@0@2 array@2@1 array@3@0 + array@0@1 array@2@2 array@3@0 + array@0@2 array@2@0 array@3@1 - array@0@0 array@2@2 array@3@1 - array@0@1 array@2@0 array@3@2 + array@0@0 array@2@1 array@3@2, array@0@2 array@1@1 array@3@0 - array@0@1 array@1@2 array@3@0 - array@0@2 array@1@0 array@3@1 + array@0@0 array@1@2 array@3@1 + array@0@1 array@1@0 array@3@2 - array@0@0 array@1@1 array@3@2, -array@0@2 array@1@1 array@2@0 + array@0@1 array@1@2 array@2@0 + array@0@2 array@1@0 array@2@1 - array@0@0 array@1@2 array@2@1 - array@0@1 array@1@0 array@2@2 + array@0@0 array@1@1 array@2@2]]]).multiplyByScalar[invdet] } // Otherwise use the adjugate function. return adjugate[].multiplyByScalar[1/det[]] } /** Returns the Kronecker product of a and b. https://en.wikipedia.org/wiki/Kronecker_product */ class KroneckerProduct[a is Matrix, b is Matrix] := { m = a.rows n = a.cols p = b.rows q = b.cols newRows = m p newCols = n q n = new array[[newRows, newCols], 0] for i=0 to newRows-1 for j=0 to newCols-1 n@i@j = (a.array)@(i div p)@(j div q) * (b.array)@(i mod p)@(j mod q) return new Matrix[n] } /** Return the Kronecker product of this matrix and b. */ KroneckerProduct[b is Matrix] := KroneckerProduct[this, b] /** Returns true if both matrices are equal (that is, they have the same dimensions and all array elements are equal.) */ equals[other is Matrix] := { return rows==other.rows and cols==other.cols and array==other.array } /** Returns a new Matrix where each element is the complex conjugate of the corresponding element in original Matrix. */ conjugate[] := { c = makeArray[[rows, cols], {|r,c,data| conjugate[data@r@c]}, array] return new Matrix[c] } /** Returns a new Matrix where the Matrix is first transposed and then each element is the complex conjugate of the corresponding element in the transposed Matrix. */ conjugateTranspose[] := { t = array.transpose[] c = makeArray[[rows, cols], {|r,c,data| conjugate[data@r@c]}, t] return new Matrix[c] } /** Returns true if this matrix is "Hermitian", or "self-adjoint", that is, the matrix must be square and equal to its own conjugate transpose. In other words, for all elements, array@i@j == conjugate[array@j@i] where conjugate is the complex conjugate of a number. */ isHermitian[] := { if rows != cols return false for i=0 to rows-1 for j = i+1 to rows-1 if array@i@j != conjugate[array@j@i] return false return true } /** Performs the QR decomposition of a matrix. This is based on the (real-only) implementation at: https://github.com/fiji/Jama/blob/master/src/main/java/Jama/QRDecomposition.java This is an internal implementation, primarily for use by the leastSquares method, that returns raw arrays [QR, rdiag]. If you're solving for least squares, use that routine directly instead of this. */ QRDecomposeInternal[] := { QR = deepCopy[array] // Copy the original matrix's array rdiag = new array[[cols], 0] // Stores the diagonal for k = 0 to cols-1 { // Compute 2-norm of k-th column nrm = 0 QR@0@k // Make units work right for i = k to rows-1 nrm = sqrt[nrm^2 + (QR@i@k)^2] if nrm != 0 QR@0@k // Make units work right { // Form k-th Householder vector if isNegativeUnit[QR@k@k] nrm = -nrm for i = k to rows-1 QR@i@k = QR@i@k / nrm QR@k@k = QR@k@k + 1 // Apply transformation to remaining columns for j = k+1 to cols-1 { s = undef for i=k to rows-1 { if s == undef s = QR@i@k * QR@i@j else s = s + QR@i@k * QR@i@j } s = -s / QR@k@k for i = k to rows-1 QR@i@j = QR@i@j + s * QR@i@k } } rdiag@k = -nrm } return [QR, rdiag] } /** This performs a QR decomposition of this matrix and returns [Q, R] as new matrices. Q is the orthogonal factor. R is the upper triangular factor. If you are doing a least-squares solve, call leastSquares directly instead. THINK ABOUT: Should this return the Householder vectors and rdiag? */ QRDecompose[] := { // The Q and R matrices are packed into the results of the following. [QR, rdiag] = QRDecomposeInternal[] // Extract the Q matrix Q = new array[[rows, cols], 0] for k = cols-1 to 0 step -1 { for i = 0 to rows-1 Q@i@k = 0 Q@k@k = 1 for j = k to cols-1 { if QR@k@k != 0 QR@k@k // Make units work right { s = 0 (QR@0@k) * Q@0@j // Make units work right for i = k to rows-1 s = s + QR@i@k * Q@i@j s = -s / QR@k@k for i = k to rows-1 Q@i@j = Q@i@j + s * QR@i@k } } } // Extract the R matrix R = new array[[cols, cols], 0] for i = 0 to cols-1 { for j = 0 to cols-1 { if i < j R@i@j = QR@i@j else if i == j R@i@j = rdiag@i else R@i@j = 0 } } return [new Matrix[Q], new Matrix[R]] } /** Return the least-squares solution (called X) of the system A * X = B where this Matrix is A and parameter B is B. This can be performed on an overdetermined system, where there are more measurements than equations. TODO: Try to solve exactly using x = A.inverse[] * b when the system is not overdetermined? See MatrixQRTest.frink for examples. This is the best discussion I've seen of least-squares fitting: https://www.aleksandrhovhannisyan.com/blog/the-method-of-least-squares/ https://www.aleksandrhovhannisyan.com/blog/least-squares-fitting/ */ leastSquares[B is Matrix] := { if rows != B.rows { println["Matrix.leastSquares: Matrix row dimensions must agree."] return undef } [QR, rdiag] = QRDecomposeInternal[] // Check if the matrix is "full rank," that is, that each row is // independent (not a multiple of) other rows. for j = 0 to cols-1 { if rdiag@j == 0 rdiag@j { println["Matrix.leastSquares: Row entries are not all independent."] return undef } } X = deepCopy[B.array] // Copy B's array // Compute Y = Q.transpose * B for k = 0 to cols-1 { for j = 0 to B.cols-1 { s = undef for i = k to rows-1 if s == undef s = QR@i@k * X@i@j else s = s + QR@i@k * X@i@j s = -s / QR@k@k for i = k to rows-1 X@i@j = X@i@j + s * QR@i@k } } // Solve R * X = y for k = cols-1 to 0 step -1 { for j = 0 to B.cols-1 X@k@j = X@k@j / rdiag@k for i = 0 to k-1 for j = 0 to B.cols-1 X@i@j = X@i@j - X@k@j * QR@i@k } X1 = new Matrix[X] return X1.getSubMatrix[0, cols-1, 0, B.cols-1] } /** Gets a submatrix of this matrix. */ getSubMatrix[fromRow, toRow, fromCol, toCol] := { a = new array[[toRow-fromRow+1, toCol-fromCol+1], 0] for row = fromRow to toRow for col = fromCol to toCol a@(row-fromRow)@(col-fromCol) = array@row@col // println["a is $a"] return new Matrix[a] } /** Rounds values to the nearest integer if they are less than the specified relative error away from an integer and returns a new array. */ roundToInt[relerror = 1e-14] := { ret = new array[[rows, cols]] for r = 0 to rows-1 for c = 0 to cols-1 { v = array@r@c rnd = round[v] if (rnd==0 and abs[rnd-v] <= relerror) or abs[(rnd-v)/rnd] <= relerror ret@r@c = rnd else ret@r@c = v } return new Matrix[ret] } /** Adds the multiple of row1 to row2, modifying row2 in place. */ addMultipleOfRow[row1, mulBy, row2] := { r1 = row1-1 r2 = row2-1 array@r2 = add[array@r2, mul[array@r1, mulBy]] } /** Static method to count the number of leading zeroes in a row (represented as a one-dimensional array.) */ class countLeadingZeroes[a] := { len = length[a] count = 0 while count < len and a@count == 0 count = count + 1 return count } /** Sorts rows by leading zeros, with elements with more leading zeroes at the bottom. This modifies the matrix in place.*/ sortByLeadingZeroes[] := { sort[array, {|a,b| Matrix.countLeadingZeroes[a] <=> Matrix.countLeadingZeroes[b]}] } /** Solves a matrix equation (this * x = B) for x. The result is a one-dimensional array of values for x. This matrix should have the same number of rows as matrix B. This is a dispatcher method that calls the appropriate method. This even works with units of measure. THINK ABOUT: What should this return for an inconsistent system? THINK ABOUT: If a row is all zeroes, then don't return that row? */ solve[B] := { if ! B.hasUnits[] return solveDimensionless[B] else return solveWithUnits[B] } /** Solves a matrix equation (this * x = B) for x. This uses Gauss-Jordan reduction to solve the equation. The result is a column Matrix of values for x. */ solveToMatrix[B] := { tr = solve[B].transpose[] return new Matrix[noEval[tr]] } /** Primary internal matrix solver. Solves a matrix equation (this * x = B) for x. This uses Gauss-Jordan reduction to solve the equation. The result is a one-dimensional array of values for x. This matrix should have the same number of rows as matrix B. THINK ABOUT: What should this return for an inconsistent system? THINK ABOUT: If a row is all zeroes, then don't return that row? */ solveDimensionless[B] := { augment = this.augment[B] // Make augmented matrix augment.reduceRows[] // Do the actual solve augment.simplifySymbolic[] // Clean up symbolic solutions // Check for inconsistent system. if augment.getPivotRow[augment.cols] == augment.rows println["Matrix.solve: Warning: inconsistent system."] col = augment.getColumnAsArray[augment.getCols[]] // Extract solution col return col } /** Specialized internal solver to work with units of measure. This factors out units of measure and then derives the required units of measure and multiplies them back in. The result is a one-dimensional array of values for x. */ solveWithUnits[B] := { [A1, A2] = factorUnits[] [B1, B2] = B.factorUnits[] x = A1.solveDimensionless[B1] // Derive necessary units of measure and multiply them back in len = length[x] for i = 0 to len-1 x@i = x@i * (B2.array@i@0 / A2.array@i@0) return x } /** This turns a Matrix into reduced row-echelon form, modifying this matrix in place. This uses Gauss-Jordan reduction. This can be used to solve a system of equations in "augmented matrix" form. See the solve function. See: https://www.cfm.brown.edu/people/dobrush/cs52/Mathematica/Part7/Part1/LinearEquations_rref.html https://chrisj.math.gatech.edu/24s/1553/materials/1_2-1_3-full.pdf */ reduceRows[] := { r = 0 for j = 1 to cols { i = r+1 while i <= rows and get[i,j] == 0 i = i + 1 if i < rows+1 { r = r + 1 exchangeRows[i,r] scaleToOne[r,j] ROW: for k = 1 to rows { if k != r and get[k,j] != 0 makeZero[k, j, r] } } } } /** Scales a cell to 1 by multiplying the specified row by the reciprocal of the cell. This modifies the matrix in place. */ scaleToOne[row, col] := { n = get[row, col] if n == 0 { println["Cannot scale [$row, $col] to 1 because it's zero"] return } if n != 1 multiplyRowByScalar[row, recip[n]] // Force the cell to exactly 1 in case there's some floating-point // or symbolic issue set[row, col, 1] } /** Scales a cell to 0 by adding another scaled row to its row. This modifies the matrix in place. */ makeZero[row, col, otherRow] := { n = get[row,col] other = get[otherRow, col] if n != 0 // Should check if other != 0 ? addMultipleOfRow[otherRow, negate[n/other], row] // Force the cell to exactly 0 in case there's some floating-point // or symbolic issue set[row, col, 0] } /** Augments a matrix with another column matrix, returning a new matrix. That is, the column matrix is appended to the right of this matrix. The matrices should have the same number of rows. */ augment[colMatrix is Matrix] := { a1 = deepCopy[array] for c = 0 to colMatrix.cols - 1 a1.setColumn[colMatrix.array.getColumn[c],cols+c] return new Matrix[a1] } /** Given a Matrix in reduced row-echelon form (i.e. one that comes from calling reduceRows), returns the index of the pivot row for the specified column (1-based). The to have a pivot in the column, the value in that row has to be 1 and all values to the left of that 1 must be zero. If the column does not have a pivot, this returns undef. When solving a system of equations, if a column does not have a pivot row, it is a free variable. An augmented matrix corresponds to an inconsistent sytem of equations if and only if the last (i.e. the augmented) column is a pivot column. */ getPivotRow[col] := { for r = rows to 1 step -1 if get[r, col] == 1 { for c = 1 to col-1 if get[r,c] != 0 return undef return r } return undef } /** When solving a symbolic matrix with ReduceRows, the reduced matrix may contain terms like 1 + 0 r^-1 s which has unreduced multiplicative terms. This eliminates those terms, modifying the matrix in place. */ simplifySymbolic[] := { // Check to see if array has unresolved symbols. This is a relatively // quick check. if length[getSymbols[array]] > 0 { for r = 1 to rows { for c = 1 to cols { before = get[r,c] after = substituteExpression[before, noEval[0 _x], 0] after = substituteExpression[after, noEval[0 + _x], noEval[_x]] if after != before { // println["$before simplified to $after"] set[r, c, after] } } } } } /** Calculate the trace of a matrix. The trace is the sum of the entries on the diagonal. */ trace[] := { if ! isSquare[] { println["Matrix.trace: matrix is not square!"] return undef } sum = array@0@0 // Make dimensions come out right for i = 1 to rows-1 sum = sum + array@i@i return sum } /** Returns true if any element has units of measure, that is, if any element is not dimensionless. */ hasUnits[] := { for r = 0 to rows-1 { row = array@r for c = 0 to cols-1 if ! ((row@c) conforms 0) return true } return false } /** Factors a matrix into two matrices [A, B] where A contains only dimensionless numbers and B contains dimensions. This is basically calling factorUnits on each element. returns [A, B] */ factorUnits[] := { A = new Matrix[rows, cols] B = new Matrix[rows, cols] for r = 1 to rows for c = 1 to cols { v = get[r,c] if isUnit[v] [a, b] = factorUnits[get[r,c]] else [a, b] = [v, 1] A.set[r,c,a] B.set[r,c,b] } return [A, B] } } "class Matrix loaded OK."