summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--doc/maintainers.rdoc6
-rw-r--r--doc/standard_library.rdoc2
-rw-r--r--gems/bundled_gems1
-rw-r--r--lib/matrix.rb2493
-rw-r--r--lib/matrix/eigenvalue_decomposition.rb882
-rw-r--r--lib/matrix/lup_decomposition.rb219
-rw-r--r--lib/matrix/matrix.gemspec26
-rw-r--r--lib/matrix/version.rb5
-rwxr-xr-xmisc/expand_tabs.rb1
-rw-r--r--test/matrix/test_matrix.rb888
-rw-r--r--test/matrix/test_vector.rb335
-rw-r--r--tool/sync_default_gems.rb1
12 files changed, 4 insertions, 4855 deletions
diff --git a/doc/maintainers.rdoc b/doc/maintainers.rdoc
index 805fbc227b..79c30b8732 100644
--- a/doc/maintainers.rdoc
+++ b/doc/maintainers.rdoc
@@ -150,10 +150,6 @@ Yukihiro Matsumoto (matz)
Naotoshi Seo (sonots)
https://github.com/ruby/logger
https://rubygems.org/gems/logger
-[lib/matrix.rb]
- Marc-André Lafortune (marcandre)
- https://github.com/ruby/matrix
- https://rubygems.org/gems/matrix
[lib/mutex_m.rb]
Keiju ISHITSUKA (keiju)
https://github.com/ruby/mutex_m
@@ -395,6 +391,8 @@ Yukihiro Matsumoto (matz)
https://github.com/ruby/rexml
[rss]
https://github.com/ruby/rss
+[matrix]
+ https://github.com/ruby/matrix
[prime]
https://github.com/ruby/prime
[rbs]
diff --git a/doc/standard_library.rdoc b/doc/standard_library.rdoc
index 28dedcb24f..e05eb6f62d 100644
--- a/doc/standard_library.rdoc
+++ b/doc/standard_library.rdoc
@@ -45,7 +45,6 @@ IPAddr:: Provides methods to manipulate IPv4 and IPv6 IP addresses
IRB:: Interactive Ruby command-line tool for REPL (Read Eval Print Loop)
OptionParser:: Ruby-oriented class for command-line option analysis
Logger:: Provides a simple logging utility for outputting messages
-Matrix:: Represents a mathematical matrix.
Mutex_m:: Mixin to extend objects to be handled like a Mutex
Net::FTP:: Support for the File Transfer Protocol
Net::HTTP:: HTTP client api for Ruby
@@ -110,6 +109,7 @@ Rake:: Ruby build program with capabilities similar to make
Test::Unit:: A compatibility layer for MiniTest
REXML:: An XML toolkit for Ruby
RSS:: Family of libraries that support various formats of XML "feeds"
+Matrix:: Represents a mathematical matrix.
Prime:: Prime numbers and factorization library
RBS:: RBS is a language to describe the structure of Ruby programs
TypeProf:: A type analysis tool for Ruby code based on abstract interpretation
diff --git a/gems/bundled_gems b/gems/bundled_gems
index 982faeb8f0..b42c6794d3 100644
--- a/gems/bundled_gems
+++ b/gems/bundled_gems
@@ -5,6 +5,7 @@ rake 13.0.3 https://github.com/ruby/rake
test-unit 3.4.1 https://github.com/test-unit/test-unit 3.4.1
rexml 3.2.5 https://github.com/ruby/rexml
rss 0.2.9 https://github.com/ruby/rss 0.2.9
+matrix 0.4.1 https://github.com/ruby/matrix
prime 0.1.2 https://github.com/ruby/prime
typeprof 0.14.1 https://github.com/ruby/typeprof
rbs 1.2.0 https://github.com/ruby/rbs
diff --git a/lib/matrix.rb b/lib/matrix.rb
deleted file mode 100644
index 26cda55a4a..0000000000
--- a/lib/matrix.rb
+++ /dev/null
@@ -1,2493 +0,0 @@
-# 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 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
diff --git a/lib/matrix/eigenvalue_decomposition.rb b/lib/matrix/eigenvalue_decomposition.rb
deleted file mode 100644
index bf6637635a..0000000000
--- a/lib/matrix/eigenvalue_decomposition.rb
+++ /dev/null
@@ -1,882 +0,0 @@
-# frozen_string_literal: false
-class Matrix
- # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
-
- # Eigenvalues and eigenvectors of a real matrix.
- #
- # Computes the eigenvalues and eigenvectors of a matrix A.
- #
- # If A is diagonalizable, this provides matrices V and D
- # such that A = V*D*V.inv, where D is the diagonal matrix with entries
- # equal to the eigenvalues and V is formed by the eigenvectors.
- #
- # If A is symmetric, then V is orthogonal and thus A = V*D*V.t
-
- class EigenvalueDecomposition
-
- # Constructs the eigenvalue decomposition for a square matrix +A+
- #
- def initialize(a)
- # @d, @e: Arrays for internal storage of eigenvalues.
- # @v: Array for internal storage of eigenvectors.
- # @h: Array for internal storage of nonsymmetric Hessenberg form.
- raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
- @size = a.row_count
- @d = Array.new(@size, 0)
- @e = Array.new(@size, 0)
-
- if (@symmetric = a.symmetric?)
- @v = a.to_a
- tridiagonalize
- diagonalize
- else
- @v = Array.new(@size) { Array.new(@size, 0) }
- @h = a.to_a
- @ort = Array.new(@size, 0)
- reduce_to_hessenberg
- hessenberg_to_real_schur
- end
- end
-
- # Returns the eigenvector matrix +V+
- #
- def eigenvector_matrix
- Matrix.send(:new, build_eigenvectors.transpose)
- end
- alias_method :v, :eigenvector_matrix
-
- # Returns the inverse of the eigenvector matrix +V+
- #
- def eigenvector_matrix_inv
- r = Matrix.send(:new, build_eigenvectors)
- r = r.transpose.inverse unless @symmetric
- r
- end
- alias_method :v_inv, :eigenvector_matrix_inv
-
- # Returns the eigenvalues in an array
- #
- def eigenvalues
- values = @d.dup
- @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0}
- values
- end
-
- # Returns an array of the eigenvectors
- #
- def eigenvectors
- build_eigenvectors.map{|ev| Vector.send(:new, ev)}
- end
-
- # Returns the block diagonal eigenvalue matrix +D+
- #
- def eigenvalue_matrix
- Matrix.diagonal(*eigenvalues)
- end
- alias_method :d, :eigenvalue_matrix
-
- # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv]
- #
- def to_ary
- [v, d, v_inv]
- end
- alias_method :to_a, :to_ary
-
-
- private def build_eigenvectors
- # JAMA stores complex eigenvectors in a strange way
- # See http://web.archive.org/web/20111016032731/http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html
- @e.each_with_index.map do |imag, i|
- if imag == 0
- Array.new(@size){|j| @v[j][i]}
- elsif imag > 0
- Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])}
- else
- Array.new(@size){|j| Complex(@v[j][i-1], -@v[j][i])}
- end
- end
- end
-
- # Complex scalar division.
-
- private def cdiv(xr, xi, yr, yi)
- if (yr.abs > yi.abs)
- r = yi/yr
- d = yr + r*yi
- [(xr + r*xi)/d, (xi - r*xr)/d]
- else
- r = yr/yi
- d = yi + r*yr
- [(r*xr + xi)/d, (r*xi - xr)/d]
- end
- end
-
-
- # Symmetric Householder reduction to tridiagonal form.
-
- private def tridiagonalize
-
- # This is derived from the Algol procedures tred2 by
- # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
- # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
- # Fortran subroutine in EISPACK.
-
- @size.times do |j|
- @d[j] = @v[@size-1][j]
- end
-
- # Householder reduction to tridiagonal form.
-
- (@size-1).downto(0+1) do |i|
-
- # Scale to avoid under/overflow.
-
- scale = 0.0
- h = 0.0
- i.times do |k|
- scale = scale + @d[k].abs
- end
- if (scale == 0.0)
- @e[i] = @d[i-1]
- i.times do |j|
- @d[j] = @v[i-1][j]
- @v[i][j] = 0.0
- @v[j][i] = 0.0
- end
- else
-
- # Generate Householder vector.
-
- i.times do |k|
- @d[k] /= scale
- h += @d[k] * @d[k]
- end
- f = @d[i-1]
- g = Math.sqrt(h)
- if (f > 0)
- g = -g
- end
- @e[i] = scale * g
- h -= f * g
- @d[i-1] = f - g
- i.times do |j|
- @e[j] = 0.0
- end
-
- # Apply similarity transformation to remaining columns.
-
- i.times do |j|
- f = @d[j]
- @v[j][i] = f
- g = @e[j] + @v[j][j] * f
- (j+1).upto(i-1) do |k|
- g += @v[k][j] * @d[k]
- @e[k] += @v[k][j] * f
- end
- @e[j] = g
- end
- f = 0.0
- i.times do |j|
- @e[j] /= h
- f += @e[j] * @d[j]
- end
- hh = f / (h + h)
- i.times do |j|
- @e[j] -= hh * @d[j]
- end
- i.times do |j|
- f = @d[j]
- g = @e[j]
- j.upto(i-1) do |k|
- @v[k][j] -= (f * @e[k] + g * @d[k])
- end
- @d[j] = @v[i-1][j]
- @v[i][j] = 0.0
- end
- end
- @d[i] = h
- end
-
- # Accumulate transformations.
-
- 0.upto(@size-1-1) do |i|
- @v[@size-1][i] = @v[i][i]
- @v[i][i] = 1.0
- h = @d[i+1]
- if (h != 0.0)
- 0.upto(i) do |k|
- @d[k] = @v[k][i+1] / h
- end
- 0.upto(i) do |j|
- g = 0.0
- 0.upto(i) do |k|
- g += @v[k][i+1] * @v[k][j]
- end
- 0.upto(i) do |k|
- @v[k][j] -= g * @d[k]
- end
- end
- end
- 0.upto(i) do |k|
- @v[k][i+1] = 0.0
- end
- end
- @size.times do |j|
- @d[j] = @v[@size-1][j]
- @v[@size-1][j] = 0.0
- end
- @v[@size-1][@size-1] = 1.0
- @e[0] = 0.0
- end
-
-
- # Symmetric tridiagonal QL algorithm.
-
- private def diagonalize
- # This is derived from the Algol procedures tql2, by
- # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
- # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
- # Fortran subroutine in EISPACK.
-
- 1.upto(@size-1) do |i|
- @e[i-1] = @e[i]
- end
- @e[@size-1] = 0.0
-
- f = 0.0
- tst1 = 0.0
- eps = Float::EPSILON
- @size.times do |l|
-
- # Find small subdiagonal element
-
- tst1 = [tst1, @d[l].abs + @e[l].abs].max
- m = l
- while (m < @size) do
- if (@e[m].abs <= eps*tst1)
- break
- end
- m+=1
- end
-
- # If m == l, @d[l] is an eigenvalue,
- # otherwise, iterate.
-
- if (m > l)
- iter = 0
- begin
- iter = iter + 1 # (Could check iteration count here.)
-
- # Compute implicit shift
-
- g = @d[l]
- p = (@d[l+1] - g) / (2.0 * @e[l])
- r = Math.hypot(p, 1.0)
- if (p < 0)
- r = -r
- end
- @d[l] = @e[l] / (p + r)
- @d[l+1] = @e[l] * (p + r)
- dl1 = @d[l+1]
- h = g - @d[l]
- (l+2).upto(@size-1) do |i|
- @d[i] -= h
- end
- f += h
-
- # Implicit QL transformation.
-
- p = @d[m]
- c = 1.0
- c2 = c
- c3 = c
- el1 = @e[l+1]
- s = 0.0
- s2 = 0.0
- (m-1).downto(l) do |i|
- c3 = c2
- c2 = c
- s2 = s
- g = c * @e[i]
- h = c * p
- r = Math.hypot(p, @e[i])
- @e[i+1] = s * r
- s = @e[i] / r
- c = p / r
- p = c * @d[i] - s * g
- @d[i+1] = h + s * (c * g + s * @d[i])
-
- # Accumulate transformation.
-
- @size.times do |k|
- h = @v[k][i+1]
- @v[k][i+1] = s * @v[k][i] + c * h
- @v[k][i] = c * @v[k][i] - s * h
- end
- end
- p = -s * s2 * c3 * el1 * @e[l] / dl1
- @e[l] = s * p
- @d[l] = c * p
-
- # Check for convergence.
-
- end while (@e[l].abs > eps*tst1)
- end
- @d[l] = @d[l] + f
- @e[l] = 0.0
- end
-
- # Sort eigenvalues and corresponding vectors.
-
- 0.upto(@size-2) do |i|
- k = i
- p = @d[i]
- (i+1).upto(@size-1) do |j|
- if (@d[j] < p)
- k = j
- p = @d[j]
- end
- end
- if (k != i)
- @d[k] = @d[i]
- @d[i] = p
- @size.times do |j|
- p = @v[j][i]
- @v[j][i] = @v[j][k]
- @v[j][k] = p
- end
- end
- end
- end
-
- # Nonsymmetric reduction to Hessenberg form.
-
- private def reduce_to_hessenberg
- # This is derived from the Algol procedures orthes and ortran,
- # by Martin and Wilkinson, Handbook for Auto. Comp.,
- # Vol.ii-Linear Algebra, and the corresponding
- # Fortran subroutines in EISPACK.
-
- low = 0
- high = @size-1
-
- (low+1).upto(high-1) do |m|
-
- # Scale column.
-
- scale = 0.0
- m.upto(high) do |i|
- scale = scale + @h[i][m-1].abs
- end
- if (scale != 0.0)
-
- # Compute Householder transformation.
-
- h = 0.0
- high.downto(m) do |i|
- @ort[i] = @h[i][m-1]/scale
- h += @ort[i] * @ort[i]
- end
- g = Math.sqrt(h)
- if (@ort[m] > 0)
- g = -g
- end
- h -= @ort[m] * g
- @ort[m] = @ort[m] - g
-
- # Apply Householder similarity transformation
- # @h = (I-u*u'/h)*@h*(I-u*u')/h)
-
- m.upto(@size-1) do |j|
- f = 0.0
- high.downto(m) do |i|
- f += @ort[i]*@h[i][j]
- end
- f = f/h
- m.upto(high) do |i|
- @h[i][j] -= f*@ort[i]
- end
- end
-
- 0.upto(high) do |i|
- f = 0.0
- high.downto(m) do |j|
- f += @ort[j]*@h[i][j]
- end
- f = f/h
- m.upto(high) do |j|
- @h[i][j] -= f*@ort[j]
- end
- end
- @ort[m] = scale*@ort[m]
- @h[m][m-1] = scale*g
- end
- end
-
- # Accumulate transformations (Algol's ortran).
-
- @size.times do |i|
- @size.times do |j|
- @v[i][j] = (i == j ? 1.0 : 0.0)
- end
- end
-
- (high-1).downto(low+1) do |m|
- if (@h[m][m-1] != 0.0)
- (m+1).upto(high) do |i|
- @ort[i] = @h[i][m-1]
- end
- m.upto(high) do |j|
- g = 0.0
- m.upto(high) do |i|
- g += @ort[i] * @v[i][j]
- end
- # Double division avoids possible underflow
- g = (g / @ort[m]) / @h[m][m-1]
- m.upto(high) do |i|
- @v[i][j] += g * @ort[i]
- end
- end
- end
- end
- end
-
- # Nonsymmetric reduction from Hessenberg to real Schur form.
-
- private def hessenberg_to_real_schur
-
- # This is derived from the Algol procedure hqr2,
- # by Martin and Wilkinson, Handbook for Auto. Comp.,
- # Vol.ii-Linear Algebra, and the corresponding
- # Fortran subroutine in EISPACK.
-
- # Initialize
-
- nn = @size
- n = nn-1
- low = 0
- high = nn-1
- eps = Float::EPSILON
- exshift = 0.0
- p = q = r = s = z = 0
-
- # Store roots isolated by balanc and compute matrix norm
-
- norm = 0.0
- nn.times do |i|
- if (i < low || i > high)
- @d[i] = @h[i][i]
- @e[i] = 0.0
- end
- ([i-1, 0].max).upto(nn-1) do |j|
- norm = norm + @h[i][j].abs
- end
- end
-
- # Outer loop over eigenvalue index
-
- iter = 0
- while (n >= low) do
-
- # Look for single small sub-diagonal element
-
- l = n
- while (l > low) do
- s = @h[l-1][l-1].abs + @h[l][l].abs
- if (s == 0.0)
- s = norm
- end
- if (@h[l][l-1].abs < eps * s)
- break
- end
- l-=1
- end
-
- # Check for convergence
- # One root found
-
- if (l == n)
- @h[n][n] = @h[n][n] + exshift
- @d[n] = @h[n][n]
- @e[n] = 0.0
- n-=1
- iter = 0
-
- # Two roots found
-
- elsif (l == n-1)
- w = @h[n][n-1] * @h[n-1][n]
- p = (@h[n-1][n-1] - @h[n][n]) / 2.0
- q = p * p + w
- z = Math.sqrt(q.abs)
- @h[n][n] = @h[n][n] + exshift
- @h[n-1][n-1] = @h[n-1][n-1] + exshift
- x = @h[n][n]
-
- # Real pair
-
- if (q >= 0)
- if (p >= 0)
- z = p + z
- else
- z = p - z
- end
- @d[n-1] = x + z
- @d[n] = @d[n-1]
- if (z != 0.0)
- @d[n] = x - w / z
- end
- @e[n-1] = 0.0
- @e[n] = 0.0
- x = @h[n][n-1]
- s = x.abs + z.abs
- p = x / s
- q = z / s
- r = Math.sqrt(p * p+q * q)
- p /= r
- q /= r
-
- # Row modification
-
- (n-1).upto(nn-1) do |j|
- z = @h[n-1][j]
- @h[n-1][j] = q * z + p * @h[n][j]
- @h[n][j] = q * @h[n][j] - p * z
- end
-
- # Column modification
-
- 0.upto(n) do |i|
- z = @h[i][n-1]
- @h[i][n-1] = q * z + p * @h[i][n]
- @h[i][n] = q * @h[i][n] - p * z
- end
-
- # Accumulate transformations
-
- low.upto(high) do |i|
- z = @v[i][n-1]
- @v[i][n-1] = q * z + p * @v[i][n]
- @v[i][n] = q * @v[i][n] - p * z
- end
-
- # Complex pair
-
- else
- @d[n-1] = x + p
- @d[n] = x + p
- @e[n-1] = z
- @e[n] = -z
- end
- n -= 2
- iter = 0
-
- # No convergence yet
-
- else
-
- # Form shift
-
- x = @h[n][n]
- y = 0.0
- w = 0.0
- if (l < n)
- y = @h[n-1][n-1]
- w = @h[n][n-1] * @h[n-1][n]
- end
-
- # Wilkinson's original ad hoc shift
-
- if (iter == 10)
- exshift += x
- low.upto(n) do |i|
- @h[i][i] -= x
- end
- s = @h[n][n-1].abs + @h[n-1][n-2].abs
- x = y = 0.75 * s
- w = -0.4375 * s * s
- end
-
- # MATLAB's new ad hoc shift
-
- if (iter == 30)
- s = (y - x) / 2.0
- s *= s + w
- if (s > 0)
- s = Math.sqrt(s)
- if (y < x)
- s = -s
- end
- s = x - w / ((y - x) / 2.0 + s)
- low.upto(n) do |i|
- @h[i][i] -= s
- end
- exshift += s
- x = y = w = 0.964
- end
- end
-
- iter = iter + 1 # (Could check iteration count here.)
-
- # Look for two consecutive small sub-diagonal elements
-
- m = n-2
- while (m >= l) do
- z = @h[m][m]
- r = x - z
- s = y - z
- p = (r * s - w) / @h[m+1][m] + @h[m][m+1]
- q = @h[m+1][m+1] - z - r - s
- r = @h[m+2][m+1]
- s = p.abs + q.abs + r.abs
- p /= s
- q /= s
- r /= s
- if (m == l)
- break
- end
- if (@h[m][m-1].abs * (q.abs + r.abs) <
- eps * (p.abs * (@h[m-1][m-1].abs + z.abs +
- @h[m+1][m+1].abs)))
- break
- end
- m-=1
- end
-
- (m+2).upto(n) do |i|
- @h[i][i-2] = 0.0
- if (i > m+2)
- @h[i][i-3] = 0.0
- end
- end
-
- # Double QR step involving rows l:n and columns m:n
-
- m.upto(n-1) do |k|
- notlast = (k != n-1)
- if (k != m)
- p = @h[k][k-1]
- q = @h[k+1][k-1]
- r = (notlast ? @h[k+2][k-1] : 0.0)
- x = p.abs + q.abs + r.abs
- next if x == 0
- p /= x
- q /= x
- r /= x
- end
- s = Math.sqrt(p * p + q * q + r * r)
- if (p < 0)
- s = -s
- end
- if (s != 0)
- if (k != m)
- @h[k][k-1] = -s * x
- elsif (l != m)
- @h[k][k-1] = -@h[k][k-1]
- end
- p += s
- x = p / s
- y = q / s
- z = r / s
- q /= p
- r /= p
-
- # Row modification
-
- k.upto(nn-1) do |j|
- p = @h[k][j] + q * @h[k+1][j]
- if (notlast)
- p += r * @h[k+2][j]
- @h[k+2][j] = @h[k+2][j] - p * z
- end
- @h[k][j] = @h[k][j] - p * x
- @h[k+1][j] = @h[k+1][j] - p * y
- end
-
- # Column modification
-
- 0.upto([n, k+3].min) do |i|
- p = x * @h[i][k] + y * @h[i][k+1]
- if (notlast)
- p += z * @h[i][k+2]
- @h[i][k+2] = @h[i][k+2] - p * r
- end
- @h[i][k] = @h[i][k] - p
- @h[i][k+1] = @h[i][k+1] - p * q
- end
-
- # Accumulate transformations
-
- low.upto(high) do |i|
- p = x * @v[i][k] + y * @v[i][k+1]
- if (notlast)
- p += z * @v[i][k+2]
- @v[i][k+2] = @v[i][k+2] - p * r
- end
- @v[i][k] = @v[i][k] - p
- @v[i][k+1] = @v[i][k+1] - p * q
- end
- end # (s != 0)
- end # k loop
- end # check convergence
- end # while (n >= low)
-
- # Backsubstitute to find vectors of upper triangular form
-
- if (norm == 0.0)
- return
- end
-
- (nn-1).downto(0) do |k|
- p = @d[k]
- q = @e[k]
-
- # Real vector
-
- if (q == 0)
- l = k
- @h[k][k] = 1.0
- (k-1).downto(0) do |i|
- w = @h[i][i] - p
- r = 0.0
- l.upto(k) do |j|
- r += @h[i][j] * @h[j][k]
- end
- if (@e[i] < 0.0)
- z = w
- s = r
- else
- l = i
- if (@e[i] == 0.0)
- if (w != 0.0)
- @h[i][k] = -r / w
- else
- @h[i][k] = -r / (eps * norm)
- end
-
- # Solve real equations
-
- else
- x = @h[i][i+1]
- y = @h[i+1][i]
- q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i]
- t = (x * s - z * r) / q
- @h[i][k] = t
- if (x.abs > z.abs)
- @h[i+1][k] = (-r - w * t) / x
- else
- @h[i+1][k] = (-s - y * t) / z
- end
- end
-
- # Overflow control
-
- t = @h[i][k].abs
- if ((eps * t) * t > 1)
- i.upto(k) do |j|
- @h[j][k] = @h[j][k] / t
- end
- end
- end
- end
-
- # Complex vector
-
- elsif (q < 0)
- l = n-1
-
- # Last vector component imaginary so matrix is triangular
-
- if (@h[n][n-1].abs > @h[n-1][n].abs)
- @h[n-1][n-1] = q / @h[n][n-1]
- @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1]
- else
- cdivr, cdivi = cdiv(0.0, -@h[n-1][n], @h[n-1][n-1]-p, q)
- @h[n-1][n-1] = cdivr
- @h[n-1][n] = cdivi
- end
- @h[n][n-1] = 0.0
- @h[n][n] = 1.0
- (n-2).downto(0) do |i|
- ra = 0.0
- sa = 0.0
- l.upto(n) do |j|
- ra = ra + @h[i][j] * @h[j][n-1]
- sa = sa + @h[i][j] * @h[j][n]
- end
- w = @h[i][i] - p
-
- if (@e[i] < 0.0)
- z = w
- r = ra
- s = sa
- else
- l = i
- if (@e[i] == 0)
- cdivr, cdivi = cdiv(-ra, -sa, w, q)
- @h[i][n-1] = cdivr
- @h[i][n] = cdivi
- else
-
- # Solve complex equations
-
- x = @h[i][i+1]
- y = @h[i+1][i]
- vr = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - q * q
- vi = (@d[i] - p) * 2.0 * q
- if (vr == 0.0 && vi == 0.0)
- vr = eps * norm * (w.abs + q.abs +
- x.abs + y.abs + z.abs)
- end
- cdivr, cdivi = cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi)
- @h[i][n-1] = cdivr
- @h[i][n] = cdivi
- if (x.abs > (z.abs + q.abs))
- @h[i+1][n-1] = (-ra - w * @h[i][n-1] + q * @h[i][n]) / x
- @h[i+1][n] = (-sa - w * @h[i][n] - q * @h[i][n-1]) / x
- else
- cdivr, cdivi = cdiv(-r-y*@h[i][n-1], -s-y*@h[i][n], z, q)
- @h[i+1][n-1] = cdivr
- @h[i+1][n] = cdivi
- end
- end
-
- # Overflow control
-
- t = [@h[i][n-1].abs, @h[i][n].abs].max
- if ((eps * t) * t > 1)
- i.upto(n) do |j|
- @h[j][n-1] = @h[j][n-1] / t
- @h[j][n] = @h[j][n] / t
- end
- end
- end
- end
- end
- end
-
- # Vectors of isolated roots
-
- nn.times do |i|
- if (i < low || i > high)
- i.upto(nn-1) do |j|
- @v[i][j] = @h[i][j]
- end
- end
- end
-
- # Back transformation to get eigenvectors of original matrix
-
- (nn-1).downto(low) do |j|
- low.upto(high) do |i|
- z = 0.0
- low.upto([j, high].min) do |k|
- z += @v[i][k] * @h[k][j]
- end
- @v[i][j] = z
- end
- end
- end
-
- end
-end
diff --git a/lib/matrix/lup_decomposition.rb b/lib/matrix/lup_decomposition.rb
deleted file mode 100644
index e37def75f6..0000000000
--- a/lib/matrix/lup_decomposition.rb
+++ /dev/null
@@ -1,219 +0,0 @@
-# frozen_string_literal: false
-class Matrix
- # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
-
- #
- # For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
- # unit lower triangular matrix L, an n-by-n upper triangular matrix U,
- # and a m-by-m permutation matrix P so that L*U = P*A.
- # If m < n, then L is m-by-m and U is m-by-n.
- #
- # The LUP decomposition with pivoting always exists, even if the matrix is
- # singular, so the constructor will never fail. The primary use of the
- # LU decomposition is in the solution of square systems of simultaneous
- # linear equations. This will fail if singular? returns true.
- #
-
- class LUPDecomposition
- # Returns the lower triangular factor +L+
-
- include Matrix::ConversionHelper
-
- def l
- Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j|
- if (i > j)
- @lu[i][j]
- elsif (i == j)
- 1
- else
- 0
- end
- end
- end
-
- # Returns the upper triangular factor +U+
-
- def u
- Matrix.build([@column_count, @row_count].min, @column_count) do |i, j|
- if (i <= j)
- @lu[i][j]
- else
- 0
- end
- end
- end
-
- # Returns the permutation matrix +P+
-
- def p
- rows = Array.new(@row_count){Array.new(@row_count, 0)}
- @pivots.each_with_index{|p, i| rows[i][p] = 1}
- Matrix.send :new, rows, @row_count
- end
-
- # Returns +L+, +U+, +P+ in an array
-
- def to_ary
- [l, u, p]
- end
- alias_method :to_a, :to_ary
-
- # Returns the pivoting indices
-
- attr_reader :pivots
-
- # Returns +true+ if +U+, and hence +A+, is singular.
-
- def singular?
- @column_count.times do |j|
- if (@lu[j][j] == 0)
- return true
- end
- end
- false
- end
-
- # Returns the determinant of +A+, calculated efficiently
- # from the factorization.
-
- def det
- if (@row_count != @column_count)
- raise Matrix::ErrDimensionMismatch
- end
- d = @pivot_sign
- @column_count.times do |j|
- d *= @lu[j][j]
- end
- d
- end
- alias_method :determinant, :det
-
- # Returns +m+ so that <tt>A*m = b</tt>,
- # or equivalently so that <tt>L*U*m = P*b</tt>
- # +b+ can be a Matrix or a Vector
-
- def solve b
- if (singular?)
- raise Matrix::ErrNotRegular, "Matrix is singular."
- end
- if b.is_a? Matrix
- if (b.row_count != @row_count)
- raise Matrix::ErrDimensionMismatch
- end
-
- # Copy right hand side with pivoting
- nx = b.column_count
- m = @pivots.map{|row| b.row(row).to_a}
-
- # Solve L*Y = P*b
- @column_count.times do |k|
- (k+1).upto(@column_count-1) do |i|
- nx.times do |j|
- m[i][j] -= m[k][j]*@lu[i][k]
- end
- end
- end
- # Solve U*m = Y
- (@column_count-1).downto(0) do |k|
- nx.times do |j|
- m[k][j] = m[k][j].quo(@lu[k][k])
- end
- k.times do |i|
- nx.times do |j|
- m[i][j] -= m[k][j]*@lu[i][k]
- end
- end
- end
- Matrix.send :new, m, nx
- else # same algorithm, specialized for simpler case of a vector
- b = convert_to_array(b)
- if (b.size != @row_count)
- raise Matrix::ErrDimensionMismatch
- end
-
- # Copy right hand side with pivoting
- m = b.values_at(*@pivots)
-
- # Solve L*Y = P*b
- @column_count.times do |k|
- (k+1).upto(@column_count-1) do |i|
- m[i] -= m[k]*@lu[i][k]
- end
- end
- # Solve U*m = Y
- (@column_count-1).downto(0) do |k|
- m[k] = m[k].quo(@lu[k][k])
- k.times do |i|
- m[i] -= m[k]*@lu[i][k]
- end
- end
- Vector.elements(m, false)
- end
- end
-
- def initialize a
- raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
- # Use a "left-looking", dot-product, Crout/Doolittle algorithm.
- @lu = a.to_a
- @row_count = a.row_count
- @column_count = a.column_count
- @pivots = Array.new(@row_count)
- @row_count.times do |i|
- @pivots[i] = i
- end
- @pivot_sign = 1
- lu_col_j = Array.new(@row_count)
-
- # Outer loop.
-
- @column_count.times do |j|
-
- # Make a copy of the j-th column to localize references.
-
- @row_count.times do |i|
- lu_col_j[i] = @lu[i][j]
- end
-
- # Apply previous transformations.
-
- @row_count.times do |i|
- lu_row_i = @lu[i]
-
- # Most of the time is spent in the following dot product.
-
- kmax = [i, j].min
- s = 0
- kmax.times do |k|
- s += lu_row_i[k]*lu_col_j[k]
- end
-
- lu_row_i[j] = lu_col_j[i] -= s
- end
-
- # Find pivot and exchange if necessary.
-
- p = j
- (j+1).upto(@row_count-1) do |i|
- if (lu_col_j[i].abs > lu_col_j[p].abs)
- p = i
- end
- end
- if (p != j)
- @column_count.times do |k|
- t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
- end
- k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
- @pivot_sign = -@pivot_sign
- end
-
- # Compute multipliers.
-
- if (j < @row_count && @lu[j][j] != 0)
- (j+1).upto(@row_count-1) do |i|
- @lu[i][j] = @lu[i][j].quo(@lu[j][j])
- end
- end
- end
- end
- end
-end
diff --git a/lib/matrix/matrix.gemspec b/lib/matrix/matrix.gemspec
deleted file mode 100644
index 6f68cbd816..0000000000
--- a/lib/matrix/matrix.gemspec
+++ /dev/null
@@ -1,26 +0,0 @@
-# frozen_string_literal: true
-
-begin
- require_relative "lib/matrix/version"
-rescue LoadError
- # for Ruby core repository
- require_relative "version"
-end
-
-Gem::Specification.new do |spec|
- spec.name = "matrix"
- spec.version = Matrix::VERSION
- spec.authors = ["Marc-Andre Lafortune"]
- spec.email = ["ruby-core@marc-andre.ca"]
-
- spec.summary = %q{An implementation of Matrix and Vector classes.}
- spec.description = %q{An implementation of Matrix and Vector classes.}
- spec.homepage = "https://github.com/ruby/matrix"
- spec.licenses = ["Ruby", "BSD-2-Clause"]
- spec.required_ruby_version = ">= 2.5.0"
-
- spec.files = [".gitignore", "Gemfile", "LICENSE.txt", "README.md", "Rakefile", "bin/console", "bin/setup", "lib/matrix.rb", "lib/matrix/eigenvalue_decomposition.rb", "lib/matrix/lup_decomposition.rb", "lib/matrix/version.rb", "matrix.gemspec"]
- spec.bindir = "exe"
- spec.executables = spec.files.grep(%r{^exe/}) { |f| File.basename(f) }
- spec.require_paths = ["lib"]
-end
diff --git a/lib/matrix/version.rb b/lib/matrix/version.rb
deleted file mode 100644
index 82b835f038..0000000000
--- a/lib/matrix/version.rb
+++ /dev/null
@@ -1,5 +0,0 @@
-# frozen_string_literal: true
-
-class Matrix
- VERSION = "0.4.1"
-end
diff --git a/misc/expand_tabs.rb b/misc/expand_tabs.rb
index 284999a2e9..630e012996 100755
--- a/misc/expand_tabs.rb
+++ b/misc/expand_tabs.rb
@@ -78,7 +78,6 @@ DEFAULT_GEM_LIBS = %w[
ipaddr
irb
logger
- matrix
mutex_m
net-ftp
net-http
diff --git a/test/matrix/test_matrix.rb b/test/matrix/test_matrix.rb
deleted file mode 100644
index 25ad7b6b95..0000000000
--- a/test/matrix/test_matrix.rb
+++ /dev/null
@@ -1,888 +0,0 @@
-# frozen_string_literal: false
-require 'test/unit'
-require 'matrix'
-
-class SubMatrix < Matrix
-end
-
-class TestMatrix < Test::Unit::TestCase
- def setup
- @m1 = Matrix[[1,2,3], [4,5,6]]
- @m2 = Matrix[[1,2,3], [4,5,6]]
- @m3 = @m1.clone
- @m4 = Matrix[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
- @n1 = Matrix[[2,3,4], [5,6,7]]
- @c1 = Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
- @e1 = Matrix.empty(2,0)
- @e2 = Matrix.empty(0,3)
- @a3 = Matrix[[4, 1, -3], [0, 3, 7], [11, -4, 2]]
- @a5 = Matrix[[2, 0, 9, 3, 9], [8, 7, 0, 1, 9], [7, 5, 6, 6, 5], [0, 7, 8, 3, 0], [7, 8, 2, 3, 1]]
- @b3 = Matrix[[-7, 7, -10], [9, -3, -2], [-1, 3, 9]]
- @rot = Matrix[[0, -1, 0], [1, 0, 0], [0, 0, -1]]
- end
-
- def test_matrix
- assert_equal(1, @m1[0, 0])
- assert_equal(2, @m1[0, 1])
- assert_equal(3, @m1[0, 2])
- assert_equal(4, @m1[1, 0])
- assert_equal(5, @m1[1, 1])
- assert_equal(6, @m1[1, 2])
- end
-
- def test_identity
- assert_same @m1, @m1
- assert_not_same @m1, @m2
- assert_not_same @m1, @m3
- assert_not_same @m1, @m4
- assert_not_same @m1, @n1
- end
-
- def test_equality
- assert_equal @m1, @m1
- assert_equal @m1, @m2
- assert_equal @m1, @m3
- assert_equal @m1, @m4
- assert_not_equal @m1, @n1
- end
-
- def test_hash_equality
- assert @m1.eql?(@m1)
- assert @m1.eql?(@m2)
- assert @m1.eql?(@m3)
- assert !@m1.eql?(@m4)
- assert !@m1.eql?(@n1)
-
- hash = { @m1 => :value }
- assert hash.key?(@m1)
- assert hash.key?(@m2)
- assert hash.key?(@m3)
- assert !hash.key?(@m4)
- assert !hash.key?(@n1)
- end
-
- def test_hash
- assert_equal @m1.hash, @m1.hash
- assert_equal @m1.hash, @m2.hash
- assert_equal @m1.hash, @m3.hash
- end
-
- def test_uplus
- assert_equal(@m1, +@m1)
- end
-
- def test_negate
- assert_equal(Matrix[[-1, -2, -3], [-4, -5, -6]], -@m1)
- assert_equal(@m1, -(-@m1))
- end
-
- def test_rank
- [
- [[0]],
- [[0], [0]],
- [[0, 0], [0, 0]],
- [[0, 0], [0, 0], [0, 0]],
- [[0, 0, 0]],
- [[0, 0, 0], [0, 0, 0]],
- [[0, 0, 0], [0, 0, 0], [0, 0, 0]],
- [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]],
- ].each do |rows|
- assert_equal 0, Matrix[*rows].rank
- end
-
- [
- [[1], [0]],
- [[1, 0], [0, 0]],
- [[1, 0], [1, 0]],
- [[0, 0], [1, 0]],
- [[1, 0], [0, 0], [0, 0]],
- [[0, 0], [1, 0], [0, 0]],
- [[0, 0], [0, 0], [1, 0]],
- [[1, 0], [1, 0], [0, 0]],
- [[0, 0], [1, 0], [1, 0]],
- [[1, 0], [1, 0], [1, 0]],
- [[1, 0, 0]],
- [[1, 0, 0], [0, 0, 0]],
- [[0, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [0, 0, 0]],
- [[0, 0, 0], [1, 0, 0], [0, 0, 0]],
- [[0, 0, 0], [0, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0], [0, 0, 0]],
- [[0, 0, 0], [1, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]],
- [[1, 0, 0], [1, 0, 0], [0, 0, 0], [0, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [1, 0, 0], [0, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [0, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0], [1, 0, 0], [0, 0, 0]],
- [[1, 0, 0], [0, 0, 0], [1, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0], [0, 0, 0], [1, 0, 0]],
- [[1, 0, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0]],
-
- [[1]],
- [[1], [1]],
- [[1, 1]],
- [[1, 1], [1, 1]],
- [[1, 1], [1, 1], [1, 1]],
- [[1, 1, 1]],
- [[1, 1, 1], [1, 1, 1]],
- [[1, 1, 1], [1, 1, 1], [1, 1, 1]],
- [[1, 1, 1], [1, 1, 1], [1, 1, 1], [1, 1, 1]],
- ].each do |rows|
- matrix = Matrix[*rows]
- assert_equal 1, matrix.rank
- assert_equal 1, matrix.transpose.rank
- end
-
- [
- [[1, 0], [0, 1]],
- [[1, 0], [0, 1], [0, 0]],
- [[1, 0], [0, 1], [0, 1]],
- [[1, 0], [0, 1], [1, 1]],
- [[1, 0, 0], [0, 1, 0]],
- [[1, 0, 0], [0, 0, 1]],
- [[1, 0, 0], [0, 1, 0], [0, 0, 0]],
- [[1, 0, 0], [0, 0, 1], [0, 0, 0]],
-
- [[1, 0, 0], [0, 0, 0], [0, 1, 0]],
- [[1, 0, 0], [0, 0, 0], [0, 0, 1]],
-
- [[1, 0], [1, 1]],
- [[1, 2], [1, 1]],
- [[1, 2], [0, 1], [1, 1]],
- ].each do |rows|
- m = Matrix[*rows]
- assert_equal 2, m.rank
- assert_equal 2, m.transpose.rank
- end
-
- [
- [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
- [[1, 1, 0], [0, 1, 1], [1, 0, 1]],
- [[1, 1, 0], [0, 1, 1], [1, 0, 1]],
- [[1, 1, 0], [0, 1, 1], [1, 0, 1], [0, 0, 0]],
- [[1, 1, 0], [0, 1, 1], [1, 0, 1], [1, 1, 1]],
- [[1, 1, 1], [1, 1, 2], [1, 3, 1], [4, 1, 1]],
- ].each do |rows|
- m = Matrix[*rows]
- assert_equal 3, m.rank
- assert_equal 3, m.transpose.rank
- end
- end
-
- def test_inverse
- assert_equal(Matrix.empty(0, 0), Matrix.empty.inverse)
- assert_equal(Matrix[[-1, 1], [0, -1]], Matrix[[-1, -1], [0, -1]].inverse)
- assert_raise(ExceptionForMatrix::ErrDimensionMismatch) { @m1.inverse }
- end
-
- def test_determinant
- assert_equal(0, Matrix[[0,0],[0,0]].determinant)
- assert_equal(45, Matrix[[7,6], [3,9]].determinant)
- assert_equal(-18, Matrix[[2,0,1],[0,-2,2],[1,2,3]].determinant)
- assert_equal(-7, Matrix[[0,0,1],[0,7,6],[1,3,9]].determinant)
- assert_equal(42, Matrix[[7,0,1,0,12],[8,1,1,9,1],[4,0,0,-7,17],[-1,0,0,-4,8],[10,1,1,8,6]].determinant)
- end
-
- def test_new_matrix
- assert_raise(TypeError) { Matrix[Object.new] }
- o = Object.new
- def o.to_ary; [1,2,3]; end
- assert_equal(@m1, Matrix[o, [4,5,6]])
- end
-
- def test_round
- a = Matrix[[1.0111, 2.32320, 3.04343], [4.81, 5.0, 6.997]]
- b = Matrix[[1.01, 2.32, 3.04], [4.81, 5.0, 7.0]]
- assert_equal(a.round(2), b)
- end
-
- def test_rows
- assert_equal(@m1, Matrix.rows([[1, 2, 3], [4, 5, 6]]))
- end
-
- def test_rows_copy
- rows1 = [[1], [1]]
- rows2 = [[1], [1]]
-
- m1 = Matrix.rows(rows1, copy = false)
- m2 = Matrix.rows(rows2, copy = true)
-
- rows1.uniq!
- rows2.uniq!
-
- assert_equal([[1]], m1.to_a)
- assert_equal([[1], [1]], m2.to_a)
- end
-
- def test_to_matrix
- assert @m1.equal? @m1.to_matrix
- end
-
- def test_columns
- assert_equal(@m1, Matrix.columns([[1, 4], [2, 5], [3, 6]]))
- end
-
- def test_diagonal
- assert_equal(Matrix.empty(0, 0), Matrix.diagonal( ))
- assert_equal(Matrix[[3,0,0],[0,2,0],[0,0,1]], Matrix.diagonal(3, 2, 1))
- assert_equal(Matrix[[4,0,0,0],[0,3,0,0],[0,0,2,0],[0,0,0,1]], Matrix.diagonal(4, 3, 2, 1))
- end
-
- def test_scalar
- assert_equal(Matrix.empty(0, 0), Matrix.scalar(0, 1))
- assert_equal(Matrix[[2,0,0],[0,2,0],[0,0,2]], Matrix.scalar(3, 2))
- assert_equal(Matrix[[2,0,0,0],[0,2,0,0],[0,0,2,0],[0,0,0,2]], Matrix.scalar(4, 2))
- end
-
- def test_identity2
- assert_equal(Matrix[[1,0,0],[0,1,0],[0,0,1]], Matrix.identity(3))
- assert_equal(Matrix[[1,0,0],[0,1,0],[0,0,1]], Matrix.unit(3))
- assert_equal(Matrix[[1,0,0],[0,1,0],[0,0,1]], Matrix.I(3))
- assert_equal(Matrix[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]], Matrix.identity(4))
- end
-
- def test_zero
- assert_equal(Matrix[[0,0,0],[0,0,0],[0,0,0]], Matrix.zero(3))
- assert_equal(Matrix[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], Matrix.zero(4))
- assert_equal(Matrix[[0]], Matrix.zero(1))
- end
-
- def test_row_vector
- assert_equal(Matrix[[1,2,3,4]], Matrix.row_vector([1,2,3,4]))
- end
-
- def test_column_vector
- assert_equal(Matrix[[1],[2],[3],[4]], Matrix.column_vector([1,2,3,4]))
- end
-
- def test_empty
- m = Matrix.empty(2, 0)
- assert_equal(Matrix[ [], [] ], m)
- n = Matrix.empty(0, 3)
- assert_equal(Matrix.columns([ [], [], [] ]), n)
- assert_equal(Matrix[[0, 0, 0], [0, 0, 0]], m * n)
- end
-
- def test_row
- assert_equal(Vector[1, 2, 3], @m1.row(0))
- assert_equal(Vector[4, 5, 6], @m1.row(1))
- a = []; @m1.row(0) {|x| a << x }
- assert_equal([1, 2, 3], a)
- end
-
- def test_column
- assert_equal(Vector[1, 4], @m1.column(0))
- assert_equal(Vector[2, 5], @m1.column(1))
- assert_equal(Vector[3, 6], @m1.column(2))
- a = []; @m1.column(0) {|x| a << x }
- assert_equal([1, 4], a)
- end
-
- def test_collect
- m1 = Matrix.zero(2,2)
- m2 = Matrix.build(3,4){|row, col| 1}
-
- assert_equal(Matrix[[5, 5, 5, 5], [5, 5, 5, 5], [5, 5, 5, 5]], m2.collect{|e| e * 5})
- assert_equal(Matrix[[7, 0],[0, 7]], m1.collect(:diagonal){|e| e + 7})
- assert_equal(Matrix[[0, 5],[5, 0]], m1.collect(:off_diagonal){|e| e + 5})
- assert_equal(Matrix[[8, 1, 1, 1], [8, 8, 1, 1], [8, 8, 8, 1]], m2.collect(:lower){|e| e + 7})
- assert_equal(Matrix[[1, 1, 1, 1], [-11, 1, 1, 1], [-11, -11, 1, 1]], m2.collect(:strict_lower){|e| e - 12})
- assert_equal(Matrix[[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]], m2.collect(:strict_upper){|e| e ** 2})
- assert_equal(Matrix[[-1, -1, -1, -1], [1, -1, -1, -1], [1, 1, -1, -1]], m2.collect(:upper){|e| -e})
- assert_raise(ArgumentError) {m1.collect(:test){|e| e + 7}}
- assert_not_equal(m2, m2.collect {|e| e * 2 })
- end
-
- def test_minor
- assert_equal(Matrix[[1, 2], [4, 5]], @m1.minor(0..1, 0..1))
- assert_equal(Matrix[[2], [5]], @m1.minor(0..1, 1..1))
- assert_equal(Matrix[[4, 5]], @m1.minor(1..1, 0..1))
- assert_equal(Matrix[[1, 2], [4, 5]], @m1.minor(0, 2, 0, 2))
- assert_equal(Matrix[[4, 5]], @m1.minor(1, 1, 0, 2))
- assert_equal(Matrix[[2], [5]], @m1.minor(0, 2, 1, 1))
- assert_raise(ArgumentError) { @m1.minor(0) }
- end
-
- def test_first_minor
- assert_equal(Matrix.empty(0, 0), Matrix[[1]].first_minor(0, 0))
- assert_equal(Matrix.empty(0, 2), Matrix[[1, 4, 2]].first_minor(0, 1))
- assert_equal(Matrix[[1, 3]], @m1.first_minor(1, 1))
- assert_equal(Matrix[[4, 6]], @m1.first_minor(0, 1))
- assert_equal(Matrix[[1, 2]], @m1.first_minor(1, 2))
- assert_raise(RuntimeError) { Matrix.empty(0, 0).first_minor(0, 0) }
- assert_raise(ArgumentError) { @m1.first_minor(4, 0) }
- assert_raise(ArgumentError) { @m1.first_minor(0, -1) }
- assert_raise(ArgumentError) { @m1.first_minor(-1, 4) }
- end
-
- def test_cofactor
- assert_equal(1, Matrix[[1]].cofactor(0, 0))
- assert_equal(9, Matrix[[7,6],[3,9]].cofactor(0, 0))
- assert_equal(0, Matrix[[0,0],[0,0]].cofactor(0, 0))
- assert_equal(3, Matrix[[0,0,1],[0,7,6],[1,3,9]].cofactor(1, 0))
- assert_equal(-21, Matrix[[7,0,1,0,12],[8,1,1,9,1],[4,0,0,-7,17],[-1,0,0,-4,8],[10,1,1,8,6]].cofactor(2, 3))
- assert_raise(RuntimeError) { Matrix.empty(0, 0).cofactor(0, 0) }
- assert_raise(ArgumentError) { Matrix[[0,0],[0,0]].cofactor(-1, 4) }
- assert_raise(ExceptionForMatrix::ErrDimensionMismatch) { Matrix[[2,0,1],[0,-2,2]].cofactor(0, 0) }
- end
-
- def test_adjugate
- assert_equal(Matrix.empty, Matrix.empty.adjugate)
- assert_equal(Matrix[[1]], Matrix[[5]].adjugate)
- assert_equal(Matrix[[9,-6],[-3,7]], Matrix[[7,6],[3,9]].adjugate)
- assert_equal(Matrix[[45,3,-7],[6,-1,0],[-7,0,0]], Matrix[[0,0,1],[0,7,6],[1,3,9]].adjugate)
- assert_equal(Matrix.identity(5), (@a5.adjugate * @a5) / @a5.det)
- assert_equal(Matrix.I(3), Matrix.I(3).adjugate)
- assert_equal((@a3 * @b3).adjugate, @b3.adjugate * @a3.adjugate)
- assert_equal(4**(@a3.row_count-1) * @a3.adjugate, (4 * @a3).adjugate)
- assert_raise(ExceptionForMatrix::ErrDimensionMismatch) { @m1.adjugate }
- end
-
- def test_laplace_expansion
- assert_equal(1, Matrix[[1]].laplace_expansion(row: 0))
- assert_equal(45, Matrix[[7,6], [3,9]].laplace_expansion(row: 1))
- assert_equal(0, Matrix[[0,0],[0,0]].laplace_expansion(column: 0))
- assert_equal(-7, Matrix[[0,0,1],[0,7,6],[1,3,9]].laplace_expansion(column: 2))
-
- assert_equal(Vector[3, -2], Matrix[[Vector[1, 0], Vector[0, 1]], [2, 3]].laplace_expansion(row: 0))
-
- assert_raise(ExceptionForMatrix::ErrDimensionMismatch) { @m1.laplace_expansion(row: 1) }
- assert_raise(ArgumentError) { Matrix[[7,6], [3,9]].laplace_expansion() }
- assert_raise(ArgumentError) { Matrix[[7,6], [3,9]].laplace_expansion(foo: 1) }
- assert_raise(ArgumentError) { Matrix[[7,6], [3,9]].laplace_expansion(row: 1, column: 1) }
- assert_raise(ArgumentError) { Matrix[[7,6], [3,9]].laplace_expansion(row: 2) }
- assert_raise(ArgumentError) { Matrix[[0,0,1],[0,7,6],[1,3,9]].laplace_expansion(column: -1) }
-
- assert_raise(RuntimeError) { Matrix.empty(0, 0).laplace_expansion(row: 0) }
- end
-
- def test_regular?
- assert(Matrix[[1, 0], [0, 1]].regular?)
- assert(Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].regular?)
- assert(!Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].regular?)
- end
-
- def test_singular?
- assert(!Matrix[[1, 0], [0, 1]].singular?)
- assert(!Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].singular?)
- assert(Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].singular?)
- end
-
- def test_square?
- assert(Matrix[[1, 0], [0, 1]].square?)
- assert(Matrix[[1, 0, 0], [0, 1, 0], [0, 0, 1]].square?)
- assert(Matrix[[1, 0, 0], [0, 0, 1], [0, 0, 1]].square?)
- assert(!Matrix[[1, 0, 0], [0, 1, 0]].square?)
- end
-
- def test_mul
- assert_equal(Matrix[[2,4],[6,8]], Matrix[[2,4],[6,8]] * Matrix.I(2))
- assert_equal(Matrix[[4,8],[12,16]], Matrix[[2,4],[6,8]] * 2)
- assert_equal(Matrix[[4,8],[12,16]], 2 * Matrix[[2,4],[6,8]])
- assert_equal(Matrix[[14,32],[32,77]], @m1 * @m1.transpose)
- assert_equal(Matrix[[17,22,27],[22,29,36],[27,36,45]], @m1.transpose * @m1)
- assert_equal(Vector[14,32], @m1 * Vector[1,2,3])
- o = Object.new
- def o.coerce(m)
- [m, m.transpose]
- end
- assert_equal(Matrix[[14,32],[32,77]], @m1 * o)
- end
-
- def test_add
- assert_equal(Matrix[[6,0],[-4,12]], Matrix.scalar(2,5) + Matrix[[1,0],[-4,7]])
- assert_equal(Matrix[[3,5,7],[9,11,13]], @m1 + @n1)
- assert_equal(Matrix[[3,5,7],[9,11,13]], @n1 + @m1)
- assert_equal(Matrix[[2],[4],[6]], Matrix[[1],[2],[3]] + Vector[1,2,3])
- assert_raise(Matrix::ErrOperationNotDefined) { @m1 + 1 }
- o = Object.new
- def o.coerce(m)
- [m, m]
- end
- assert_equal(Matrix[[2,4,6],[8,10,12]], @m1 + o)
- end
-
- def test_sub
- assert_equal(Matrix[[4,0],[4,-2]], Matrix.scalar(2,5) - Matrix[[1,0],[-4,7]])
- assert_equal(Matrix[[-1,-1,-1],[-1,-1,-1]], @m1 - @n1)
- assert_equal(Matrix[[1,1,1],[1,1,1]], @n1 - @m1)
- assert_equal(Matrix[[0],[0],[0]], Matrix[[1],[2],[3]] - Vector[1,2,3])
- assert_raise(Matrix::ErrOperationNotDefined) { @m1 - 1 }
- o = Object.new
- def o.coerce(m)
- [m, m]
- end
- assert_equal(Matrix[[0,0,0],[0,0,0]], @m1 - o)
- end
-
- def test_div
- assert_equal(Matrix[[0,1,1],[2,2,3]], @m1 / 2)
- assert_equal(Matrix[[1,1],[1,1]], Matrix[[2,2],[2,2]] / Matrix.scalar(2,2))
- o = Object.new
- def o.coerce(m)
- [m, Matrix.scalar(2,2)]
- end
- assert_equal(Matrix[[1,1],[1,1]], Matrix[[2,2],[2,2]] / o)
- end
-
- def test_hadamard_product
- assert_equal(Matrix[[1,4], [9,16]], Matrix[[1,2], [3,4]].hadamard_product(Matrix[[1,2], [3,4]]))
- assert_equal(Matrix[[2, 6, 12], [20, 30, 42]], @m1.hadamard_product(@n1))
- o = Object.new
- def o.to_matrix
- Matrix[[1, 2, 3], [-1, 0, 1]]
- end
- assert_equal(Matrix[[1, 4, 9], [-4, 0, 6]], @m1.hadamard_product(o))
- e = Matrix.empty(3, 0)
- assert_equal(e, e.hadamard_product(e))
- e = Matrix.empty(0, 3)
- assert_equal(e, e.hadamard_product(e))
- end
-
- def test_exp
- assert_equal(Matrix[[67,96],[48,99]], Matrix[[7,6],[3,9]] ** 2)
- assert_equal(Matrix.I(5), Matrix.I(5) ** -1)
- assert_raise(Matrix::ErrOperationNotDefined) { Matrix.I(5) ** Object.new }
-
- m = Matrix[[0,2],[1,0]]
- exp = 0b11101000
- assert_equal(Matrix.scalar(2, 1 << (exp/2)), m ** exp)
- exp = 0b11101001
- assert_equal(Matrix[[0, 2 << (exp/2)], [1 << (exp/2), 0]], m ** exp)
- end
-
- def test_det
- assert_equal(Matrix.instance_method(:determinant), Matrix.instance_method(:det))
- end
-
- def test_rank2
- assert_equal(2, Matrix[[7,6],[3,9]].rank)
- assert_equal(0, Matrix[[0,0],[0,0]].rank)
- assert_equal(3, Matrix[[0,0,1],[0,7,6],[1,3,9]].rank)
- assert_equal(1, Matrix[[0,1],[0,1],[0,1]].rank)
- assert_equal(2, @m1.rank)
- end
-
- def test_trace
- assert_equal(1+5+9, Matrix[[1,2,3],[4,5,6],[7,8,9]].trace)
- end
-
- def test_transpose
- assert_equal(Matrix[[1,4],[2,5],[3,6]], @m1.transpose)
- end
-
- def test_conjugate
- assert_equal(Matrix[[Complex(1,-2), Complex(0,-1), 0], [1, 2, 3]], @c1.conjugate)
- end
-
- def test_eigensystem
- m = Matrix[[1, 2], [3, 4]]
- v, d, v_inv = m.eigensystem
- assert(d.diagonal?)
- assert_equal(v.inv, v_inv)
- assert_equal((v * d * v_inv).round(5), m)
- end
-
- def test_imaginary
- assert_equal(Matrix[[2, 1, 0], [0, 0, 0]], @c1.imaginary)
- end
-
- def test_lup
- m = Matrix[[1, 2], [3, 4]]
- l, u, p = m.lup
- assert(l.lower_triangular?)
- assert(u.upper_triangular?)
- assert(p.permutation?)
- assert(l * u == p * m)
- assert_equal(m.lup.solve([2, 5]), Vector[1, Rational(1,2)])
- end
-
- def test_real
- assert_equal(Matrix[[1, 0, 0], [1, 2, 3]], @c1.real)
- end
-
- def test_rect
- assert_equal([Matrix[[1, 0, 0], [1, 2, 3]], Matrix[[2, 1, 0], [0, 0, 0]]], @c1.rect)
- end
-
- def test_row_vectors
- assert_equal([Vector[1,2,3], Vector[4,5,6]], @m1.row_vectors)
- end
-
- def test_column_vectors
- assert_equal([Vector[1,4], Vector[2,5], Vector[3,6]], @m1.column_vectors)
- end
-
- def test_to_s
- assert_equal("Matrix[[1, 2, 3], [4, 5, 6]]", @m1.to_s)
- assert_equal("Matrix.empty(0, 0)", Matrix[].to_s)
- assert_equal("Matrix.empty(1, 0)", Matrix[[]].to_s)
- end
-
- def test_inspect
- assert_equal("Matrix[[1, 2, 3], [4, 5, 6]]", @m1.inspect)
- assert_equal("Matrix.empty(0, 0)", Matrix[].inspect)
- assert_equal("Matrix.empty(1, 0)", Matrix[[]].inspect)
- end
-
- def test_scalar_add
- s1 = @m1.coerce(1).first
- assert_equal(Matrix[[1]], (s1 + 0) * Matrix[[1]])
- assert_raise(Matrix::ErrOperationNotDefined) { s1 + Vector[0] }
- assert_raise(Matrix::ErrOperationNotDefined) { s1 + Matrix[[0]] }
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(2, s1 + o)
- end
-
- def test_scalar_sub
- s1 = @m1.coerce(1).first
- assert_equal(Matrix[[1]], (s1 - 0) * Matrix[[1]])
- assert_raise(Matrix::ErrOperationNotDefined) { s1 - Vector[0] }
- assert_raise(Matrix::ErrOperationNotDefined) { s1 - Matrix[[0]] }
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(0, s1 - o)
- end
-
- def test_scalar_mul
- s1 = @m1.coerce(1).first
- assert_equal(Matrix[[1]], (s1 * 1) * Matrix[[1]])
- assert_equal(Vector[2], s1 * Vector[2])
- assert_equal(Matrix[[2]], s1 * Matrix[[2]])
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(1, s1 * o)
- end
-
- def test_scalar_div
- s1 = @m1.coerce(1).first
- assert_equal(Matrix[[1]], (s1 / 1) * Matrix[[1]])
- assert_raise(Matrix::ErrOperationNotDefined) { s1 / Vector[0] }
- assert_equal(Matrix[[Rational(1,2)]], s1 / Matrix[[2]])
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(1, s1 / o)
- end
-
- def test_scalar_pow
- s1 = @m1.coerce(1).first
- assert_equal(Matrix[[1]], (s1 ** 1) * Matrix[[1]])
- assert_raise(Matrix::ErrOperationNotDefined) { s1 ** Vector[0] }
- assert_raise(Matrix::ErrOperationNotImplemented) { s1 ** Matrix[[1]] }
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(1, s1 ** o)
- end
-
- def test_abs
- s1 = @a3.abs
- assert_equal(s1, Matrix[[4, 1, 3], [0, 3, 7], [11, 4, 2]])
- end
-
- def test_hstack
- assert_equal Matrix[[1,2,3,2,3,4,1,2,3], [4,5,6,5,6,7,4,5,6]],
- @m1.hstack(@n1, @m1)
- # Error checking:
- assert_raise(TypeError) { @m1.hstack(42) }
- assert_raise(TypeError) { Matrix.hstack(42, @m1) }
- assert_raise(Matrix::ErrDimensionMismatch) { @m1.hstack(Matrix.identity(3)) }
- assert_raise(Matrix::ErrDimensionMismatch) { @e1.hstack(@e2) }
- # Corner cases:
- assert_equal @m1, @m1.hstack
- assert_equal @e1, @e1.hstack(@e1)
- assert_equal Matrix.empty(0,6), @e2.hstack(@e2)
- assert_equal SubMatrix, SubMatrix.hstack(@e1).class
- # From Vectors:
- assert_equal Matrix[[1, 3],[2, 4]], Matrix.hstack(Vector[1,2], Vector[3, 4])
- end
-
- def test_vstack
- assert_equal Matrix[[1,2,3], [4,5,6], [2,3,4], [5,6,7], [1,2,3], [4,5,6]],
- @m1.vstack(@n1, @m1)
- # Error checking:
- assert_raise(TypeError) { @m1.vstack(42) }
- assert_raise(TypeError) { Matrix.vstack(42, @m1) }
- assert_raise(Matrix::ErrDimensionMismatch) { @m1.vstack(Matrix.identity(2)) }
- assert_raise(Matrix::ErrDimensionMismatch) { @e1.vstack(@e2) }
- # Corner cases:
- assert_equal @m1, @m1.vstack
- assert_equal Matrix.empty(4,0), @e1.vstack(@e1)
- assert_equal @e2, @e2.vstack(@e2)
- assert_equal SubMatrix, SubMatrix.vstack(@e1).class
- # From Vectors:
- assert_equal Matrix[[1],[2],[3]], Matrix.vstack(Vector[1,2], Vector[3])
- end
-
- def test_combine
- x = Matrix[[6, 6], [4, 4]]
- y = Matrix[[1, 2], [3, 4]]
- assert_equal Matrix[[5, 4], [1, 0]], Matrix.combine(x, y) {|a, b| a - b}
- assert_equal Matrix[[5, 4], [1, 0]], x.combine(y) {|a, b| a - b}
- # Without block
- assert_equal Matrix[[5, 4], [1, 0]], Matrix.combine(x, y).each {|a, b| a - b}
- # With vectors
- assert_equal Matrix[[111], [222]], Matrix.combine(Matrix[[1], [2]], Vector[10,20], Vector[100,200], &:sum)
- # Basic checks
- assert_raise(Matrix::ErrDimensionMismatch) { @m1.combine(x) { raise } }
- # Edge cases
- assert_equal Matrix.empty, Matrix.combine{ raise }
- assert_equal Matrix.empty(3,0), Matrix.combine(Matrix.empty(3,0), Matrix.empty(3,0)) { raise }
- assert_equal Matrix.empty(0,3), Matrix.combine(Matrix.empty(0,3), Matrix.empty(0,3)) { raise }
- end
-
- def test_set_element
- src = Matrix[
- [1, 2, 3, 4],
- [5, 6, 7, 8],
- [9, 10, 11, 12],
- ]
- rows = {
- range: [1..2, 1...3, 1..-1, -2..2, 1.., 1..., -2.., -2...],
- int: [2, -1],
- invalid: [-4, 4, -4..2, 2..-4, 0...0, 2..0, -4..],
- }
- columns = {
- range: [2..3, 2...4, 2..-1, -2..3, 2.., 2..., -2..., -2..],
- int: [3, -1],
- invalid: [-5, 5, -5..2, 2..-5, 0...0, -5..],
- }
- values = {
- element: 42,
- matrix: Matrix[[20, 21], [22, 23]],
- vector: Vector[30, 31],
- row: Matrix[[60, 61]],
- column: Matrix[[50], [51]],
- mismatched_matrix: Matrix.identity(3),
- mismatched_vector: Vector[0, 1, 2, 3],
- }
- solutions = {
- [:int, :int] => {
- element: Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 42]],
- },
- [:range , :int] => {
- element: Matrix[[1, 2, 3, 4], [5, 6, 7, 42], [9, 10, 11, 42]],
- column: Matrix[[1, 2, 3, 4], [5, 6, 7, 50], [9, 10, 11, 51]],
- vector: Matrix[[1, 2, 3, 4], [5, 6, 7, 30], [9, 10, 11, 31]],
- },
- [:int, :range] => {
- element: Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 42, 42]],
- row: Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 60, 61]],
- vector: Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 30, 31]],
- },
- [:range , :range] => {
- element: Matrix[[1, 2, 3, 4], [5, 6, 42, 42], [9, 10, 42, 42]],
- matrix: Matrix[[1, 2, 3, 4], [5, 6, 20, 21], [9, 10, 22, 23]],
- },
- }
- solutions.default = Hash.new(IndexError)
-
- rows.each do |row_style, row_arguments|
- row_arguments.each do |row_argument|
- columns.each do |column_style, column_arguments|
- column_arguments.each do |column_argument|
- values.each do |value_type, value|
- expected = solutions[[row_style, column_style]][value_type] || Matrix::ErrDimensionMismatch
-
- result = src.clone
- begin
- result[row_argument, column_argument] = value
- assert_equal expected, result,
- "m[#{row_argument.inspect}][#{column_argument.inspect}] = #{value.inspect} failed"
- rescue Exception => e
- raise if e.class != expected
- end
- end
- end
- end
- end
- end
- end
-
- def test_map!
- m1 = Matrix.zero(2,2)
- m2 = Matrix.build(3,4){|row, col| 1}
- m3 = Matrix.zero(3,5).freeze
- m4 = Matrix.empty.freeze
-
- assert_equal Matrix[[5, 5, 5, 5], [5, 5, 5, 5], [5, 5, 5, 5]], m2.map!{|e| e * 5}
- assert_equal Matrix[[7, 0],[0, 7]], m1.map!(:diagonal){|e| e + 7}
- assert_equal Matrix[[7, 5],[5, 7]], m1.map!(:off_diagonal){|e| e + 5}
- assert_equal Matrix[[12, 5, 5, 5], [12, 12, 5, 5], [12, 12, 12, 5]], m2.map!(:lower){|e| e + 7}
- assert_equal Matrix[[12, 5, 5, 5], [0, 12, 5, 5], [0, 0, 12, 5]], m2.map!(:strict_lower){|e| e - 12}
- assert_equal Matrix[[12, 25, 25, 25], [0, 12, 25, 25], [0, 0, 12, 25]], m2.map!(:strict_upper){|e| e ** 2}
- assert_equal Matrix[[-12, -25, -25, -25], [0, -12, -25, -25], [0, 0, -12, -25]], m2.map!(:upper){|e| -e}
- assert_equal m1, m1.map!{|e| e ** 2 }
- assert_equal m2, m2.map!(:lower){ |e| e - 3 }
- assert_raise(ArgumentError) {m1.map!(:test){|e| e + 7}}
- assert_raise(FrozenError) { m3.map!{|e| e * 2} }
- assert_raise(FrozenError) { m4.map!{} }
- end
-
- def test_freeze
- m = Matrix[[1, 2, 3],[4, 5, 6]]
- f = m.freeze
- assert_equal true, f.frozen?
- assert m.equal?(f)
- assert m.equal?(f.freeze)
- assert_raise(FrozenError){ m[0, 1] = 56 }
- assert_equal m.dup, m
- end
-
- def test_clone
- a = Matrix[[4]]
- def a.foo
- 42
- end
-
- m = a.clone
- m[0, 0] = 2
- assert_equal a, m * 2
- assert_equal 42, m.foo
-
- a.freeze
- m = a.clone
- assert m.frozen?
- assert_equal 42, m.foo
- end
-
- def test_dup
- a = Matrix[[4]]
- def a.foo
- 42
- end
- a.freeze
-
- m = a.dup
- m[0, 0] = 2
- assert_equal a, m * 2
- assert !m.respond_to?(:foo)
- end
-
- def test_eigenvalues_and_eigenvectors_symmetric
- m = Matrix[
- [8, 1],
- [1, 8]
- ]
- values = m.eigensystem.eigenvalues
- assert_in_epsilon(7.0, values[0])
- assert_in_epsilon(9.0, values[1])
- vectors = m.eigensystem.eigenvectors
- assert_in_epsilon(-vectors[0][0], vectors[0][1])
- assert_in_epsilon(vectors[1][0], vectors[1][1])
- end
-
- def test_eigenvalues_and_eigenvectors_nonsymmetric
- m = Matrix[
- [8, 1],
- [4, 5]
- ]
- values = m.eigensystem.eigenvalues
- assert_in_epsilon(9.0, values[0])
- assert_in_epsilon(4.0, values[1])
- vectors = m.eigensystem.eigenvectors
- assert_in_epsilon(vectors[0][0], vectors[0][1])
- assert_in_epsilon(-4 * vectors[1][0], vectors[1][1])
- end
-
- def test_unitary?
- assert_equal true, @rot.unitary?
- assert_equal true, ((0+1i) * @rot).unitary?
- assert_equal false, @a3.unitary?
- assert_raise(Matrix::ErrDimensionMismatch) { @m1.unitary? }
- end
-
- def test_orthogonal
- assert_equal true, @rot.orthogonal?
- assert_equal false, ((0+1i) * @rot).orthogonal?
- assert_equal false, @a3.orthogonal?
- assert_raise(Matrix::ErrDimensionMismatch) { @m1.orthogonal? }
- end
-
- def test_adjoint
- assert_equal(Matrix[[(1-2i), 1], [(0-1i), 2], [0, 3]], @c1.adjoint)
- assert_equal(Matrix.empty(0,2), @e1.adjoint)
- end
-
- def test_ractor
- assert_ractor(<<~RUBY, require: 'matrix')
- obj1 = Matrix[[1, 2], [3, 4]].freeze
-
- obj2 = Ractor.new obj1 do |obj|
- obj
- end.take
- assert_same obj1, obj2
- RUBY
- end if defined?(Ractor)
-
- def test_rotate_with_symbol
- assert_equal(Matrix[[4, 1], [5, 2], [6, 3]], @m1.rotate_entries)
- assert_equal(@m1.rotate_entries, @m1.rotate_entries(:clockwise))
- assert_equal(Matrix[[4, 1], [5, 2], [6, 3]],
- @m1.rotate_entries(:clockwise))
- assert_equal(Matrix[[3, 6], [2, 5], [1, 4]],
- @m1.rotate_entries(:counter_clockwise))
- assert_equal(Matrix[[6, 5, 4], [3, 2, 1]],
- @m1.rotate_entries(:half_turn))
- assert_equal(Matrix[[6, 5, 4], [3, 2, 1]],
- @m1.rotate_entries(:half_turn))
- assert_equal(Matrix.empty(0,2),
- @e1.rotate_entries(:clockwise))
- assert_equal(Matrix.empty(0,2),
- @e1.rotate_entries(:counter_clockwise))
- assert_equal(Matrix.empty(2,0),
- @e1.rotate_entries(:half_turn))
- assert_equal(Matrix.empty(0,3),
- @e2.rotate_entries(:half_turn))
- end
-
- def test_rotate_with_numeric
- assert_equal(Matrix[[4, 1], [5, 2], [6, 3]],
- @m1.rotate_entries(1))
- assert_equal(@m2.rotate_entries(:half_turn),
- @m2.rotate_entries(2))
- assert_equal(@m2.rotate_entries(:half_turn),
- @m1.rotate_entries(2))
- assert_equal(@m1.rotate_entries(:counter_clockwise),
- @m1.rotate_entries(3))
- assert_equal(@m1,
- @m1.rotate_entries(4))
- assert_equal(@m1,
- @m1.rotate_entries(4))
- assert_not_same(@m1,
- @m1.rotate_entries(4))
- assert_equal(@m1.rotate_entries(:clockwise),
- @m1.rotate_entries(5))
- assert_equal(Matrix.empty(0,2),
- @e1.rotate_entries(1))
- assert_equal(@e2,
- @e2.rotate_entries(2))
- assert_equal(@e2.rotate_entries(1),
- @e2.rotate_entries(3))
- assert_equal(@e2.rotate_entries(:counter_clockwise),
- @e2.rotate_entries(-1))
- assert_equal(@m1.rotate_entries(:counter_clockwise),
- @m1.rotate_entries(-1))
- assert_equal(Matrix[[6, 5, 4], [3, 2, 1]],
- @m1.rotate_entries(-2))
- assert_equal(@m1,
- @m1.rotate_entries(-4))
- assert_equal(@m1.rotate_entries(-1),
- @m1.rotate_entries(-5))
- end
-end
diff --git a/test/matrix/test_vector.rb b/test/matrix/test_vector.rb
deleted file mode 100644
index 8b771efee3..0000000000
--- a/test/matrix/test_vector.rb
+++ /dev/null
@@ -1,335 +0,0 @@
-# frozen_string_literal: false
-require 'test/unit'
-require 'matrix'
-
-class TestVector < Test::Unit::TestCase
- def setup
- @v1 = Vector[1,2,3]
- @v2 = Vector[1,2,3]
- @v3 = @v1.clone
- @v4 = Vector[1.0, 2.0, 3.0]
- @w1 = Vector[2,3,4]
- end
-
- def test_zero
- assert_equal Vector[0, 0, 0, 0], Vector.zero(4)
- assert_equal Vector[], Vector.zero(0)
- assert_raise(ArgumentError) { Vector.zero(-1) }
- assert Vector[0, 0, 0, 0].zero?
- end
-
- def test_basis
- assert_equal(Vector[1, 0, 0], Vector.basis(size: 3, index: 0))
- assert_raise(ArgumentError) { Vector.basis(size: -1, index: 2) }
- assert_raise(ArgumentError) { Vector.basis(size: 4, index: -1) }
- assert_raise(ArgumentError) { Vector.basis(size: 3, index: 3) }
- assert_raise(ArgumentError) { Vector.basis(size: 3) }
- assert_raise(ArgumentError) { Vector.basis(index: 3) }
- end
-
- def test_get_element
- assert_equal(@v1[0..], [1, 2, 3])
- assert_equal(@v1[0..1], [1, 2])
- assert_equal(@v1[2], 3)
- assert_equal(@v1[4], nil)
- end
-
- def test_set_element
-
- assert_block do
- v = Vector[5, 6, 7, 8, 9]
- v[1..2] = Vector[1, 2]
- v == Vector[5, 1, 2, 8, 9]
- end
-
- assert_block do
- v = Vector[6, 7, 8]
- v[1..2] = Matrix[[1, 3]]
- v == Vector[6, 1, 3]
- end
-
- assert_block do
- v = Vector[1, 2, 3, 4, 5, 6]
- v[0..2] = 8
- v == Vector[8, 8, 8, 4, 5, 6]
- end
-
- assert_block do
- v = Vector[1, 3, 4, 5]
- v[2] = 5
- v == Vector[1, 3, 5, 5]
- end
-
- assert_block do
- v = Vector[2, 3, 5]
- v[-2] = 13
- v == Vector[2, 13, 5]
- end
-
- assert_block do
- v = Vector[4, 8, 9, 11, 30]
- v[1..-2] = Vector[1, 2, 3]
- v == Vector[4, 1, 2, 3, 30]
- end
-
- assert_raise(IndexError) {Vector[1, 3, 4, 5][5..6] = 17}
- assert_raise(IndexError) {Vector[1, 3, 4, 5][6] = 17}
- assert_raise(Matrix::ErrDimensionMismatch) {Vector[1, 3, 4, 5][0..2] = Matrix[[1], [2], [3]]}
- assert_raise(ArgumentError) {Vector[1, 2, 3, 4, 5, 6][0..2] = Vector[1, 2, 3, 4, 5, 6]}
- assert_raise(FrozenError) { Vector[7, 8, 9].freeze[0..1] = 5}
- end
-
- def test_map!
- v1 = Vector[1, 2, 3]
- v2 = Vector[1, 3, 5].freeze
- v3 = Vector[].freeze
- assert_equal Vector[1, 4, 9], v1.map!{|e| e ** 2}
- assert_equal v1, v1.map!{|e| e - 8}
- assert_raise(FrozenError) { v2.map!{|e| e + 2 }}
- assert_raise(FrozenError){ v3.map!{} }
- end
-
- def test_freeze
- v = Vector[1,2,3]
- f = v.freeze
- assert_equal true, f.frozen?
- assert v.equal?(f)
- assert v.equal?(f.freeze)
- assert_raise(FrozenError){ v[1] = 56 }
- assert_equal v.dup, v
- end
-
- def test_clone
- a = Vector[4]
- def a.foo
- 42
- end
-
- v = a.clone
- v[0] = 2
- assert_equal a, v * 2
- assert_equal 42, v.foo
-
- a.freeze
- v = a.clone
- assert v.frozen?
- assert_equal 42, v.foo
- end
-
- def test_dup
- a = Vector[4]
- def a.foo
- 42
- end
- a.freeze
-
- v = a.dup
- v[0] = 2
- assert_equal a, v * 2
- assert !v.respond_to?(:foo)
- end
-
- def test_identity
- assert_same @v1, @v1
- assert_not_same @v1, @v2
- assert_not_same @v1, @v3
- assert_not_same @v1, @v4
- assert_not_same @v1, @w1
- end
-
- def test_equality
- assert_equal @v1, @v1
- assert_equal @v1, @v2
- assert_equal @v1, @v3
- assert_equal @v1, @v4
- assert_not_equal @v1, @w1
- end
-
- def test_hash_equality
- assert @v1.eql?(@v1)
- assert @v1.eql?(@v2)
- assert @v1.eql?(@v3)
- assert !@v1.eql?(@v4)
- assert !@v1.eql?(@w1)
-
- hash = { @v1 => :value }
- assert hash.key?(@v1)
- assert hash.key?(@v2)
- assert hash.key?(@v3)
- assert !hash.key?(@v4)
- assert !hash.key?(@w1)
- end
-
- def test_hash
- assert_equal @v1.hash, @v1.hash
- assert_equal @v1.hash, @v2.hash
- assert_equal @v1.hash, @v3.hash
- end
-
- def test_aref
- assert_equal(1, @v1[0])
- assert_equal(2, @v1[1])
- assert_equal(3, @v1[2])
- assert_equal(3, @v1[-1])
- assert_equal(nil, @v1[3])
- end
-
- def test_size
- assert_equal(3, @v1.size)
- end
-
- def test_each2
- a = []
- @v1.each2(@v4) {|x, y| a << [x, y] }
- assert_equal([[1,1.0],[2,2.0],[3,3.0]], a)
- end
-
- def test_collect
- a = @v1.collect {|x| x + 1 }
- assert_equal(Vector[2,3,4], a)
- end
-
- def test_collect2
- a = @v1.collect2(@v4) {|x, y| x + y }
- assert_equal([2.0,4.0,6.0], a)
- end
-
- def test_map2
- a = @v1.map2(@v4) {|x, y| x + y }
- assert_equal(Vector[2.0,4.0,6.0], a)
- end
-
- def test_independent?
- assert(Vector.independent?(@v1, @w1))
- assert(
- Vector.independent?(
- Vector.basis(size: 3, index: 0),
- Vector.basis(size: 3, index: 1),
- Vector.basis(size: 3, index: 2),
- )
- )
-
- refute(Vector.independent?(@v1, Vector[2,4,6]))
- refute(Vector.independent?(Vector[2,4], Vector[1,3], Vector[5,6]))
-
- assert_raise(TypeError) { Vector.independent?(@v1, 3) }
- assert_raise(Vector::ErrDimensionMismatch) { Vector.independent?(@v1, Vector[2,4]) }
-
- assert(@v1.independent?(Vector[1,2,4], Vector[1,3,4]))
- end
-
- def test_mul
- assert_equal(Vector[2,4,6], @v1 * 2)
- assert_equal(Matrix[[1, 4, 9], [2, 8, 18], [3, 12, 27]], @v1 * Matrix[[1,4,9]])
- assert_raise(Matrix::ErrOperationNotDefined) { @v1 * Vector[1,4,9] }
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(1, Vector[1, 2, 3] * o)
- end
-
- def test_add
- assert_equal(Vector[2,4,6], @v1 + @v1)
- assert_equal(Matrix[[2],[6],[12]], @v1 + Matrix[[1],[4],[9]])
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(2, Vector[1, 2, 3] + o)
- end
-
- def test_sub
- assert_equal(Vector[0,0,0], @v1 - @v1)
- assert_equal(Matrix[[0],[-2],[-6]], @v1 - Matrix[[1],[4],[9]])
- o = Object.new
- def o.coerce(x)
- [1, 1]
- end
- assert_equal(0, Vector[1, 2, 3] - o)
- end
-
- def test_uplus
- assert_equal(@v1, +@v1)
- end
-
- def test_negate
- assert_equal(Vector[-1, -2, -3], -@v1)
- assert_equal(@v1, -(-@v1))
- end
-
- def test_inner_product
- assert_equal(1+4+9, @v1.inner_product(@v1))
- assert_equal(1+4+9, @v1.dot(@v1))
- end
-
- def test_r
- assert_equal(5, Vector[3, 4].r)
- end
-
- def test_round
- assert_equal(Vector[1.234, 2.345, 3.40].round(2), Vector[1.23, 2.35, 3.4])
- end
-
- def test_covector
- assert_equal(Matrix[[1,2,3]], @v1.covector)
- end
-
- def test_to_s
- assert_equal("Vector[1, 2, 3]", @v1.to_s)
- end
-
- def test_to_matrix
- assert_equal Matrix[[1], [2], [3]], @v1.to_matrix
- end
-
- def test_inspect
- assert_equal("Vector[1, 2, 3]", @v1.inspect)
- end
-
- def test_magnitude
- assert_in_epsilon(3.7416573867739413, @v1.norm)
- assert_in_epsilon(3.7416573867739413, @v4.norm)
- end
-
- def test_complex_magnitude
- bug6966 = '[ruby-dev:46100]'
- v = Vector[Complex(0,1), 0]
- assert_equal(1.0, v.norm, bug6966)
- end
-
- def test_rational_magnitude
- v = Vector[Rational(1,2), 0]
- assert_equal(0.5, v.norm)
- end
-
- def test_cross_product
- v = Vector[1, 0, 0].cross_product Vector[0, 1, 0]
- assert_equal(Vector[0, 0, 1], v)
- v2 = Vector[1, 2].cross_product
- assert_equal(Vector[-2, 1], v2)
- v3 = Vector[3, 5, 2, 1].cross(Vector[4, 3, 1, 8], Vector[2, 9, 4, 3])
- assert_equal(Vector[16, -65, 139, -1], v3)
- assert_equal Vector[0, 0, 0, 1],
- Vector[1, 0, 0, 0].cross(Vector[0, 1, 0, 0], Vector[0, 0, 1, 0])
- assert_equal Vector[0, 0, 0, 0, 1],
- Vector[1, 0, 0, 0, 0].cross(Vector[0, 1, 0, 0, 0], Vector[0, 0, 1, 0, 0], Vector[0, 0, 0, 1, 0])
- assert_raise(Vector::ErrDimensionMismatch) { Vector[1, 2, 3].cross_product(Vector[1, 4]) }
- assert_raise(TypeError) { Vector[1, 2, 3].cross_product(42) }
- assert_raise(ArgumentError) { Vector[1, 2].cross_product(Vector[2, -1]) }
- assert_raise(Vector::ErrOperationNotDefined) { Vector[1].cross_product }
- end
-
- def test_angle_with
- assert_in_epsilon(Math::PI, Vector[1, 0].angle_with(Vector[-1, 0]))
- assert_in_epsilon(Math::PI/2, Vector[1, 0].angle_with(Vector[0, -1]))
- assert_in_epsilon(Math::PI/4, Vector[2, 2].angle_with(Vector[0, 1]))
- assert_in_delta(0.0, Vector[1, 1].angle_with(Vector[1, 1]), 0.00001)
- assert_equal(Vector[6, 6].angle_with(Vector[7, 7]), 0.0)
- assert_equal(Vector[6, 6].angle_with(Vector[-7, -7]), Math::PI)
-
- assert_raise(Vector::ZeroVectorError) { Vector[1, 1].angle_with(Vector[0, 0]) }
- assert_raise(Vector::ZeroVectorError) { Vector[0, 0].angle_with(Vector[1, 1]) }
- assert_raise(Matrix::ErrDimensionMismatch) { Vector[1, 2, 3].angle_with(Vector[0, 1]) }
- end
-end
diff --git a/tool/sync_default_gems.rb b/tool/sync_default_gems.rb
index 37737b9711..5f95a3d80b 100644
--- a/tool/sync_default_gems.rb
+++ b/tool/sync_default_gems.rb
@@ -24,7 +24,6 @@ REPOSITORIES = {
strscan: 'ruby/strscan',
ipaddr: 'ruby/ipaddr',
logger: 'ruby/logger',
- matrix: 'ruby/matrix',
ostruct: 'ruby/ostruct',
irb: 'ruby/irb',
forwardable: "ruby/forwardable",
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-23 10:29:22 +0000