Added a lot to Matrix
This commit is contained in:
Chris Watson
parent
131ab63faf
commit
1b69ec120e
|
|
@ -6,18 +6,21 @@ module Apatite
|
|||
include Indexable(Vector)
|
||||
include Comparable(Matrix)
|
||||
|
||||
@col_count : Int32
|
||||
getter column_count : Int32
|
||||
|
||||
getter row_count : Int32
|
||||
|
||||
@buffer : Pointer(Vector)
|
||||
|
||||
def initialize(rows : Array(Indexable(Number)), col_count : Int32 = rows[0].size)
|
||||
@buffer = rows.map { |r| Vector.create(r) }.to_unsafe
|
||||
@col_count = col_count
|
||||
def initialize(rows, column_count : Int32 = rows[0].size)
|
||||
@buffer = rows.map { |r| Vector.create(r) }.to_a.to_unsafe
|
||||
@row_count = rows.size
|
||||
@column_count = column_count
|
||||
end
|
||||
|
||||
# Creates a matrix where each argument is a row.
|
||||
def self.[](*rows)
|
||||
rows(rows, false)
|
||||
rows(rows)
|
||||
end
|
||||
|
||||
# Creates a matrix of size `row_count x column_count`. It fills the values by calling
|
||||
|
|
@ -39,28 +42,29 @@ module Apatite
|
|||
|
||||
# Creates a matrix using `columns` as an array of column vectors.
|
||||
def self.columns(columns)
|
||||
rows(columns, false).transpose
|
||||
rows(columns).transpose
|
||||
end
|
||||
|
||||
# Create a matrix by combining matrices entrywise, using the given block.
|
||||
def self.combine(matrices, &block)
|
||||
return Matrix.empty if matrices.empty?
|
||||
x = matrices.first
|
||||
matrices.each do |m|
|
||||
rase "Dimension mismatch" unless x.row_count == m.row_count && x.column_count == m.column_count
|
||||
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
|
||||
new rows, x.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
|
||||
|
||||
def self.combine(*matrices, &block)
|
||||
Matrix.combine(matrices, &block)
|
||||
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)
|
||||
|
|
@ -100,7 +104,7 @@ module Apatite
|
|||
|
||||
# Creates a single-row matrix where the values of that row are as given in `row`.
|
||||
def self.row_vector(row)
|
||||
Matrix.new([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.
|
||||
|
|
@ -126,5 +130,399 @@ module Apatite
|
|||
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
|
||||
col = Array(Float64).new(row_count) { |i|
|
||||
rows[i][j]
|
||||
}
|
||||
Vector.create(col)
|
||||
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
|
||||
col = Array(Float64).new(row_count) { |i|
|
||||
rows[i][j]
|
||||
}
|
||||
Vector.create(col)
|
||||
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
|
||||
Array(Vector).new(column_count) { |i|
|
||||
column(i)
|
||||
}
|
||||
end
|
||||
|
||||
# def combine(*matrices, &block)
|
||||
# Matrix.combine(self, matrices, &block)
|
||||
# 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
|
||||
}
|
||||
|
||||
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
|
||||
|
||||
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)
|
||||
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
|
||||
end
|
||||
|
||||
def row(i, &block)
|
||||
@rows.fetch(i){return self}.each(&block)
|
||||
self
|
||||
end
|
||||
|
||||
def row(i)
|
||||
Vector.create(rows.fetch(i) { [] of Float64 })
|
||||
end
|
||||
|
||||
def rows
|
||||
rows = [] of Vector
|
||||
row_count.times do |i|
|
||||
rows << self[i - 1]
|
||||
end
|
||||
rows
|
||||
end
|
||||
|
||||
def square?
|
||||
row_count == column_count
|
||||
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)
|
||||
end
|
||||
|
||||
def to_s(io)
|
||||
if empty?
|
||||
"Matrix.empty(#{row_count}, #{column_count})"
|
||||
else
|
||||
io << "Matrix["
|
||||
|
||||
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
|
||||
|
||||
@[AlwaysInline]
|
||||
def unsafe_fetch(index : Int)
|
||||
@buffer[index]
|
||||
end
|
||||
|
||||
# To be in compliance with `Indexable`
|
||||
private def size
|
||||
@row_count
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -614,7 +614,7 @@ module Apatite
|
|||
|
||||
# Transpose this vector into a 1xn `Matrix`
|
||||
def transpose
|
||||
Matrix.column_vector(self)
|
||||
Matrix.col_vector(self)
|
||||
end
|
||||
|
||||
# Returns a copy of the vector with elements set to `value` if
|
||||
|
|
@ -739,7 +739,7 @@ module Apatite
|
|||
end
|
||||
|
||||
def to_s(io)
|
||||
io << "Vector{"
|
||||
io << "{"
|
||||
join ", ", io, &.inspect(io)
|
||||
io << "}"
|
||||
end
|
||||
|
|
@ -810,9 +810,9 @@ module Apatite
|
|||
return map { |v| yield(v, value) }
|
||||
end
|
||||
|
||||
values = value.is_a?(Vector) ? value.elements : value
|
||||
values = value.is_a?(Vector) ? value : Vector.create(value)
|
||||
|
||||
unless @elements.size == values.size
|
||||
unless size == values.size
|
||||
raise "Cannot perform operations on vectors with different dimensions."
|
||||
end
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue