Skip to content

fixes the det method for obect dtypes #572

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
+182 −188
Closed

fixes the det method for obect dtypes #572

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
214 changes: 107 additions & 107 deletions lib/nmatrix/math.rb
View file
Original file line number Diff line number Diff line change
@@ -65,8 +65,8 @@ def permutation_array_for(pivot_array)
# call-seq:
# invert! -> NMatrix
#
# Use LAPACK to calculate the inverse of the matrix (in-place) if available.
# Only works on dense matrices. Alternatively uses in-place Gauss-Jordan
# Use LAPACK to calculate the inverse of the matrix (in-place) if available.
# Only works on dense matrices. Alternatively uses in-place Gauss-Jordan
# elimination.
#
# * *Raises* :
@@ -116,8 +116,8 @@ def invert
# call-seq:
# adjugate! -> NMatrix
#
# Calculate the adjugate of the matrix (in-place).
# Only works on dense matrices.
# Calculate the adjugate of the matrix (in-place).
# Only works on dense matrices.
#
# * *Raises* :
# - +StorageTypeError+ -> only implemented on dense matrices.
@@ -130,7 +130,7 @@ def adjugate!
raise(DataTypeError, "Cannot calculate adjugate of an integer matrix in-place") if self.integer_dtype?
d = self.det
self.invert!
self.map! { |e| e * d }
self.map! { |e| e * d }
self
end
alias :adjoint! :adjugate!
@@ -139,13 +139,13 @@ def adjugate!
# call-seq:
# adjugate -> NMatrix
#
# Make a copy of the matrix and calculate the adjugate of the matrix.
# Only works on dense matrices.
# Make a copy of the matrix and calculate the adjugate of the matrix.
# Only works on dense matrices.
#
# * *Returns* :
# - A dense NMatrix. Will be the same type as the input NMatrix,
# except if the input is an integral dtype, in which case it will be a
# :float64 NMatrix.
# :float64 NMatrix.
#
# * *Raises* :
# - +StorageTypeError+ -> only implemented on dense matrices.
@@ -156,8 +156,8 @@ def adjugate
raise(ShapeError, "Cannot calculate adjugate of a non-square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1]
d = self.det
mat = self.invert
mat.map! { |e| e * d }
mat
mat.map! { |e| e * d }
mat
end
alias :adjoint :adjugate

@@ -204,13 +204,13 @@ def getrf!

#
# call-seq:
# geqrf! -> shape.min x 1 NMatrix
# geqrf! -> shape.min x 1 NMatrix
#
# QR factorization of a general M-by-N matrix +A+.
# QR factorization of a general M-by-N matrix +A+.
#
# The QR factorization is A = QR, where Q is orthogonal and R is Upper Triangular
# +A+ is overwritten with the elements of R and Q with Q being represented by the
# elements below A's diagonal and an array of scalar factors in the output NMatrix.
# +A+ is overwritten with the elements of R and Q with Q being represented by the
# elements below A's diagonal and an array of scalar factors in the output NMatrix.
#
# The matrix Q is represented as a product of elementary reflectors
# Q = H(1) H(2) . . . H(k), where k = min(m,n).
@@ -220,7 +220,7 @@ def getrf!
# H(i) = I - tau * v * v'
#
# http://www.netlib.org/lapack/explore-html/d3/d69/dgeqrf_8f.html
#
#
# Only works for dense matrices.
#
# * *Returns* :
@@ -232,29 +232,29 @@ def geqrf!
# The real implementation is in lib/nmatrix/lapacke.rb
raise(NotImplementedError, "geqrf! requires the nmatrix-lapacke gem")
end

#
# call-seq:
# ormqr(tau) -> NMatrix
# ormqr(tau, side, transpose, c) -> NMatrix
#
# Returns the product Q * c or c * Q after a call to geqrf! used in QR factorization.
# +c+ is overwritten with the elements of the result NMatrix if supplied. Q is the orthogonal matrix
# Returns the product Q * c or c * Q after a call to geqrf! used in QR factorization.
# +c+ is overwritten with the elements of the result NMatrix if supplied. Q is the orthogonal matrix
# represented by tau and the calling NMatrix
#
#
# Only works on float types, use unmqr for complex types.
#
# == Arguments
#
# * +tau+ - vector containing scalar factors of elementary reflectors
# * +side+ - direction of multiplication [:left, :right]
# * +transpose+ - apply Q with or without transpose [false, :transpose]
# * +transpose+ - apply Q with or without transpose [false, :transpose]
# * +c+ - NMatrix multplication argument that is overwritten, no argument assumes c = identity
#
# * *Returns* :
#
# - Q * c or c * Q Where Q may be transposed before multiplication.
#
# - Q * c or c * Q Where Q may be transposed before multiplication.
#
#
# * *Raises* :
# - +StorageTypeError+ -> LAPACK functions only work on dense matrices.
@@ -264,31 +264,31 @@ def geqrf!
def ormqr(tau, side=:left, transpose=false, c=nil)
# The real implementation is in lib/nmatrix/lapacke.rb
raise(NotImplementedError, "ormqr requires the nmatrix-lapacke gem")

end

#
# call-seq:
# unmqr(tau) -> NMatrix
# unmqr(tau, side, transpose, c) -> NMatrix
#
# Returns the product Q * c or c * Q after a call to geqrf! used in QR factorization.
# +c+ is overwritten with the elements of the result NMatrix if it is supplied. Q is the orthogonal matrix
# Returns the product Q * c or c * Q after a call to geqrf! used in QR factorization.
# +c+ is overwritten with the elements of the result NMatrix if it is supplied. Q is the orthogonal matrix
# represented by tau and the calling NMatrix
#
#
# Only works on complex types, use ormqr for float types.
#
# == Arguments
#
# * +tau+ - vector containing scalar factors of elementary reflectors
# * +side+ - direction of multiplication [:left, :right]
# * +transpose+ - apply Q as Q or its complex conjugate [false, :complex_conjugate]
# * +transpose+ - apply Q as Q or its complex conjugate [false, :complex_conjugate]
# * +c+ - NMatrix multplication argument that is overwritten, no argument assumes c = identity
#
# * *Returns* :
#
# - Q * c or c * Q Where Q may be transformed to its complex conjugate before multiplication.
#
# - Q * c or c * Q Where Q may be transformed to its complex conjugate before multiplication.
#
#
# * *Raises* :
# - +StorageTypeError+ -> LAPACK functions only work on dense matrices.
@@ -361,13 +361,13 @@ def factorize_cholesky
#
# LU factorization of a matrix. Optionally return the permutation matrix.
# Note that computing the permutation matrix will introduce a slight memory
# and time overhead.
#
# and time overhead.
#
# == Arguments
#
# +with_permutation_matrix+ - If set to *true* will return the permutation
#
# +with_permutation_matrix+ - If set to *true* will return the permutation
# matrix alongwith the LU factorization as a second return value.
#
#
def factorize_lu with_permutation_matrix=nil
raise(NotImplementedError, "only implemented for dense storage") unless self.stype == :dense
raise(NotImplementedError, "matrix is not 2-dimensional") unless self.dimensions == 2
@@ -383,9 +383,9 @@ def factorize_lu with_permutation_matrix=nil
# call-seq:
# factorize_qr -> [Q,R]
#
# QR factorization of a matrix without column pivoting.
# Q is orthogonal and R is upper triangular if input is square or upper trapezoidal if
# input is rectangular.
# QR factorization of a matrix without column pivoting.
# Q is orthogonal and R is upper triangular if input is square or upper trapezoidal if
# input is rectangular.
#
# Only works for dense matrices.
#
@@ -403,11 +403,11 @@ def factorize_qr
rows, columns = self.shape
r = self.clone
tau = r.geqrf!
#Obtain Q

#Obtain Q
q = self.complex_dtype? ? r.unmqr(tau) : r.ormqr(tau)
#Obtain R

#Obtain R
if rows <= columns
r.upper_triangle!
#Need to account for upper trapezoidal structure if R is a tall rectangle (rows > columns)
@@ -420,7 +420,7 @@ def factorize_qr
end

# Reduce self to upper hessenberg form using householder transforms.
#
#
# == References
#
# * http://en.wikipedia.org/wiki/Hessenberg_matrix
@@ -431,12 +431,12 @@ def hessenberg

# Destructive version of #hessenberg
def hessenberg!
raise ShapeError, "Trying to reduce non 2D matrix to hessenberg form" if
raise ShapeError, "Trying to reduce non 2D matrix to hessenberg form" if
shape.size != 2
raise ShapeError, "Trying to reduce non-square matrix to hessenberg form" if
raise ShapeError, "Trying to reduce non-square matrix to hessenberg form" if
shape[0] != shape[1]
raise StorageTypeError, "Matrix must be dense" if stype != :dense
raise TypeError, "Works with float matrices only" unless
raise TypeError, "Works with float matrices only" unless
[:float64,:float32].include?(dtype)

__hessenberg__(self)
@@ -446,28 +446,28 @@ def hessenberg!
# Solve the matrix equation AX = B, where A is +self+, B is the first
# argument, and X is returned. A must be a nxn square matrix, while B must be
# nxm. Only works with dense matrices and non-integer, non-object data types.
#
#
# == Arguments
#
# * +b+ - the right hand side
#
# == Options
#
# * +form+ - Signifies the form of the matrix A in the linear system AX=B.
# If not set then it defaults to +:general+, which uses an LU solver.
# If not set then it defaults to +:general+, which uses an LU solver.
# Other possible values are +:lower_tri+, +:upper_tri+ and +:pos_def+ (alternatively,
# non-abbreviated symbols +:lower_triangular+, +:upper_triangular+,
# and +:positive_definite+ can be used.
# If +:lower_tri+ or +:upper_tri+ is set, then a specialized linear solver for linear
# systems AX=B with a lower or upper triangular matrix A is used. If +:pos_def+ is chosen,
# non-abbreviated symbols +:lower_triangular+, +:upper_triangular+,
# and +:positive_definite+ can be used.
# If +:lower_tri+ or +:upper_tri+ is set, then a specialized linear solver for linear
# systems AX=B with a lower or upper triangular matrix A is used. If +:pos_def+ is chosen,
# then the linear system is solved via the Cholesky factorization.
# Note that when +:lower_tri+ or +:upper_tri+ is used, then the algorithm just assumes that
# all entries in the lower/upper triangle of the matrix are zeros without checking (which
# can be useful in certain applications).
#
#
#
# == Usage
#
#
# a = NMatrix.new [2,2], [3,1,1,2], dtype: dtype
# b = NMatrix.new [2,1], [9,8], dtype: dtype
# a.solve(b)
@@ -488,18 +488,18 @@ def hessenberg!
#
def solve(b, opts = {})
raise(ShapeError, "Must be called on square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1]
raise(ShapeError, "number of rows of b must equal number of cols of self") if
raise(ShapeError, "number of rows of b must equal number of cols of self") if
self.shape[1] != b.shape[0]
raise(ArgumentError, "only works with dense matrices") if self.stype != :dense
raise(ArgumentError, "only works for non-integer, non-object dtypes") if
raise(ArgumentError, "only works for non-integer, non-object dtypes") if
integer_dtype? or object_dtype? or b.integer_dtype? or b.object_dtype?

opts = { form: :general }.merge(opts)
x = b.clone
n = self.shape[0]
nrhs = b.shape[1]

case opts[:form]
case opts[:form]
when :general
clone = self.clone
ipiv = NMatrix::LAPACK.clapack_getrf(:row, n, n, clone, n)
@@ -536,15 +536,15 @@ def solve(b, opts = {})
# least_squares(b) -> NMatrix
# least_squares(b, tolerance: 10e-10) -> NMatrix
#
# Provides the linear least squares approximation of an under-determined system
# Provides the linear least squares approximation of an under-determined system
# using QR factorization provided that the matrix is not rank-deficient.
#
# Only works for dense matrices.
#
# * *Arguments* :
# - +b+ -> The solution column vector NMatrix of A * X = b.
# - +tolerance:+ -> Absolute tolerance to check if a diagonal element in A = QR is near 0
#
# - +tolerance:+ -> Absolute tolerance to check if a diagonal element in A = QR is near 0
#
# * *Returns* :
# - NMatrix that is a column vector with the LLS solution
#
@@ -554,8 +554,8 @@ def solve(b, opts = {})
#
# Examples :-
#
# a = NMatrix.new([3,2], [2.0, 0, -1, 1, 0, 2])
#
# a = NMatrix.new([3,2], [2.0, 0, -1, 1, 0, 2])
#
# b = NMatrix.new([3,1], [1.0, 0, -1])
#
# a.least_squares(b)
@@ -564,30 +564,30 @@ def solve(b, opts = {})
# [ -0.3333333333333334 ]
# ]
#
def least_squares(b, tolerance: 10e-6)
raise(ArgumentError, "least squares approximation only works for non-complex types") if
def least_squares(b, tolerance: 10e-6)
raise(ArgumentError, "least squares approximation only works for non-complex types") if
self.complex_dtype?

rows, columns = self.shape

raise(ShapeError, "system must be under-determined ( rows > columns )") unless
raise(ShapeError, "system must be under-determined ( rows > columns )") unless
rows > columns

#Perform economical QR factorization
r = self.clone
tau = r.geqrf!
q_transpose_b = r.ormqr(tau, :left, :transpose, b)

#Obtain R from geqrf! intermediate
r[0...columns, 0...columns].upper_triangle!
r[columns...rows, 0...columns] = 0

diagonal = r.diagonal

raise(ArgumentError, "rank deficient matrix") if diagonal.any? { |x| x == 0 }

if diagonal.any? { |x| x.abs < tolerance }
warn "warning: A diagonal element of R in A = QR is close to zero ;" <<
warn "warning: A diagonal element of R in A = QR is close to zero ;" <<
" indicates a possible loss of precision"
end

@@ -669,28 +669,28 @@ def gesdd(workspace_size=nil)
#
# In-place permute the columns of a dense matrix using LASWP according to the order given as an array +ary+.
#
# If +:convention+ is +:lapack+, then +ary+ represents a sequence of pair-wise permutations which are
# performed successively. That is, the i'th entry of +ary+ is the index of the column to swap
# the i'th column with, having already applied all earlier swaps.
# If +:convention+ is +:lapack+, then +ary+ represents a sequence of pair-wise permutations which are
# performed successively. That is, the i'th entry of +ary+ is the index of the column to swap
# the i'th column with, having already applied all earlier swaps.
#
# If +:convention+ is +:intuitive+, then +ary+ represents the order of columns after the permutation.
# That is, the i'th entry of +ary+ is the index of the column that will be in position i after the
# If +:convention+ is +:intuitive+, then +ary+ represents the order of columns after the permutation.
# That is, the i'th entry of +ary+ is the index of the column that will be in position i after the
# reordering (Matlab-like behaviour). This is the default.
#
# Not yet implemented for yale or list.
# Not yet implemented for yale or list.
#
# == Arguments
#
# * +ary+ - An Array specifying the order of the columns. See above for details.
#
#
# == Options
#
#
# * +:covention+ - Possible values are +:lapack+ and +:intuitive+. Default is +:intuitive+. See above for details.
#
def laswp!(ary, opts={})
raise(StorageTypeError, "ATLAS functions only work on dense matrices") unless self.dense?
opts = { convention: :intuitive }.merge(opts)

if opts[:convention] == :intuitive
if ary.length != ary.uniq.length
raise(ArgumentError, "No duplicated entries in the order array are allowed under convention :intuitive")
@@ -716,22 +716,22 @@ def laswp!(ary, opts={})
#
# Permute the columns of a dense matrix using LASWP according to the order given in an array +ary+.
#
# If +:convention+ is +:lapack+, then +ary+ represents a sequence of pair-wise permutations which are
# performed successively. That is, the i'th entry of +ary+ is the index of the column to swap
# If +:convention+ is +:lapack+, then +ary+ represents a sequence of pair-wise permutations which are
# performed successively. That is, the i'th entry of +ary+ is the index of the column to swap
# the i'th column with, having already applied all earlier swaps. This is the default.
#
# If +:convention+ is +:intuitive+, then +ary+ represents the order of columns after the permutation.
# That is, the i'th entry of +ary+ is the index of the column that will be in position i after the
# reordering (Matlab-like behaviour).
# If +:convention+ is +:intuitive+, then +ary+ represents the order of columns after the permutation.
# That is, the i'th entry of +ary+ is the index of the column that will be in position i after the
# reordering (Matlab-like behaviour).
#
# Not yet implemented for yale or list.
# Not yet implemented for yale or list.
#
# == Arguments
#
# * +ary+ - An Array specifying the order of the columns. See above for details.
#
#
# == Options
#
#
# * +:covention+ - Possible values are +:lapack+ and +:intuitive+. Default is +:lapack+. See above for details.
#
def laswp(ary, opts={})
@@ -767,7 +767,7 @@ def det
raise(ShapeError, "determinant can be calculated only for square matrices") unless self.dim == 2 && self.shape[0] == self.shape[1]

# Cast to a dtype for which getrf is implemented
new_dtype = self.integer_dtype? ? :float64 : self.dtype
new_dtype = self.object_dtype? || self.integer_dtype? ? :float64 : self.dtype
copy = self.cast(:dense, new_dtype)

# Need to know the number of permutations. We'll add up the diagonals of
@@ -810,23 +810,23 @@ def complex_conjugate(new_stype = self.stype)
end

# Calculate the variance co-variance matrix
#
#
# == Options
#
#
# * +:for_sample_data+ - Default true. If set to false will consider the denominator for
# population data (i.e. N, as opposed to N-1 for sample data).
#
#
# == References
#
#
# * http://stattrek.com/matrix-algebra/covariance-matrix.aspx
def cov(opts={})
raise TypeError, "Only works for non-integer dtypes" if integer_dtype?
opts = {
for_sample_data: true
}.merge(opts)

denominator = opts[:for_sample_data] ? rows - 1 : rows
ones = NMatrix.ones [rows,1]
ones = NMatrix.ones [rows,1]
deviation_scores = self - ones.dot(ones.transpose).dot(self) / rows
deviation_scores.transpose.dot(deviation_scores) / denominator
end
@@ -840,24 +840,24 @@ def corr

# Raise a square matrix to a power. Be careful of numeric overflows!
# In case *n* is 0, an identity matrix of the same dimension is returned. In case
# of negative *n*, the matrix is inverted and the absolute value of *n* taken
# of negative *n*, the matrix is inverted and the absolute value of *n* taken
# for computing the power.
#
#
# == Arguments
#
#
# * +n+ - Integer to which self is to be raised.
#
#
# == References
#
# * R.G Dromey - How to Solve it by Computer. Link -
#
# * R.G Dromey - How to Solve it by Computer. Link -
# http://www.amazon.com/Solve-Computer-Prentice-Hall-International-Science/dp/0134340019/ref=sr_1_1?ie=UTF8&qid=1422605572&sr=8-1&keywords=how+to+solve+it+by+computer
def pow n
raise ShapeError, "Only works with 2D square matrices." if
raise ShapeError, "Only works with 2D square matrices." if
shape[0] != shape[1] or shape.size != 2
raise TypeError, "Only works with integer powers" unless n.is_a?(Integer)

sequence = (integer_dtype? ? self.cast(dtype: :int64) : self).clone
product = NMatrix.eye shape[0], dtype: sequence.dtype, stype: sequence.stype
product = NMatrix.eye shape[0], dtype: sequence.dtype, stype: sequence.stype

if n == 0
return NMatrix.eye(shape, dtype: dtype, stype: stype)
@@ -885,8 +885,8 @@ def pow n
#
# * +mat+ - A 2D NMatrix object
#
# === Usage
#
# === Usage
#
# a = NMatrix.new([2,2],[1,2,
# 3,4])
# b = NMatrix.new([2,3],[1,1,1,
@@ -895,7 +895,7 @@ def pow n
# [1.0, 1.0, 1.0, 2.0, 2.0, 2.0]
# [3.0, 3.0, 3.0, 4.0, 4.0, 4.0]
# [3.0, 3.0, 3.0, 4.0, 4.0, 4.0] ]
#
#
def kron_prod(mat)
unless self.dimensions==2 and mat.dimensions==2
raise ShapeError, "Implemented for 2D NMatrix objects only."
@@ -920,7 +920,7 @@ def kron_prod(mat)
end
end

NMatrix.new([n,m], kron_prod_array)
NMatrix.new([n,m], kron_prod_array)
end

#
@@ -1169,7 +1169,7 @@ def scale!(alpha, incx=1, n=nil)
# - +n+ -> Number of elements of +vector+.
#
# Return the scaling result of the matrix. BLAS scal will be invoked if provided.

def scale(alpha, incx=1, n=nil)
return self.clone.scale!(alpha, incx, n)
end
156 changes: 75 additions & 81 deletions spec/math_spec.rb
View file
Original file line number Diff line number Diff line change
@@ -99,14 +99,14 @@
[:asin, :acos, :atan, :atanh].each do |atf|

it "should correctly apply elementwise #{atf}" do
expect(@m.send(atf)).to eq N.new(@size,
expect(@m.send(atf)).to eq N.new(@size,
@a.map{ |e| Math.send(atf, e) },
dtype: :float64, stype: stype)
end
end

it "should correctly apply elementtwise atan2" do
expect(@m.atan2(@m*0+1)).to eq N.new(@size,
expect(@m.atan2(@m*0+1)).to eq N.new(@size,
@a.map { |e| Math.send(:atan2, e, 1) }, dtype: :float64, stype: stype)
end

@@ -120,8 +120,8 @@
end
end
end
context "Floor and ceil for #{stype}" do

context "Floor and ceil for #{stype}" do

[:floor, :ceil].each do |meth|
ALL_DTYPES.each do |dtype|
@@ -134,7 +134,7 @@

if dtype.to_s.match(/int/) or [:byte, :object].include?(dtype)
it "should return #{dtype} for #{dtype}" do

expect(@m.send(meth)).to eq N.new(@size, @a.map { |e| e.send(meth) }, dtype: dtype, stype: stype)

if dtype == :object
@@ -147,10 +147,10 @@
it "should return dtype int64 for #{dtype}" do

expect(@m.send(meth)).to eq N.new(@size, @a.map { |e| e.send(meth) }, dtype: dtype, stype: stype)

expect(@m.send(meth).dtype).to eq :int64
end
elsif dtype.to_s.match(/complex/)
elsif dtype.to_s.match(/complex/)
it "should properly calculate #{meth} for #{dtype}" do

expect(@m.send(meth)).to eq N.new(@size, @a.map { |e| e = Complex(e.real.send(meth), e.imag.send(meth)) }, dtype: dtype, stype: stype)
@@ -169,35 +169,35 @@
context dtype do
before :each do
@size = [2,2]
@mat = NMatrix.new @size, [1.33334, 0.9998, 1.9999, -8.9999],
@mat = NMatrix.new @size, [1.33334, 0.9998, 1.9999, -8.9999],
dtype: dtype, stype: stype
@ans = @mat.to_a.flatten
end

it "rounds" do
expect(@mat.round).to eq(N.new(@size, @ans.map { |a| a.round},
expect(@mat.round).to eq(N.new(@size, @ans.map { |a| a.round},
dtype: dtype, stype: stype))
end unless(/complex/ =~ dtype)

it "rounds with args" do
expect(@mat.round(2)).to eq(N.new(@size, @ans.map { |a| a.round(2)},
expect(@mat.round(2)).to eq(N.new(@size, @ans.map { |a| a.round(2)},
dtype: dtype, stype: stype))
end unless(/complex/ =~ dtype)

it "rounds complex with args" do
puts @mat.round(2)
expect(@mat.round(2)).to be_within(0.0001).of(N.new [2,2], @ans.map {|a|
expect(@mat.round(2)).to be_within(0.0001).of(N.new [2,2], @ans.map {|a|
Complex(a.real.round(2), a.imag.round(2))},dtype: dtype, stype: stype)
end if(/complex/ =~ dtype)

it "rounds complex" do
expect(@mat.round).to eq(N.new [2,2], @ans.map {|a|
expect(@mat.round).to eq(N.new [2,2], @ans.map {|a|
Complex(a.real.round, a.imag.round)},dtype: dtype, stype: stype)
end if(/complex/ =~ dtype)
end
end
end

end
end
end
@@ -288,17 +288,17 @@
NON_INTEGER_DTYPES.each do |dtype|
next if dtype == :object
context dtype do

it "calculates QR decomposition using factorize_qr for a square matrix" do

a = NMatrix.new(3, [12.0, -51.0, 4.0,
6.0, 167.0, -68.0,
a = NMatrix.new(3, [12.0, -51.0, 4.0,
6.0, 167.0, -68.0,
-4.0, 24.0, -41.0] , dtype: dtype)

q_solution = NMatrix.new([3,3], Q_SOLUTION_ARRAY_2, dtype: dtype)

r_solution = NMatrix.new([3,3], [-14.0, -21.0, 14,
0.0, -175, 70,
r_solution = NMatrix.new([3,3], [-14.0, -21.0, 14,
0.0, -175, 70,
0.0, 0.0, -35] , dtype: dtype)

err = case dtype
@@ -307,30 +307,30 @@
when :float64, :complex128
1e-13
end

begin
q,r = a.factorize_qr

expect(q).to be_within(err).of(q_solution)
expect(r).to be_within(err).of(r_solution)

rescue NotImplementedError
pending "Suppressing a NotImplementedError when the lapacke plugin is not available"
end
end

it "calculates QR decomposition using factorize_qr for a tall and narrow rectangular matrix" do

a = NMatrix.new([4,2], [34.0, 21.0,
23.0, 53.0,
26.0, 346.0,
a = NMatrix.new([4,2], [34.0, 21.0,
23.0, 53.0,
26.0, 346.0,
23.0, 121.0] , dtype: dtype)

q_solution = NMatrix.new([4,4], Q_SOLUTION_ARRAY_1, dtype: dtype)

r_solution = NMatrix.new([4,2], [-53.75872022286244, -255.06559574252242,
0.0, 269.34836526051555,
0.0, 0.0,
r_solution = NMatrix.new([4,2], [-53.75872022286244, -255.06559574252242,
0.0, 269.34836526051555,
0.0, 0.0,
0.0, 0.0] , dtype: dtype)

err = case dtype
@@ -339,22 +339,22 @@
when :float64, :complex128
1e-13
end

begin
q,r = a.factorize_qr

expect(q).to be_within(err).of(q_solution)
expect(r).to be_within(err).of(r_solution)

rescue NotImplementedError
pending "Suppressing a NotImplementedError when the lapacke plugin is not available"
end
end

it "calculates QR decomposition using factorize_qr for a short and wide rectangular matrix" do

a = NMatrix.new([3,4], [123,31,57,81,92,14,17,36,42,34,11,28], dtype: dtype)

q_solution = NMatrix.new([3,3], Q_SOLUTION_ARRAY_3, dtype: dtype)

r_solution = NMatrix.new([3,4], R_SOLUTION_ARRAY, dtype: dtype)
@@ -365,37 +365,37 @@
when :float64, :complex128
1e-13
end

begin
q,r = a.factorize_qr

expect(q).to be_within(err).of(q_solution)
expect(r).to be_within(err).of(r_solution)

rescue NotImplementedError
pending "Suppressing a NotImplementedError when the lapacke plugin is not available"
end
end

it "calculates QR decomposition such that A - QR ~ 0" do

a = NMatrix.new([3,3], [ 9.0, 0.0, 26.0,
12.0, 0.0, -7.0,
a = NMatrix.new([3,3], [ 9.0, 0.0, 26.0,
12.0, 0.0, -7.0,
0.0, 4.0, 0.0] , dtype: dtype)

err = case dtype
when :float32, :complex64
1e-4
when :float64, :complex128
1e-13
end

begin
q,r = a.factorize_qr
a_expected = q.dot(r)
a_expected = q.dot(r)

expect(a_expected).to be_within(err).of(a)

rescue NotImplementedError
pending "Suppressing a NotImplementedError when the lapacke plugin is not available"
end
@@ -412,10 +412,10 @@
when :float64, :complex128
1e-13
end

begin
q,r = a.factorize_qr

#Q is orthogonal if Q x Q.transpose = I
product = q.dot(q.transpose)

@@ -444,9 +444,9 @@
end

it "should correctly invert a matrix in place (bang)" do
a = NMatrix.new(:dense, 5, [1, 8,-9, 7, 5,
0, 1, 0, 4, 4,
0, 0, 1, 2, 5,
a = NMatrix.new(:dense, 5, [1, 8,-9, 7, 5,
0, 1, 0, 4, 4,
0, 0, 1, 2, 5,
0, 0, 0, 1,-5,
0, 0, 0, 0, 1 ], dtype)
b = NMatrix.new(:dense, 5, [1,-8, 9, 7, 17,
@@ -613,7 +613,7 @@
ALL_DTYPES.each do |dtype|
next if integer_dtype?(dtype)
context "#cov dtype #{dtype}" do
before do
before do
@n = NMatrix.new( [5,3], [4.0,2.0,0.60,
4.2,2.1,0.59,
3.9,2.0,0.58,
@@ -622,15 +622,15 @@
end

it "calculates variance co-variance matrix (sample)" do
expect(@n.cov).to be_within(0.0001).of(NMatrix.new([3,3],
expect(@n.cov).to be_within(0.0001).of(NMatrix.new([3,3],
[0.025 , 0.0075, 0.00175,
0.0075, 0.007 , 0.00135,
0.00175, 0.00135 , 0.00043 ], dtype: dtype)
)
end

it "calculates variance co-variance matrix (population)" do
expect(@n.cov(for_sample_data: false)).to be_within(0.0001).of(NMatrix.new([3,3],
expect(@n.cov(for_sample_data: false)).to be_within(0.0001).of(NMatrix.new([3,3],
[2.0000e-02, 6.0000e-03, 1.4000e-03,
6.0000e-03, 5.6000e-03, 1.0800e-03,
1.4000e-03, 1.0800e-03, 3.4400e-04], dtype: dtype)
@@ -645,15 +645,15 @@
3.9,2.0,0.58,
4.3,2.1,0.62,
4.1,2.2,0.63], dtype: dtype)
expect(n.corr).to be_within(0.001).of(NMatrix.new([3,3],
expect(n.corr).to be_within(0.001).of(NMatrix.new([3,3],
[1.00000, 0.56695, 0.53374,
0.56695, 1.00000, 0.77813,
0.53374, 0.77813, 1.00000], dtype: dtype))
end unless dtype =~ /complex/
end

context "#symmetric? for #{dtype}" do
it "should return true for symmetric matrix" do
it "should return true for symmetric matrix" do
n = NMatrix.new([3,3], [1.00000, 0.56695, 0.53374,
0.56695, 1.00000, 0.77813,
0.53374, 0.77813, 1.00000], dtype: dtype)
@@ -662,7 +662,7 @@
end

context "#hermitian? for #{dtype}" do
it "should return true for complex hermitian or non-complex symmetric matrix" do
it "should return true for complex hermitian or non-complex symmetric matrix" do
n = NMatrix.new([3,3], [1.00000, 0.56695, 0.53374,
0.56695, 1.00000, 0.77813,
0.53374, 0.77813, 1.00000], dtype: dtype) unless dtype =~ /complex/
@@ -835,7 +835,7 @@
"Suppressing a NotImplementedError when the lapacke or atlas plugin is not available"
end
end

it "solves a linear system A * X = B with positive definite A and matrix B" do
b = NMatrix.new([3,2], [8,3,14,13,14,19], dtype: dtype)
begin
@@ -851,37 +851,37 @@

context "#least_squares" do
it "finds the least squares approximation to the equation A * X = B" do
a = NMatrix.new([3,2], [2.0, 0, -1, 1, 0, 2])
a = NMatrix.new([3,2], [2.0, 0, -1, 1, 0, 2])
b = NMatrix.new([3,1], [1.0, 0, -1])
solution = NMatrix.new([2,1], [1.0 / 3 , -1.0 / 3], dtype: :float64)

begin
least_squares = a.least_squares(b)
expect(least_squares).to be_within(0.0001).of solution
rescue NotImplementedError
"Suppressing a NotImplementedError when the lapacke or atlas plugin is not available"
end
end

it "finds the least squares approximation to the equation A * X = B with high tolerance" do
a = NMatrix.new([4,2], [1.0, 1, 1, 2, 1, 3,1,4])
a = NMatrix.new([4,2], [1.0, 1, 1, 2, 1, 3,1,4])
b = NMatrix.new([4,1], [6.0, 5, 7, 10])
solution = NMatrix.new([2,1], [3.5 , 1.4], dtype: :float64)

begin
least_squares = a.least_squares(b, tolerance: 10e-5)
expect(least_squares).to be_within(0.0001).of solution
rescue NotImplementedError
"Suppressing a NotImplementedError when the lapacke or atlas plugin is not available"
end
end
end
end

context "#hessenberg" do
FLOAT_DTYPES.each do |dtype|
context dtype do
before do
@n = NMatrix.new [5,5],
@n = NMatrix.new [5,5],
[0, 2, 0, 1, 1,
2, 2, 3, 2, 2,
4,-3, 0, 1, 3,
@@ -890,7 +890,7 @@
end

it "transforms a matrix to Hessenberg form" do
expect(@n.hessenberg).to be_within(0.0001).of(NMatrix.new([5,5],
expect(@n.hessenberg).to be_within(0.0001).of(NMatrix.new([5,5],
[0.00000,-1.66667, 0.79432,-0.45191,-1.54501,
-9.00000, 2.95062,-6.89312, 3.22250,-0.19012,
0.00000,-8.21682,-0.57379, 5.26966,-1.69976,
@@ -905,20 +905,20 @@
[:dense, :yale].each do |stype|
answer_dtype = integer_dtype?(dtype) ? :int64 : dtype
next if dtype == :byte

context "#pow #{dtype} #{stype}" do
before do
before do
@n = NMatrix.new [4,4], [0, 2, 0, 1,
2, 2, 3, 2,
4,-3, 0, 1,
6, 1,-6,-5], dtype: dtype, stype: stype
end

it "raises a square matrix to even power" do
expect(@n.pow(4)).to eq(NMatrix.new([4,4], [292, 28,-63, -42,
expect(@n.pow(4)).to eq(NMatrix.new([4,4], [292, 28,-63, -42,
360, 96, 51, -14,
448,-231,-24,-87,
-1168, 595,234, 523],
-1168, 595,234, 523],
dtype: answer_dtype,
stype: stype))
end
@@ -934,14 +934,14 @@
it "raises a sqaure matrix to negative power" do
expect(@n.pow(-3)).to be_within(0.00001).of (NMatrix.new([4,4],
[1.0647e-02, 4.2239e-04,-6.2281e-05, 2.7680e-03,
-1.6415e-02, 2.1296e-02, 1.0718e-02, 4.8589e-03,
-1.6415e-02, 2.1296e-02, 1.0718e-02, 4.8589e-03,
8.6956e-03,-8.6569e-03, 2.8993e-02, 7.2015e-03,
5.0034e-02,-1.7500e-02,-3.6777e-02,-1.2128e-02], dtype: answer_dtype,
stype: stype))
5.0034e-02,-1.7500e-02,-3.6777e-02,-1.2128e-02], dtype: answer_dtype,
stype: stype))
end unless stype =~ /yale/ or dtype == :object or ALL_DTYPES.grep(/int/).include? dtype

it "raises a square matrix to zero" do
expect(@n.pow(0)).to eq(NMatrix.eye([4,4], dtype: answer_dtype,
expect(@n.pow(0)).to eq(NMatrix.eye([4,4], dtype: answer_dtype,
stype: stype))
end

@@ -955,15 +955,15 @@
ALL_DTYPES.each do |dtype|
[:dense, :yale].each do |stype|
context "#kron_prod #{dtype} #{stype}" do
before do
before do
@a = NMatrix.new([2,2], [1,2,
3,4], dtype: dtype, stype: stype)
@b = NMatrix.new([2,3], [1,1,1,
1,1,1], dtype: dtype, stype: stype)
@c = NMatrix.new([4,6], [1, 1, 1, 2, 2, 2,
1, 1, 1, 2, 2, 2,
3, 3, 3, 4, 4, 4,
3, 3, 3, 4, 4, 4], dtype: dtype, stype: stype)
3, 3, 3, 4, 4, 4], dtype: dtype, stype: stype)
end
it "Compute the Kronecker product of two NMatrix" do
expect(@a.kron_prod(@b)).to eq(@c)
@@ -995,19 +995,13 @@
end
end
it "computes the determinant of 2x2 matrix" do
if dtype != :object
expect(@a.det).to be_within(@err).of(-2)
end
end
it "computes the determinant of 3x3 matrix" do
if dtype != :object
expect(@b.det).to be_within(@err).of(-8)
end
end
it "computes the determinant of 4x4 matrix" do
if dtype != :object
expect(@c.det).to be_within(@err).of(-18)
end
end
it "computes the exact determinant of 2x2 matrix" do
if dtype == :byte