Matrix updates
This commit is contained in:
Chris
parent
3354307427
commit
77b2bc62c7
|
|
@ -1,6 +1,6 @@
|
||||||
require "./linear_algebra/error"
|
require "./linear_algebra/error"
|
||||||
# require "./linear_algebra/ndarray"
|
# require "./linear_algebra/ndarray"
|
||||||
# require "./linear_algebra/matrix"
|
require "./linear_algebra/matrix"
|
||||||
require "./linear_algebra/vector"
|
require "./linear_algebra/vector"
|
||||||
|
|
||||||
module Apatite
|
module Apatite
|
||||||
|
|
|
||||||
|
|
@ -2,569 +2,302 @@ require "json"
|
||||||
require "./vector"
|
require "./vector"
|
||||||
|
|
||||||
module Apatite::LinearAlgebra
|
module Apatite::LinearAlgebra
|
||||||
class Matrix
|
class Matrix(T)
|
||||||
include Enumerable(Vector)
|
include Enumerable(Vector)
|
||||||
include Indexable(Vector)
|
include Indexable(Vector)
|
||||||
include Comparable(Matrix)
|
include Comparable(Matrix)
|
||||||
|
|
||||||
getter column_count : Int32
|
protected getter rows : Array(Array(T))
|
||||||
|
|
||||||
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
|
|
||||||
|
|
||||||
# Creates a matrix where each argument is a row.
|
# Creates a matrix where each argument is a row.
|
||||||
|
#
|
||||||
|
# ```
|
||||||
|
# Matrix[[25, 93], [-1, 66]]
|
||||||
|
# # => [ 25, 93,
|
||||||
|
# # -1, 66 ]
|
||||||
|
# ```
|
||||||
def self.[](*rows)
|
def self.[](*rows)
|
||||||
rows(rows)
|
rows(rows, false)
|
||||||
end
|
end
|
||||||
|
|
||||||
# Creates a matrix of size `row_count x column_count`. It fills the values by calling
|
# Creates a matrix where +rows+ is an array of arrays, each of which is a row
|
||||||
# the given block, passing the current row and column.
|
# 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)
|
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
|
raise ArgumentError.new if row_count < 0 || column_count < 0
|
||||||
rows = Array.new(row_count) do |i|
|
rows = Array(T).new(row_count) do |i|
|
||||||
Vector.new(column_count) do |j|
|
Array(T).new(column_count) do |j|
|
||||||
yield i, j
|
yield i, j
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
Matrix.new(rows, column_count)
|
new(rows, column_count)
|
||||||
end
|
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`.
|
# 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
|
size = values.size
|
||||||
return Matrix.empty if size == 0
|
return Matrix(T).empty if size == 0
|
||||||
|
rows = Array(Array(T)).new(size) do |j|
|
||||||
rows = Array.new(size) do |j|
|
row = Array(T).new(size, T.new(0))
|
||||||
row = Array.new(size, 0)
|
|
||||||
row[j] = values[j]
|
row[j] = values[j]
|
||||||
row
|
row
|
||||||
end
|
end
|
||||||
|
new(rows)
|
||||||
new rows
|
|
||||||
end
|
end
|
||||||
|
|
||||||
# :ditto:
|
def self.diagonal(*values : T)
|
||||||
def self.diagonal(*values)
|
diagonal(values, nil)
|
||||||
Matrix.diagonal(values)
|
|
||||||
end
|
end
|
||||||
|
|
||||||
# Creates a empty matrix of `row_count x column_count`. At least one of
|
# Creates an +n+ by +n+ diagonal matrix where each diagonal element is
|
||||||
# `row_count` or `column_count` must be 0.
|
# `value`.
|
||||||
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.scalar(2, 5)
|
||||||
Matrix.new(([] of Vector) * row_count, column_count)
|
# # => [ 5, 0,
|
||||||
|
# # 0, 5 ]
|
||||||
|
# ```
|
||||||
|
def self.scalar(n, value : T)
|
||||||
|
diagonal(Array(T).new(n, value))
|
||||||
end
|
end
|
||||||
|
|
||||||
# Creates a new diagonal matrix of size `n` with ones in the diagonal
|
# Creates an `n` by `n` identity matrix.
|
||||||
# and zeros elsewhere.
|
#
|
||||||
def self.eye(n)
|
# ```
|
||||||
Matrix.diagonal([1] * n)
|
# Matrix.identity(2)
|
||||||
|
# # => [ 1, 0,
|
||||||
|
# # 0, 1 ]
|
||||||
|
# ```
|
||||||
|
def self.identity(n : T)
|
||||||
|
scalar(n, T.new(1))
|
||||||
end
|
end
|
||||||
|
|
||||||
# TODO
|
# ditto
|
||||||
def self.hstack(x, *matrices)
|
def self.unit(n : T)
|
||||||
end
|
identity(n)
|
||||||
|
|
||||||
# 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)
|
|
||||||
end
|
end
|
||||||
|
|
||||||
# Creates a zero matrix.
|
# 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
|
# Matrix.zero(2)
|
||||||
last = size - 1
|
# # => [ 0, 0,
|
||||||
a = to_a
|
# # 0, 0 ]
|
||||||
no_pivot = Proc(Int32).new { return 0 }
|
# ```
|
||||||
sign = +1
|
def self.zero(row_count, column_count = row_count)
|
||||||
pivot = 1
|
rows = Array(T).new(row_count) { Array(T).new(column_count, T.new(0)) }
|
||||||
size.times do |k|
|
new(rows, column_count)
|
||||||
previous_pivot = pivot
|
|
||||||
if (pivot = a[k][k]) == 0
|
|
||||||
switch = (k + 1...size).find(0) { |row|
|
|
||||||
a[row][k] != 0
|
|
||||||
}
|
|
||||||
|
|
||||||
a[switch], a[k] = a[k], a[switch]
|
|
||||||
pivot = a[k][k]
|
|
||||||
sign = -sign
|
|
||||||
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
|
|
||||||
end
|
|
||||||
end
|
|
||||||
sign * pivot
|
|
||||||
end
|
end
|
||||||
|
|
||||||
def first_minor(row, column)
|
# Creates a single-row matrix where the values of that row are as given in
|
||||||
raise "first_minor of empty matrix is not defined" if empty?
|
# `row`.
|
||||||
|
#
|
||||||
unless 0 <= row && row < row_count
|
# ```
|
||||||
raise ArgumentError.new("invalid row (#{row.inspect} for 0..#{row_count - 1})")
|
# Matrix.row_vector([4,5,6])
|
||||||
|
# # => [ 4, 5, 6 ]
|
||||||
|
# ```
|
||||||
|
def self.row_vector(row)
|
||||||
|
row = row.to_a
|
||||||
|
new([row])
|
||||||
end
|
end
|
||||||
|
|
||||||
unless 0 <= column && column < column_count
|
# Creates a single-column matrix where the values of that column are as given
|
||||||
raise ArgumentError.new("invalid column (#{column.inspect} for 0..#{column_count - 1})")
|
# 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
|
end
|
||||||
|
|
||||||
arrays = to_a.map(&.to_a)
|
# Creates a empty matrix of `row_count` x `column_count`.
|
||||||
arrays.delete_at(row)
|
# At least one of `row_count` or `column_count` must be 0.
|
||||||
arrays.each do |array|
|
#
|
||||||
array.delete_at(column)
|
# ```
|
||||||
|
# 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
|
end
|
||||||
|
|
||||||
Matrix.new(arrays, column_count - 1)
|
# 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
|
||||||
|
result.concat(m.rows)
|
||||||
|
end
|
||||||
|
new(result, x.column_count)
|
||||||
end
|
end
|
||||||
|
|
||||||
# Returns the Laplace expansion along given row or column.
|
# Create a matrix by stacking matrices horizontally
|
||||||
def laplace_expansion(*, row = nil, column = nil)
|
#
|
||||||
num = row || column
|
# ```
|
||||||
|
# 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
|
||||||
|
|
||||||
if !num || (row && column)
|
matrices.each do |m|
|
||||||
raise ArgumentError.new("exactly one the row or column arguments must be specified")
|
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
|
end
|
||||||
|
|
||||||
raise "Dimention mismatch: `Matrix#determinant` requires a square matrix." unless square?
|
result.each_with_index do |row, i|
|
||||||
raise "laplace_expansion of empty matrix is not defined" if empty?
|
row.concat m.rows[i]
|
||||||
|
|
||||||
unless 0 <= num && num < row_count
|
|
||||||
raise ArgumentError.new("invalid num (#{num.inspect} for 0..#{row_count - 1})")
|
|
||||||
end
|
end
|
||||||
|
|
||||||
if row
|
total_column_count += m.column_count
|
||||||
row(num).map_with_index { |e, k|
|
end
|
||||||
e * cofactor(num, k)
|
|
||||||
}.reduce(&.+)
|
new(result, total_column_count)
|
||||||
else
|
end
|
||||||
column(num).map_with_index { |e, k|
|
|
||||||
e * cofactor(k, num)
|
# Create a matrix by combining matrices entrywise, using the given block
|
||||||
}.reduce(&.+)
|
#
|
||||||
|
# ```
|
||||||
|
# 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
|
||||||
end
|
end
|
||||||
|
|
||||||
def row(i)
|
new(rows, x.column_count)
|
||||||
self[i - 1].dup
|
|
||||||
end
|
end
|
||||||
|
|
||||||
def rows(start = 0, stop = row_count)
|
# ditto
|
||||||
rows = [] of Vector
|
def combine(*matrices, &block)
|
||||||
start.upto(stop) do |i|
|
Matrix.combine(self, *matrices, &block)
|
||||||
rows << self[i - 1]
|
|
||||||
end
|
|
||||||
rows
|
|
||||||
end
|
end
|
||||||
|
|
||||||
def square?
|
private def initialize(rows : Array(Array(T)), column_count = nil)
|
||||||
row_count == column_count
|
# 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
|
end
|
||||||
|
|
||||||
def transpose
|
# Returns element (`i`, `j`) of the matrix. That is: row `i`, column `j`.
|
||||||
return Matrix.empty(column_count, 0) if row_count.zero?
|
# Throws if either index is not found.
|
||||||
transposed = rows.map { |v| v.to_a }.transpose
|
def [](i, j)
|
||||||
Matrix.new(transposed, row_count)
|
@rows[i][j]
|
||||||
end
|
end
|
||||||
|
|
||||||
def to_s(io)
|
# Returns element (`i`, `j`) of the matrix. That is: row `i`, column `j`.
|
||||||
if empty?
|
# Returns nil if either index is not found.
|
||||||
"Matrix.empty(#{row_count}, #{column_count})"
|
def []?(i, j)
|
||||||
else
|
@rows[i]?.try &.[j]?
|
||||||
io << "Matrix["
|
|
||||||
|
|
||||||
io << map { |row|
|
|
||||||
"{" + row.to_a.map { |e| e.to_s }.join(", ") + "}"
|
|
||||||
}.join(", ")
|
|
||||||
|
|
||||||
io << "]"
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
def pretty_print(pp) : Nil
|
# Set the value at index (`i`, `j`). That is: row `i`, column `j`.
|
||||||
pp.list("[", self, "]") do |vec|
|
def []=(i, j, v : T)
|
||||||
pp.group do
|
@rows[i][j] = v
|
||||||
vec.to_a.pretty_print(pp)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
def to_json(json : JSON::Builder)
|
# Returns the number of rows.
|
||||||
json.array do
|
def row_count
|
||||||
each &.to_json(json)
|
@rows.size
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
def to_unsafe
|
# Returns the number of columns.
|
||||||
@buffer
|
getter column_count : Int32
|
||||||
end
|
|
||||||
|
|
||||||
@[AlwaysInline]
|
def unsafe_fetch(i)
|
||||||
def unsafe_fetch(index : Int)
|
@rows.unsafe_fetch(i)
|
||||||
@buffer[index]
|
|
||||||
end
|
|
||||||
|
|
||||||
# To be in compliance with `Indexable`
|
|
||||||
private def size
|
|
||||||
@row_count
|
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue