mirror of
https://github.com/ruby/ruby.git
synced 2022-11-09 12:17:21 -05:00
871b4612b0
iff means if and only if, but readers without that knowledge might assume this to be a spelling mistake. To me, this seems like exclusionary language that is unnecessary. Simply using "whether" when describing return values should suffice.
2493 lines
63 KiB
Ruby
2493 lines
63 KiB
Ruby
# encoding: utf-8
|
|
# frozen_string_literal: false
|
|
#
|
|
# = matrix.rb
|
|
#
|
|
# An implementation of Matrix and Vector classes.
|
|
#
|
|
# See classes Matrix and Vector for documentation.
|
|
#
|
|
# Current Maintainer:: Marc-André Lafortune
|
|
# Original Author:: Keiju ISHITSUKA
|
|
# Original Documentation:: Gavin Sinclair (sourced from <i>Ruby in a Nutshell</i> (Matsumoto, O'Reilly))
|
|
##
|
|
|
|
require_relative "matrix/version"
|
|
|
|
module ExceptionForMatrix # :nodoc:
|
|
class ErrDimensionMismatch < StandardError
|
|
def initialize(val = nil)
|
|
if val
|
|
super(val)
|
|
else
|
|
super("Dimension mismatch")
|
|
end
|
|
end
|
|
end
|
|
|
|
class ErrNotRegular < StandardError
|
|
def initialize(val = nil)
|
|
if val
|
|
super(val)
|
|
else
|
|
super("Not Regular Matrix")
|
|
end
|
|
end
|
|
end
|
|
|
|
class ErrOperationNotDefined < StandardError
|
|
def initialize(vals)
|
|
if vals.is_a?(Array)
|
|
super("Operation(#{vals[0]}) can't be defined: #{vals[1]} op #{vals[2]}")
|
|
else
|
|
super(vals)
|
|
end
|
|
end
|
|
end
|
|
|
|
class ErrOperationNotImplemented < StandardError
|
|
def initialize(vals)
|
|
super("Sorry, Operation(#{vals[0]}) not implemented: #{vals[1]} op #{vals[2]}")
|
|
end
|
|
end
|
|
end
|
|
|
|
#
|
|
# The +Matrix+ class represents a mathematical matrix. It provides methods for creating
|
|
# matrices, operating on them arithmetically and algebraically,
|
|
# and determining their mathematical properties such as trace, rank, inverse, determinant,
|
|
# or eigensystem.
|
|
#
|
|
class Matrix
|
|
include Enumerable
|
|
include ExceptionForMatrix
|
|
autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition"
|
|
autoload :LUPDecomposition, "matrix/lup_decomposition"
|
|
|
|
# instance creations
|
|
private_class_method :new
|
|
attr_reader :rows
|
|
protected :rows
|
|
|
|
#
|
|
# Creates a matrix where each argument is a row.
|
|
# Matrix[ [25, 93], [-1, 66] ]
|
|
# # => 25 93
|
|
# # -1 66
|
|
#
|
|
def Matrix.[](*rows)
|
|
rows(rows, false)
|
|
end
|
|
|
|
#
|
|
# 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 Matrix.rows(rows, copy = true)
|
|
rows = convert_to_array(rows, copy)
|
|
rows.map! do |row|
|
|
convert_to_array(row, copy)
|
|
end
|
|
size = (rows[0] || []).size
|
|
rows.each do |row|
|
|
raise ErrDimensionMismatch, "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 Matrix.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 Matrix.build(row_count, column_count = row_count)
|
|
row_count = CoercionHelper.coerce_to_int(row_count)
|
|
column_count = CoercionHelper.coerce_to_int(column_count)
|
|
raise ArgumentError if row_count < 0 || column_count < 0
|
|
return to_enum :build, row_count, column_count unless block_given?
|
|
rows = Array.new(row_count) do |i|
|
|
Array.new(column_count) do |j|
|
|
yield i, j
|
|
end
|
|
end
|
|
new rows, column_count
|
|
end
|
|
|
|
#
|
|
# Creates a matrix where the diagonal elements are composed of +values+.
|
|
# Matrix.diagonal(9, 5, -3)
|
|
# # => 9 0 0
|
|
# # 0 5 0
|
|
# # 0 0 -3
|
|
#
|
|
def Matrix.diagonal(*values)
|
|
size = values.size
|
|
return Matrix.empty if size == 0
|
|
rows = Array.new(size) {|j|
|
|
row = Array.new(size, 0)
|
|
row[j] = values[j]
|
|
row
|
|
}
|
|
new rows
|
|
end
|
|
|
|
#
|
|
# Creates an +n+ by +n+ diagonal matrix where each diagonal element is
|
|
# +value+.
|
|
# Matrix.scalar(2, 5)
|
|
# # => 5 0
|
|
# # 0 5
|
|
#
|
|
def Matrix.scalar(n, value)
|
|
diagonal(*Array.new(n, value))
|
|
end
|
|
|
|
#
|
|
# Creates an +n+ by +n+ identity matrix.
|
|
# Matrix.identity(2)
|
|
# # => 1 0
|
|
# # 0 1
|
|
#
|
|
def Matrix.identity(n)
|
|
scalar(n, 1)
|
|
end
|
|
class << Matrix
|
|
alias_method :unit, :identity
|
|
alias_method :I, :identity
|
|
end
|
|
|
|
#
|
|
# Creates a zero matrix.
|
|
# Matrix.zero(2)
|
|
# # => 0 0
|
|
# # 0 0
|
|
#
|
|
def Matrix.zero(row_count, column_count = row_count)
|
|
rows = Array.new(row_count){Array.new(column_count, 0)}
|
|
new rows, column_count
|
|
end
|
|
|
|
#
|
|
# 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 Matrix.row_vector(row)
|
|
row = convert_to_array(row)
|
|
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 Matrix.column_vector(column)
|
|
column = convert_to_array(column)
|
|
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.empty(2, 0)
|
|
# m == Matrix[ [], [] ]
|
|
# # => true
|
|
# n = Matrix.empty(0, 3)
|
|
# n == Matrix.columns([ [], [], [] ])
|
|
# # => true
|
|
# m * n
|
|
# # => Matrix[[0, 0, 0], [0, 0, 0]]
|
|
#
|
|
def Matrix.empty(row_count = 0, column_count = 0)
|
|
raise ArgumentError, "One size must be 0" if column_count != 0 && row_count != 0
|
|
raise ArgumentError, "Negative size" if column_count < 0 || row_count < 0
|
|
|
|
new([[]]*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)
|
|
x = CoercionHelper.coerce_to_matrix(x)
|
|
result = x.send(:rows).map(&:dup)
|
|
matrices.each do |m|
|
|
m = CoercionHelper.coerce_to_matrix(m)
|
|
if m.column_count != x.column_count
|
|
raise ErrDimensionMismatch, "The given matrices must have #{x.column_count} columns, but one has #{m.column_count}"
|
|
end
|
|
result.concat(m.send(: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)
|
|
x = CoercionHelper.coerce_to_matrix(x)
|
|
result = x.send(:rows).map(&:dup)
|
|
total_column_count = x.column_count
|
|
matrices.each do |m|
|
|
m = CoercionHelper.coerce_to_matrix(m)
|
|
if m.row_count != x.row_count
|
|
raise ErrDimensionMismatch, "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.send(:rows)[i]
|
|
end
|
|
total_column_count += m.column_count
|
|
end
|
|
new result, total_column_count
|
|
end
|
|
|
|
# :call-seq:
|
|
# Matrix.combine(*matrices) { |*elements| ... }
|
|
#
|
|
# 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 Matrix.combine(*matrices)
|
|
return to_enum(__method__, *matrices) unless block_given?
|
|
|
|
return Matrix.empty if matrices.empty?
|
|
matrices.map!(&CoercionHelper.method(:coerce_to_matrix))
|
|
x = matrices.first
|
|
matrices.each do |m|
|
|
raise ErrDimensionMismatch unless x.row_count == m.row_count && x.column_count == m.column_count
|
|
end
|
|
|
|
rows = Array.new(x.row_count) do |i|
|
|
Array.new(x.column_count) do |j|
|
|
yield matrices.map{|m| m[i,j]}
|
|
end
|
|
end
|
|
new rows, x.column_count
|
|
end
|
|
|
|
# :call-seq:
|
|
# combine(*other_matrices) { |*elements| ... }
|
|
#
|
|
# Creates new matrix by combining with <i>other_matrices</i> entrywise,
|
|
# using the given block.
|
|
#
|
|
# x = Matrix[[6, 6], [4, 4]]
|
|
# y = Matrix[[1, 2], [3, 4]]
|
|
# x.combine(y) {|a, b| a - b} # => Matrix[[5, 4], [1, 0]]
|
|
def combine(*matrices, &block)
|
|
Matrix.combine(self, *matrices, &block)
|
|
end
|
|
|
|
#
|
|
# Matrix.new is private; use ::rows, ::columns, ::[], etc... to create.
|
|
#
|
|
def initialize(rows, column_count = rows[0].size)
|
|
# 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
|
|
end
|
|
|
|
private def new_matrix(rows, column_count = rows[0].size) # :nodoc:
|
|
self.class.send(:new, rows, column_count) # bypass privacy of Matrix.new
|
|
end
|
|
|
|
#
|
|
# Returns element (+i+,+j+) of the matrix. That is: row +i+, column +j+.
|
|
#
|
|
def [](i, j)
|
|
@rows.fetch(i){return nil}[j]
|
|
end
|
|
alias element []
|
|
alias component []
|
|
|
|
#
|
|
# :call-seq:
|
|
# matrix[range, range] = matrix/element
|
|
# matrix[range, integer] = vector/column_matrix/element
|
|
# matrix[integer, range] = vector/row_matrix/element
|
|
# matrix[integer, integer] = element
|
|
#
|
|
# Set element or elements of matrix.
|
|
def []=(i, j, v)
|
|
raise FrozenError, "can't modify frozen Matrix" if frozen?
|
|
rows = check_range(i, :row) or row = check_int(i, :row)
|
|
columns = check_range(j, :column) or column = check_int(j, :column)
|
|
if rows && columns
|
|
set_row_and_col_range(rows, columns, v)
|
|
elsif rows
|
|
set_row_range(rows, column, v)
|
|
elsif columns
|
|
set_col_range(row, columns, v)
|
|
else
|
|
set_value(row, column, v)
|
|
end
|
|
end
|
|
alias set_element []=
|
|
alias set_component []=
|
|
private :set_element, :set_component
|
|
|
|
# Returns range or nil
|
|
private def check_range(val, direction)
|
|
return unless val.is_a?(Range)
|
|
count = direction == :row ? row_count : column_count
|
|
CoercionHelper.check_range(val, count, direction)
|
|
end
|
|
|
|
private def check_int(val, direction)
|
|
count = direction == :row ? row_count : column_count
|
|
CoercionHelper.check_int(val, count, direction)
|
|
end
|
|
|
|
private def set_value(row, col, value)
|
|
raise ErrDimensionMismatch, "Expected a a value, got a #{value.class}" if value.respond_to?(:to_matrix)
|
|
|
|
@rows[row][col] = value
|
|
end
|
|
|
|
private def set_row_and_col_range(row_range, col_range, value)
|
|
if value.is_a?(Matrix)
|
|
if row_range.size != value.row_count || col_range.size != value.column_count
|
|
raise ErrDimensionMismatch, [
|
|
'Expected a Matrix of dimensions',
|
|
"#{row_range.size}x#{col_range.size}",
|
|
'got',
|
|
"#{value.row_count}x#{value.column_count}",
|
|
].join(' ')
|
|
end
|
|
source = value.instance_variable_get :@rows
|
|
row_range.each_with_index do |row, i|
|
|
@rows[row][col_range] = source[i]
|
|
end
|
|
elsif value.is_a?(Vector)
|
|
raise ErrDimensionMismatch, 'Expected a Matrix or a value, got a Vector'
|
|
else
|
|
value_to_set = Array.new(col_range.size, value)
|
|
row_range.each do |i|
|
|
@rows[i][col_range] = value_to_set
|
|
end
|
|
end
|
|
end
|
|
|
|
private def set_row_range(row_range, col, value)
|
|
if value.is_a?(Vector)
|
|
raise ErrDimensionMismatch unless row_range.size == value.size
|
|
set_column_vector(row_range, col, value)
|
|
elsif value.is_a?(Matrix)
|
|
raise ErrDimensionMismatch unless value.column_count == 1
|
|
value = value.column(0)
|
|
raise ErrDimensionMismatch unless row_range.size == value.size
|
|
set_column_vector(row_range, col, value)
|
|
else
|
|
@rows[row_range].each{|e| e[col] = value }
|
|
end
|
|
end
|
|
|
|
private def set_column_vector(row_range, col, value)
|
|
value.each_with_index do |e, index|
|
|
r = row_range.begin + index
|
|
@rows[r][col] = e
|
|
end
|
|
end
|
|
|
|
private def set_col_range(row, col_range, value)
|
|
value = if value.is_a?(Vector)
|
|
value.to_a
|
|
elsif value.is_a?(Matrix)
|
|
raise ErrDimensionMismatch unless value.row_count == 1
|
|
value.row(0).to_a
|
|
else
|
|
Array.new(col_range.size, value)
|
|
end
|
|
raise ErrDimensionMismatch unless col_range.size == value.size
|
|
@rows[row][col_range] = value
|
|
end
|
|
|
|
#
|
|
# Returns the number of rows.
|
|
#
|
|
def row_count
|
|
@rows.size
|
|
end
|
|
|
|
alias_method :row_size, :row_count
|
|
#
|
|
# Returns the number of columns.
|
|
#
|
|
attr_reader :column_count
|
|
alias_method :column_size, :column_count
|
|
|
|
#
|
|
# Returns row vector number +i+ of the matrix as a Vector (starting at 0 like
|
|
# an array). When a block is given, the elements of that vector are iterated.
|
|
#
|
|
def row(i, &block) # :yield: e
|
|
if block_given?
|
|
@rows.fetch(i){return self}.each(&block)
|
|
self
|
|
else
|
|
Vector.elements(@rows.fetch(i){return nil})
|
|
end
|
|
end
|
|
|
|
#
|
|
# Returns column vector number +j+ of the matrix as a Vector (starting at 0
|
|
# like an array). When a block is given, the elements of that vector are
|
|
# iterated.
|
|
#
|
|
def column(j) # :yield: e
|
|
if block_given?
|
|
return self if j >= column_count || j < -column_count
|
|
row_count.times do |i|
|
|
yield @rows[i][j]
|
|
end
|
|
self
|
|
else
|
|
return nil if j >= column_count || j < -column_count
|
|
col = Array.new(row_count) {|i|
|
|
@rows[i][j]
|
|
}
|
|
Vector.elements(col, false)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Returns a matrix that is the result of iteration of the given block over all
|
|
# elements of the matrix.
|
|
# Elements can be restricted by passing an argument:
|
|
# * :all (default): yields all elements
|
|
# * :diagonal: yields only elements on the diagonal
|
|
# * :off_diagonal: yields all elements except on the diagonal
|
|
# * :lower: yields only elements on or below the diagonal
|
|
# * :strict_lower: yields only elements below the diagonal
|
|
# * :strict_upper: yields only elements above the diagonal
|
|
# * :upper: yields only elements on or above the diagonal
|
|
# Matrix[ [1,2], [3,4] ].collect { |e| e**2 }
|
|
# # => 1 4
|
|
# # 9 16
|
|
#
|
|
def collect(which = :all, &block) # :yield: e
|
|
return to_enum(:collect, which) unless block_given?
|
|
dup.collect!(which, &block)
|
|
end
|
|
alias_method :map, :collect
|
|
|
|
#
|
|
# Invokes the given block for each element of matrix, replacing the element with the value
|
|
# returned by the block.
|
|
# Elements can be restricted by passing an argument:
|
|
# * :all (default): yields all elements
|
|
# * :diagonal: yields only elements on the diagonal
|
|
# * :off_diagonal: yields all elements except on the diagonal
|
|
# * :lower: yields only elements on or below the diagonal
|
|
# * :strict_lower: yields only elements below the diagonal
|
|
# * :strict_upper: yields only elements above the diagonal
|
|
# * :upper: yields only elements on or above the diagonal
|
|
#
|
|
def collect!(which = :all)
|
|
return to_enum(:collect!, which) unless block_given?
|
|
raise FrozenError, "can't modify frozen Matrix" if frozen?
|
|
each_with_index(which){ |e, row_index, col_index| @rows[row_index][col_index] = yield e }
|
|
end
|
|
|
|
alias map! collect!
|
|
|
|
def freeze
|
|
@rows.each(&:freeze).freeze
|
|
|
|
super
|
|
end
|
|
|
|
#
|
|
# Yields all elements of the matrix, starting with those of the first row,
|
|
# or returns an Enumerator if no block given.
|
|
# Elements can be restricted by passing an argument:
|
|
# * :all (default): yields all elements
|
|
# * :diagonal: yields only elements on the diagonal
|
|
# * :off_diagonal: yields all elements except on the diagonal
|
|
# * :lower: yields only elements on or below the diagonal
|
|
# * :strict_lower: yields only elements below the diagonal
|
|
# * :strict_upper: yields only elements above the diagonal
|
|
# * :upper: yields only elements on or above the diagonal
|
|
#
|
|
# Matrix[ [1,2], [3,4] ].each { |e| puts e }
|
|
# # => prints the numbers 1 to 4
|
|
# Matrix[ [1,2], [3,4] ].each(:strict_lower).to_a # => [3]
|
|
#
|
|
def each(which = :all, &block) # :yield: e
|
|
return to_enum :each, which unless block_given?
|
|
last = column_count - 1
|
|
case which
|
|
when :all
|
|
@rows.each do |row|
|
|
row.each(&block)
|
|
end
|
|
when :diagonal
|
|
@rows.each_with_index do |row, row_index|
|
|
yield row.fetch(row_index){return self}
|
|
end
|
|
when :off_diagonal
|
|
@rows.each_with_index do |row, row_index|
|
|
column_count.times do |col_index|
|
|
yield row[col_index] unless row_index == col_index
|
|
end
|
|
end
|
|
when :lower
|
|
@rows.each_with_index do |row, row_index|
|
|
0.upto([row_index, last].min) do |col_index|
|
|
yield row[col_index]
|
|
end
|
|
end
|
|
when :strict_lower
|
|
@rows.each_with_index do |row, row_index|
|
|
[row_index, column_count].min.times do |col_index|
|
|
yield row[col_index]
|
|
end
|
|
end
|
|
when :strict_upper
|
|
@rows.each_with_index do |row, row_index|
|
|
(row_index+1).upto(last) do |col_index|
|
|
yield row[col_index]
|
|
end
|
|
end
|
|
when :upper
|
|
@rows.each_with_index do |row, row_index|
|
|
row_index.upto(last) do |col_index|
|
|
yield row[col_index]
|
|
end
|
|
end
|
|
else
|
|
raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
|
|
end
|
|
self
|
|
end
|
|
|
|
#
|
|
# Same as #each, but the row index and column index in addition to the element
|
|
#
|
|
# Matrix[ [1,2], [3,4] ].each_with_index do |e, row, col|
|
|
# puts "#{e} at #{row}, #{col}"
|
|
# end
|
|
# # => Prints:
|
|
# # 1 at 0, 0
|
|
# # 2 at 0, 1
|
|
# # 3 at 1, 0
|
|
# # 4 at 1, 1
|
|
#
|
|
def each_with_index(which = :all) # :yield: e, row, column
|
|
return to_enum :each_with_index, which unless block_given?
|
|
last = column_count - 1
|
|
case which
|
|
when :all
|
|
@rows.each_with_index do |row, row_index|
|
|
row.each_with_index do |e, col_index|
|
|
yield e, row_index, col_index
|
|
end
|
|
end
|
|
when :diagonal
|
|
@rows.each_with_index do |row, row_index|
|
|
yield row.fetch(row_index){return self}, row_index, row_index
|
|
end
|
|
when :off_diagonal
|
|
@rows.each_with_index do |row, row_index|
|
|
column_count.times do |col_index|
|
|
yield row[col_index], row_index, col_index unless row_index == col_index
|
|
end
|
|
end
|
|
when :lower
|
|
@rows.each_with_index do |row, row_index|
|
|
0.upto([row_index, last].min) do |col_index|
|
|
yield row[col_index], row_index, col_index
|
|
end
|
|
end
|
|
when :strict_lower
|
|
@rows.each_with_index do |row, row_index|
|
|
[row_index, column_count].min.times do |col_index|
|
|
yield row[col_index], row_index, col_index
|
|
end
|
|
end
|
|
when :strict_upper
|
|
@rows.each_with_index do |row, row_index|
|
|
(row_index+1).upto(last) do |col_index|
|
|
yield row[col_index], row_index, col_index
|
|
end
|
|
end
|
|
when :upper
|
|
@rows.each_with_index do |row, row_index|
|
|
row_index.upto(last) do |col_index|
|
|
yield row[col_index], row_index, col_index
|
|
end
|
|
end
|
|
else
|
|
raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
|
|
end
|
|
self
|
|
end
|
|
|
|
SELECTORS = {all: true, diagonal: true, off_diagonal: true, lower: true, strict_lower: true, strict_upper: true, upper: true}.freeze
|
|
#
|
|
# :call-seq:
|
|
# index(value, selector = :all) -> [row, column]
|
|
# index(selector = :all){ block } -> [row, column]
|
|
# index(selector = :all) -> an_enumerator
|
|
#
|
|
# The index method is specialized to return the index as [row, column]
|
|
# It also accepts an optional +selector+ argument, see #each for details.
|
|
#
|
|
# Matrix[ [1,2], [3,4] ].index(&:even?) # => [0, 1]
|
|
# Matrix[ [1,1], [1,1] ].index(1, :strict_lower) # => [1, 0]
|
|
#
|
|
def index(*args)
|
|
raise ArgumentError, "wrong number of arguments(#{args.size} for 0-2)" if args.size > 2
|
|
which = (args.size == 2 || SELECTORS.include?(args.last)) ? args.pop : :all
|
|
return to_enum :find_index, which, *args unless block_given? || args.size == 1
|
|
if args.size == 1
|
|
value = args.first
|
|
each_with_index(which) do |e, row_index, col_index|
|
|
return row_index, col_index if e == value
|
|
end
|
|
else
|
|
each_with_index(which) do |e, row_index, col_index|
|
|
return row_index, col_index if yield e
|
|
end
|
|
end
|
|
nil
|
|
end
|
|
alias_method :find_index, :index
|
|
|
|
#
|
|
# Returns a section of the matrix. The parameters are either:
|
|
# * start_row, nrows, start_col, ncols; OR
|
|
# * row_range, col_range
|
|
#
|
|
# Matrix.diagonal(9, 5, -3).minor(0..1, 0..2)
|
|
# # => 9 0 0
|
|
# # 0 5 0
|
|
#
|
|
# Like Array#[], negative indices count backward from the end of the
|
|
# row or column (-1 is the last element). Returns nil if the starting
|
|
# row or column is greater than row_count or column_count respectively.
|
|
#
|
|
def minor(*param)
|
|
case param.size
|
|
when 2
|
|
row_range, col_range = param
|
|
from_row = row_range.first
|
|
from_row += row_count if from_row < 0
|
|
to_row = row_range.end
|
|
to_row += row_count if to_row < 0
|
|
to_row += 1 unless row_range.exclude_end?
|
|
size_row = to_row - from_row
|
|
|
|
from_col = col_range.first
|
|
from_col += column_count if from_col < 0
|
|
to_col = col_range.end
|
|
to_col += column_count if to_col < 0
|
|
to_col += 1 unless col_range.exclude_end?
|
|
size_col = to_col - from_col
|
|
when 4
|
|
from_row, size_row, from_col, size_col = param
|
|
return nil if size_row < 0 || size_col < 0
|
|
from_row += row_count if from_row < 0
|
|
from_col += column_count if from_col < 0
|
|
else
|
|
raise ArgumentError, param.inspect
|
|
end
|
|
|
|
return nil if from_row > row_count || from_col > column_count || from_row < 0 || from_col < 0
|
|
rows = @rows[from_row, size_row].collect{|row|
|
|
row[from_col, size_col]
|
|
}
|
|
new_matrix rows, [column_count - from_col, size_col].min
|
|
end
|
|
|
|
#
|
|
# Returns the submatrix obtained by deleting the specified row and column.
|
|
#
|
|
# Matrix.diagonal(9, 5, -3, 4).first_minor(1, 2)
|
|
# # => 9 0 0
|
|
# # 0 0 0
|
|
# # 0 0 4
|
|
#
|
|
def first_minor(row, column)
|
|
raise RuntimeError, "first_minor of empty matrix is not defined" if empty?
|
|
|
|
unless 0 <= row && row < row_count
|
|
raise ArgumentError, "invalid row (#{row.inspect} for 0..#{row_count - 1})"
|
|
end
|
|
|
|
unless 0 <= column && column < column_count
|
|
raise ArgumentError, "invalid column (#{column.inspect} for 0..#{column_count - 1})"
|
|
end
|
|
|
|
arrays = to_a
|
|
arrays.delete_at(row)
|
|
arrays.each do |array|
|
|
array.delete_at(column)
|
|
end
|
|
|
|
new_matrix arrays, column_count - 1
|
|
end
|
|
|
|
#
|
|
# Returns the (row, column) cofactor which is obtained by multiplying
|
|
# the first minor by (-1)**(row + column).
|
|
#
|
|
# Matrix.diagonal(9, 5, -3, 4).cofactor(1, 1)
|
|
# # => -108
|
|
#
|
|
def cofactor(row, column)
|
|
raise RuntimeError, "cofactor of empty matrix is not defined" if empty?
|
|
raise ErrDimensionMismatch unless square?
|
|
|
|
det_of_minor = first_minor(row, column).determinant
|
|
det_of_minor * (-1) ** (row + column)
|
|
end
|
|
|
|
#
|
|
# Returns the adjugate of the matrix.
|
|
#
|
|
# Matrix[ [7,6],[3,9] ].adjugate
|
|
# # => 9 -6
|
|
# # -3 7
|
|
#
|
|
def adjugate
|
|
raise ErrDimensionMismatch unless square?
|
|
Matrix.build(row_count, column_count) do |row, column|
|
|
cofactor(column, row)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Returns the Laplace expansion along given row or column.
|
|
#
|
|
# Matrix[[7,6], [3,9]].laplace_expansion(column: 1)
|
|
# # => 45
|
|
#
|
|
# Matrix[[Vector[1, 0], Vector[0, 1]], [2, 3]].laplace_expansion(row: 0)
|
|
# # => Vector[3, -2]
|
|
#
|
|
#
|
|
def laplace_expansion(row: nil, column: nil)
|
|
num = row || column
|
|
|
|
if !num || (row && column)
|
|
raise ArgumentError, "exactly one the row or column arguments must be specified"
|
|
end
|
|
|
|
raise ErrDimensionMismatch unless square?
|
|
raise RuntimeError, "laplace_expansion of empty matrix is not defined" if empty?
|
|
|
|
unless 0 <= num && num < row_count
|
|
raise ArgumentError, "invalid num (#{num.inspect} for 0..#{row_count - 1})"
|
|
end
|
|
|
|
send(row ? :row : :column, num).map.with_index { |e, k|
|
|
e * cofactor(*(row ? [num, k] : [k,num]))
|
|
}.inject(:+)
|
|
end
|
|
alias_method :cofactor_expansion, :laplace_expansion
|
|
|
|
|
|
#--
|
|
# TESTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Returns +true+ if this is a diagonal matrix.
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def diagonal?
|
|
raise ErrDimensionMismatch unless square?
|
|
each(:off_diagonal).all?(&:zero?)
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is an empty matrix, i.e. if the number of rows
|
|
# or the number of columns is 0.
|
|
#
|
|
def empty?
|
|
column_count == 0 || row_count == 0
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is an hermitian matrix.
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def hermitian?
|
|
raise ErrDimensionMismatch unless square?
|
|
each_with_index(:upper).all? do |e, row, col|
|
|
e == rows[col][row].conj
|
|
end
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a lower triangular matrix.
|
|
#
|
|
def lower_triangular?
|
|
each(:strict_upper).all?(&:zero?)
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a normal matrix.
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def normal?
|
|
raise ErrDimensionMismatch unless square?
|
|
rows.each_with_index do |row_i, i|
|
|
rows.each_with_index do |row_j, j|
|
|
s = 0
|
|
rows.each_with_index do |row_k, k|
|
|
s += row_i[k] * row_j[k].conj - row_k[i].conj * row_k[j]
|
|
end
|
|
return false unless s == 0
|
|
end
|
|
end
|
|
true
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is an orthogonal matrix
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def orthogonal?
|
|
raise ErrDimensionMismatch unless square?
|
|
|
|
rows.each_with_index do |row_i, i|
|
|
rows.each_with_index do |row_j, j|
|
|
s = 0
|
|
row_count.times do |k|
|
|
s += row_i[k] * row_j[k]
|
|
end
|
|
return false unless s == (i == j ? 1 : 0)
|
|
end
|
|
end
|
|
true
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a permutation matrix
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def permutation?
|
|
raise ErrDimensionMismatch unless square?
|
|
cols = Array.new(column_count)
|
|
rows.each_with_index do |row, i|
|
|
found = false
|
|
row.each_with_index do |e, j|
|
|
if e == 1
|
|
return false if found || cols[j]
|
|
found = cols[j] = true
|
|
elsif e != 0
|
|
return false
|
|
end
|
|
end
|
|
return false unless found
|
|
end
|
|
true
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if all entries of the matrix are real.
|
|
#
|
|
def real?
|
|
all?(&:real?)
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a regular (i.e. non-singular) matrix.
|
|
#
|
|
def regular?
|
|
not singular?
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a singular matrix.
|
|
#
|
|
def singular?
|
|
determinant == 0
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a square matrix.
|
|
#
|
|
def square?
|
|
column_count == row_count
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a symmetric matrix.
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def symmetric?
|
|
raise ErrDimensionMismatch unless square?
|
|
each_with_index(:strict_upper) do |e, row, col|
|
|
return false if e != rows[col][row]
|
|
end
|
|
true
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is an antisymmetric matrix.
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def antisymmetric?
|
|
raise ErrDimensionMismatch unless square?
|
|
each_with_index(:upper) do |e, row, col|
|
|
return false unless e == -rows[col][row]
|
|
end
|
|
true
|
|
end
|
|
alias_method :skew_symmetric?, :antisymmetric?
|
|
|
|
#
|
|
# Returns +true+ if this is a unitary matrix
|
|
# Raises an error if matrix is not square.
|
|
#
|
|
def unitary?
|
|
raise ErrDimensionMismatch unless square?
|
|
rows.each_with_index do |row_i, i|
|
|
rows.each_with_index do |row_j, j|
|
|
s = 0
|
|
row_count.times do |k|
|
|
s += row_i[k].conj * row_j[k]
|
|
end
|
|
return false unless s == (i == j ? 1 : 0)
|
|
end
|
|
end
|
|
true
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is an upper triangular matrix.
|
|
#
|
|
def upper_triangular?
|
|
each(:strict_lower).all?(&:zero?)
|
|
end
|
|
|
|
#
|
|
# Returns +true+ if this is a matrix with only zero elements
|
|
#
|
|
def zero?
|
|
all?(&:zero?)
|
|
end
|
|
|
|
#--
|
|
# OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Returns whether the two matrices contain equal elements.
|
|
#
|
|
def ==(other)
|
|
return false unless Matrix === other &&
|
|
column_count == other.column_count # necessary for empty matrices
|
|
rows == other.rows
|
|
end
|
|
|
|
def eql?(other)
|
|
return false unless Matrix === other &&
|
|
column_count == other.column_count # necessary for empty matrices
|
|
rows.eql? other.rows
|
|
end
|
|
|
|
#
|
|
# Called for dup & clone.
|
|
#
|
|
private def initialize_copy(m)
|
|
super
|
|
@rows = @rows.map(&:dup) unless frozen?
|
|
end
|
|
|
|
#
|
|
# Returns a hash-code for the matrix.
|
|
#
|
|
def hash
|
|
@rows.hash
|
|
end
|
|
|
|
#--
|
|
# ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Matrix multiplication.
|
|
# Matrix[[2,4], [6,8]] * Matrix.identity(2)
|
|
# # => 2 4
|
|
# # 6 8
|
|
#
|
|
def *(m) # m is matrix or vector or number
|
|
case(m)
|
|
when Numeric
|
|
new_rows = @rows.collect {|row|
|
|
row.collect {|e| e * m }
|
|
}
|
|
return new_matrix new_rows, column_count
|
|
when Vector
|
|
m = self.class.column_vector(m)
|
|
r = self * m
|
|
return r.column(0)
|
|
when Matrix
|
|
raise ErrDimensionMismatch if column_count != m.row_count
|
|
m_rows = m.rows
|
|
new_rows = rows.map do |row_i|
|
|
Array.new(m.column_count) do |j|
|
|
vij = 0
|
|
column_count.times do |k|
|
|
vij += row_i[k] * m_rows[k][j]
|
|
end
|
|
vij
|
|
end
|
|
end
|
|
return new_matrix new_rows, m.column_count
|
|
else
|
|
return apply_through_coercion(m, __method__)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Matrix addition.
|
|
# Matrix.scalar(2,5) + Matrix[[1,0], [-4,7]]
|
|
# # => 6 0
|
|
# # -4 12
|
|
#
|
|
def +(m)
|
|
case m
|
|
when Numeric
|
|
raise ErrOperationNotDefined, ["+", self.class, m.class]
|
|
when Vector
|
|
m = self.class.column_vector(m)
|
|
when Matrix
|
|
else
|
|
return apply_through_coercion(m, __method__)
|
|
end
|
|
|
|
raise ErrDimensionMismatch unless row_count == m.row_count && column_count == m.column_count
|
|
|
|
rows = Array.new(row_count) {|i|
|
|
Array.new(column_count) {|j|
|
|
self[i, j] + m[i, j]
|
|
}
|
|
}
|
|
new_matrix rows, column_count
|
|
end
|
|
|
|
#
|
|
# Matrix subtraction.
|
|
# Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]]
|
|
# # => -8 2
|
|
# # 8 1
|
|
#
|
|
def -(m)
|
|
case m
|
|
when Numeric
|
|
raise ErrOperationNotDefined, ["-", self.class, m.class]
|
|
when Vector
|
|
m = self.class.column_vector(m)
|
|
when Matrix
|
|
else
|
|
return apply_through_coercion(m, __method__)
|
|
end
|
|
|
|
raise ErrDimensionMismatch unless row_count == m.row_count && column_count == m.column_count
|
|
|
|
rows = Array.new(row_count) {|i|
|
|
Array.new(column_count) {|j|
|
|
self[i, j] - m[i, j]
|
|
}
|
|
}
|
|
new_matrix rows, column_count
|
|
end
|
|
|
|
#
|
|
# Matrix division (multiplication by the inverse).
|
|
# Matrix[[7,6], [3,9]] / Matrix[[2,9], [3,1]]
|
|
# # => -7 1
|
|
# # -3 -6
|
|
#
|
|
def /(other)
|
|
case other
|
|
when Numeric
|
|
rows = @rows.collect {|row|
|
|
row.collect {|e| e / other }
|
|
}
|
|
return new_matrix rows, column_count
|
|
when Matrix
|
|
return self * other.inverse
|
|
else
|
|
return apply_through_coercion(other, __method__)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Hadamard product
|
|
# Matrix[[1,2], [3,4]].hadamard_product(Matrix[[1,2], [3,2]])
|
|
# # => 1 4
|
|
# # 9 8
|
|
#
|
|
def hadamard_product(m)
|
|
combine(m){|a, b| a * b}
|
|
end
|
|
alias_method :entrywise_product, :hadamard_product
|
|
|
|
#
|
|
# Returns the inverse of the matrix.
|
|
# Matrix[[-1, -1], [0, -1]].inverse
|
|
# # => -1 1
|
|
# # 0 -1
|
|
#
|
|
def inverse
|
|
raise ErrDimensionMismatch unless square?
|
|
self.class.I(row_count).send(:inverse_from, self)
|
|
end
|
|
alias_method :inv, :inverse
|
|
|
|
private def inverse_from(src) # :nodoc:
|
|
last = row_count - 1
|
|
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 ErrNotRegular 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].quo(akk)
|
|
a[ii][k] = 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].quo(akk)
|
|
end
|
|
0.upto(last) do |j|
|
|
@rows[k][j] = @rows[k][j].quo(akk)
|
|
end
|
|
end
|
|
self
|
|
end
|
|
|
|
#
|
|
# Matrix exponentiation.
|
|
# Equivalent to multiplying the matrix by itself N times.
|
|
# Non integer exponents will be handled by diagonalizing the matrix.
|
|
#
|
|
# Matrix[[7,6], [3,9]] ** 2
|
|
# # => 67 96
|
|
# # 48 99
|
|
#
|
|
def **(exp)
|
|
case exp
|
|
when Integer
|
|
case
|
|
when exp == 0
|
|
raise ErrDimensionMismatch unless square?
|
|
self.class.identity(column_count)
|
|
when exp < 0
|
|
inverse.power_int(-exp)
|
|
else
|
|
power_int(exp)
|
|
end
|
|
when Numeric
|
|
v, d, v_inv = eigensystem
|
|
v * self.class.diagonal(*d.each(:diagonal).map{|e| e ** exp}) * v_inv
|
|
else
|
|
raise ErrOperationNotDefined, ["**", self.class, exp.class]
|
|
end
|
|
end
|
|
|
|
protected def power_int(exp)
|
|
# assumes `exp` is an Integer > 0
|
|
#
|
|
# Previous algorithm:
|
|
# build M**2, M**4 = (M**2)**2, M**8, ... and multiplying those you need
|
|
# e.g. M**0b1011 = M**11 = M * M**2 * M**8
|
|
# ^ ^
|
|
# (highlighted the 2 out of 5 multiplications involving `M * x`)
|
|
#
|
|
# Current algorithm has same number of multiplications but with lower exponents:
|
|
# M**11 = M * (M * M**4)**2
|
|
# ^ ^ ^
|
|
# (highlighted the 3 out of 5 multiplications involving `M * x`)
|
|
#
|
|
# This should be faster for all (non nil-potent) matrices.
|
|
case
|
|
when exp == 1
|
|
self
|
|
when exp.odd?
|
|
self * power_int(exp - 1)
|
|
else
|
|
sqrt = power_int(exp / 2)
|
|
sqrt * sqrt
|
|
end
|
|
end
|
|
|
|
def +@
|
|
self
|
|
end
|
|
|
|
# Unary matrix negation.
|
|
#
|
|
# -Matrix[[1,5], [4,2]]
|
|
# # => -1 -5
|
|
# # -4 -2
|
|
def -@
|
|
collect {|e| -e }
|
|
end
|
|
|
|
#
|
|
# Returns the absolute value elementwise
|
|
#
|
|
def abs
|
|
collect(&:abs)
|
|
end
|
|
|
|
#--
|
|
# MATRIX FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Returns the determinant of the matrix.
|
|
#
|
|
# Beware that using Float values can yield erroneous results
|
|
# because of their lack of precision.
|
|
# Consider using exact types like Rational or BigDecimal instead.
|
|
#
|
|
# Matrix[[7,6], [3,9]].determinant
|
|
# # => 45
|
|
#
|
|
def determinant
|
|
raise ErrDimensionMismatch 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
|
|
alias_method :det, :determinant
|
|
|
|
#
|
|
# Private. Use Matrix#determinant
|
|
#
|
|
# 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.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(no_pivot) {|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
|
|
|
|
#
|
|
# deprecated; use Matrix#determinant
|
|
#
|
|
def determinant_e
|
|
warn "Matrix#determinant_e is deprecated; use #determinant", uplevel: 1
|
|
determinant
|
|
end
|
|
alias_method :det_e, :determinant_e
|
|
|
|
#
|
|
# Returns a new matrix resulting by stacking horizontally
|
|
# the receiver with the given matrices
|
|
#
|
|
# x = Matrix[[1, 2], [3, 4]]
|
|
# y = Matrix[[5, 6], [7, 8]]
|
|
# x.hstack(y) # => Matrix[[1, 2, 5, 6], [3, 4, 7, 8]]
|
|
#
|
|
def hstack(*matrices)
|
|
self.class.hstack(self, *matrices)
|
|
end
|
|
|
|
#
|
|
# Returns the rank of the matrix.
|
|
# Beware that using Float values can yield erroneous results
|
|
# because of their lack of precision.
|
|
# Consider using exact types like Rational or BigDecimal instead.
|
|
#
|
|
# Matrix[[7,6], [3,9]].rank
|
|
# # => 2
|
|
#
|
|
def rank
|
|
# We currently use Bareiss' multistep integer-preserving gaussian elimination
|
|
# (see comments on determinant)
|
|
a = to_a
|
|
last_column = column_count - 1
|
|
last_row = row_count - 1
|
|
pivot_row = 0
|
|
previous_pivot = 1
|
|
0.upto(last_column) do |k|
|
|
switch_row = (pivot_row .. last_row).find {|row|
|
|
a[row][k] != 0
|
|
}
|
|
if switch_row
|
|
a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
|
|
pivot = a[pivot_row][k]
|
|
(pivot_row+1).upto(last_row) do |i|
|
|
ai = a[i]
|
|
(k+1).upto(last_column) do |j|
|
|
ai[j] = (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
|
|
end
|
|
end
|
|
pivot_row += 1
|
|
previous_pivot = pivot
|
|
end
|
|
end
|
|
pivot_row
|
|
end
|
|
|
|
#
|
|
# deprecated; use Matrix#rank
|
|
#
|
|
def rank_e
|
|
warn "Matrix#rank_e is deprecated; use #rank", uplevel: 1
|
|
rank
|
|
end
|
|
|
|
#
|
|
# Returns a new matrix with rotated elements.
|
|
# The argument specifies the rotation (defaults to `:clockwise`):
|
|
# * :clockwise, 1, -3, etc.: "turn right" - first row becomes last column
|
|
# * :half_turn, 2, -2, etc.: first row becomes last row, elements in reverse order
|
|
# * :counter_clockwise, -1, 3: "turn left" - first row becomes first column
|
|
# (but with elements in reverse order)
|
|
#
|
|
# m = Matrix[ [1, 2], [3, 4] ]
|
|
# r = m.rotate_entries(:clockwise)
|
|
# # => Matrix[[3, 1], [4, 2]]
|
|
#
|
|
def rotate_entries(rotation = :clockwise)
|
|
rotation %= 4 if rotation.respond_to? :to_int
|
|
|
|
case rotation
|
|
when 0
|
|
dup
|
|
when 1, :clockwise
|
|
new_matrix @rows.transpose.each(&:reverse!), row_count
|
|
when 2, :half_turn
|
|
new_matrix @rows.map(&:reverse).reverse!, column_count
|
|
when 3, :counter_clockwise
|
|
new_matrix @rows.transpose.reverse!, row_count
|
|
else
|
|
raise ArgumentError, "expected #{rotation.inspect} to be one of :clockwise, :counter_clockwise, :half_turn or an integer"
|
|
end
|
|
end
|
|
|
|
# Returns a matrix with entries rounded to the given precision
|
|
# (see Float#round)
|
|
#
|
|
def round(ndigits=0)
|
|
map{|e| e.round(ndigits)}
|
|
end
|
|
|
|
#
|
|
# Returns the trace (sum of diagonal elements) of the matrix.
|
|
# Matrix[[7,6], [3,9]].trace
|
|
# # => 16
|
|
#
|
|
def trace
|
|
raise ErrDimensionMismatch unless square?
|
|
(0...column_count).inject(0) do |tr, i|
|
|
tr + @rows[i][i]
|
|
end
|
|
end
|
|
alias_method :tr, :trace
|
|
|
|
#
|
|
# Returns the transpose of the matrix.
|
|
# Matrix[[1,2], [3,4], [5,6]]
|
|
# # => 1 2
|
|
# # 3 4
|
|
# # 5 6
|
|
# Matrix[[1,2], [3,4], [5,6]].transpose
|
|
# # => 1 3 5
|
|
# # 2 4 6
|
|
#
|
|
def transpose
|
|
return self.class.empty(column_count, 0) if row_count.zero?
|
|
new_matrix @rows.transpose, row_count
|
|
end
|
|
alias_method :t, :transpose
|
|
|
|
#
|
|
# Returns a new matrix resulting by stacking vertically
|
|
# the receiver with the given matrices
|
|
#
|
|
# x = Matrix[[1, 2], [3, 4]]
|
|
# y = Matrix[[5, 6], [7, 8]]
|
|
# x.vstack(y) # => Matrix[[1, 2], [3, 4], [5, 6], [7, 8]]
|
|
#
|
|
def vstack(*matrices)
|
|
self.class.vstack(self, *matrices)
|
|
end
|
|
|
|
#--
|
|
# DECOMPOSITIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
|
|
#++
|
|
|
|
#
|
|
# Returns the Eigensystem of the matrix; see +EigenvalueDecomposition+.
|
|
# m = Matrix[[1, 2], [3, 4]]
|
|
# v, d, v_inv = m.eigensystem
|
|
# d.diagonal? # => true
|
|
# v.inv == v_inv # => true
|
|
# (v * d * v_inv).round(5) == m # => true
|
|
#
|
|
def eigensystem
|
|
EigenvalueDecomposition.new(self)
|
|
end
|
|
alias_method :eigen, :eigensystem
|
|
|
|
#
|
|
# Returns the LUP decomposition of the matrix; see +LUPDecomposition+.
|
|
# a = Matrix[[1, 2], [3, 4]]
|
|
# l, u, p = a.lup
|
|
# l.lower_triangular? # => true
|
|
# u.upper_triangular? # => true
|
|
# p.permutation? # => true
|
|
# l * u == p * a # => true
|
|
# a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)]
|
|
#
|
|
def lup
|
|
LUPDecomposition.new(self)
|
|
end
|
|
alias_method :lup_decomposition, :lup
|
|
|
|
#--
|
|
# COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
|
|
#++
|
|
|
|
#
|
|
# Returns the conjugate of the matrix.
|
|
# Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
|
|
# # => 1+2i i 0
|
|
# # 1 2 3
|
|
# Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].conjugate
|
|
# # => 1-2i -i 0
|
|
# # 1 2 3
|
|
#
|
|
def conjugate
|
|
collect(&:conjugate)
|
|
end
|
|
alias_method :conj, :conjugate
|
|
|
|
#
|
|
# Returns the adjoint of the matrix.
|
|
#
|
|
# Matrix[ [i,1],[2,-i] ].adjoint
|
|
# # => -i 2
|
|
# # 1 i
|
|
#
|
|
def adjoint
|
|
conjugate.transpose
|
|
end
|
|
|
|
#
|
|
# Returns the imaginary part of the matrix.
|
|
# Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
|
|
# # => 1+2i i 0
|
|
# # 1 2 3
|
|
# Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].imaginary
|
|
# # => 2i i 0
|
|
# # 0 0 0
|
|
#
|
|
def imaginary
|
|
collect(&:imaginary)
|
|
end
|
|
alias_method :imag, :imaginary
|
|
|
|
#
|
|
# Returns the real part of the matrix.
|
|
# Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
|
|
# # => 1+2i i 0
|
|
# # 1 2 3
|
|
# Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].real
|
|
# # => 1 0 0
|
|
# # 1 2 3
|
|
#
|
|
def real
|
|
collect(&:real)
|
|
end
|
|
|
|
#
|
|
# Returns an array containing matrices corresponding to the real and imaginary
|
|
# parts of the matrix
|
|
#
|
|
# m.rect == [m.real, m.imag] # ==> true for all matrices m
|
|
#
|
|
def rect
|
|
[real, imag]
|
|
end
|
|
alias_method :rectangular, :rect
|
|
|
|
#--
|
|
# CONVERTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# The coerce method provides support for Ruby type coercion.
|
|
# This coercion mechanism is used by Ruby to handle mixed-type
|
|
# numeric operations: it is intended to find a compatible common
|
|
# type between the two operands of the operator.
|
|
# See also Numeric#coerce.
|
|
#
|
|
def coerce(other)
|
|
case other
|
|
when Numeric
|
|
return Scalar.new(other), self
|
|
else
|
|
raise TypeError, "#{self.class} can't be coerced into #{other.class}"
|
|
end
|
|
end
|
|
|
|
#
|
|
# Returns an array of the row vectors of the matrix. See Vector.
|
|
#
|
|
def row_vectors
|
|
Array.new(row_count) {|i|
|
|
row(i)
|
|
}
|
|
end
|
|
|
|
#
|
|
# Returns an array of the column vectors of the matrix. See Vector.
|
|
#
|
|
def column_vectors
|
|
Array.new(column_count) {|i|
|
|
column(i)
|
|
}
|
|
end
|
|
|
|
#
|
|
# Explicit conversion to a Matrix. Returns self
|
|
#
|
|
def to_matrix
|
|
self
|
|
end
|
|
|
|
#
|
|
# Returns an array of arrays that describe the rows of the matrix.
|
|
#
|
|
def to_a
|
|
@rows.collect(&:dup)
|
|
end
|
|
|
|
# Deprecated.
|
|
#
|
|
# Use <code>map(&:to_f)</code>
|
|
def elements_to_f
|
|
warn "Matrix#elements_to_f is deprecated, use map(&:to_f)", uplevel: 1
|
|
map(&:to_f)
|
|
end
|
|
|
|
# Deprecated.
|
|
#
|
|
# Use <code>map(&:to_i)</code>
|
|
def elements_to_i
|
|
warn "Matrix#elements_to_i is deprecated, use map(&:to_i)", uplevel: 1
|
|
map(&:to_i)
|
|
end
|
|
|
|
# Deprecated.
|
|
#
|
|
# Use <code>map(&:to_r)</code>
|
|
def elements_to_r
|
|
warn "Matrix#elements_to_r is deprecated, use map(&:to_r)", uplevel: 1
|
|
map(&:to_r)
|
|
end
|
|
|
|
#--
|
|
# PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Overrides Object#to_s
|
|
#
|
|
def to_s
|
|
if empty?
|
|
"#{self.class}.empty(#{row_count}, #{column_count})"
|
|
else
|
|
"#{self.class}[" + @rows.collect{|row|
|
|
"[" + row.collect{|e| e.to_s}.join(", ") + "]"
|
|
}.join(", ")+"]"
|
|
end
|
|
end
|
|
|
|
#
|
|
# Overrides Object#inspect
|
|
#
|
|
def inspect
|
|
if empty?
|
|
"#{self.class}.empty(#{row_count}, #{column_count})"
|
|
else
|
|
"#{self.class}#{@rows.inspect}"
|
|
end
|
|
end
|
|
|
|
# Private helper modules
|
|
|
|
module ConversionHelper # :nodoc:
|
|
#
|
|
# Converts the obj to an Array. If copy is set to true
|
|
# a copy of obj will be made if necessary.
|
|
#
|
|
private def convert_to_array(obj, copy = false) # :nodoc:
|
|
case obj
|
|
when Array
|
|
copy ? obj.dup : obj
|
|
when Vector
|
|
obj.to_a
|
|
else
|
|
begin
|
|
converted = obj.to_ary
|
|
rescue Exception => e
|
|
raise TypeError, "can't convert #{obj.class} into an Array (#{e.message})"
|
|
end
|
|
raise TypeError, "#{obj.class}#to_ary should return an Array" unless converted.is_a? Array
|
|
converted
|
|
end
|
|
end
|
|
end
|
|
|
|
extend ConversionHelper
|
|
|
|
module CoercionHelper # :nodoc:
|
|
#
|
|
# Applies the operator +oper+ with argument +obj+
|
|
# through coercion of +obj+
|
|
#
|
|
private def apply_through_coercion(obj, oper)
|
|
coercion = obj.coerce(self)
|
|
raise TypeError unless coercion.is_a?(Array) && coercion.length == 2
|
|
coercion[0].public_send(oper, coercion[1])
|
|
rescue
|
|
raise TypeError, "#{obj.inspect} can't be coerced into #{self.class}"
|
|
end
|
|
|
|
#
|
|
# Helper method to coerce a value into a specific class.
|
|
# Raises a TypeError if the coercion fails or the returned value
|
|
# is not of the right class.
|
|
# (from Rubinius)
|
|
#
|
|
def self.coerce_to(obj, cls, meth) # :nodoc:
|
|
return obj if obj.kind_of?(cls)
|
|
raise TypeError, "Expected a #{cls} but got a #{obj.class}" unless obj.respond_to? meth
|
|
begin
|
|
ret = obj.__send__(meth)
|
|
rescue Exception => e
|
|
raise TypeError, "Coercion error: #{obj.inspect}.#{meth} => #{cls} failed:\n" \
|
|
"(#{e.message})"
|
|
end
|
|
raise TypeError, "Coercion error: obj.#{meth} did NOT return a #{cls} (was #{ret.class})" unless ret.kind_of? cls
|
|
ret
|
|
end
|
|
|
|
def self.coerce_to_int(obj)
|
|
coerce_to(obj, Integer, :to_int)
|
|
end
|
|
|
|
def self.coerce_to_matrix(obj)
|
|
coerce_to(obj, Matrix, :to_matrix)
|
|
end
|
|
|
|
# Returns `nil` for non Ranges
|
|
# Checks range validity, return canonical range with 0 <= begin <= end < count
|
|
def self.check_range(val, count, kind)
|
|
canonical = (val.begin + (val.begin < 0 ? count : 0))..
|
|
(val.end ? val.end + (val.end < 0 ? count : 0) - (val.exclude_end? ? 1 : 0)
|
|
: count - 1)
|
|
unless 0 <= canonical.begin && canonical.begin <= canonical.end && canonical.end < count
|
|
raise IndexError, "given range #{val} is outside of #{kind} dimensions: 0...#{count}"
|
|
end
|
|
canonical
|
|
end
|
|
|
|
def self.check_int(val, count, kind)
|
|
val = CoercionHelper.coerce_to_int(val)
|
|
if val >= count || val < -count
|
|
raise IndexError, "given #{kind} #{val} is outside of #{-count}...#{count}"
|
|
end
|
|
val
|
|
end
|
|
end
|
|
|
|
include CoercionHelper
|
|
|
|
# Private CLASS
|
|
|
|
class Scalar < Numeric # :nodoc:
|
|
include ExceptionForMatrix
|
|
include CoercionHelper
|
|
|
|
def initialize(value)
|
|
@value = value
|
|
end
|
|
|
|
# ARITHMETIC
|
|
def +(other)
|
|
case other
|
|
when Numeric
|
|
Scalar.new(@value + other)
|
|
when Vector, Matrix
|
|
raise ErrOperationNotDefined, ["+", @value.class, other.class]
|
|
else
|
|
apply_through_coercion(other, __method__)
|
|
end
|
|
end
|
|
|
|
def -(other)
|
|
case other
|
|
when Numeric
|
|
Scalar.new(@value - other)
|
|
when Vector, Matrix
|
|
raise ErrOperationNotDefined, ["-", @value.class, other.class]
|
|
else
|
|
apply_through_coercion(other, __method__)
|
|
end
|
|
end
|
|
|
|
def *(other)
|
|
case other
|
|
when Numeric
|
|
Scalar.new(@value * other)
|
|
when Vector, Matrix
|
|
other.collect{|e| @value * e}
|
|
else
|
|
apply_through_coercion(other, __method__)
|
|
end
|
|
end
|
|
|
|
def /(other)
|
|
case other
|
|
when Numeric
|
|
Scalar.new(@value / other)
|
|
when Vector
|
|
raise ErrOperationNotDefined, ["/", @value.class, other.class]
|
|
when Matrix
|
|
self * other.inverse
|
|
else
|
|
apply_through_coercion(other, __method__)
|
|
end
|
|
end
|
|
|
|
def **(other)
|
|
case other
|
|
when Numeric
|
|
Scalar.new(@value ** other)
|
|
when Vector
|
|
raise ErrOperationNotDefined, ["**", @value.class, other.class]
|
|
when Matrix
|
|
#other.powered_by(self)
|
|
raise ErrOperationNotImplemented, ["**", @value.class, other.class]
|
|
else
|
|
apply_through_coercion(other, __method__)
|
|
end
|
|
end
|
|
end
|
|
|
|
end
|
|
|
|
|
|
#
|
|
# The +Vector+ class represents a mathematical vector, which is useful in its own right, and
|
|
# also constitutes a row or column of a Matrix.
|
|
#
|
|
# == Method Catalogue
|
|
#
|
|
# To create a Vector:
|
|
# * Vector.[](*array)
|
|
# * Vector.elements(array, copy = true)
|
|
# * Vector.basis(size: n, index: k)
|
|
# * Vector.zero(n)
|
|
#
|
|
# To access elements:
|
|
# * #[](i)
|
|
#
|
|
# To set elements:
|
|
# * #[]=(i, v)
|
|
#
|
|
# To enumerate the elements:
|
|
# * #each2(v)
|
|
# * #collect2(v)
|
|
#
|
|
# Properties of vectors:
|
|
# * #angle_with(v)
|
|
# * Vector.independent?(*vs)
|
|
# * #independent?(*vs)
|
|
# * #zero?
|
|
#
|
|
# Vector arithmetic:
|
|
# * #*(x) "is matrix or number"
|
|
# * #+(v)
|
|
# * #-(v)
|
|
# * #/(v)
|
|
# * #+@
|
|
# * #-@
|
|
#
|
|
# Vector functions:
|
|
# * #inner_product(v), #dot(v)
|
|
# * #cross_product(v), #cross(v)
|
|
# * #collect
|
|
# * #collect!
|
|
# * #magnitude
|
|
# * #map
|
|
# * #map!
|
|
# * #map2(v)
|
|
# * #norm
|
|
# * #normalize
|
|
# * #r
|
|
# * #round
|
|
# * #size
|
|
#
|
|
# Conversion to other data types:
|
|
# * #covector
|
|
# * #to_a
|
|
# * #coerce(other)
|
|
#
|
|
# String representations:
|
|
# * #to_s
|
|
# * #inspect
|
|
#
|
|
class Vector
|
|
include ExceptionForMatrix
|
|
include Enumerable
|
|
include Matrix::CoercionHelper
|
|
extend Matrix::ConversionHelper
|
|
#INSTANCE CREATION
|
|
|
|
private_class_method :new
|
|
attr_reader :elements
|
|
protected :elements
|
|
|
|
#
|
|
# Creates a Vector from a list of elements.
|
|
# Vector[7, 4, ...]
|
|
#
|
|
def Vector.[](*array)
|
|
new convert_to_array(array, false)
|
|
end
|
|
|
|
#
|
|
# Creates a vector from an Array. The optional second argument specifies
|
|
# whether the array itself or a copy is used internally.
|
|
#
|
|
def Vector.elements(array, copy = true)
|
|
new convert_to_array(array, copy)
|
|
end
|
|
|
|
#
|
|
# Returns a standard basis +n+-vector, where k is the index.
|
|
#
|
|
# Vector.basis(size:, index:) # => Vector[0, 1, 0]
|
|
#
|
|
def Vector.basis(size:, index:)
|
|
raise ArgumentError, "invalid size (#{size} for 1..)" if size < 1
|
|
raise ArgumentError, "invalid index (#{index} for 0...#{size})" unless 0 <= index && index < size
|
|
array = Array.new(size, 0)
|
|
array[index] = 1
|
|
new convert_to_array(array, false)
|
|
end
|
|
|
|
#
|
|
# Return a zero vector.
|
|
#
|
|
# Vector.zero(3) # => Vector[0, 0, 0]
|
|
#
|
|
def Vector.zero(size)
|
|
raise ArgumentError, "invalid size (#{size} for 0..)" if size < 0
|
|
array = Array.new(size, 0)
|
|
new convert_to_array(array, false)
|
|
end
|
|
|
|
#
|
|
# Vector.new is private; use Vector[] or Vector.elements to create.
|
|
#
|
|
def initialize(array)
|
|
# No checking is done at this point.
|
|
@elements = array
|
|
end
|
|
|
|
# ACCESSING
|
|
|
|
#
|
|
# :call-seq:
|
|
# vector[range]
|
|
# vector[integer]
|
|
#
|
|
# Returns element or elements of the vector.
|
|
#
|
|
def [](i)
|
|
@elements[i]
|
|
end
|
|
alias element []
|
|
alias component []
|
|
|
|
#
|
|
# :call-seq:
|
|
# vector[range] = new_vector
|
|
# vector[range] = row_matrix
|
|
# vector[range] = new_element
|
|
# vector[integer] = new_element
|
|
#
|
|
# Set element or elements of vector.
|
|
#
|
|
def []=(i, v)
|
|
raise FrozenError, "can't modify frozen Vector" if frozen?
|
|
if i.is_a?(Range)
|
|
range = Matrix::CoercionHelper.check_range(i, size, :vector)
|
|
set_range(range, v)
|
|
else
|
|
index = Matrix::CoercionHelper.check_int(i, size, :index)
|
|
set_value(index, v)
|
|
end
|
|
end
|
|
alias set_element []=
|
|
alias set_component []=
|
|
private :set_element, :set_component
|
|
|
|
private def set_value(index, value)
|
|
@elements[index] = value
|
|
end
|
|
|
|
private def set_range(range, value)
|
|
if value.is_a?(Vector)
|
|
raise ArgumentError, "vector to be set has wrong size" unless range.size == value.size
|
|
@elements[range] = value.elements
|
|
elsif value.is_a?(Matrix)
|
|
raise ErrDimensionMismatch unless value.row_count == 1
|
|
@elements[range] = value.row(0).elements
|
|
else
|
|
@elements[range] = Array.new(range.size, value)
|
|
end
|
|
end
|
|
|
|
# Returns a vector with entries rounded to the given precision
|
|
# (see Float#round)
|
|
#
|
|
def round(ndigits=0)
|
|
map{|e| e.round(ndigits)}
|
|
end
|
|
|
|
#
|
|
# Returns the number of elements in the vector.
|
|
#
|
|
def size
|
|
@elements.size
|
|
end
|
|
|
|
#--
|
|
# ENUMERATIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Iterate over the elements of this vector
|
|
#
|
|
def each(&block)
|
|
return to_enum(:each) unless block_given?
|
|
@elements.each(&block)
|
|
self
|
|
end
|
|
|
|
#
|
|
# Iterate over the elements of this vector and +v+ in conjunction.
|
|
#
|
|
def each2(v) # :yield: e1, e2
|
|
raise TypeError, "Integer is not like Vector" if v.kind_of?(Integer)
|
|
raise ErrDimensionMismatch if size != v.size
|
|
return to_enum(:each2, v) unless block_given?
|
|
size.times do |i|
|
|
yield @elements[i], v[i]
|
|
end
|
|
self
|
|
end
|
|
|
|
#
|
|
# Collects (as in Enumerable#collect) over the elements of this vector and +v+
|
|
# in conjunction.
|
|
#
|
|
def collect2(v) # :yield: e1, e2
|
|
raise TypeError, "Integer is not like Vector" if v.kind_of?(Integer)
|
|
raise ErrDimensionMismatch if size != v.size
|
|
return to_enum(:collect2, v) unless block_given?
|
|
Array.new(size) do |i|
|
|
yield @elements[i], v[i]
|
|
end
|
|
end
|
|
|
|
#--
|
|
# PROPERTIES -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Returns whether all of vectors are linearly independent.
|
|
#
|
|
# Vector.independent?(Vector[1,0], Vector[0,1])
|
|
# # => true
|
|
#
|
|
# Vector.independent?(Vector[1,2], Vector[2,4])
|
|
# # => false
|
|
#
|
|
def Vector.independent?(*vs)
|
|
vs.each do |v|
|
|
raise TypeError, "expected Vector, got #{v.class}" unless v.is_a?(Vector)
|
|
raise ErrDimensionMismatch unless v.size == vs.first.size
|
|
end
|
|
return false if vs.count > vs.first.size
|
|
Matrix[*vs].rank.eql?(vs.count)
|
|
end
|
|
|
|
#
|
|
# Returns whether all of vectors are linearly independent.
|
|
#
|
|
# Vector[1,0].independent?(Vector[0,1])
|
|
# # => true
|
|
#
|
|
# Vector[1,2].independent?(Vector[2,4])
|
|
# # => false
|
|
#
|
|
def independent?(*vs)
|
|
self.class.independent?(self, *vs)
|
|
end
|
|
|
|
#
|
|
# Returns whether all elements are zero.
|
|
#
|
|
def zero?
|
|
all?(&:zero?)
|
|
end
|
|
|
|
#
|
|
# Makes the matrix frozen and Ractor-shareable
|
|
#
|
|
def freeze
|
|
@elements.freeze
|
|
super
|
|
end
|
|
|
|
#
|
|
# Called for dup & clone.
|
|
#
|
|
private def initialize_copy(v)
|
|
super
|
|
@elements = @elements.dup unless frozen?
|
|
end
|
|
|
|
|
|
#--
|
|
# COMPARING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Returns whether the two vectors have the same elements in the same order.
|
|
#
|
|
def ==(other)
|
|
return false unless Vector === other
|
|
@elements == other.elements
|
|
end
|
|
|
|
def eql?(other)
|
|
return false unless Vector === other
|
|
@elements.eql? other.elements
|
|
end
|
|
|
|
#
|
|
# Returns a hash-code for the vector.
|
|
#
|
|
def hash
|
|
@elements.hash
|
|
end
|
|
|
|
#--
|
|
# ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Multiplies the vector by +x+, where +x+ is a number or a matrix.
|
|
#
|
|
def *(x)
|
|
case x
|
|
when Numeric
|
|
els = @elements.collect{|e| e * x}
|
|
self.class.elements(els, false)
|
|
when Matrix
|
|
Matrix.column_vector(self) * x
|
|
when Vector
|
|
raise ErrOperationNotDefined, ["*", self.class, x.class]
|
|
else
|
|
apply_through_coercion(x, __method__)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Vector addition.
|
|
#
|
|
def +(v)
|
|
case v
|
|
when Vector
|
|
raise ErrDimensionMismatch if size != v.size
|
|
els = collect2(v) {|v1, v2|
|
|
v1 + v2
|
|
}
|
|
self.class.elements(els, false)
|
|
when Matrix
|
|
Matrix.column_vector(self) + v
|
|
else
|
|
apply_through_coercion(v, __method__)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Vector subtraction.
|
|
#
|
|
def -(v)
|
|
case v
|
|
when Vector
|
|
raise ErrDimensionMismatch if size != v.size
|
|
els = collect2(v) {|v1, v2|
|
|
v1 - v2
|
|
}
|
|
self.class.elements(els, false)
|
|
when Matrix
|
|
Matrix.column_vector(self) - v
|
|
else
|
|
apply_through_coercion(v, __method__)
|
|
end
|
|
end
|
|
|
|
#
|
|
# Vector division.
|
|
#
|
|
def /(x)
|
|
case x
|
|
when Numeric
|
|
els = @elements.collect{|e| e / x}
|
|
self.class.elements(els, false)
|
|
when Matrix, Vector
|
|
raise ErrOperationNotDefined, ["/", self.class, x.class]
|
|
else
|
|
apply_through_coercion(x, __method__)
|
|
end
|
|
end
|
|
|
|
def +@
|
|
self
|
|
end
|
|
|
|
def -@
|
|
collect {|e| -e }
|
|
end
|
|
|
|
#--
|
|
# VECTOR FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Returns the inner product of this vector with the other.
|
|
# Vector[4,7].inner_product Vector[10,1] # => 47
|
|
#
|
|
def inner_product(v)
|
|
raise ErrDimensionMismatch if size != v.size
|
|
|
|
p = 0
|
|
each2(v) {|v1, v2|
|
|
p += v1 * v2.conj
|
|
}
|
|
p
|
|
end
|
|
alias_method :dot, :inner_product
|
|
|
|
#
|
|
# Returns the cross product of this vector with the others.
|
|
# Vector[1, 0, 0].cross_product Vector[0, 1, 0] # => Vector[0, 0, 1]
|
|
#
|
|
# It is generalized to other dimensions to return a vector perpendicular
|
|
# to the arguments.
|
|
# Vector[1, 2].cross_product # => Vector[-2, 1]
|
|
# Vector[1, 0, 0, 0].cross_product(
|
|
# Vector[0, 1, 0, 0],
|
|
# Vector[0, 0, 1, 0]
|
|
# ) #=> Vector[0, 0, 0, 1]
|
|
#
|
|
def cross_product(*vs)
|
|
raise ErrOperationNotDefined, "cross product is not defined on vectors of dimension #{size}" unless size >= 2
|
|
raise ArgumentError, "wrong number of arguments (#{vs.size} for #{size - 2})" unless vs.size == size - 2
|
|
vs.each do |v|
|
|
raise TypeError, "expected Vector, got #{v.class}" unless v.is_a? Vector
|
|
raise ErrDimensionMismatch unless v.size == size
|
|
end
|
|
case size
|
|
when 2
|
|
Vector[-@elements[1], @elements[0]]
|
|
when 3
|
|
v = vs[0]
|
|
Vector[ v[2]*@elements[1] - v[1]*@elements[2],
|
|
v[0]*@elements[2] - v[2]*@elements[0],
|
|
v[1]*@elements[0] - v[0]*@elements[1] ]
|
|
else
|
|
rows = self, *vs, Array.new(size) {|i| Vector.basis(size: size, index: i) }
|
|
Matrix.rows(rows).laplace_expansion(row: size - 1)
|
|
end
|
|
end
|
|
alias_method :cross, :cross_product
|
|
|
|
#
|
|
# Like Array#collect.
|
|
#
|
|
def collect(&block) # :yield: e
|
|
return to_enum(:collect) unless block_given?
|
|
els = @elements.collect(&block)
|
|
self.class.elements(els, false)
|
|
end
|
|
alias_method :map, :collect
|
|
|
|
#
|
|
# Like Array#collect!
|
|
#
|
|
def collect!(&block)
|
|
return to_enum(:collect!) unless block_given?
|
|
raise FrozenError, "can't modify frozen Vector" if frozen?
|
|
@elements.collect!(&block)
|
|
self
|
|
end
|
|
alias map! collect!
|
|
|
|
#
|
|
# Returns the modulus (Pythagorean distance) of the vector.
|
|
# Vector[5,8,2].r # => 9.643650761
|
|
#
|
|
def magnitude
|
|
Math.sqrt(@elements.inject(0) {|v, e| v + e.abs2})
|
|
end
|
|
alias_method :r, :magnitude
|
|
alias_method :norm, :magnitude
|
|
|
|
#
|
|
# Like Vector#collect2, but returns a Vector instead of an Array.
|
|
#
|
|
def map2(v, &block) # :yield: e1, e2
|
|
return to_enum(:map2, v) unless block_given?
|
|
els = collect2(v, &block)
|
|
self.class.elements(els, false)
|
|
end
|
|
|
|
class ZeroVectorError < StandardError
|
|
end
|
|
#
|
|
# Returns a new vector with the same direction but with norm 1.
|
|
# v = Vector[5,8,2].normalize
|
|
# # => Vector[0.5184758473652127, 0.8295613557843402, 0.20739033894608505]
|
|
# v.norm # => 1.0
|
|
#
|
|
def normalize
|
|
n = magnitude
|
|
raise ZeroVectorError, "Zero vectors can not be normalized" if n == 0
|
|
self / n
|
|
end
|
|
|
|
#
|
|
# Returns an angle with another vector. Result is within the [0..Math::PI].
|
|
# Vector[1,0].angle_with(Vector[0,1])
|
|
# # => Math::PI / 2
|
|
#
|
|
def angle_with(v)
|
|
raise TypeError, "Expected a Vector, got a #{v.class}" unless v.is_a?(Vector)
|
|
raise ErrDimensionMismatch if size != v.size
|
|
prod = magnitude * v.magnitude
|
|
raise ZeroVectorError, "Can't get angle of zero vector" if prod == 0
|
|
dot = inner_product(v)
|
|
if dot.abs >= prod
|
|
dot.positive? ? 0 : Math::PI
|
|
else
|
|
Math.acos(dot / prod)
|
|
end
|
|
end
|
|
|
|
#--
|
|
# CONVERTING
|
|
#++
|
|
|
|
#
|
|
# Creates a single-row matrix from this vector.
|
|
#
|
|
def covector
|
|
Matrix.row_vector(self)
|
|
end
|
|
|
|
#
|
|
# Returns the elements of the vector in an array.
|
|
#
|
|
def to_a
|
|
@elements.dup
|
|
end
|
|
|
|
#
|
|
# Return a single-column matrix from this vector
|
|
#
|
|
def to_matrix
|
|
Matrix.column_vector(self)
|
|
end
|
|
|
|
def elements_to_f
|
|
warn "Vector#elements_to_f is deprecated", uplevel: 1
|
|
map(&:to_f)
|
|
end
|
|
|
|
def elements_to_i
|
|
warn "Vector#elements_to_i is deprecated", uplevel: 1
|
|
map(&:to_i)
|
|
end
|
|
|
|
def elements_to_r
|
|
warn "Vector#elements_to_r is deprecated", uplevel: 1
|
|
map(&:to_r)
|
|
end
|
|
|
|
#
|
|
# The coerce method provides support for Ruby type coercion.
|
|
# This coercion mechanism is used by Ruby to handle mixed-type
|
|
# numeric operations: it is intended to find a compatible common
|
|
# type between the two operands of the operator.
|
|
# See also Numeric#coerce.
|
|
#
|
|
def coerce(other)
|
|
case other
|
|
when Numeric
|
|
return Matrix::Scalar.new(other), self
|
|
else
|
|
raise TypeError, "#{self.class} can't be coerced into #{other.class}"
|
|
end
|
|
end
|
|
|
|
#--
|
|
# PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
|
|
#++
|
|
|
|
#
|
|
# Overrides Object#to_s
|
|
#
|
|
def to_s
|
|
"Vector[" + @elements.join(", ") + "]"
|
|
end
|
|
|
|
#
|
|
# Overrides Object#inspect
|
|
#
|
|
def inspect
|
|
"Vector" + @elements.inspect
|
|
end
|
|
end
|