Index: lib/matrix.rb
===================================================================
--- lib/matrix.rb (revision 23601)
+++ lib/matrix.rb (working copy)
@@ -1,7 +1,7 @@
#--
# matrix.rb -
# $Release Version: 1.0$
-# $Revision: 1.13 $
+# $Revision: 1.14 $
# Original Version from Smalltalk-80 version
# on July 23, 1985 at 8:37:17 am
# by Keiju ISHITSUKA
@@ -35,10 +35,9 @@
# arithmetically and algebraically, and determining their mathematical properties (trace, rank,
# inverse, determinant).
#
-# Note that although matrices should theoretically be rectangular, this is not
-# enforced by the class.
+# Note that matrices must be rectangular.
#
-# Also note that the determinant of integer matrices may be incorrectly calculated unless you
+# Also note that the determinant of integer matrices may be approximated unless you
# also require 'mathn'. This may be fixed in the future.
#
# == Method Catalogue
@@ -50,7 +49,6 @@
# * Matrix.columns(columns)
# * Matrix.diagonal(*values)
# * Matrix.scalar(n, value)
-# * Matrix.scalar(n, value)
# * Matrix.identity(n)
# * Matrix.unit(n)
# * Matrix.I(n)
@@ -102,13 +100,14 @@
# * #inspect
#
class Matrix
- @RCS_ID='-$Id: matrix.rb,v 1.13 2001/12/09 14:22:23 keiju Exp keiju $-'
+ @RCS_ID='-$Id: matrix.rb,v 1.14 2009/05/28 00:00:00 keiju Exp keiju $-'
-# extend Exception2MessageMapper
include ExceptionForMatrix
# instance creations
private_class_method :new
+ attr_reader :rows
+ protected :rows
#
# Creates a matrix where each argument is a row.
@@ -117,18 +116,27 @@
# -1 66
#
def Matrix.[](*rows)
- new(:init_rows, rows, false)
+ Matrix.rows(rows, false)
end
#
# Creates a matrix where +rows+ is an array of arrays, each of which is a row
- # to the matrix. If the optional argument +copy+ is false, use the given
+ # 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)
- new(:init_rows, rows, copy)
+ rows = Matrix.convert_to_array(rows)
+ rows.map! do |row|
+ Matrix.convert_to_array(row, copy)
+ end
+ size = (rows[0] || []).size
+ rows.each do |row|
+ Matrix.Raise ErrDimensionMismatch, "element size differs (#{row.size} should be #{size})" unless row.size == size
+ end
+ new rows, size
end
#
@@ -137,14 +145,8 @@
# => 25 -1
# 93 66
#
- #
def Matrix.columns(columns)
- rows = (0 .. columns[0].size - 1).collect {|i|
- (0 .. columns.size - 1).collect {|j|
- columns[j][i]
- }
- }
- Matrix.rows(rows, false)
+ Matrix.rows(columns, false).transpose
end
#
@@ -156,12 +158,12 @@
#
def Matrix.diagonal(*values)
size = values.size
- rows = (0 .. size - 1).collect {|j|
+ rows = (0 ... size).collect {|j|
row = Array.new(size).fill(0, 0, size)
row[j] = values[j]
row
}
- rows(rows, false)
+ new rows
end
#
@@ -206,14 +208,8 @@
# => 4 5 6
#
def Matrix.row_vector(row)
- case row
- when Vector
- Matrix.rows([row.to_a], false)
- when Array
- Matrix.rows([row.dup], false)
- else
- Matrix.rows([[row]], false)
- end
+ row = Matrix.convert_to_array(row)
+ new [row]
end
#
@@ -225,39 +221,34 @@
# 6
#
def Matrix.column_vector(column)
- case column
- when Vector
- Matrix.columns([column.to_a])
- when Array
- Matrix.columns([column])
- else
- Matrix.columns([[column]])
- end
+ column = Matrix.convert_to_array(column)
+ new column.map{|elem| [elem]}, 1
end
#
# This method is used by the other methods that create matrices, and is of no
# use to general users.
+ # No checking is done at this point. rows must be an Array of Arrays.
+ # column_size must be passed if rows is empty
#
- def initialize(init_method, *argv)
- self.send(init_method, *argv)
+ def initialize(rows, column_size = rows[0].size)
+ @rows = rows
+ @column_size = column_size
end
- def init_rows(rows, copy)
- if copy
- @rows = rows.collect{|row| row.dup}
- else
- @rows = rows
- end
- self
+ #
+ # Use to bypass privacy of Matrix.new
+ #
+ def new(rows, column_size = rows[0].size)
+ Matrix.send(:new, rows, column_size)
end
- private :init_rows
+ private :new
#
# Returns element (+i+,+j+) of the matrix. That is: row +i+, column +j+.
#
def [](i, j)
- @rows[i][j]
+ (@rows[i] || [])[j]
end
alias element []
alias component []
@@ -277,26 +268,27 @@
end
#
- # Returns the number of columns. Note that it is possible to construct a
- # matrix with uneven columns (e.g. Matrix[ [1,2,3], [4,5] ]), but this is
- # mathematically unsound. This method uses the first row to determine the
- # result.
+ # Returns the number of columns.
#
- def column_size
- @rows[0].size
- end
+ 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) # :yield: e
- if block_given?
- for e in @rows[i]
- yield e
+ def row(i, &block) # :yield: e
+ if row = @rows[i]
+ if block_given?
+ row.each(&block)
+ else
+ Vector.elements(row, false)
end
else
- Vector.elements(@rows[i])
+ if block_given?
+ self
+ else
+ nil
+ end
end
end
@@ -307,11 +299,13 @@
#
def column(j) # :yield: e
if block_given?
- 0.upto(row_size - 1) do |i|
+ row_size.times do |i|
yield @rows[i][j]
- end
+ end unless j >= column_size
+ self
else
- col = (0 .. row_size - 1).collect {|i|
+ return nil if j >= column_size
+ col = (0 ... row_size).collect {|i|
@rows[i][j]
}
Vector.elements(col, false)
@@ -325,9 +319,10 @@
# => 1 4
# 9 16
#
- def collect # :yield: e
- rows = @rows.collect{|row| row.collect{|e| yield e}}
- Matrix.rows(rows, false)
+ def collect(&block) # :yield: e
+ return to_enum(:collect) unless block_given?
+ rows = @rows.collect{|row| row.collect(&block)}
+ new(rows)
end
alias map collect
@@ -350,18 +345,15 @@
size_col = param[1].end - from_col
size_col += 1 unless param[1].exclude_end?
when 4
- from_row = param[0]
- size_row = param[1]
- from_col = param[2]
- size_col = param[3]
+ from_row, size_row, from_col, size_col = param
else
Matrix.Raise ArgumentError, param.inspect
end
-
+ return nil if from_row > row_size || from_col > column_size
rows = @rows[from_row, size_row].collect{|row|
row[from_col, size_col]
}
- Matrix.rows(rows, false)
+ new(rows, column_size - from_col)
end
#--
@@ -383,8 +375,7 @@
end
#
- # Returns +true+ is this is a square matrix. See note in column_size about this
- # being unreliable, though.
+ # Returns +true+ is this is a square matrix.
#
def square?
column_size == row_size
@@ -399,46 +390,28 @@
#
def ==(other)
return false unless Matrix === other
-
- other.compare_by_row_vectors(@rows)
+ @rows == other.rows
end
+
def eql?(other)
return false unless Matrix === other
-
- other.compare_by_row_vectors(@rows, :eql?)
+ @rows.eql? other.rows
end
#
- # Not really intended for general consumption.
- #
- def compare_by_row_vectors(rows, comparison = :==)
- return false unless @rows.size == rows.size
-
- 0.upto(@rows.size - 1) do |i|
- return false unless @rows[i].send(comparison, rows[i])
- end
- true
- 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
- Matrix.rows(@rows)
+ new(@rows.map{|row| row.dup}, column_size)
end
#
# Returns a hash-code for the matrix.
#
def hash
- value = 0
- for row in @rows
- for e in row
- value ^= e.hash
- end
- end
- return value
+ @rows.hash
end
#--
@@ -459,7 +432,7 @@
e * m
}
}
- return Matrix.rows(rows, false)
+ return new(rows, column_size)
when Vector
m = Matrix.column_vector(m)
r = self * m
@@ -467,16 +440,16 @@
when Matrix
Matrix.Raise ErrDimensionMismatch if column_size != m.row_size
- rows = (0 .. row_size - 1).collect {|i|
- (0 .. m.column_size - 1).collect {|j|
+ rows = (0 ... row_size).collect {|i|
+ (0 ... m.column_size).collect {|j|
vij = 0
- 0.upto(column_size - 1) do |k|
+ column_size.times do |k|
vij += self[i, k] * m[k, j]
end
vij
}
}
- return Matrix.rows(rows, false)
+ return new(rows, m.column_size)
else
x, y = m.coerce(self)
return x * y
@@ -503,12 +476,12 @@
Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
- rows = (0 .. row_size - 1).collect {|i|
- (0 .. column_size - 1).collect {|j|
+ rows = (0 ... row_size).collect {|i|
+ (0 ... column_size).collect {|j|
self[i, j] + m[i, j]
}
}
- Matrix.rows(rows, false)
+ new(rows, column_size)
end
#
@@ -531,12 +504,12 @@
Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
- rows = (0 .. row_size - 1).collect {|i|
- (0 .. column_size - 1).collect {|j|
+ rows = (0 ... row_size).collect {|i|
+ (0 ... column_size).collect {|j|
self[i, j] - m[i, j]
}
}
- Matrix.rows(rows, false)
+ new(rows, column_size)
end
#
@@ -553,7 +526,7 @@
e / other
}
}
- return Matrix.rows(rows, false)
+ return new(rows, column_size)
when Matrix
return self * other.inverse
else
@@ -578,13 +551,13 @@
# Not for public consumption?
#
def inverse_from(src)
- size = row_size - 1
+ size = row_size
a = src.to_a
- for k in 0..size
+ for k in 0...size
i = k
akk = a[k][k].abs
- ((k+1)..size).each do |j|
+ ((k+1)...size).each do |j|
v = a[j][k].abs
if v > akk
i = j
@@ -598,23 +571,23 @@
end
akk = a[k][k]
- for i in 0 .. size
+ for i in 0 ... size
next if i == k
q = a[i][k].quo(akk)
a[i][k] = 0
- for j in (k + 1).. size
+ for j in (k + 1)... size
a[i][j] -= a[k][j] * q
end
- for j in 0..size
+ for j in 0...size
@rows[i][j] -= @rows[k][j] * q
end
end
- for j in (k + 1).. size
+ for j in (k + 1)...size
a[k][j] = a[k][j].quo(akk)
end
- for j in 0..size
+ for j in 0...size
@rows[k][j] = @rows[k][j].quo(akk)
end
end
@@ -630,11 +603,12 @@
# 48 99
#
def ** (other)
+ Matrix.Raise ErrDimensionMismatch unless square?
if other.kind_of?(Integer)
x = self
if other <= 0
- x = self.inverse
return Matrix.identity(self.column_size) if other == 0
+ x = self.inverse
other = -other
end
z = x
@@ -662,27 +636,26 @@
#
# Returns the determinant of the matrix. If the matrix is not square, the
- # result is 0. This method's algorism is Gaussian elimination method
+ # result is 0. This method's algorithm is Gaussian elimination method
# and using Numeric#quo(). Beware that using Float values, with their
# usual lack of precision, can affect the value returned by this method. Use
# Rational values or Matrix#det_e instead if this is important to you.
#
# Matrix[[7,6], [3,9]].determinant
- # => 63.0
+ # => 45.0
#
def determinant
return 0 unless square?
- size = row_size - 1
+ size = row_size
a = to_a
det = 1
- k = 0
- loop do
+ for k in 0...size do
if (akk = a[k][k]) == 0
i = k
loop do
- return 0 if (ii += 1) > size
+ return 0 if (i += 1) >= size
break unless a[i][k] == 0
end
a[i], a[k] = a[k], a[i]
@@ -690,14 +663,13 @@
det *= -1
end
- for i in k + 1 .. size
+ for i in k + 1 ... size
q = a[i][k].quo(akk)
- (k + 1).upto(size) do |j|
+ for j in k + 1 ... size
a[i][j] -= a[k][j] * q
end
end
det *= akk
- break unless (k += 1) <= size
end
det
end
@@ -705,8 +677,8 @@
#
# Returns the determinant of the matrix. If the matrix is not square, the
- # result is 0. This method's algorism is Gaussian elimination method.
- # This method uses Euclidean algorism. If all elements are integer,
+ # result is 0. This method's algorithm is Gaussian elimination method.
+ # This method uses Euclidean algorithm. If all elements are integer,
# really exact value. But, if an element is a float, can't return
# exact value.
#
@@ -716,25 +688,24 @@
def determinant_e
return 0 unless square?
- size = row_size - 1
+ size = row_size
a = to_a
det = 1
- k = 0
- loop do
+ for k in 0...size
if a[k][k].zero?
i = k
loop do
- return 0 if (i += 1) > size
+ return 0 if (i += 1) >= size
break unless a[i][k].zero?
end
a[i], a[k] = a[k], a[i]
det *= -1
end
- for i in (k + 1)..size
+ for i in (k + 1)...size
q = a[i][k].quo(a[k][k])
- k.upto(size) do |j|
+ for j in k...size do
a[i][j] -= a[k][j] * q
end
unless a[i][k].zero?
@@ -744,7 +715,6 @@
end
end
det *= a[k][k]
- break unless (k += 1) <= size
end
det
end
@@ -769,13 +739,12 @@
a_row_size = row_size
end
rank = 0
- k = 0
- loop do
+ for k in 0...a_column_size
if (akk = a[k][k]) == 0
i = k
exists = true
loop do
- if (i += 1) > a_column_size - 1
+ if (i += 1) >= a_column_size
exists = false
break
end
@@ -788,14 +757,14 @@
i = k
exists = true
loop do
- if (i += 1) > a_row_size - 1
+ if (i += 1) >= a_row_size
exists = false
break
end
break unless a[k][i] == 0
end
if exists
- k.upto(a_column_size - 1) do |j|
+ for j in k...a_column_size do
a[j][k], a[j][i] = a[j][i], a[j][k]
end
akk = a[k][k]
@@ -805,21 +774,20 @@
end
end
- for i in (k + 1)..(a_row_size - 1)
+ for i in (k + 1)...a_row_size
q = a[i][k].quo(akk)
- for j in (k + 1)..(a_column_size - 1)
+ for j in (k + 1)...a_column_size
a[i][j] -= a[k][j] * q
end
end
rank += 1
- break unless (k += 1) <= a_column_size - 1
end
return rank
end
#
# Returns the rank of the matrix. This method uses Euclidean
- # algorism. If all elements are integer, really exact value. But, if
+ # algorithm. If all elements are integer, really exact value. But, if
# an element is a float, can't return exact value.
#
# Matrix[[7,6], [3,9]].rank
@@ -830,7 +798,7 @@
a_column_size = column_size
a_row_size = row_size
pi = 0
- (0 ... a_column_size).each do |j|
+ a_column_size.times do |j|
if i = (pi ... a_row_size).find{|i0| !a[i0][j].zero?}
if i != pi
a[pi], a[i] = a[i], a[pi]
@@ -858,8 +826,9 @@
# => 16
#
def trace
+ Matrix.Raise ErrDimensionMismatch unless square?
tr = 0
- 0.upto(column_size - 1) do |i|
+ column_size.times do |i|
tr += @rows[i][i]
end
tr
@@ -877,7 +846,7 @@
# 2 4 6
#
def transpose
- Matrix.columns(@rows)
+ new @rows.transpose, row_size
end
alias t transpose
@@ -901,7 +870,7 @@
# Returns an array of the row vectors of the matrix. See Vector.
#
def row_vectors
- rows = (0 .. row_size - 1).collect {|i|
+ rows = (0 ... row_size).collect {|i|
row(i)
}
rows
@@ -911,7 +880,7 @@
# Returns an array of the column vectors of the matrix. See Vector.
#
def column_vectors
- columns = (0 .. column_size - 1).collect {|i|
+ columns = (0 ... column_size).collect {|i|
column(i)
}
columns
@@ -921,7 +890,7 @@
# Returns an array of arrays that describe the rows of the matrix.
#
def to_a
- @rows.collect{|row| row.collect{|e| e}}
+ @rows.collect{|row| row.dup}
end
def elements_to_f
@@ -944,18 +913,44 @@
# Overrides Object#to_s
#
def to_s
- "Matrix[" + @rows.collect{|row|
- "[" + row.collect{|e| e.to_s}.join(", ") + "]"
- }.join(", ")+"]"
+ if row_size == 0 && column_size > 0
+ "Matrix.columns([" + (["[]"] * column_size).join(", ") + "])"
+ else
+ "Matrix[" + @rows.collect{|row|
+ "[" + row.collect{|e| e.to_s}.join(", ") + "]"
+ }.join(", ")+"]"
+ end
end
#
# Overrides Object#inspect
#
def inspect
+ return to_s if row_size == 0 && column_size > 0
"Matrix"+@rows.inspect
end
+ #
+ # Converts the obj to an Array. If copy is set to true
+ # a copy of obj will be made if necessary.
+ #
+ def Matrix.convert_to_array(obj, copy = false)
+ 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 CLASS
class Scalar < Numeric # :nodoc:
@@ -1013,7 +1008,7 @@
when Vector
Scalar.Raise WrongArgType, other.class, "Numeric or Scalar or Matrix"
when Matrix
- self * other.inverse
+ self * other.inverse
else
x, y = other.coerce(self)
x.quo(y)
@@ -1082,13 +1077,14 @@
#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(:init_elements, array, copy = false)
+ new array
end
#
@@ -1096,27 +1092,16 @@
# whether the array itself or a copy is used internally.
#
def Vector.elements(array, copy = true)
- new(:init_elements, array, copy)
+ new Matrix.convert_to_array(array, copy)
end
#
# For internal use.
#
- def initialize(method, array, copy)
- self.send(method, array, copy)
+ def initialize(array)
+ @elements = array
end
- #
- # For internal use.
- #
- def init_elements(array, copy)
- if copy
- @elements = array.dup
- else
- @elements = array
- end
- end
-
# ACCESSING
#
@@ -1150,8 +1135,9 @@
# Iterate over the elements of this vector and +v+ in conjunction.
#
def each2(v) # :yield: e1, e2
+ return to_enum(:each2, v) unless block_given?
Vector.Raise ErrDimensionMismatch if size != v.size
- 0.upto(size - 1) do |i|
+ size.times do |i|
yield @elements[i], v[i]
end
end
@@ -1161,8 +1147,9 @@
# in conjunction.
#
def collect2(v) # :yield: e1, e2
+ return to_enum(:collect2, v) unless block_given?
Vector.Raise ErrDimensionMismatch if size != v.size
- (0 .. size - 1).collect do |i|
+ (0 ... size).collect do |i|
yield @elements[i], v[i]
end
end
@@ -1176,23 +1163,14 @@
#
def ==(other)
return false unless Vector === other
-
- other.compare_by(@elements)
+ @elements == other.elements
end
def eql?(other)
return false unless Vector === other
-
- other.compare_by(@elements, :eql?)
+ @elements.eql? other.elements
end
#
- # For internal use.
- #
- def compare_by(elements, comparison = :==)
- @elements.send(comparison, elements)
- end
-
- #
# Return a copy of the vector.
#
def clone
@@ -1285,10 +1263,9 @@
#
# Like Array#collect.
#
- def collect # :yield: e
- els = @elements.collect {|v|
- yield v
- }
+ def collect(&block) # :yield: e
+ return to_enum(:collect) unless block_given?
+ els = @elements.collect(&block)
Vector.elements(els, false)
end
alias map collect
@@ -1296,10 +1273,9 @@
#
# Like Vector#collect2, but returns a Vector instead of an Array.
#
- def map2(v) # :yield: e1, e2
- els = collect2(v) {|v1, v2|
- yield v1, v2
- }
+ def map2(v, &block) # :yield: e1, e2
+ return to_enum(:map2, v) unless block_given?
+ els = collect2(v, &block)
Vector.elements(els, false)
end
@@ -1376,6 +1352,5 @@
end
end
-
# Documentation comments:
# - Matrix#coerce and Vector#coerce need to be documented