1berry/ruby
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Promote matrix to the bundled gems

Browse Source
This commit is contained in:
Hiroshi SHIBATA 2021-05-26 15:36:16 +09:00
parent c9178c1127
commit 454a36794f
Notes: git 2021-05-27 14:42:39 +09:00 
Merged: https://github.com/ruby/ruby/pull/4530
12 changed files with 4 additions and 4855 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
doc/maintainers.rdoc
View File

@ -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]

2
doc/standard_library.rdoc
View File

@ -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

1
gems/bundled_gems
View File

@ -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

2493
lib/matrix.rb
View File

File diff suppressed because it is too large Load Diff

882
lib/matrix/eigenvalue_decomposition.rb
View File

@ -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

219
lib/matrix/lup_decomposition.rb
View File

@ -1,219 +0,0 @@
# frozen_string_literal: false
class Matrix
# Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
#
# For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
# unit lower triangular matrix L, an n-by-n upper triangular matrix U,
# and a m-by-m permutation matrix P so that L*U = P*A.
# If m < n, then L is m-by-m and U is m-by-n.
#
# The LUP decomposition with pivoting always exists, even if the matrix is
# singular, so the constructor will never fail. The primary use of the
# LU decomposition is in the solution of square systems of simultaneous
# linear equations. This will fail if singular? returns true.
#
class LUPDecomposition
# Returns the lower triangular factor +L+
include Matrix::ConversionHelper
def l
Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j|
if (i > j)
@lu[i][j]
elsif (i == j)
1
else
0
end
end
end
# Returns the upper triangular factor +U+
def u
Matrix.build([@column_count, @row_count].min, @column_count) do |i, j|
if (i <= j)
@lu[i][j]
else
0
end
end
end
# Returns the permutation matrix +P+
def p
rows = Array.new(@row_count){Array.new(@row_count, 0)}
@pivots.each_with_index{|p, i| rows[i][p] = 1}
Matrix.send :new, rows, @row_count
end
# Returns +L+, +U+, +P+ in an array
def to_ary
[l, u, p]
end
alias_method :to_a, :to_ary
# Returns the pivoting indices
attr_reader :pivots
# Returns +true+ if +U+, and hence +A+, is singular.
def singular?
@column_count.times do |j|
if (@lu[j][j] == 0)
return true
end
end
false
end
# Returns the determinant of +A+, calculated efficiently
# from the factorization.
def det
if (@row_count != @column_count)
raise Matrix::ErrDimensionMismatch
end
d = @pivot_sign
@column_count.times do |j|
d *= @lu[j][j]
end
d
end
alias_method :determinant, :det
# Returns +m+ so that <tt>A*m = b</tt>,
# or equivalently so that <tt>L*U*m = P*b</tt>
# +b+ can be a Matrix or a Vector
def solve b
if (singular?)
raise Matrix::ErrNotRegular, "Matrix is singular."
end
if b.is_a? Matrix
if (b.row_count != @row_count)
raise Matrix::ErrDimensionMismatch
end
# Copy right hand side with pivoting
nx = b.column_count
m = @pivots.map{|row| b.row(row).to_a}
# Solve L*Y = P*b
@column_count.times do |k|
(k+1).upto(@column_count-1) do |i|
nx.times do |j|
m[i][j] -= m[k][j]*@lu[i][k]
end
end
end
# Solve U*m = Y
(@column_count-1).downto(0) do |k|
nx.times do |j|
m[k][j] = m[k][j].quo(@lu[k][k])
end
k.times do |i|
nx.times do |j|
m[i][j] -= m[k][j]*@lu[i][k]
end
end
end
Matrix.send :new, m, nx
else # same algorithm, specialized for simpler case of a vector
b = convert_to_array(b)
if (b.size != @row_count)
raise Matrix::ErrDimensionMismatch
end
# Copy right hand side with pivoting
m = b.values_at(*@pivots)
# Solve L*Y = P*b
@column_count.times do |k|
(k+1).upto(@column_count-1) do |i|
m[i] -= m[k]*@lu[i][k]
end
end
# Solve U*m = Y
(@column_count-1).downto(0) do |k|
m[k] = m[k].quo(@lu[k][k])
k.times do |i|
m[i] -= m[k]*@lu[i][k]
end
end
Vector.elements(m, false)
end
end
def initialize a
raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
# Use a "left-looking", dot-product, Crout/Doolittle algorithm.
@lu = a.to_a
@row_count = a.row_count
@column_count = a.column_count
@pivots = Array.new(@row_count)
@row_count.times do |i|
@pivots[i] = i
end
@pivot_sign = 1
lu_col_j = Array.new(@row_count)
# Outer loop.
@column_count.times do |j|
# Make a copy of the j-th column to localize references.
@row_count.times do |i|
lu_col_j[i] = @lu[i][j]
end
# Apply previous transformations.
@row_count.times do |i|
lu_row_i = @lu[i]
# Most of the time is spent in the following dot product.
kmax = [i, j].min
s = 0
kmax.times do |k|
s += lu_row_i[k]*lu_col_j[k]
end
lu_row_i[j] = lu_col_j[i] -= s
end
# Find pivot and exchange if necessary.
p = j
(j+1).upto(@row_count-1) do |i|
if (lu_col_j[i].abs > lu_col_j[p].abs)
p = i
end
end
if (p != j)
@column_count.times do |k|
t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
end
k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
@pivot_sign = -@pivot_sign
end
# Compute multipliers.
if (j < @row_count && @lu[j][j] != 0)
(j+1).upto(@row_count-1) do |i|
@lu[i][j] = @lu[i][j].quo(@lu[j][j])
end
end
end
end
end
end

26
lib/matrix/matrix.gemspec
View File

@ -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

5
lib/matrix/version.rb
View File

@ -1,5 +0,0 @@
# frozen_string_literal: true
class Matrix
VERSION = "0.4.1"
end

1
misc/expand_tabs.rb
View File

@ -78,7 +78,6 @@ DEFAULT_GEM_LIBS = %w[
ipaddr
irb
logger
matrix
mutex_m
net-ftp
net-http

888
test/matrix/test_matrix.rb
View File

@ -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

335
test/matrix/test_vector.rb
View File

@ -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

1
tool/sync_default_gems.rb
View File

@ -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",