diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index 778b4486..b84c2abf 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* : @@ -214,13 +214,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 @@ -233,7 +233,7 @@ def factorize_lu with_permutation_matrix=nil end # Reduce self to upper hessenberg form using householder transforms. - # + # # == References # # * http://en.wikipedia.org/wiki/Hessenberg_matrix @@ -244,12 +244,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) @@ -260,18 +260,18 @@ def hessenberg! # 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. - # + # # == Usage - # + # # a = NMatrix.new [2,2], [3,1,1,2], dtype: dtype # b = NMatrix.new [2,1], [9,8], dtype: dtype # a.solve(b) def solve b 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? x = b.clone @@ -357,28 +357,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") @@ -404,22 +404,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={}) @@ -498,23 +498,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 @@ -528,24 +528,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) @@ -573,8 +573,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, @@ -583,7 +583,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." @@ -608,7 +608,7 @@ def kron_prod(mat) end end - NMatrix.new([n,m], kron_prod_array) + NMatrix.new([n,m], kron_prod_array) end # @@ -793,6 +793,47 @@ def abs end + # + # call-seq: + # norm -> Numeric + # + # Calculates the selected norm (defaults to 2-norm) of a 2D matrix. + # + # This should be used for small or medium sized matrices. + # For greater matrices, there should be a separate implementation where + # the norm is estimated rather than computed, for the sake of computation speed. + # + # Currently implemented norms are 1-norm, 2-norm, Frobenius, Infinity. + # A minus on the 1, 2 and inf norms returns the minimum instead of the maximum value. + # + # Tested mainly with dense matrices. Further checks and modifications might + # be necessary for sparse matrices. + # + # * *Returns* : + # - The selected norm of the matrix. + # * *Raises* : + # - +NotImplementedError+ -> norm can be calculated only for 2D matrices + # - +ArgumentError+ -> unrecognized norm + # + def norm type = 2 + raise(NotImplementedError, "norm can be calculated only for 2D matrices") unless self.dim == 2 + raise(NotImplementedError, "norm only implemented for dense storage") unless self.stype == :dense + + case type + when nil, 2, -2 + return self.two_norm (type == -2) + when 1, -1 + return self.one_norm (type == -1) + when :frobenius, :fro + return self.fro_norm + when :infinity, :inf, :'-inf', :'-infinity' + return self.inf_norm (type == :'-inf' || type == :'-infinity') + else + raise ArgumentError.new("argument must be a valid integer or symbol") + end + end + + # # call-seq: # absolute_sum -> Numeric @@ -1045,4 +1086,53 @@ def dtype_for_floor_or_ceil self.__dense_map__ { |l| l.send(op,rhs) } end end + + # Norm calculation methods + # Frobenius norm: the Euclidean norm of the matrix, treated as if it were a vector + def fro_norm + #float64 has to be used in any case, since nrm2 will not yield correct result for float32 + self_cast = self.cast(:dtype => :float64) unless self.dtype == :float64 + + column_vector = self_cast.reshape([self.size, 1]) + + return column_vector.nrm2 + end + + # 2-norm: the largest/smallest singular value of the matrix + def two_norm minus = false + #self.dtype == :int32 ? self_cast = self.cast(:dtype => :float32) : self_cast = self.cast(:dtype => :float64) + self_cast = self.cast(:dtype => :float64) + + #TODO: confirm if this is the desired svd calculation + svd = self_cast.gesvd + return svd[1][0, 0] unless minus + return svd[1][svd[1].rows-1, svd[1].cols-1] + end + + # 1-norm: the maximum/minimum absolute column sum of the matrix + def one_norm minus = false + #TODO: change traversing method for sparse matrices + number_of_columns = self.cols + col_sums = [] + + number_of_columns.times do |i| + col_sums << self.col(i).inject(0) { |sum, number| sum += number.abs} + end + + return col_sums.sort!.last unless minus + return col_sums.sort!.first + end + + # Infinity norm: the maximum/minimum absolute row sum of the matrix + def inf_norm minus = false + number_of_rows = self.rows + row_sums = [] + + number_of_rows.times do |i| + row_sums << self.row(i).inject(0) { |sum, number| sum += number.abs} + end + + return row_sums.sort!.last unless minus + return row_sums.sort!.first + end end diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 89c27c63..16e9afd7 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 @@ -298,9 +298,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, @@ -424,7 +424,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, @@ -433,7 +433,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) @@ -441,7 +441,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) @@ -456,7 +456,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)) @@ -464,7 +464,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) @@ -473,7 +473,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/ @@ -587,7 +587,7 @@ 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, @@ -596,7 +596,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, @@ -611,9 +611,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, @@ -621,11 +621,11 @@ 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], - dtype: answer_dtype, + -1168, 595,234, 523], + dtype: answer_dtype, stype: stype)) end @@ -640,14 +640,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 @@ -661,7 +661,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, @@ -669,7 +669,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) diff --git a/spec/plugins/atlas/atlas_spec.rb b/spec/plugins/atlas/atlas_spec.rb index 3448603a..c9a8ee88 100644 --- a/spec/plugins/atlas/atlas_spec.rb +++ b/spec/plugins/atlas/atlas_spec.rb @@ -239,4 +239,39 @@ end end end + + context "math#norm" do + NORM_TOLERANCE = 1.0e-10 + + it "should default to 2-norm" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm).to be_within(NORM_TOLERANCE).of(7.348469228349535) + end + + it "should reject invalid arguments" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + + expect{n.norm(0.5)}.to raise_error(ArgumentError) + end + + it "should calculate 1 and 2 norms correctly" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm(1)).to eq(7) + #FIXME: change to the correct value when overflow issue is resolved + #expect(n.norm(-2)).to eq(1.8628605857884395e-07) + expect(n.norm(-2)).to be_within(NORM_TOLERANCE).of(0.0) + expect(n.norm(-1)).to eq(6) + end + + it "should calculate infinity norms correctly" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm(:inf)).to eq(9) + expect(n.norm(:'-inf')).to eq(2) + end + + it "should calculate frobenius norms correctly" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm(:fro)).to be_within(NORM_TOLERANCE).of(7.745966692414834) + end + end end diff --git a/spec/plugins/lapacke/lapacke_spec.rb b/spec/plugins/lapacke/lapacke_spec.rb index fe4fc8fe..a05a7ff3 100644 --- a/spec/plugins/lapacke/lapacke_spec.rb +++ b/spec/plugins/lapacke/lapacke_spec.rb @@ -300,4 +300,39 @@ end end end + + context "math#norm" do + NORM_TOLERANCE = 1.0e-10 + + it "should default to 2-norm" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm).to be_within(NORM_TOLERANCE).of(7.348469228349535) + end + + it "should reject invalid arguments" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + + expect{n.norm(0.5)}.to raise_error(ArgumentError) + end + + it "should calculate 1 and 2 norms correctly" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm(1)).to eq(7) + #FIXME: change to the correct value when overflow issue is resolved + #expect(n.norm(-2)).to eq(1.8628605857884395e-07) + expect(n.norm(-2)).to be_within(NORM_TOLERANCE).of(0.0) + expect(n.norm(-1)).to eq(6) + end + + it "should calculate infinity norms correctly" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm(:inf)).to eq(9) + expect(n.norm(:'-inf')).to eq(2) + end + + it "should calculate frobenius norms correctly" do + n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4]) + expect(n.norm(:fro)).to be_within(NORM_TOLERANCE).of(7.745966692414834) + end + end end