Project

General

Profile

Download (21.2 KB)

Bug #1532 » matrix.patch

marcandre (Marc-Andre Lafortune), 05/29/2009 04:18 AM

View differences:

lib/matrix.rb (working copy)
#--
# 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
......
# 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 <tt>require 'mathn'</tt>. This may be fixed in the future.
#
# == Method Catalogue
......
# * <tt> Matrix.columns(columns) </tt>
# * <tt> Matrix.diagonal(*values) </tt>
# * <tt> Matrix.scalar(n, value) </tt>
# * <tt> Matrix.scalar(n, value) </tt>
# * <tt> Matrix.identity(n) </tt>
# * <tt> Matrix.unit(n) </tt>
# * <tt> Matrix.I(n) </tt>
......
# * <tt> #inspect </tt>
#
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.
......
# -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
#
......
# => 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
#
......
#
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
#
......
# => 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
#
......
# 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 []
......
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
......
#
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)
......
# => 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
......
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
#--
......
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
......
#
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
#--
......
e * m
}
}
return Matrix.rows(rows, false)
return new(rows, column_size)
when Vector
m = Matrix.column_vector(m)
r = self * m
......
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
......
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
#
......
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
#
......
e / other
}
}
return Matrix.rows(rows, false)
return new(rows, column_size)
when Matrix
return self * other.inverse
else
......
# 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
......
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
......
# 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
......
#
# 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]
......
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
......
#
# 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.
#
......
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?
......
end
end
det *= a[k][k]
break unless (k += 1) <= size
end
det
end
......
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
......
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]
......
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
......
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]
......
# => 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
......
# 2 4 6
#
def transpose
Matrix.columns(@rows)
new @rows.transpose, row_size
end
alias t transpose
......
# 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
......
# 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
......
# 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
......
# 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:
......
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)
......
#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
#
......
# 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
#
......
# 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
......
# 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
......
#
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
......
#
# 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
......
#
# 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
......
end
end
# Documentation comments:
# - Matrix#coerce and Vector#coerce need to be documented
(1-1/7)