From 454a36794f83395d0827a9e2e85ac8f0d9e53e16 Mon Sep 17 00:00:00 2001 From: Hiroshi SHIBATA Date: Wed, 26 May 2021 15:36:16 +0900 Subject: Promote matrix to the bundled gems --- doc/maintainers.rdoc | 6 +- doc/standard_library.rdoc | 2 +- gems/bundled_gems | 1 + lib/matrix.rb | 2493 -------------------------------- lib/matrix/eigenvalue_decomposition.rb | 882 ----------- lib/matrix/lup_decomposition.rb | 219 --- lib/matrix/matrix.gemspec | 26 - lib/matrix/version.rb | 5 - misc/expand_tabs.rb | 1 - test/matrix/test_matrix.rb | 888 ------------ test/matrix/test_vector.rb | 335 ----- tool/sync_default_gems.rb | 1 - 12 files changed, 4 insertions(+), 4855 deletions(-) delete mode 100644 lib/matrix.rb delete mode 100644 lib/matrix/eigenvalue_decomposition.rb delete mode 100644 lib/matrix/lup_decomposition.rb delete mode 100644 lib/matrix/matrix.gemspec delete mode 100644 lib/matrix/version.rb delete mode 100644 test/matrix/test_matrix.rb delete mode 100644 test/matrix/test_vector.rb 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 Ruby in a Nutshell (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 other_matrices 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 map(&:to_f) - def elements_to_f - warn "Matrix#elements_to_f is deprecated, use map(&:to_f)", uplevel: 1 - map(&:to_f) - end - - # Deprecated. - # - # Use map(&:to_i) - def elements_to_i - warn "Matrix#elements_to_i is deprecated, use map(&:to_i)", uplevel: 1 - map(&:to_i) - end - - # Deprecated. - # - # Use map(&:to_r) - 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 A*m = b, - # or equivalently so that L*U*m = P*b - # +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", -- cgit v1.2.3