From 77b2bc62c77ec9101c5a7db441eae6165477a327 Mon Sep 17 00:00:00 2001 From: Chris Date: Tue, 9 Jul 2019 00:56:15 -0700 Subject: [PATCH] Matrix updates --- src/apatite/linear_algebra.cr | 2 +- src/apatite/linear_algebra/matrix.cr | 755 +++++++++------------------ 2 files changed, 245 insertions(+), 512 deletions(-) diff --git a/src/apatite/linear_algebra.cr b/src/apatite/linear_algebra.cr index 36b242a..42c553a 100644 --- a/src/apatite/linear_algebra.cr +++ b/src/apatite/linear_algebra.cr @@ -1,6 +1,6 @@ require "./linear_algebra/error" # require "./linear_algebra/ndarray" -# require "./linear_algebra/matrix" +require "./linear_algebra/matrix" require "./linear_algebra/vector" module Apatite diff --git a/src/apatite/linear_algebra/matrix.cr b/src/apatite/linear_algebra/matrix.cr index 3856af5..aed0114 100644 --- a/src/apatite/linear_algebra/matrix.cr +++ b/src/apatite/linear_algebra/matrix.cr @@ -2,569 +2,302 @@ require "json" require "./vector" module Apatite::LinearAlgebra - class Matrix + class Matrix(T) include Enumerable(Vector) include Indexable(Vector) include Comparable(Matrix) - getter column_count : Int32 - - getter row_count : Int32 - - @buffer : Pointer(Vector) - - def initialize(rows, column_count : Int32? = nil) - @buffer = rows.map { |r| Vector.create(r) }.to_a.to_unsafe - @row_count = rows.size - @column_count = column_count || rows[0].size - end - - # Creates a new `Vector` instance from a `JSON::PullParser` - def self.new(pull : JSON::PullParser) - arr = [] of Vector - new(pull) do |element| - arr << element - end - rows(arr) - end - - def self.new(pull : JSON::PullParser, &block) - pull.read_array do - yield Vector.new(pull) - end - end + protected getter rows : Array(Array(T)) # Creates a matrix where each argument is a row. + # + # ``` + # Matrix[[25, 93], [-1, 66]] + # # => [ 25, 93, + # # -1, 66 ] + # ``` def self.[](*rows) - rows(rows) + rows(rows, false) end - # Creates a matrix of size `row_count x column_count`. It fills the values by calling - # the given block, passing the current row and column. + # Creates a matrix where +rows+ is an array of arrays, each of which is a row + # of the matrix. If the optional argument +copy+ is false, use the given + # arrays as the internal structure of the matrix without copying. + # + # ``` + # Matrix.rows([[25, 93], [-1, 66]]) + # # => [ 25, 93, + # # -1, 66 ] + # ``` + def self.rows(rows : Indexable(Array(T)), copy = true) + rows = rows.dup if copy + rows = rows.to_a + rows.map! do |row| + row = row.dup if row + row.to_a + end + size = (rows[0] || [] of T).size + rows.each do |row| + raise ErrDimensionMismatch.new("row size differs (#{row.size} should be #{size})") unless row.size == size + end + new(rows, size) + end + + # Creates a matrix using +columns+ as an array of column vectors. + # + # ``` + # Matrix.columns([[25, 93], [-1, 66]]) + # # => [ 25, -1, + # # 93, 66 ] + # ``` + def self.columns(columns) + rows(columns, false).transpose + end + + # Creates a matrix of size +row_count+ x +column_count+. + # It fills the values by calling the given block, + # passing the current row and column. + # Returns an enumerator if no block is given. + # + # ``` + # m = Matrix.build(2, 4) { |row, col| col - row } + # # => Matrix[[0, 1, 2, 3], [-1, 0, 1, 2]] + # m = Matrix.build(3) { rand } + # # => a 3x3 matrix with random elements + # ``` def self.build(row_count, column_count = row_count, &block) + row_count = row_count.to_i + column_count = column_count.to_i raise ArgumentError.new if row_count < 0 || column_count < 0 - rows = Array.new(row_count) do |i| - Vector.new(column_count) do |j| + rows = Array(T).new(row_count) do |i| + Array(T).new(column_count) do |j| yield i, j end end - Matrix.new(rows, column_count) + new(rows, column_count) end - # Creates a single-column matrix where the values of that column are as given in column. - def self.col_vector(column) - Matrix.new([column].transpose, 1) - end - - # Creates a matrix using `columns` as an array of column vectors. - def self.columns(columns) - rows(columns).transpose - end - - # # Create a matrix by combining matrices entrywise, using the given block. - # def self.combine(matrices : Array(Matrix), &block : Matrix -> Matrix -> Matrix) - # return Matrix.empty if matrices.empty? - # x = matrices.first - # matrices.each do |m| - # raise "Dimension mismatch" unless x.row_count == m.row_count && x.column_count == m.column_count - # end - - # rows = Array.new(x.row_count) do |i| - # Vector.new(x.column_count) do |j| - # yield matrices.map { |m| m[i, j] } - # end - # end - - # Matrix.new(rows, x.column_count) - # end - - # def self.combine(*matrices, &block : Matrix -> Matrix) - # Matrix.combine(matrices, &block) - # end - # Creates a matrix where the diagonal elements are composed of `values`. - def self.diagonal(values) + # + # ``` + # Matrix.diagonal(9, 5, -3) + # # => [ 9, 0, 0, + # # 0, 5, 0, + # # 0, 0, -3 ] + # ``` + def self.diagonal(values : Indexable(T), dummy = nil) size = values.size - return Matrix.empty if size == 0 - - rows = Array.new(size) do |j| - row = Array.new(size, 0) + return Matrix(T).empty if size == 0 + rows = Array(Array(T)).new(size) do |j| + row = Array(T).new(size, T.new(0)) row[j] = values[j] row end - - new rows + new(rows) end - # :ditto: - def self.diagonal(*values) - Matrix.diagonal(values) + def self.diagonal(*values : T) + diagonal(values, nil) end - # Creates a empty matrix of `row_count x column_count`. At least one of - # `row_count` or `column_count` must be 0. - def self.empty(row_count = 0, column_count = 0) - raise ArgumentError.new("One size must be 0") if column_count != 0 && row_count != 0 - raise ArgumentError.new("Negative size") if column_count < 0 || row_count < 0 - Matrix.new(([] of Vector) * row_count, column_count) + # Creates an +n+ by +n+ diagonal matrix where each diagonal element is + # `value`. + # + # ``` + # Matrix.scalar(2, 5) + # # => [ 5, 0, + # # 0, 5 ] + # ``` + def self.scalar(n, value : T) + diagonal(Array(T).new(n, value)) end - # Creates a new diagonal matrix of size `n` with ones in the diagonal - # and zeros elsewhere. - def self.eye(n) - Matrix.diagonal([1] * n) + # Creates an `n` by `n` identity matrix. + # + # ``` + # Matrix.identity(2) + # # => [ 1, 0, + # # 0, 1 ] + # ``` + def self.identity(n : T) + scalar(n, T.new(1)) end - # TODO - def self.hstack(x, *matrices) - end - - # Creates a `n x n` identity matrix. - def self.identity(n) - scalar(n, 1) - end - - # Creates a matrix of the given shape with random vectors. - def self.random(num_rows, num_columns, range = nil) - Matrix.build(num_rows, num_columns) do |i, j| - rand(range || -1e+1..1e+1) - end - end - - # Creates a single-row matrix where the values of that row are as given in `row`. - def self.row_vector(row) - Matrix.new([row], 0) - end - - # Creates a matrix where rows is an array of arrays, each of which is a row of the matrix. - def self.rows(rows) - size = rows[0]? ? rows[0].size : 0 - rows.each do |row| - raise "Dimension mismatch: row size differs (#{row.size} should be #{size})" unless row.size == size - end - Matrix.new(rows, size) - end - - # Creates an `n` by `n` diagonal matrix where each diagonal element is value. - def self.scalar(n, value) - Matrix.diagonal(Array.new(n, value)) - end - - # TODO - def self.vstack(x, y) + # ditto + def self.unit(n : T) + identity(n) end # Creates a zero matrix. - def self.zero(row_count, column_count = row_count) - rows = Array.new(row_count) { Vector.new(column_count) } - Matrix.new(rows, column_count) - end - - def *(other : Matrix) - raise "Dimension mismatch" if column_count != other.column_count - - rows = Array.new(row_count) do |i| - Vector.new(other.column_count) do |j| - (0...column_count).reduce(0.0) do |vij, k| - vij + self[i, k] * other[k, j] - end - end - end - - return Matrix.new(rows, other.column_count) - end - - def *(int : Int) - rows = self.rows.map do |row| - row.map { |e| e * int } - end - Matrix.new(rows, column_count) - end - - def *(ind : Indexable) - m = column_vector - r = self * m - r.column(0) - end - - def **(int) - raise "Number can not the less than 1" unless int >= 1 - mat = self - (int - 1).times do - mat = mat * self - end - mat - end - - def +(other : Matrix) - raise "Dimension mismatch" if column_count != other.column_count - - rows = Array.new(row_count) do |i| - Vector.new(other.column_count) do |j| - self[i, j] + other[i, j] - end - end - - return Matrix.new(rows, other.column_count) - end - - def +(vec : Indexable) - vec = vec.is_a?(Vector) ? vec : Vector.create(vec) - self + column_vector(vec) - end - - def -(other : Matrix) - raise "Dimension mismatch" if column_count != other.column_count - - rows = Array.new(row_count) do |i| - Vector.new(other.column_count) do |j| - self[i, j] - other[i, j] - end - end - - return Matrix.new(rows, other.column_count) - end - - def -(vec : Indexable) - vec = vec.is_a?(Vector) ? vec : Vector.create(vec) - self + column_vector(vec) - end - - def /(other : Matrix) - self * other.inverse - end - - def /(vec : Indexable) - rows = self.rows.map { |row| - row.map { |e| e / other } - } - return new_matrix rows, column_count - end - - def ==(other : Matrix) - return false unless Matrix === other && - column_count == other.column_count # necessary for empty matrices - rows == other.rows - end - - # Returns element `(row, col)` of the matrix. Throws error on index error. - def [](row : Int, col : Int) - self[row][col] - end - - # Returns element `(row, col)` of the matrix, or nil if the index is not found. - def []?(row : Int, col : Int) - v = fetch(row) { nil } - v[col]? unless v.nil? - end - - # Returns the adjugate of the matrix. - def adjugate - raise "Dimention mismatch: `Matrix#adjugate` requires a square matrix." unless square? - Matrix.build(row_count, column_count) do |row, column| - cofactor(column, row) - end - end - - # Returns the (row, column) cofactor which is obtained by multiplying the first minor by (-1)**(row + column) - def cofactor(row, column) - raise "cofactor of empty matrix is not defined" if empty? - raise "Dimention mismatch: `Matrix#cofactor` requires a square matrix." unless square? - - det_of_minor = first_minor(row, column).determinant - det_of_minor * (-1.0) ** (row + column) - end - - # Returns column vector number `j` of the matrix as a `Vector` (starting at 0 like an array). - def column?(j) - return nil if j >= column_count || j < -column_count - Vector.new(row_count) { |i| rows[i][j] } - end - - # Returns column vector number `j` of the matrix as a `Vector` (starting at 0 like an array). - def column(j) - raise "Index out of range" if j >= column_count || j < -column_count - Vector.new(row_count) { |i| rows[i][j] } - end - - # Iterates over the specified column in the matrix, returning the Vector's items. - def column(j, &block) - return self if j >= column_count || j < -column_count - row_count.times do |i| - yield rows[i][j] - end - self - end - - # Returns an array of the column vectors of the matrix. See `Vector`. - def column_vectors - Matrix.new(column_count) { |i| - column(i) - } - end - - # def combine(*matrices, &block) - # Matrix.combine(self, matrices, &block) - # end - - # Iterates over each column, yielding the column - def each_column(&block) - vectors = column_vectors.map { |vec| yield(vec) } - @buffer = Matrix.columns(vectors).to_unsafe - vectors - end - - # Iterates over each row, yielding the row - def each_row(&block) - vectors = rows.map { |vec| yield(vec) } - @buffer = Matrix.rows(vectors).to_unsafe - vectors - end - - # # Hadamard product - # def hadamard_product(m) - # combine(m){|a, b| a * b} - # end - - # # Returns a new matrix resulting by stacking horizontally the receiver with the given matrices - # def hstack(*matrices) - # Matrix.hstack(self, *matrices) - # end - - # Returns the inverse of the matrix. - def inverse - raise "Dimention mismatch: `Matrix#inverse` requires a square matrix." unless square? unless square? - Matrix.identity(row_count).inverse_from(self) - end - - # :nodoc: - def inverse_from(src) - last = row_count - 1.0 - a = src.to_a - - 0.upto(last) do |k| - i = k - akk = a[k][k].abs - (k + 1).upto(last) do |j| - v = a[j][k].abs - if v > akk - i = j - akk = v - end - end - raise "Not regular" if akk == 0 - if i != k - a[i], a[k] = a[k], a[i] - rows[i], rows[k] = rows[k], rows[i] - end - akk = a[k][k] - - 0.upto(last) do |ii| - next if ii == k - q = a[ii][k] / akk - a[ii][k] = 0.0 - - (k + 1).upto(last) do |j| - a[ii][j] -= a[k][j] * q - end - 0.upto(last) do |j| - rows[ii][j] -= rows[k][j] * q - end - end - - (k + 1).upto(last) do |j| - a[k][j] = a[k][j] / akk - end - 0.upto(last) do |j| - rows[k][j] = rows[k][j] / akk - end - end - self - end - - def determinant - raise "Dimention mismatch: `Matrix#determinant` requires a square matrix." unless square? - m = rows - case row_count - # Up to 4x4, give result using Laplacian expansion by minors. - # This will typically be faster, as well as giving good results - # in case of Floats - when 0 - +1 - when 1 - +m[0][0] - when 2 - +m[0][0] * m[1][1] - m[0][1] * m[1][0] - when 3 - m0, m1, m2 = m - +m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \ - - m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \ - + m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0] - when 4 - m0, m1, m2, m3 = m - +m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \ - - m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \ - + m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \ - - m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \ - + m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \ - - m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \ - + m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \ - - m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \ - + m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \ - - m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \ - + m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \ - - m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0] - else - # For bigger matrices, use an efficient and general algorithm. - # Currently, we use the Gauss-Bareiss algorithm - determinant_bareiss - end - end - - # Returns the determinant of the matrix, using - # Bareiss' multistep integer-preserving gaussian elimination. - # It has the same computational cost order O(n^3) as standard Gaussian elimination. - # Intermediate results are fraction free and of lower complexity. - # A matrix of Integers will have thus intermediate results that are also Integers, - # with smaller bignums (if any), while a matrix of Float will usually have - # intermediate results with better precision. # - private def determinant_bareiss - size = row_count - last = size - 1 - a = to_a - no_pivot = Proc(Int32).new { return 0 } - sign = +1 - pivot = 1 - size.times do |k| - previous_pivot = pivot - if (pivot = a[k][k]) == 0 - switch = (k + 1...size).find(0) { |row| - a[row][k] != 0 - } + # ``` + # Matrix.zero(2) + # # => [ 0, 0, + # # 0, 0 ] + # ``` + def self.zero(row_count, column_count = row_count) + rows = Array(T).new(row_count) { Array(T).new(column_count, T.new(0)) } + new(rows, column_count) + end - a[switch], a[k] = a[k], a[switch] - pivot = a[k][k] - sign = -sign + # Creates a single-row matrix where the values of that row are as given in + # `row`. + # + # ``` + # Matrix.row_vector([4,5,6]) + # # => [ 4, 5, 6 ] + # ``` + def self.row_vector(row) + row = row.to_a + new([row]) + end + + # Creates a single-column matrix where the values of that column are as given + # in `column`. + # + # ``` + # Matrix.column_vector([4,5,6]) + # # => [ 4, + # # 5, + # # 6 ] + # ``` + def self.column_vector(column) + column = column.to_a + new([column].transpose, 1) + end + + # Creates a empty matrix of `row_count` x `column_count`. + # At least one of `row_count` or `column_count` must be 0. + # + # ``` + # m = Matrix(Int32).empty(2, 0) + # m == Matrix[ [], [] ] + # # => true + # n = Matrix(Int32).empty(0, 3) + # m * n + # # => Matrix[[0, 0, 0], [0, 0, 0]] + # ``` + def self.empty(row_count = 0, column_count = 0) + raise ArgumentError.new("One size must be 0") if column_count != 0 && row_count != 0 + raise ArgumentError.new("Negative size") if column_count < 0 || row_count < 0 + + new([[] of T] * row_count, column_count) + end + + # Create a matrix by stacking matrices vertically + # + # ``` + # x = Matrix[[1, 2], [3, 4]] + # y = Matrix[[5, 6], [7, 8]] + # Matrix.vstack(x, y) + # # => Matrix[[1, 2], [3, 4], [5, 6], [7, 8]] + # ``` + def Matrix.vstack(x, *matrices) + result = x.rows + matrices.each do |m| + m = m.is_a?(Matrix) ? m : rows(m) + if m.column_count != x.column_count + raise ErrDimensionMismatch.new("The given matrices must have #{x.column_count} columns, but one has #{m.column_count}") end - (k + 1).upto(last) do |i| - ai = a[i] - (k + 1).upto(last) do |j| - ai[j] = (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot - end + result.concat(m.rows) + end + new(result, x.column_count) + end + + # Create a matrix by stacking matrices horizontally + # + # ``` + # x = Matrix[[1, 2], [3, 4]] + # y = Matrix[[5, 6], [7, 8]] + # Matrix.hstack(x, y) + # # => Matrix[[1, 2, 5, 6], [3, 4, 7, 8]] + # ``` + def Matrix.hstack(x, *matrices) + result = x.rows + total_column_count = x.column_count + + matrices.each do |m| + m = m.is_a?(Matrix) ? m : rows(m) + if m.row_count != x.row_count + raise ErrDimensionMismatch.new("The given matrices must have #{x.row_count} rows, but one has #{m.row_count}") + end + + result.each_with_index do |row, i| + row.concat m.rows[i] + end + + total_column_count += m.column_count + end + + new(result, total_column_count) + end + + # Create a matrix by combining matrices entrywise, using the given block + # + # ``` + # x = Matrix[[6, 6], [4, 4]] + # y = Matrix[[1, 2], [3, 4]] + # Matrix.combine(x, y) {|a, b| a - b} + # # => Matrix[[5, 4], [1, 0]] + # ``` + def self.combine(*matrices, &block) + return Matrix.empty if matrices.empty? + + matrices = matrices.map { |m| m = m.is_a?(Matrix) ? m : rows(m) } + x = matrices.first + matrices.each do |m| + raise ErrDimensionMismatch.new unless x.row_count == m.row_count && x.column_count == m.column_count + end + + rows = Array(T).new(x.row_count) do |i| + Array(T).new(x.column_count) do |j| + yield matrices.map{ |m| m[i,j] } end end - sign * pivot + + new(rows, x.column_count) end - def first_minor(row, column) - raise "first_minor of empty matrix is not defined" if empty? - - unless 0 <= row && row < row_count - raise ArgumentError.new("invalid row (#{row.inspect} for 0..#{row_count - 1})") - end - - unless 0 <= column && column < column_count - raise ArgumentError.new("invalid column (#{column.inspect} for 0..#{column_count - 1})") - end - - arrays = to_a.map(&.to_a) - arrays.delete_at(row) - arrays.each do |array| - array.delete_at(column) - end - - Matrix.new(arrays, column_count - 1) + # ditto + def combine(*matrices, &block) + Matrix.combine(self, *matrices, &block) end - # Returns the Laplace expansion along given row or column. - def laplace_expansion(*, row = nil, column = nil) - num = row || column - - if !num || (row && column) - raise ArgumentError.new("exactly one the row or column arguments must be specified") - end - - raise "Dimention mismatch: `Matrix#determinant` requires a square matrix." unless square? - raise "laplace_expansion of empty matrix is not defined" if empty? - - unless 0 <= num && num < row_count - raise ArgumentError.new("invalid num (#{num.inspect} for 0..#{row_count - 1})") - end - - if row - row(num).map_with_index { |e, k| - e * cofactor(num, k) - }.reduce(&.+) - else - column(num).map_with_index { |e, k| - e * cofactor(k, num) - }.reduce(&.+) - end + private def initialize(rows : Array(Array(T)), column_count = nil) + # No checking is done at this point. rows must be an Array of Arrays. + # column_count must be the size of the first row, if there is one, + # otherwise it *must* be specified and can be any integer >= 0 + @rows = rows + @column_count = column_count || rows[0].try &.size end - def row(i) - self[i - 1].dup + # Returns element (`i`, `j`) of the matrix. That is: row `i`, column `j`. + # Throws if either index is not found. + def [](i, j) + @rows[i][j] end - def rows(start = 0, stop = row_count) - rows = [] of Vector - start.upto(stop) do |i| - rows << self[i - 1] - end - rows + # Returns element (`i`, `j`) of the matrix. That is: row `i`, column `j`. + # Returns nil if either index is not found. + def []?(i, j) + @rows[i]?.try &.[j]? end - def square? - row_count == column_count + # Set the value at index (`i`, `j`). That is: row `i`, column `j`. + def []=(i, j, v : T) + @rows[i][j] = v end - def transpose - return Matrix.empty(column_count, 0) if row_count.zero? - transposed = rows.map { |v| v.to_a }.transpose - Matrix.new(transposed, row_count) + # Returns the number of rows. + def row_count + @rows.size end - def to_s(io) - if empty? - "Matrix.empty(#{row_count}, #{column_count})" - else - io << "Matrix[" + # Returns the number of columns. + getter column_count : Int32 - io << map { |row| - "{" + row.to_a.map { |e| e.to_s }.join(", ") + "}" - }.join(", ") - - io << "]" - end - end - - def pretty_print(pp) : Nil - pp.list("[", self, "]") do |vec| - pp.group do - vec.to_a.pretty_print(pp) - end - end - end - - def to_json(json : JSON::Builder) - json.array do - each &.to_json(json) - end - end - - def to_unsafe - @buffer - end - - @[AlwaysInline] - def unsafe_fetch(index : Int) - @buffer[index] - end - - # To be in compliance with `Indexable` - private def size - @row_count + def unsafe_fetch(i) + @rows.unsafe_fetch(i) end end end