diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index c3e7fbb9..c4683483 100644 --- a/lib/nmatrix/math.rb +++ b/lib/nmatrix/math.rb @@ -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,7 +264,7 @@ 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 # @@ -272,23 +272,23 @@ def ormqr(tau, side=:left, transpose=false, c=nil) # 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,7 +446,7 @@ 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 @@ -454,20 +454,20 @@ def hessenberg! # == 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,10 +488,10 @@ 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) @@ -499,7 +499,7 @@ def solve(b, opts = {}) 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 diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 66783c7e..b179aa8f 100644 --- a/spec/math_spec.rb +++ b/spec/math_spec.rb @@ -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,13 +307,13 @@ 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 @@ -321,16 +321,16 @@ 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,7 +622,7 @@ 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) @@ -630,7 +630,7 @@ 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,7 +645,7 @@ 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)) @@ -653,7 +653,7 @@ 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,10 +851,10 @@ 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 @@ -862,26 +862,26 @@ "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,9 +905,9 @@ [: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, @@ -915,10 +915,10 @@ 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,7 +955,7 @@ 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, @@ -963,7 +963,7 @@ @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