diff options
Diffstat (limited to 'ruby_1_9_3/lib/matrix.rb')
| -rw-r--r-- | ruby_1_9_3/lib/matrix.rb | 1866 |
1 files changed, 1866 insertions, 0 deletions
diff --git a/ruby_1_9_3/lib/matrix.rb b/ruby_1_9_3/lib/matrix.rb new file mode 100644 index 0000000000..3c078c0c06 --- /dev/null +++ b/ruby_1_9_3/lib/matrix.rb @@ -0,0 +1,1866 @@ +# encoding: utf-8 +# +# = matrix.rb +# +# An implementation of Matrix and Vector classes. +# +# See classes Matrix and Vector for documentation. +# +# Current Maintainer:: Marc-André Lafortune +# Original Author:: Keiju ISHITSUKA +# Original Documentation:: Gavin Sinclair (sourced from <i>Ruby in a Nutshell</i> (Matsumoto, O'Reilly)) +## + +require "e2mmap.rb" + +module ExceptionForMatrix # :nodoc: + extend Exception2MessageMapper + def_e2message(TypeError, "wrong argument type %s (expected %s)") + def_e2message(ArgumentError, "Wrong # of arguments(%d for %d)") + + def_exception("ErrDimensionMismatch", "\#{self.name} dimension mismatch") + def_exception("ErrNotRegular", "Not Regular Matrix") + def_exception("ErrOperationNotDefined", "Operation(%s) can\\'t be defined: %s op %s") + def_exception("ErrOperationNotImplemented", "Sorry, Operation(%s) not implemented: %s op %s") +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 (trace, rank, inverse, determinant). +# +# == Method Catalogue +# +# To create a matrix: +# * <tt> Matrix[*rows] </tt> +# * <tt> Matrix.[](*rows) </tt> +# * <tt> Matrix.rows(rows, copy = true) </tt> +# * <tt> Matrix.columns(columns) </tt> +# * <tt> Matrix.build(row_size, column_size, &block) </tt> +# * <tt> Matrix.diagonal(*values) </tt> +# * <tt> Matrix.scalar(n, value) </tt> +# * <tt> Matrix.identity(n) </tt> +# * <tt> Matrix.unit(n) </tt> +# * <tt> Matrix.I(n) </tt> +# * <tt> Matrix.zero(n) </tt> +# * <tt> Matrix.row_vector(row) </tt> +# * <tt> Matrix.column_vector(column) </tt> +# +# To access Matrix elements/columns/rows/submatrices/properties: +# * <tt> [](i, j) </tt> +# * <tt> #row_size </tt> +# * <tt> #column_size </tt> +# * <tt> #row(i) </tt> +# * <tt> #column(j) </tt> +# * <tt> #collect </tt> +# * <tt> #map </tt> +# * <tt> #each </tt> +# * <tt> #each_with_index </tt> +# * <tt> #find_index </tt> +# * <tt> #minor(*param) </tt> +# +# Properties of a matrix: +# * <tt> #diagonal? </tt> +# * <tt> #empty? </tt> +# * <tt> #hermitian? </tt> +# * <tt> #lower_triangular? </tt> +# * <tt> #normal? </tt> +# * <tt> #orthogonal? </tt> +# * <tt> #permutation? </tt> +# * <tt> #real? </tt> +# * <tt> #regular? </tt> +# * <tt> #singular? </tt> +# * <tt> #square? </tt> +# * <tt> #symmetric? </tt> +# * <tt> #unitary? </tt> +# * <tt> #upper_triangular? </tt> +# * <tt> #zero? </tt> +# +# Matrix arithmetic: +# * <tt> *(m) </tt> +# * <tt> +(m) </tt> +# * <tt> -(m) </tt> +# * <tt> #/(m) </tt> +# * <tt> #inverse </tt> +# * <tt> #inv </tt> +# * <tt> ** </tt> +# +# Matrix functions: +# * <tt> #determinant </tt> +# * <tt> #det </tt> +# * <tt> #rank </tt> +# * <tt> #round </tt> +# * <tt> #trace </tt> +# * <tt> #tr </tt> +# * <tt> #transpose </tt> +# * <tt> #t </tt> +# +# Matrix decompositions: +# * <tt> #eigen </tt> +# * <tt> #eigensystem </tt> +# * <tt> #lup </tt> +# * <tt> #lup_decomposition </tt> +# +# Complex arithmetic: +# * <tt> conj </tt> +# * <tt> conjugate </tt> +# * <tt> imag </tt> +# * <tt> imaginary </tt> +# * <tt> real </tt> +# * <tt> rect </tt> +# * <tt> rectangular </tt> +# +# Conversion to other data types: +# * <tt> #coerce(other) </tt> +# * <tt> #row_vectors </tt> +# * <tt> #column_vectors </tt> +# * <tt> #to_a </tt> +# +# String representations: +# * <tt> #to_s </tt> +# * <tt> #inspect </tt> +# +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) + Matrix.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) + rows.map! do |row| + convert_to_array(row, copy) + end + size = (rows[0] || []).size + rows.each do |row| + Matrix.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) + Matrix.rows(columns, false).transpose + end + + # + # Creates a matrix of size +row_size+ x +column_size+. + # 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_size, column_size = row_size) + row_size = CoercionHelper.coerce_to_int(row_size) + column_size = CoercionHelper.coerce_to_int(column_size) + raise ArgumentError if row_size < 0 || column_size < 0 + return to_enum :build, row_size, column_size unless block_given? + rows = Array.new(row_size) do |i| + Array.new(column_size) do |j| + yield i, j + end + end + new rows, column_size + 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 + 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) + Matrix.diagonal(*Array.new(n, value)) + end + + # + # Creates an +n+ by +n+ identity matrix. + # Matrix.identity(2) + # => 1 0 + # 0 1 + # + def Matrix.identity(n) + Matrix.scalar(n, 1) + end + class << Matrix + alias unit identity + alias I identity + end + + # + # Creates a zero matrix. + # Matrix.zero(2) + # => 0 0 + # 0 0 + # + def Matrix.zero(row_size, column_size = row_size) + rows = Array.new(row_size){Array.new(column_size, 0)} + new rows, column_size + 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_size+ x +column_size+. + # At least one of +row_size+ or +column_size+ 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_size = 0, column_size = 0) + Matrix.Raise ArgumentError, "One size must be 0" if column_size != 0 && row_size != 0 + Matrix.Raise ArgumentError, "Negative size" if column_size < 0 || row_size < 0 + + new([[]]*row_size, column_size) + end + + # + # Matrix.new is private; use Matrix.rows, columns, [], etc... to create. + # + def initialize(rows, column_size = rows[0].size) + # No checking is done at this point. rows must be an Array of Arrays. + # column_size 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_size = column_size + end + + def new_matrix(rows, column_size = rows[0].size) # :nodoc: + Matrix.send(:new, rows, column_size) # bypass privacy of Matrix.new + end + private :new_matrix + + # + # 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 [] + + def []=(i, j, v) + @rows[i][j] = v + end + alias set_element []= + alias set_component []= + private :[]=, :set_element, :set_component + + # + # Returns the number of rows. + # + def row_size + @rows.size + end + + # + # Returns the number of columns. + # + attr_reader :column_size + + # + # 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_size || j < -column_size + row_size.times do |i| + yield @rows[i][j] + end + self + else + return nil if j >= column_size || j < -column_size + col = Array.new(row_size) {|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. + # Matrix[ [1,2], [3,4] ].collect { |e| e**2 } + # => 1 4 + # 9 16 + # + def collect(&block) # :yield: e + return to_enum(:collect) unless block_given? + rows = @rows.collect{|row| row.collect(&block)} + new_matrix rows, column_size + end + alias map collect + + # + # Yields all elements of the matrix, starting with those of the first row, + # or returns an Enumerator is 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) # :yield: e + return to_enum :each, which unless block_given? + last = column_size - 1 + case which + when :all + block = Proc.new + @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_size.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_size].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 + Matrix.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_size - 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_size.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_size].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 + Matrix.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_size or column_size respectively. + # + def minor(*param) + case param.size + when 2 + row_range, col_range = param + from_row = row_range.first + from_row += row_size if from_row < 0 + to_row = row_range.end + to_row += row_size 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_size if from_col < 0 + to_col = col_range.end + to_col += column_size 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_size if from_row < 0 + from_col += column_size if from_col < 0 + else + Matrix.Raise ArgumentError, param.inspect + end + + return nil if from_row > row_size || from_col > column_size || from_row < 0 || from_col < 0 + rows = @rows[from_row, size_row].collect{|row| + row[from_col, size_col] + } + new_matrix rows, [column_size - from_col, size_col].min + end + + #-- + # TESTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + #++ + + # + # Returns +true+ is this is a diagonal matrix. + # Raises an error if matrix is not square. + # + def diagonal? + Matrix.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_size == 0 || row_size == 0 + end + + # + # Returns +true+ is this is an hermitian matrix. + # Raises an error if matrix is not square. + # + def hermitian? + Matrix.Raise ErrDimensionMismatch unless square? + each_with_index(:strict_upper).all? do |e, row, col| + e == rows[col][row].conj + end + end + + # + # Returns +true+ is this is a lower triangular matrix. + # + def lower_triangular? + each(:strict_upper).all?(&:zero?) + end + + # + # Returns +true+ is this is a normal matrix. + # Raises an error if matrix is not square. + # + def normal? + Matrix.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+ is this is an orthogonal matrix + # Raises an error if matrix is not square. + # + def orthogonal? + Matrix.Raise ErrDimensionMismatch unless square? + rows.each_with_index do |row, i| + column_size.times do |j| + s = 0 + row_size.times do |k| + s += row[k] * rows[k][j] + end + return false unless s == (i == j ? 1 : 0) + end + end + true + end + + # + # Returns +true+ is this is a permutation matrix + # Raises an error if matrix is not square. + # + def permutation? + Matrix.Raise ErrDimensionMismatch unless square? + cols = Array.new(column_size) + 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+ is this is a singular matrix. + # + def singular? + determinant == 0 + end + + # + # Returns +true+ is this is a square matrix. + # + def square? + column_size == row_size + end + + # + # Returns +true+ is this is a symmetric matrix. + # Raises an error if matrix is not square. + # + def symmetric? + Matrix.Raise ErrDimensionMismatch unless square? + each_with_index(:strict_upper).all? do |e, row, col| + e == rows[col][row] + end + end + + # + # Returns +true+ is this is a unitary matrix + # Raises an error if matrix is not square. + # + def unitary? + Matrix.Raise ErrDimensionMismatch unless square? + rows.each_with_index do |row, i| + column_size.times do |j| + s = 0 + row_size.times do |k| + s += row[k].conj * rows[k][j] + end + return false unless s == (i == j ? 1 : 0) + end + end + true + end + + # + # Returns +true+ is this is an upper triangular matrix. + # + def upper_triangular? + each(:strict_lower).all?(&:zero?) + end + + # + # Returns +true+ is this is a matrix with only zero elements + # + def zero? + all?(&:zero?) + end + + #-- + # OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + #++ + + # + # Returns +true+ if and only if the two matrices contain equal elements. + # + def ==(other) + return false unless Matrix === other && + column_size == other.column_size # necessary for empty matrices + rows == other.rows + end + + def eql?(other) + return false unless Matrix === other && + column_size == other.column_size # necessary for empty matrices + rows.eql? other.rows + end + + # + # Returns a clone of the matrix, so that the contents of each do not reference + # identical objects. + # There should be no good reason to do this since Matrices are immutable. + # + def clone + new_matrix @rows.map(&:dup), column_size + 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 + rows = @rows.collect {|row| + row.collect {|e| e * m } + } + return new_matrix rows, column_size + when Vector + m = Matrix.column_vector(m) + r = self * m + return r.column(0) + when Matrix + Matrix.Raise ErrDimensionMismatch if column_size != m.row_size + + rows = Array.new(row_size) {|i| + Array.new(m.column_size) {|j| + (0 ... column_size).inject(0) do |vij, k| + vij + self[i, k] * m[k, j] + end + } + } + return new_matrix rows, m.column_size + 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 + Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class + when Vector + m = Matrix.column_vector(m) + when Matrix + else + return apply_through_coercion(m, __method__) + end + + Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size + + rows = Array.new(row_size) {|i| + Array.new(column_size) {|j| + self[i, j] + m[i, j] + } + } + new_matrix rows, column_size + end + + # + # Matrix subtraction. + # Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]] + # => -8 2 + # 8 1 + # + def -(m) + case m + when Numeric + Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class + when Vector + m = Matrix.column_vector(m) + when Matrix + else + return apply_through_coercion(m, __method__) + end + + Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size + + rows = Array.new(row_size) {|i| + Array.new(column_size) {|j| + self[i, j] - m[i, j] + } + } + new_matrix rows, column_size + 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_size + when Matrix + return self * other.inverse + else + return apply_through_coercion(other, __method__) + end + end + + # + # Returns the inverse of the matrix. + # Matrix[[-1, -1], [0, -1]].inverse + # => -1 1 + # 0 -1 + # + def inverse + Matrix.Raise ErrDimensionMismatch unless square? + Matrix.I(row_size).send(:inverse_from, self) + end + alias inv inverse + + def inverse_from(src) # :nodoc: + last = row_size - 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 + Matrix.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 + private :inverse_from + + # + # 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 ** (other) + case other + when Integer + x = self + if other <= 0 + x = self.inverse + return Matrix.identity(self.column_size) if other == 0 + other = -other + end + z = nil + loop do + z = z ? z * x : x if other[0] == 1 + return z if (other >>= 1).zero? + x *= x + end + when Numeric + v, d, v_inv = eigensystem + v * Matrix.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv + else + Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class + end + 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 + Matrix.Raise ErrDimensionMismatch unless square? + m = @rows + case row_size + # 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. + # + def determinant_bareiss + size = row_size + 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 + private :determinant_bareiss + + # + # deprecated; use Matrix#determinant + # + def determinant_e + warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant" + rank + end + alias det_e determinant_e + + # + # 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_size - 1 + last_row = row_size - 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 "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank" + rank + 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 + Matrix.Raise ErrDimensionMismatch unless square? + (0...column_size).inject(0) do |tr, i| + tr + @rows[i][i] + end + end + alias 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 Matrix.empty(column_size, 0) if row_size.zero? + new_matrix @rows.transpose, row_size + end + alias t transpose + + #-- + # 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 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 == a * p # => true + # a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)] + # + def lup + LUPDecomposition.new(self) + end + alias 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 conj conjugate + + # + # 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 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 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_size) {|i| + row(i) + } + end + + # + # Returns an array of the column vectors of the matrix. See Vector. + # + def column_vectors + Array.new(column_size) {|i| + column(i) + } + end + + # + # Returns an array of arrays that describe the rows of the matrix. + # + def to_a + @rows.collect(&:dup) + end + + def elements_to_f + warn "#{caller(1)[0]}: warning: Matrix#elements_to_f is deprecated, use map(&:to_f)" + map(&:to_f) + end + + def elements_to_i + warn "#{caller(1)[0]}: warning: Matrix#elements_to_i is deprecated, use map(&:to_i)" + map(&:to_i) + end + + def elements_to_r + warn "#{caller(1)[0]}: warning: Matrix#elements_to_r is deprecated, use map(&:to_r)" + map(&:to_r) + end + + #-- + # PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + #++ + + # + # Overrides Object#to_s + # + def to_s + if empty? + "Matrix.empty(#{row_size}, #{column_size})" + else + "Matrix[" + @rows.collect{|row| + "[" + row.collect{|e| e.to_s}.join(", ") + "]" + }.join(", ")+"]" + end + end + + # + # Overrides Object#inspect + # + def inspect + if empty? + "Matrix.empty(#{row_size}, #{column_size})" + else + "Matrix#{@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. + # + 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 + private :convert_to_array + end + + extend ConversionHelper + + module CoercionHelper # :nodoc: + # + # Applies the operator +oper+ with argument +obj+ + # through coercion of +obj+ + # + 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 + private :apply_through_coercion + + # + # 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) + + 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 + 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 + Scalar.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 + Scalar.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 + Scalar.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 + Scalar.Raise ErrOperationNotDefined, "**", @value.class, other.class + when Matrix + #other.powered_by(self) + Scalar.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: +# * <tt> Vector.[](*array) </tt> +# * <tt> Vector.elements(array, copy = true) </tt> +# +# To access elements: +# * <tt> [](i) </tt> +# +# To enumerate the elements: +# * <tt> #each2(v) </tt> +# * <tt> #collect2(v) </tt> +# +# Vector arithmetic: +# * <tt> *(x) "is matrix or number" </tt> +# * <tt> +(v) </tt> +# * <tt> -(v) </tt> +# +# Vector functions: +# * <tt> #inner_product(v) </tt> +# * <tt> #collect </tt> +# * <tt> #magnitude </tt> +# * <tt> #map </tt> +# * <tt> #map2(v) </tt> +# * <tt> #norm </tt> +# * <tt> #normalize </tt> +# * <tt> #r </tt> +# * <tt> #size </tt> +# +# Conversion to other data types: +# * <tt> #covector </tt> +# * <tt> #to_a </tt> +# * <tt> #coerce(other) </tt> +# +# String representations: +# * <tt> #to_s </tt> +# * <tt> #inspect </tt> +# +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 + + # + # 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 + + # + # Returns element number +i+ (starting at zero) of the vector. + # + def [](i) + @elements[i] + end + alias element [] + alias component [] + + def []=(i, v) + @elements[i]= v + end + alias set_element []= + alias set_component []= + private :[]=, :set_element, :set_component + + # + # 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) + Vector.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) + Vector.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 + + #-- + # COMPARING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + #++ + + # + # Returns +true+ iff 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 + + # + # Return a copy of the vector. + # + def clone + Vector.elements(@elements) + end + + # + # Return a hash-code for the vector. + # + def hash + @elements.hash + end + + #-- + # ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + #++ + + # + # Multiplies the vector by +x+, where +x+ is a number or another vector. + # + def *(x) + case x + when Numeric + els = @elements.collect{|e| e * x} + Vector.elements(els, false) + when Matrix + Matrix.column_vector(self) * x + when Vector + Vector.Raise ErrOperationNotDefined, "*", self.class, x.class + else + apply_through_coercion(x, __method__) + end + end + + # + # Vector addition. + # + def +(v) + case v + when Vector + Vector.Raise ErrDimensionMismatch if size != v.size + els = collect2(v) {|v1, v2| + v1 + v2 + } + Vector.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 + Vector.Raise ErrDimensionMismatch if size != v.size + els = collect2(v) {|v1, v2| + v1 - v2 + } + Vector.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} + Vector.elements(els, false) + when Matrix, Vector + Vector.Raise ErrOperationNotDefined, "/", self.class, x.class + else + apply_through_coercion(x, __method__) + end + 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) + Vector.Raise ErrDimensionMismatch if size != v.size + + p = 0 + each2(v) {|v1, v2| + p += v1 * v2 + } + p + end + + # + # Like Array#collect. + # + def collect(&block) # :yield: e + return to_enum(:collect) unless block_given? + els = @elements.collect(&block) + Vector.elements(els, false) + 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*e}) + end + alias r magnitude + alias 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) + Vector.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 + + #-- + # 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 + + def elements_to_f + warn "#{caller(1)[0]}: warning: Vector#elements_to_f is deprecated" + map(&:to_f) + end + + def elements_to_i + warn "#{caller(1)[0]}: warning: Vector#elements_to_i is deprecated" + map(&:to_i) + end + + def elements_to_r + warn "#{caller(1)[0]}: warning: Vector#elements_to_r is deprecated" + 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 |